Journal of Combustion
Volume 2020 (2020), Article ID 2764523, 23 pages
https://doi.org/10.1155/2020/2764523

Research Article

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

Adrien Chatelier,¹ Benoît Fiorina,¹ Vincent Moureau,² and Nicolas Bertier³

¹Université Paris-Saclay, CNRS, CentraleSupélec, Laboratoire EM2C, 91190 Gif-sur-Yvette, France
²Normandie Univ, INSA Rouen, UNIROUEN, CNRS, CORIA, 76000 Rouen, France
³ONERA—The French Aerospace Lab—Centre de Châtillon, BP 72, 92322 Chatillon Cedex, France

Correspondence should be addressed to Benoît Fiorina; benoit.fiorina@centralesupelec.fr

Received 10 September 2019; Revised 9 January 2020; Accepted 20 January 2020; Published 19 March 2020

Academic Editor: Benjamin Shaw

Copyright © 2020 Adrien Chatelier et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

This work presents Large Eddy Simulations of the unconfined CORIA Rouen Spray Burner, fed with liquid n-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.

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

Liquid n-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 and a mixture fraction Z as follows:

where TAB indicates that the variable φ is stored in a look-up table. The progress variable is defined as , where is the weight associated to species mass fraction . φ may include chemical reaction rates, species mass fractions, density but also thermodynamical and transport properties such as the heat capacity 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 , , and in a chemical look-up table such as

where and 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 temperature and the density are tabulated in the filtered chemical look-up table given by equation (2), as any other thermochemical variables [28]. and are solutions of the two following balance equations:

where ρ is the density, is the turbulent viscosity, is the turbulent Schmidt number, is the source term of mixture fraction due to the evaporation of the spray, is the subgrid scale flame wrinkling, is the progress variable diffusion factor, is the density in fresh gases, is the diffusion coefficient in fresh gases, is the progress variable unresolved convective fluxes due to thermal expansion, and is the progress variable reaction rate.

The functions , , and 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:

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 :

where is the adiabatic consumption speed of a freely propagating laminar premixed flame. The model for is modeled using the Charlette et al. formulation [29]:

where and are the subgrid scale Reynolds number and velocity fluctuations, respectively, while 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 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:

where is the particle position, is the particle velocity, u is the gas velocity, is the particle mass, is the drag coefficient, ν is the kinematic viscosity, is the particle density, and 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

where is the particle diameter, is the diffusion coefficient, is the Sherwood number, and is the Spalding mass number. The term in the mixture fraction equation is obtained by adding the mass transfer contribution of all the droplets around each node of the mesh:

where is the volume around the node. The other droplet parameters are derived by integrating either the droplet mass or energy equations. Droplet temperature and diameter are obtained by solving the following set of equations:

where is the thermal characteristic time of the Spalding model, is the gas temperature in the far field, is the latent heat of vaporization of the fuel, is the Spalding thermal number, is the heat capacity at a classical reference state assuming a one-third/two-third equilibrium between the far field and the droplet surface, is the dynamic viscosity at the same reference state, and is the Schmidt number.

3. Experimental Configuration

The experimental configuration is an n-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. 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. 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 Figure 3 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.

(a)

(b)

Figure 3: Instantaneous (a) and mean (b) OH-PLIF shots, from [19].

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 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 . Consequently, by assuming the flamelet regime, the flame wrinkling is also fully resolved on the LES grid and . At the opposite, the flamelet library is filtered in Case B by using a filter width 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) 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 at Z = 10, 20, and 40 mm. Blue lines: axial profiles at r = 0, 5, 10, and 15 mm.

5.1. Flow Topology and Gas Velocity

Figures 6 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.

(a)

(b)

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.

Figures 8 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.

(a)

(b)

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.

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. 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 . The axial velocity in nonreacting conditions drops rapidly to 0 m/s between and while in reacting conditions, the axial velocity decreases slowly between and . 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 .

(a)

(b)

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.

The results for the radial velocity component are plotted in Figure 12 (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 , the peaks at 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.

(a)

(b)

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.

5.2. Flame Topology

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. 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.

(a)

(b)

Figure 14: Normalized OH mass fraction (Case A top, case B middle). Experiments (bottom): OH-PLIF, from [19]. (a) Instantaneous signal, (b) mean signal.

The 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) . The mean inner flame front position is well captured by the simulations and is located at up to . The mean OH-PLIF measurements show that the outer flame front extends up to . 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 mm and quickly merges with the inner flame front for .

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 Figure 10. 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, Figures 15 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.

(a)

(b)

Figure 15: Contour of instantaneous temperature colored by the instantaneous velocity magnitude for cases A-R (a) and B-R (b).

(a)

(b)

