CFD Simulations AC2-11
Delft-Jet-in-Hot-Coflow (DJHC) burner
Application Challenge AC2-11 © copyright ERCOFTAC 2023
The DJHC flames have been studied in the literature using several different approaches. For a recent overview we refer to Perpignan et al. (2018). Here we present results from a comparative study of three different turbulence-chemistry approaches also discussed in Perpignan et al. (2018). These three approaches are: Conditional Source-term Estimation (CSE) model, a new variant of the Eddy Dissipation Concept (EDC) model and a new variant of the Flamelet Generated Manifolds (FGM) / presumed PDF model.
Before reviewing all details of the model setup, the three modelling approaches for turbulence chemistry interaction are described in some detail.
Turbulence-chemistry Interaction Approaches
Eddy Dissipation Concept
The Eddy Dissipation Concept (Magnussen, 1981) has been frequently employed in simulations related to the Flameless combustion regime due to its availability in several CFD codes, the fact that it contains no assumptions related to the combustion regime, and the possibility of directly employing a chemical reaction mechanism of choice.
In many simulations of Jet-in-Hot Coflow flames, the standard EDC model was shown to predict ignition far upstream, leading to a clear peak in the radial mean temperature profile not present in the experiments. This is due to the fact that the EDC by its definition in no way represents ignition by separate local events (ignition kernels) responsible for the ignition in many JHC flames. Nevertheless, several authors attempt to improve the model by modifying the model constants (Christo & Dally, 2005; Sarras et al., 2014; De & Dongre, 2015; De et al., 2011; Aminian et al., 2012; Frassoldati et al., 2010; Mardani, 2017). Indeed, in the EDC model, two model parameters appear. The key model features are the volume fraction of the reaction zones (γ) and the residence time in these zones (t), which are dependent on two constants, Ct and Cγ respectively. For example, Aminian et al., 2012 proposed adjusting the value of Ct, while Evans et al. (2015) investigated a range of values for the model constants and recommended values for CH4/H2 and C2H4/H2 flames. Both works were based on the Adelaide burner experiments and showed improvement in relation to the standard constants.
The proposed changes of model constants in the EDC model were not only purely empirical but were also motivated by the argument that the flameless combustion regime is different from the conditions for which the EDC was formulated originally. The modification of EDC models can potentially take into account the increased volumes caused by the different shapes of the reaction zones in the FC regime. This was in corroboration with earlier findings, also pointed in the conclusions of De & Dongre, 2015.
EDC model with local parameters (EDC-LP)
Using the EDC, Bao, (2017) presented a formulation with locally determined constants (EDC-LP), which was compared to the standard EDC formulation. The laminar and turbulent flame speeds are calculated in order to evaluate the constants, leading an expression for Cγ proportional to Da*3/4, where Da* is the Damk”hler number based on the Kolmogorov scale (i.e. the inverse of the more commonly used Karlovitz number). Simulations were performed with the standard k-ε model and with a ε-equation based Reynolds stress model (RSM) using ANSYS Fluent®. with linear pressure-strain model and Lien and Leschziner model for turbulent diffusion (RSM). Simulations were performed in ANSYS Fluent®.
According to the paper of Ertesvåg & Magnussen (2000) in which the energy cascade model is discussed, it holds that the supply of mechanical energy (w) and the energy dissipation rate at a given energy cascade level (q) of the n-th eddy are given by Eqs. 1 and 2, where ω is the strain rate, u is the velocity, ν is the viscosity, and CD1 and CD2 are model constants.
By assuming quasi-steadiness, from energy conservation:
Considering the smallest eddy (*):
Ertesvåg & Magnussen (2000) showed that Therefore:
Instead of imposing the model constant values, Bao, (2017) used the ratio of the two constant (obtained by combining Eqs. 5 and 6) to find expressions for the model constants that can be locally defined, a formulation here named EDC-LP (EDC with local paramenters).
Considering the velocity to be the turbulent flame speed, which can in turn be approximated by a relation containing the laminar flame speed and the turbulent Reynolds number (), the following relation is obtained:
Having the chemical time-scale and the Damköhler number of the Kolmogorov scale :
Using the turbulent flame speed again, and the definition of the Damköhler number of the Kolmogorov scale , where is the Kolmogorov time-scale, Eq. 5 can be rewritten as:
From this an expression for the constant can be written:
Using Eqs. 9 and 11, the original EDC constant values can be locally calculated:
Conditional Source-term Estimation
The Conditional Source-term Estimation (CSE) model (Labahn et al, 2015;Labahn et al, 2016) is not as widely employed as the EDC, possibly due to its complexity, the fact that it is relatively recent, and its absence in commercial codes.
The CSE shares similarities with the Conditional Moment Closure (CMC) model in its formulation. Both CSE and CMC use conditional averages to compute the chemical source term, while the conditional fluctuations are assumed to be negligible (conditioned to mixture fraction). In the CMC, transport equations for the conditional means of species and temperature are solved, while in the CSE these conditional means are obtained by inverting the integral equation that relates the unconditional to the conditional mean in a relation involving the (presumed) PDF. The model's absence of assumptions related to the combustion regime and its representation of turbulent fluctuations make it attractive for simulations involving the Flameless regime.
CSE as applied to the DJHC burner
Labahn et al. (2015, 2016) applied the CSE to model the DJHC flames using RANS and LES. Their approach employed two mixture fractions as conditioning variables: the first was the conventional mixture fraction and the second was a variable representing at the same time the oxygen profile and temperature (or enthalpy deficit profile) at the inflow boundary. Handling oxygen and temperature profile with only one variable induces an approximation but was found to be effective. LES calculations were performed using the standard Smagorinsky SGS model. RANS calculations were performed with both the standard and the realizable k-ε models and were found to give only slightly different results. The results herein shown are for the standard k-ε model. Chemistry was tabulated with an approach called Trajectory Generated Low Dimension Manifold (TGLDM), explained in (Labahn, 2015), in which chemical states and reaction rates are stored as a function of the two mixture fractions and mass fractions of CO2 and H2O. Simulations were performed using OpenFOAM (release not mentioned).
The CSE shares similarities with the CMC (Conditional Moment Closure) approach (Klimenko & Bilger, 1999), as both employ conditional averages to calculate the chemical source term. The underlying assumption is that the fluctuation of a given scalar can be associated to the fluctuations of a key quantity. Labahn et al. (2015) represent the conditional chemical source term of species k as the following function:
In this equation, η represents the value in sample space of the conditioning variable (e.g. mixture fraction). And for any stochastic variable A, 〈A|η〉 is the expected value conditional on the conditioning variable taking the value . The conditional average of the source term depends on the conditional averages of mass fraction (Yk), temperature (T) and density (ρ). In contrast with CMC, these quantities do not have transport equations in the CSE. Instead, transport equations only for the Favre-averaged species mass fractions are included. Therefore, the mean unconditional reaction rates are determined with:
In which the spatial coordinate(xj), time (t), and Favre-averaged PDF () are present. Labahn et al. (2015) assumed a β-PDF distribution. In order to calculate , the integral that equals the Favre-averaged mass fraction (Eq. 16) is inverted to find .
The averages were conditioned to values of mixture fraction (Z) and of a variable able to capture the influence of the air stream in the DJHC problem. This variable (W) is defined by Eq. 17 and is equal to 0 at the location with the minimum amount of oxygen in the coflow, and is equal to 1 where the oxygen concentration is maximum (air).
The integral to be inverted was then:
Having the conditional mean after the inversion, is determined by consulting the chemical look-up table.
Flamelet Generated Manifolds
First presented by van Oijen & de Goey (2000), the Flamelet Generated Manifolds (FGM) approach is based on tabulated-chemistry, in which laminar flame structures are considered and tabulated in a pre-computation step, in which detailed chemistry is adopted. The important thermochemical properties and the source terms of a selected progress variable are tabulated as functions of independent variables (mixture fraction and/or progress variable, enthalpy), which are solved for in the CFD calculation. For turbulent flows usually an assumed PDF approached is used for the independent variables.
The relatively low computational costs make the FGM approach attractive and many works have been performed with the approach. Earlier works applying an FGM approach to the DJHC are (Sarras et al., 2014), (Abtahizadeh et al., 2015) and (Abtahizadeh et al., 2017). The approach presented here is a new variant presented in (Huang et al., 2017) and (Huang, 2018)
FGM with diluted air (DAFGM)
The application of FGM herein shown was developed to represent the dilution of reactants with burnt gases. Developed by Huang et al. (2017, 2018), the DAFGM (Diluted Air FGM) are generated having four control variables: a mixture fraction, a progress variable, enthalpy loss, and air dilution level. The progress variable was defined as the linear combination of the CO2, H2O, CO and H2 mass fractions. Therefore, the counterflow diffusion flamelets were calculated and tabulated based on this four-dimension set of variables. The RANS simulations were performed with the standard k-ε model and the LES used the standard Smagorinsky model. Simulations were performed using OpenFOAM, release 2.3.1, extended with a new storage and retrieval routine.
The pre-calculated flamelets were parametrized by four variables: a mixture fraction (ξ), the progress variable mentioned above (C), enthalpy loss (η) and dilution level (γ). The mixture fraction is defined by Eq. 19, where Z is the conventional mixture fraction based on fuel and air and ? represents the mass fraction of diluent (combustion products of stoichiometric combustion) in the mixture (Eq. 20).
In Eq. 20, Yd is the dilution variable, defined by the sum of products mass fractions (CO2 and H2O). The superscript Dil denotes the dilution variable value in combustion products of stoichiometric combustion, i.e. at the stoichiometric mixture fraction (Zst). The dilution level of the air in the counterflow flames underlying the FGM then is given by:
The enthalpy loss is defined by Eqs. 22, in which is the enthalpy, is the adiabatic enthalpy defined by Eq. 23, while and are the enthalpies of the diluent for maximum and minimum enthalpy loss, respectively.
In the calculation of the DJHC flame radiative heat transfer can be neglected but the temperature variation at the coflow inlet leads to an enthalpy variation that should be taken into account. Accordingly, was taken zero and from the experimental data on the temperature, was taken as 2x106 J/kg.
The four variables employed to parametrize the manifolds have values at every location in space from the transport equations of mixture fraction, progress variable, enthalpy and dilution variable that are solved throughout the domain. In the framework of RANS and LES transport equations are solved respectively for mean and for resolved variables. And for mixture fraction also the variance is calculated in order to determine the assumed β-PDF. Therefore a five-dimensional FGM lookup-table is used. For further details we refer to (Huang, 2018).
Calculation Domain and Grids
The characteristics of the different modelling approaches for which calculations are compared are listed in Table 3. The selection of models and parameter values is explained following the rows in the table, starting with domain and grid size.
The computational domain extends from the fuel nozzle tip in axial direction beyond the last measurement station and in radial direction to a distance much wider than the width of the coflow. The latter is necessary because the entrainment of laboratory air into the coflow stream was shown to be crucial for accurate predictions and has to be included in the simulation.
|Domain (mm)||2D||2D||2D||3D cylindrical||3D hexahedral|
|Axial-radial||225 x 80||225 x 80||225 x 80||225 x 80||225 x 80|
|Grid size (#)||190 x 145
|31250||300 x 215
|Closure model||k-ε or RSM||Standard k-ε||Standard k-ε||Smagorinsky||Smagorinsky|
|Kinetic scheme||GRI 1.2||GRI 2.11||GRI 3.0||GRI 2.11||GRI 3.0|
|Multistream CSE (two mixture fractions, β-pdf)||Assumed PDF (, β-pdf for mixture fraction and for scaled progress variable)||Multistream CSE (two mixture fractions, β-pdf))||Assumed PDF (, β-pdf for mixture fraction and for scaled progress variable)|
|Thermochemical Scalar equations||Mean of species
Mean of enthalpy
|Mean and variance of two mixture fractions||Resolved mixture fraction and resolved progress variable
|Resolved mixture fractions and SGS variances
|Resolved mixture fraction and resolved progress variable
|Radiation||Not included||Not included||Not included||Optically thin.
TRI not included.
|Figure 9: Two-dimensional grid employed by Bao (2017). Axial direction from left to right. Radial direction from bottom to top.|
All simulations considered three inlets: fuel inlet through the jet nozzle, coflow inlet, and air inlet. An example of a two-dimensional computation domain and grid used in RANS is given in Figure 9. The grid is refined close to the nozzle, near the centre and in the regions where a shear layer is expected (as indicated by the dashed rectangles). A grid independence study for RANS-EDC simulations has been reported in (De et al., 2011) and concluded that a grid of 180x125=22500 cells is sufficiently fine. Details of the grids used in this study are given in Table 3.
Software and solvers
Steady RANS simulations are performed using ANSYS-Fluent, version 14.5. The EDC-LP model is implemented using user-defined functions. (Bao, 2017)
OpenFOAM is used for both RANS and LES. (Labahn, 2015, 2016)
OpenFOAM (version 2.3.1) is used for both RANS and LES. The momentum equations are solved using a pressure-based solver. The transport equations for the scalar variables are solved In every iteration after the velocity prediction, and the controlling parameters needed for table lookup are derived and properties are retrieved from the table. (Huang, 2017, 2018)
The turbulence chemistry interaction approaches have been used with different models for the statistical description of turbulence. The three considered approaches have been used with RANS (Reynolds-Averaged Navier-Stokes) using the standard k-ε model (Bao, 2017) (Labahn, 2015, 2016) (Huang, 2017, 2018).
The simulations using EDC were also done using a standard Reynolds stress model (RSM) more precisely an ε-equation based RSM with linear pressure-strain model and Lien and Leschziner model for turbulent diffusion (ANSYS-Fluent, 2013) (Bao, 2017).
CSE and DAFGM have also been combined with LES (Labahn, 2015, 2016) (Huang, 2018). The standard Smagorinsky model was used for closure of the subgrid stress.
The boundary conditions at the inflow of fuel and coflow are set using information available for measurements at z=3mm, but differences between the different approaches are mentioned below.
Neumann boundary conditions are imposed at the outflow. At the side boundary on the air side a slip wall boundary condition was assumed. In the 2D RANS simulations, symmetry is imposed at the axis of symmetry. LES simulations were performed in 3D and then the symmetry axis is not a boundary.
The experimental data measured at the axial position of 3 mm are employed to determine the jet and coflow inlet boundary conditions. Equilibrium is assumed in order to estimate the species distribution at the boundary based on the available temperature and O2 profiles. Bao (2017) imposed the turbulence kinetic energy based on the available profiles of the Reynolds stresses, with the relation shown in Eq. 24, while the dissipation rate followed the relation in Eq. 25. In these equations, are the axial, radial, azimuthal, and shear Reynolds stresses, respectively. The mean axial velocity is represented by .
The outer boundary is defined as a slip wall, having zero shear stress. A constant profile for mean axial velocity equal to 0.3 m/s is imposed at the air inlet.
The experimental data measured at the axial position of 3 mm was employed to determine the jet and coflow inlet boundary conditions. Equilibrium was assumed in order to determine the thermochemical variables at the boundary based on the information on mean temperature and O2 profiles. The oxygen mass fraction profile is assumed to have an analytical relation with the mean temperature profile. See below the thermochemical scalar equations for more details. The resulting oxygen mass fraction profile is shown in Figure 2.
The fuel composition was assumed to be 85% CH4 and 15% N2.
For the LES simulations, turbulence at the inlet is generated using the synthetic turbulence method of Kornev & Hassel (2007) based on prescribed mean velocity profiles, the corresponding Reynolds stress tensor and integral length scales. Turbulent fluctuations for temperature, mixture fraction, oxidiser split or species mass fractions are assumed zero at the inlet.
The outer boundary is defined as a slip wall. A constant profile for mean axial velocity equal to 0.5 m/s is imposed at the air inlet. The choice of this value is further discussed in the Best Practice Advice section.
For the RANS simulations the inlet boundary conditions used the available experimental data at 3 mm, except for the oxygen concentration. The turbulent kinetic energy and the dissipation rate were obtained in the same manner as in the EDC simulations explained before.
Equilibrium was assumed in order to estimate the scalar variables at the boundary based on the available temperature and a modified O2 profiles. In order to exactly satisfy the constraint that the fuel/air ratio agrees with the fuel and air flow rates to the coflow burner measured by mass flow controllers, the mixture fraction boundary condition was modified and, as can be seen in Figure 10, the resulting assumed O2 – profile has a lower value in the inner region r < 25mm than measured with the sampling probe.
For the LES simulations, turbulence at the inlet is generated using the synthetic turbulence method of Kornev & Hassel (2007) based on prescribed mean velocity profiles, the corresponding Reynolds stress tensor and integral length scales. No turbulent fluctuations for enthalpy and dilution variable were included at the inlet.
|Figure 10: O2 mass % in experiment (Oldenhof et al, 2011) and in simulations (CSE and DAFGM)|
Kinetic scheme and chemistry reduction
All simulations are based on versions of the GRI-mech detailed mechanism, but they significantly differ in the chemistry reduction method.
In EDC transport equations are solved for mean values of all species. The use of the DRM19 mechanism, which is derived from GRI-mech 1.2 allows significant reduction in computational cost.
Employing GRI-mech 2.11, a lower dimensional manifold for chemical evolution is derived using the trajectory generated low dimensional manifold (TGLDM) method. GRI-mech 2.11 is a mechanism for 49 species, but a lower dimensional manifold in composition space is derived from a set of trajectories in this space. Mixing variables and a limited number of reacting coordinates determine the dimension of the manifold (Labahn, 2015). See below under scalar variables for more information.
Employing GRI-mech 3.1, a lower dimensional manifold was derived using a new variant of the Flamelet Generated Manifold (FGM) method. GRI-mech is a mechanism for 53 species, but a lower dimensional manifold for species evolotion is derived from steady laminar counterflow flames between fuel and diluted and cooled air for strain rates from low value to extinction strain rate and also including an extinguishing flamelet. The flamelets are calculated assuming unity Lewis number for all species (Huang, 2017, 2018) The method has been explained in more detail above in the section on turbulence-chemistry interaction.
Thermochemical scalar equations
In EDC transport equations are solved for mean values of species mass fractions and enthalpy. Since the DRM19 has 19 species in total 20 equations are solved.
To describe the local conditions two non-reacting scalars and two reacting scalars (H2O and CO2 mass fractions) are used.
The first non-reacting scalar is mixture fraction Z based on fuel and air. The second non-reacting scalar W represents the enthalpy variation at the coflow inlet. It is defined as a normalised temperature: for r > 2.25mm for r < 2.25mm. The value of temperature for radial distance larger than the last experimental point, i.e. from 35.1 mm to 44.6 mm, was obtained by linear extrapolation of the known profile. The mass fraction of oxygen at the inlet boundary is assumed to vary linearly with W between minimal value in the middle of the coflow temperature maximum (W=0) and maximal value in the air (W=0). Mean (or resolved) values of these variables are solved (4 equations) and equations for the variance of the non-reacting scalars. In total 6 equations are solved. In addition to that, unconditional and conditional averages of all species mass fractions are calculated through the CSE algorithm (Labahn, 2015). In the LES, the mean enthalpy equation was also solved and the radiation heat loss was taken into account via an optically thin model (emission only).
To describe the local conditions, four scalar quantities are needed: mixture fraction, progress variable, dilution variable and enthalpy. In turbulent flow equations the mean (or resolved) values of these quantities are solved (4 equations) and equations for the variance of mixture fraction. In total 5 equations are solved. Fluctuations in progress variable are not taken into account, but it was checked that including them has negligible influence on the results. (Huang, 2018)
- Abtahizadeh E, de Goey P, van Oijen J. Development of a novel flamelet-based model to include preferential diffusion effects in autoignition of CH4/H2 flames. Combustion and Flame 2015; 162:4358-4369.
- Abtahizadeh E, de Goey P, van Oijen J. LES of Delft Jet-in-Hot Coflow burner to investigate the effect of preferential diffusion on autoignition of CH4/H2 flames. Fuel 2017;191:36-45.
- Aminian J, Galletti C, Shahhosseini S, Tognotti L. Numerical investigation of a MILD combustion burner: analysis of mixing field, chemical kinetics and turbulence-chemistry interaction. Flow, Turbulence and Combustion 2012; 88:597-623.
- Ansys Fluent, Ansys fluent theory guide 15.0, Inc, Canonsburg, PA (2013).
- Bao H. Development and validation of a new Eddy Dissipation Concept (EDC) model for MILD combustion, MSc Thesis, Delft University of Technology, 2017.
- Christo FC, Dally BB. Modeling turbulent reacting jets issuing into a hot and diluted coflow. Combustion and Flame 2005;142:117-129.
- De A, Dongre A. Assessment of turbulence-chemistry interaction models in MILD combustion regime. Flow, Turbulence and Combustion 2015;94:439-478.
- De A, Oldenhof E, Sathiah P, Roekaerts DJEM. Numerical simulation of Delft-Jet-in-Hot-Coflow (DJHC) flames using the Eddy Dissipation Concept model for turbulence-chemistry interaction. Flow, Turbulence and Combustion 2011;87:537-567.
- Ertesvåg IS, Magnussen BF. The eddy dissipation turbulence energy cascade model. Combustion science and technology 2000;159:213-235.
- Evans MJ, Medwell PR, Tian ZF. Modeling lifted jet flames in a heated coflow using an optimized eddy dissipation concept model. Combustion Science and Technology 2015;187:1093-1109.
- Frassoldati A, Sharma P, Cuoci A, Faravelli T, Ranzi E. Kinetic and fluid dynamics modeling of methane/hydrogen jet flames in diluted coflow. Applied Thermal Engineering 2010;30:376-383.
- Huang X, Tummers MJ, Roekaerts DJEM. Experimental and numerical study of MILD combustion in a lab-scale furnace. Energy Procedia 2017;120:395-402.
- Huang X. PhD Thesis, Delft University of Technology, 2018
- Klimenko AY, Bilger RW. Conditional moment closure for turbulent combustion. Progress in Energy and Combustion Science 1999;25:595-687.
- Kornev N, Hassel E. Method of random spots for generation of synthetic inhomogeneous turbulent fields with prescribed autocorrelation functions. Communications in Numerical Methods in Engineering 2007;23:35-43.
- Labahn JW, Devaud CB. Large Eddy Simulations (LES) including Conditional Source-term Estimation (CSE) applied to two Delft-Jet-in-Hot-Coflow (DJHC) flames. Combustion and Flame 2016;164:68-84.
- Labahn JW, Dovizio D, Devaud CB. Numerical simulation of the Delft-Jet-in-Hot-Coflow (DJHC) flame using conditional source-term estimation. Proceedings of the Combustion Institute 2015;35:3547-3555.
- Mardani A. Optimization of the Eddy Dissipation Concept (EDC) model for turbulence chemistry interactions under hot diluted combustion of CH4/H2. Fuel 2017;191:114-129.
- Perpignan A.A.V., Gangoli Rao A., Roekaerts, D.J.E.M., Flameless Combustion and its Potential Towards Gas Turbines, Progress in Energy and Combustion Science 2018; 69: 28-62
- Sarras G, Mahmoudi Y, Arteaga Mendez LD, van Veen EH, Tummers MJ, Roekaerts, DJEM. Modeling of turbulent natural gas and biogas flames of the Delft Jet-in-Hot-Coflow burner: effects of coflow temperature, fuel temperature and fuel composition on the flame lift-off height. Flow, Turbulence and Combustion 2014;93:607-635.
- van Oijen JA, de Goey PH. Modelling of premixed laminar flames using flamelet-generated manifolds. Combustion Science and Technology 2000;161:113-137.
Contributed by: André Perpignan, Dirk Roekaerts, E. Oldenhof, E.H. van Veen, M.J. Tummers, Hesheng Bao, Xu Huang — TU Delft
Contributed by: Jeffrey W. Labahn — Stanford University
© copyright ERCOFTAC 2018