Figure 16: Contour of mean temperature colored by the mean velocity magnitude for cases A-R (a) and B-R (b).

Figures 17 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).

(a)

(b)

Figure 17: Contour of instantaneous temperature colored by the instantaneous velocity magnitude for cases A-R (a) and B-R (b).

(a)

(b)

Figure 18: Contour of mean temperature colored by the mean velocity magnitude for cases A-R (a) and B-R (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 Figure 14. 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 is while case A recovers a lift-off of and case B a lift-off of . 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 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)

Estimated relative CPU cost
—
1
10
0.5
0.1

5.3. Droplets Diameter

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. 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 , especially on the outer region.

(a)

(b)

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).

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. 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.

(a)

(b)

Figure 21: Radial profiles of droplet diameter for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case 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 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.

5.4. Spray Velocity

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.

(a)

(b)

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).

Droplet axial velocity is reported in Figure 24 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.

(a)

(b)

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.

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.

(a)

(b)

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.

5.5. Spray Temperature

Figures 26 and 27 show the particles in the central vertical plane colored by their temperature for the cold and reacting cases, respectively. The scale is for the cold case and for the reacting cases. In Figure 26, the small droplets temperature decreases rapidly to 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 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.

(a)

(b)

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).

The droplet temperature predicted by the LES is now compared with the Global Rainbow Technique (GRT) measurements, whose uncertainty is [19]. Figure 28 presents radial profiles of temperature for the cold (left) and reacting (right) configurations.

(a)

(b)

Figure 28: Radial profiles of droplet temperature for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case B.

The experimental data highlight two zones. For , 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 , is expressed as

where is the droplet density; , its diameter; , the Schmidt number; , the Sherwood number; , the heat capacity at constant pressure of the n-heptane; and , the heat capacity at constant pressure and the dynamic viscosity of the mixture according to the 1/3-2/3 rule (see Chapter 1); , the thermal Spalding number; and , the mass Spalding number. As 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 , the droplets (which are small at this radial position) temperature drops quickly to . 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.

(a)

(b)

Figure 29: Axial profiles of droplet temperature for nonreacting (a) and reacting case (b). Symbols: experiments, solid line: case A dashed line: case 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 at and , 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.

Nomenclature

(Nomenclature entries should have the units identified)

: | Reynolds filter of variable A
: | Favre filter of variable A
: | Spalding mass number
: | Spalding thermal number
: | Heat capacity at constant pressure
D: | Molecular diffusion
: | Droplet diameter
: | Latent heat of the fuel
: | Mass of the particle
: | Total number of species
: | Weight of kth-species in the progress variable definition
: | Schmitt number
: | Turbulent Schmitt number
: | Sherwood number
: | Unstretched laminar flame speed
: | Subgrid scale turbulent flame speed
: | Droplet temperature
: | Temperature in the far field away from the droplet
t: | Time
: | Diffusion velocity of species k
: | Instantaneous velocity in the ith-coordinate direction
: | Lagrangian particle velocity vector
: | Subgrid scale velocity fluctuations
: | Cartesian coordinate in the i direction
: | Lagrangian particle position vector
: | Reaction progress variable
: | Mass fraction of species k
z: | Mixture fraction
: | Progress variable diffusion factor
: | Flame filter size
δ: | Laminar flame thickness
φ: | Thermochemical quantity
λ: | Thermal conductivity
: | Turbulent viscosity
: | Subgrid scale flame wrinkling
ρ: | Density
: | Mixture fraction evaporation source term
: | Progress variable source term
: | Progress variable unresolved convective fluxes

Subscripts

0:
Relative to fresh gases

Superscripts

TAB: | Variable stored in a chemical look-up table
∗: | From a 1D unstrained planar laminar premixed freely propagating flame.

Data Availability

The numerical data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This work was performed using HPC resources from GENCI-IDRIS (Grants 2016-x20162b0164 and 2017-x2017b0164). The authors acknowledge Antoine Verdier and Bruno Renou from CORIA laboratory for sharing the experimental data and for the fruitful discussions. They also are grateful of Francis Shum-Kivan, Eléonore Riber, and Bénédicte Cuénot for sharing the geometry and mesh of the computational domain and for the helpful discussions. A preliminary version of this manuscript has been presented at the AIAA Scitech 2019 Forum, San Diego, CA, USA.

References

P. Jenny, D. Roekaerts, and N. Beishuizen, “Modeling of turbulent dilute spray combustion,” Progress in Energy and Combustion Science, vol. 38, no. 6, pp. 846–887, 2012.

S. James, J. Zhu, and M. S. Anand, “Large-Eddy simulations as a design tool for gas turbine combustion systems,” AIAA Journal, vol. 44, no. 4, pp. 674–686, 2006.

P. Moin and S. V. Apte, “Large-Eddy simulation of realistic gas turbine combustors,” AIAA Journal, vol. 44, no. 4, pp. 698–708, 2006.

W. P. Jones, A. J. Marquis, and D. Noh, “An investigation of a turbulent spray flame using large eddy simulation with a stochastic breakup model,” Combustion and Flame, vol. 186, pp. 277–298, 2017.

W. P. Jones, A. J. Marquis, and K. Vogiatzaki, “Large-eddy simulation of spray combustion in a gas turbine combustor,” Combustion and Flame, vol. 161, no. 1, pp. 222–239, 2014.

F. Shum-Kivan, J. Marrero Santiago, A. Verdier et al., “Experimental and numerical analysis of a turbulent spray flame structure,” Proceedings of the Combustion Institute, vol. 36, no. 2, pp. 2567–2575, 2017.

B. Franzelli, A. Vié, M. Boileau, B. Fiorina, and N. Darabiha, “Large Eddy simulation of swirled spray flame using detailed and tabulated chemical descriptions,” Flow, Turbulence and Combustion, vol. 98, no. 2, pp. 633–661, 2017.

C. Heye, V. Raman, and A. R. Masri, “Influence of spray/combustion interactions on auto-ignition of methanol spray flames,” Proceedings of the Combustion Institute, vol. 35, no. 2, pp. 1639–1648, 2015.

S. B. Pope, “Small scales, many species and the manifold challenges of turbulent combustion,” Proceedings of the Combustion Institute, vol. 34, no. 1, pp. 1–31, 2013.

S. Menon and N. Patel, “Subgrid modeling for simulation of spray combustion in large-scale combustors,” AIAA Journal, vol. 44, no. 4, pp. 709–723, 2006.

B. Fiorina, D. Veynante, and S. Candel, “Modeling combustion chemistry in large eddy simulation of turbulent flames,” Flow, Turbulence and Combustion, vol. 94, no. 1, pp. 3–42, 2015.

J. A. V. Oijen, A. Donini, R. J. M. Bastiaans, J. H. M. T. Boonkkamp, and L. P. H. D. Goey, “State-of-the-art in premixed combustion modeling using flamelet generated manifolds,” Progress in Energy and Combustion Science, vol. 57, pp. 30–74, 2016.

B. Fiorina, R. Vicquelin, P. Auzillon, N. Darabiha, O. Gicquel, and D. Veynante, “A filtered tabulated chemistry model for LES of premixed combustion,” Combustion and Flame, vol. 157, no. 3, pp. 465–475, 2010.

P. Auzillon, O. Gicquel, N. Darabiha, D. Veynante, and B. Fiorina, “A filtered tabulated chemistry model for LES of stratified flames,” Combustion and Flame, vol. 159, no. 8, pp. 2704–2717, 2012.

R. Mercier, P. Auzillon, V. Moureau et al., “LES modeling of the impact of heat losses and differential diffusion on turbulent stratified flame propagation: application to the TU darmstadt stratified flame,” Flow, Turbulence and Combustion, vol. 93, no. 2, pp. 349–381, 2014.

R. Mercier, T. F. Guiberti, A. Chatelier et al., “Experimental and numerical investigation of the influence of thermal boundary conditions on premixed swirling flame stabilization,” Combustion and Flame, vol. 171, pp. 42–58, 2016.

R. S. Miller and J. W. Foster, “Survey of turbulent combustion models for large-Eddy simulations of propulsive flowfields,” AIAA Journal, vol. 54, no. 10, pp. 2930–2946, 2016.

J. C. Oefelein, R. W. Schefer, and R. S. Barlow, “Toward validation of large eddy simulation for turbulent combustion,” AIAA Journal, vol. 44, no. 3, pp. 418–433, 2006.

A. Verdier, J. Marrero Santiago, A. Vandel et al., “Experimental study of local flame structures and fuel droplet properties of a spray jet flame,” Proceedings of the Combustion Institute, vol. 36, no. 2, pp. 2595–2602, 2017.

B. Franzelli, B. Fiorina, and N. Darabiha, “A tabulated chemistry method for spray combustion,” Proceedings of the Combustion Institute, vol. 34, no. 1, pp. 1659–1666, 2013.

E. Ranzi, A. Frassoldati, A. Stagni, M. Pelucchi, A. Cuoci, and T. Faravelli, “Reduced kinetic schemes of complex reaction systems: fossil and biomass-derived transportation fuels,” International Journal of Chemical Kinetics, vol. 46, no. 9, pp. 512–542, 2014.

B. Franzelli, E. Riber, M. Sanjosé, and T. Poinsot, “A two-step chemical scheme for kerosene-air premixed flames,” Combustion and Flame, vol. 157, no. 7, pp. 1364–1373, 2010.

F. Shum-kivan, “Simulation des grandes echelles de flammes de spray et modélisation de la combustion non-prémélangée,” Toulouse University, Toulouse, France, 2017, Ph.D. thesis.

P. Pepiotdesjardins and H. Pitsch, “An efficient error-propagation-based reduction method for large chemical kinetic mechanisms,” Combustion and Flame, vol. 154, no. 1-2, pp. 67–81, 2008.

K. Kumar, J. E. Freeh, C. J. Sung, and Y. Huang, “Laminar flame speeds of preheated iso-octane/O₂/N₂ and n-heptane/O₂/N₂ mixtures,” Journal of Propulsion and Power, vol. 23, no. 2, pp. 428–436, 2007.

O. Gicquel, N. Darabiha, and D. Thévenin, “Liminar premixed hydrogen/air counterflow flame simulations using flame prolongation of ILDM with differential diffusion,” Proceedings of the Combustion Institute, vol. 28, no. 2, pp. 1901–1908, 2000.

N. Darabiha, “Transient behaviour of laminar counterflow hydrogen-air diffusion flames with complex chemistry,” Combustion Science and Technology, vol. 86, no. 1–6, pp. 163–181, 1992.

R. Vicquelin, B. Fiorina, S. Payet, N. Darabiha, and O. Gicquel, “Coupling tabulated chemistry with compressible CFD solvers,” Proceedings of the Combustion Institute, vol. 33, no. 1, pp. 1481–1488, 2011.

F. Charlette, C. Meneveau, and D. Veynante, “A power-law flame wrinkling model for LES of premixed turbulent combustion part I: non-dynamic formulation and initial tests,” Combustion and Flame, vol. 131, no. 1-2, pp. 159–180, 2002.

D. B. Spalding, “The combustion of liquid fuels,” Symposium (International) on Combustion, vol. 4, no. 1, pp. 847–864, 1953.

S. Saengkaew, V. Bodoc, G. Lavergne, and G. Grehan, “Application of global rainbow technique in sprays with a dependence of the refractive index on droplet size,” Optics Communications, vol. 286, pp. 295–303, 2013.

D. Paulhiac, “Modélisation de la combustion d’un spray dans un brûleur aéronautique,” Toulouse University, Toulouse, France, 2015, Ph.D. thesis.

J.-P. Legier, T. Poinsot, and D. Veynante, “Dynamically thickened flame LES model for premixed and non-premixed turbulent combustion,” in Proceedings of the Summer Program, pp. 157–168, Center for Turbulence Research, Stanford, 2000.

O. Colin, F. Ducros, D. Veynante, and T. Poinsot, “A thickened flame model for large eddy simulations of turbulent premixed combustion,” Physics of Fluids, vol. 12, no. 7, pp. 1843–1863, 2000.

S. Gallot Lavallee, D. Noh, W. Jones et al., “Experimental and numerical study of turbulent flame structures of a spray jet flame,” in Proceedings of the European Combustion Meeting, Combustion Institute, Croatia, Europe, 2017.

V. Moureau, P. Domingo, and L. Vervisch, “Design of a massively parallel CFD code for complex geometries,” Comptes Rendus Mécanique, vol. 339, no. 2-3, pp. 141–148, 2011.

F. Nicoud, H. B. Toda, O. Cabrit, S. Bose, and J. Lee, “Using singular values to build a subgrid-scale model for large eddy simulations,” Physics of Fluids, vol. 23, no. 8, Article ID 085106, 2011.

P. González-Tello, F. Camacho, J. M. Vicaria, and P. A. González, “A modified Nukiyama-Tanasawa distribution function and a Rosin-Rammler model for the particle-size-distribution analysis,” Powder Technology, vol. 186, no. 3, pp. 278–281, 2008.

L. Guedot, “Développement de méthodes numériques pour la caractérisation des grandes structures tourbillonnaires dans les brûleurs aéronautiques: application aux systèmes d’injection multi-points,” INSA, Rouen, France, 2015, Ph.D. thesis.

Shashank, E. Knudsen, and H. Pitsch, “Spray evaporation model sensitivities,” Annual Research Briefs, vol. 2011, Centre for Turbulence Research, Stanford, CA, USA, 1998.

R. S. Miller, K. Harstad, and J. Bellan, “Evaluation of equilibrium and non-equilibrium evaporation models for many-droplet gas-liquid flow simulations,” International Journal of Multiphase Flow, vol. 24, no. 6, pp. 1025–1055, 1998.

P. Sierra Sánchez, “Modeling the dispersion and evaporation of sprays in aeronautical combustion chambers,” University of Toulouse, Toulouse, France, 2012, Ph.D. thesis.