Transported-PDF (IEM, EMST) micromixing models in a hydrogen-air nonpremixed turbulent flame

The topic of this study is the numerical simulation of a turbulent non-premixed hydrogen flame with different micromixing models in order to investigate their predictive capability. The two micromixing models are compared. The first model is the interaction by exchange with the mean (IEM) model (Dopazo and O’Brien Acta Astronaut 1(9–10):1239, 1974). The second one is the Euclidean minimum spanning tree (EMST) model (Subramaniam and Pope Combust Flame 115(4):487, 1998). The dynamic model for the mixing time-scale is handled, by computing the individual time-scales for the reactive scalars in each cell during the simulation course using the Fluent/MM-INTAS CFD code. The predictions are validated against experimental data provided by Raman and Laser Doppler anemometry measurements for a turbulent jet Hydrogen-air diffusion flame. The interaction of turbulence-chemistry is handled with TPDF method. The chemical model used here consists of 11 chemical species and 23 reactions. Comparisons with experimental data demonstrate that predictions based on the EMST model are slightly better. The EMST improves largely the precision of the results to the detriment of the RAM and the CPU performances. Overall, profile predictions of mixture fraction, flame temperature and major species are in reasonable agreement with experimental data.


Introduction
Probability density function (PDF) transport equation methods integrated within a conventional CFD flow solver represents a rational approach for the study of turbulent combustion. The major attraction of this method is that the terms associated with chemical reaction appear in closed form, leaving only 'molecular' mixing and turbulent transport terms to be modelled. There has been considerable development in the past three decades in PDF methods, and reviews can be found in Dopazo and O'Brien [1] and Subramaniam and Pope [2].
Hydrogen combustion attached much attention recently because of the need for a clean alternative energy. Indeed, unlike hydrocarbon fuels, the combustion of hydrogen does not produce harmful emissions such as CO 2 , CO, Soot and unburned hydrocarbons. The only serious pollutant is NO x , which can be minimized by reducing the combustion temperature by carefully designed lean premixed combustors. However, this technique can be dangerous because of the very flammable character of hydrogen. In order to avoid this kind of risks, it is recommended to operate in diffusion flames.
The modelling of turbulence and chemistry is very important for the numerical simulation of combustion devices. The knowledge of the fields of the major chemical species, pollutants and temperature is essential for the design of practical combustion devices such as industrial furnaces, gas turbines or internal combustion engines. This study contributes to the modelling of turbulence in reacting flows. The objective of this work is to compare two different micromixing models in order to investigate their predictive capability. The hydrogen jet flame is an attractive candidate for a detailed study because it has well defined boundary conditions, measured data of velocity, temperature, major and minor species, high Reynolds number and a visible flame length smaller than one meter.
The studies on a turbulent non-premixed hydrogen flame are concentrated on combustion modelling and turbulence modelling in hydrogen mixtures (H 2 -N 2 , H 2 -He) and pure hydrogen cases. The most used models of turbulence in the calculation of turbulent hydrogen non-premixed flames are the k-ε and the Reynolds stress (RSM) models. It is known that the RSM model is clearly superior for situations in which the anisotropy of turbulence has a dominant effect on the mean flow, and it has greater potential to give accurate predictions for complex flows. This model is used in the present work. The H 2 flame investigated in this work is studied experimentally by Barlow [3]. These measured data are utilized in this study for comparison. Barlow et al. [4] simulated numerically NO x formation in this flame under helium dilution (0, 20 and 40 %) in order to draw comparisons between calculations made with the CMC and PDF combustion models with the (k-ε) model. Fairweather and Woolley [5] examined the same flame using a CMC approach with three kinetics schemes (5, 24 and 62 reaction steps). The predictions were based on both k-ε and Reynolds stress/scalar flux turbulence closures. They showed that simulations based on Reynolds stress closures are superior to those obtained using the simpler eddy viscosity-based approach. The predictions are optimized against non-diluted measures results of Barlow and Carter [6] in the radial profiles rather than in the axial profile. In this case, only slight overpredictions of radial mixture fraction profiles in regions away from the nozzle were observed. Moreover, all results are found to be insensitive to the kinetic scheme employed. Obieglo et al. [7] applies the (k-ε) model of turbulence with pope correction to this flame. Three combustion models are compared. These are: a probabilistic Euler Lagrangian (PEUL) model, a PDF model and eddy dissipation (EDC) model with single one step reaction. It was shown that the precision of the results depends on the choice of the combustion model. Both probabilistic methods give better predictions than the standard model. In the results obtained, the calculation data for a physical space were mostly in agreement with the experiments except for temperature discrepancies that occurred in mixture fraction space. In the simulations made at EKT (Darmstadt), a comparison between the (k-ε) model and RSM is made in this H 2 flame under chemical equilibrium assumption. The corrected constants of turbulence models are calibrated to reproduce the spreading rate of the round jet of Panchapakesan and Lumley [8]. The constants values adopted for the turbulence models are respectively: k-ε model: C ε2 = 1.83, rev. RSM of Jones and Musonge [9]: C ε2 = 1.78. For the axial and radial velocity and turbulence fields, the agreement in the far field is satisfactory, but in the near field the deviations are large. However, the Jones and Musonge model showed to perform slightly better. The spreading of major concentrations in radial direction and temperature is smaller than predicted. The axial variations of species concentrations and temperature show similar behaviour.
In this study, the numerical simulations are performed for an axisymmetric turbulent jet diffusion hydrogen/air flame under atmospheric pressure, by using a detailed chemistry in a combustion model which is, here, TPDF approach. It is observed that few studies have been done with concern to the application of TPDF approach to the H 2 flame and especially the above H 2 flame. Schlatter et al. [10] studied this flame in the case of 20 % He dilution with flamelet combustion modelling including two calculations of flamelets obtained for two different strain rates: high strain rate for the section close to the nozzle x/l = 1/8 where l is a flame measured visible length and low strain rate at the flame tip and PEUL combustion model. The turbulence model applied was (k-ε) with Pope correction. This contribution was concentred on the applicability of the two combustion model to NO x modelling rather than physical parameters of the flame. It was shown that both combustion models are equally capable of predicting the measured temperatures with reasonable accuracy. The flamelet results of the temperature at position x/l = 1/8 are very close to experiment because preferential diffusion is including the flamelet model whereas the PEUL model does not. More recently, a hybrid method (RANS-TPDF) has been used by Masri [11] for the study of the H 2 /N 2 lifted flame structure. He showed that the flame is controlled by chemical kinetics. We can also quote the experimental study of Theron [12] on the effects of heat release on mixing processes and flow structure in a high-speed subsonic turbulent H 2 jet.
The simulations are performed with the sophisticated commercially available code Fluent which solves the Navier-Stokes equations with the finite volume method.

Modelling the flow field
The numerical simulation of the turbulent flow field includes the solution of the overall continuity equation, Navier-Stokes equation energy and scalar equation. These equations are written in their Favre-averaged form. The turbulent fluctuations are taken into account with RSM. The TPDF concept necessitates the solution of transport equations of PDF, withρ the mass density,ũ i the velocity vector components,p the pressure, −ρ u i u j the Reynolds stress (unclosed term),φ the scalar (species mass fractions,Ỹ α or enthalpy, h), J φ i the molecular scalar flux, −ρ u i φ the turbulent scalar flux (unclosed term), andS φ is the source term (unclosed term). The closure of this term depends on which model of turbulence we choose.
In addition to the conservation Eqs. (1)-(3), an expression for the viscous stress tensor and equations describing the thermodynamic state are needed. For a Newtonian fluid the viscous stress tensor contains only simple shear effects and can be written as: The hydrodynamic equation of state is the gas law for multi-component mixtures: Yα Wα (5) with μ eff the effective viscosity (molecular and turbulent viscosity), R the perfect gas constant, and Wα is the atomic weight of species α. With large turbulence Reynolds numbers, the molecular effects are negligible in front of the effects due to turbulent agitation.

Modelling the turbulence
First, the Reynolds tensor is closed using RSM. Diffusion species flux is given by Fick's law. Turbulent scalar flux terms are closed using a transport gradient assumption. The RSM involves calculation of the individual Reynolds stress using differential transport equations. Unfortunately, several of the terms in the exact equation are unknown, and modelling assumptions are required in order to close the equations. The transport equations for the Reynolds stress are: Here, the buoyancy production and production by system rotation terms are equal to zero. The first term on the right hand side is the production term due to mean strain, which needs no modeling: The diffusion term in Eq. (6) is modelled as [13]: where k is the turbulent kinetic energy: The pressure-strain correlation is:φ The Launder-Reece-Rodi (LRR) model [14] is used for this term. The return to isotropy term is modelled as: where a ij is the anisotropy tensor defined by: The mean strain, also called the rapid term, is modelled as: where the productionP k of turbulent kinetic energy is: According to Fu et al. [15], Launder [16] the contribution of the term (C ij − 2 3 δ ijCk ) is negligible, unless in the case of swirling flows.
The transport equation of the dynamic dissipation is: The termψ(ε) contains the effects of production and destruction of dynamic dissipation. This term is modeled by using the form proposed by Launder et al. [6]: The modeling of the diffusion termD ε is given by the following relationship: The turbulent scalar flux is modelled as:ρ in which μ t is the turbulent or eddy viscosity and σ φ is the effective Prandtl or Schmidt number for the scalar φ.
The turbulent viscosity is given by: From the foregoing we can deduce the parabolized set of equations in cylindrical coordinates where the generalized equation is: The sources terms S Φ are given in Table 1, whereP = −ρ uv ∂ũ ∂r . The model constants used in the present study are given in Table 2.
For all the Reynolds stresses, C diff = C s , and for dynamic dissipation rate C diff = C ε Table 2 The turbulent model constants

Combustion modelling
Computational fluid dynamics (CFD) is based on the Reynolds average approach, and its application for modeling turbulent reacting flows requires the addition of a micro-mixing model. In fact, when the Reynoldsaverage is applied to the transport equation of a reacting scalar, the chemical source term appears in an unclosed form and needs a special treatment. One of the most popular is the joint probability density function (JPDF) approach, which has been successfully applied in modeling of the turbulence-chemistry interaction.
In the literature many different joint PDF models can be found, for example models for the joint PDF of composition, for the joint PDF of velocity and composition or for the joint PDF of velocity, composition and turbulent frequency.
In this work, only a joint PDF of composition vector is considered. This is a one-time, one point JPDF which has the main advantage to treat chemical reactions exactly without any modelling assumptions [17]. However, the effect of molecular mixing has to be modelled.
In composition PDF methods physical scalars, including temperature and species concentrations, are treated as independent random variables. The JPDF is then a function of spatial location, time and composition space. Once the joint PDF is obtained at a certain position x and time instant t, the mean value for any function, Q, of these scalars can be evaluated exactly as: where ϕ is the vector of physical scalars, x is the corresponding random variable vector, Q is a function of ϕ only and f is the joint PDF, which represents the probability density of a compound event ϕ = ψ. For variable density flows, it is useful to consider the joint mass density function (MDF)F φ (ψ) = ρ(ψ) f φ (ψ). Density weighted averages (Favre averages) can be considered:

Scalar PDF transport equation
When the joint scalar MDF Fϕψ is considered, the following transport equation is modelled and solved [17]: where i and α are summation indices in physical space and composition space, respectively; <A/B> is the conditional mean of the event A, given that the event B occurs. The terms in the LHS (appear in closed form) describe, respectively, evolution of probability in time, convection of probability in the physical space with the mean velocity, and transport in composition space due to chemical reaction. Terms on the RHS are unclosed. The first term represents the turbulent transport of probability in physical space (turbulent scalar flux) and is commonly modeled using the gradient diffusion hypothesis, the second term describes the transport of probability in composition space due to molecular fluxes and is further referred to as micro-mixing term. Two models, IEM (interaction by exchange with the mean) [18] and EMST (Euclidean minimum spanning tree) [19] are tested and used to close this term.

Hybrid solution method
The numerical model constructed for this study is based on the geometry and dimensions of the experimental flame configuration. The flow and mixing fields are computed by the solution of the 2D, axisymmetric forms of the density-weighted fluid flow equations, supplemented with a Reynolds stress closure. The RSM involves calculation of the individual Reynolds stress terms using differential transport equations.
In the model adopted in this study, turbulence production by buoyancy and system rotation are taken into account. The Pope correction for round jets and the buoyancy contributions are added to the turbulence dissipation rate for RSM turbulence models. Solution of the transport equations is achieved using the Fluent CFD code. For the turbulent Schmidt number (Sc t ) in the mixture fraction transport equation, a value of 0.85 is taken. A transport equation for the mean enthalpy is solved with the radiation source term which is modeled using the assumption of optically thin transfer between the hot combustion gases and the cold surroundings [6]. The radiating species considered is H 2 O. Some modifications have been made to introduce Pope correction and flame radiation.
Equation (23) is solved using the consistent hybrid Finite-Volume/Monte Carlo method presented in [20]. The finite volume method solves for the mean velocityŨ , the turbulence kinetic energy k, the turbulent dissipation rate ε and the mean pressure < P >. These values and the turbulent time scale t = k/ε are then passed to the Monte Carlo method. In the Monte Carlo method, the particles are moved through the domain by a spatially second order accurate Lagrangian method. The particles evolve due to different processes such as convection, diffusion, mixing and reaction. For the numerical solution these processes are applied in different, so-called fractional steps. In a first convection step, particles are advanced to a new position, where Xi is the position vector of the ith particle and t is the particle time step. For unsteady flows, the particle time step is the physical time step. For steady-state flows, local time steps are calculated for each cell as: The next fractional step is the mixing step. Two models are applied, the IEM and EMST model: IEM model: It is also clear that the time-scale ratio C φ = ϕ 2 ε ϕ ε k can be a function of space and time, and that it generally is a fluctuating quantity. In the past, it has been shown that simulation results are often sensitive to this model constant C φ , and different values have been used. Therefore, we will here use a transport model to evaluate directly C φ , thereby removing the need to specify the value of the mixing constant. The model for the mixing time-scale is obtained from the transport models for the sub-mesh varianceφ 2 and the scalar dissipation rate ε ϕ . This approach was detailed in our former work [21,22].The commonly used models for the conserved scalar variance and dissipation rate are employed. The scalar mixing time-scale can be computed as where the mixing time-scale for all scalars is approximated by that of the mixture-fraction [21].
It is noted that using this formulation, the individual time-scales for the reactive scalars may also be computed. The scalar-dissipation rate model is based on the assumption that production equals dissipation in the scalar variance equation. For reactive scalars, this relation could be modified to include the co-variance of scalar and scalar production rate appearing in the reactive scalar variance equation. This term is closed in the context of a transported PDF model. In the numerical implementation, the sub-mesh variance is computed directly from the particle composition vector. Hence, the dynamic model for the variance is not required in the PDF context. In the present work, it was found, that C φ varies slightly around the value 3. The computer program MM-INTAS [23] was provided as a tool to implement the scalar-dissipation rate model. Combined with the FLUENT/CFD code, this coupled algorithm was sufficient to carry out the computations discussed in this work.
EMST model: where the vth edge of the tree connects the particle pair (m v , n v ), and δ represents the Kronecker delta. The determination of B v and γ is discussed in Ref. [19]. After the mixing step, the reaction fractional step is applied in which the chemical source term is integrated over the local time step t: In practice, the integral must be evaluated using a stiff ODE solver. Because transported PDF simulations are typically used for reacting flows with complex chemistry, the chemical reaction step will often dominate the overall computational costs. The SEMI/ISAT (in situ adaptive tabulation) method is adopted to calculate the integral in (28). Integration of the chemical source term is a significant hurdle in RANS/PDF simulations. Recent advances in computational methods have somewhat alleviated this problem [24,25].
In this work, we resort to direct integration using a stiff ODE solver, but employ certain simplification procedures. First, we limit direct integration only to particles with mixture-fraction in the range of 0.03-0.95. Composition vectors outside this range are advanced using a flamelet look-up table. It was found that using wider ranges did not affect the final results.
In addition, we tabulate the direct integration-based results in a single time-step and re-use the computations for all particles with compositions that are close to these tabulated points in composition space. This follows the ISAT scheme [24] in the basic formulation.
Finally, the second convection step is calculated as: The increment of the Wiener process dWi (t) is Gaussian distributed: where ζ i are the Gaussian random variables with zero mean and unity variance.

Computational domain and boundary conditions
For the simple jet flame calculations, we use a cylindrical coordinate system with the origin at the center of the fuel jet, Fig. 1. In the computations, velocity and length are normalized respectively by the centerline velocity at the inlet U j and the jet radius R j (for radial evolution) and visible flame length L [26] (for axial evolution). The central jet velocity is 296 m/s, and the Reynolds number is Re = 10,500. The co-flow velocity is fixed at 1 m/s. The boundary condition of the inlet velocity profile is given by the so-called power law: [27]. The turbulence properties follow from k = 3 2 I 2 u 2 and ε = where the turbulence intensity is I ≈ 0.05 and the mixing length is given by L m = 0.07 R j [25,26]. The turbulence is assumed to be isotropic such that: uv = 0 and uu = vv = ww = 2 3 k. The simulations were performed with 100 particles per computational cell.
An orthogonal mesh was initially generated. The mesh is dynamically refined during numerical iterations using user-specified gradient and curvature boundaries. From the nozzle exit, the mesh characteristics are respectively 200 nodes in the axial direction and 165 nodes in the radial direction. The numerical accuracy is checked by comparing the predicted results calculated using the grid mentioned above with that obtained using a coarser grid with 160 nodes in axial direction and 120 nodes in radial direction. It is found that the two sets of results are very close to each other and therefore may be regarded as grid independent. The mass flow rate, total temperature, turbulence intensity and hydraulic diameter are specified for the inlet boundary. At the outlet region, outflow condition is assumed and the symmetry condition on the side boundary.
The measured velocity profile at the inlet stream is not used because with the present mesh resolution it was not possible to achieve the same profile. This is mainly due to Fluent's interpolation scheme.

Results
The mathematical model is discretized using the finite volume technique on a cylindrical staggered grid. Central differences are employed for the evaluation of the diffusion terms, while a first order upwind scheme is used for the evaluation of the convective one. A time-marching SIMPLE algorithm is employed to couple velocity-pressure fields. Discretized equations are solved in a segregated manner using a multigrid solver. The convergence of the time-marching iterative procedure is truncated once normalized residuals are below 10 −8 .
Calculations are compared to the experimental data in physical space. The comparison, in axial direction, includes mean values of velocity, mixture fraction, temperature and the main species. The different sections are identified where comparisons in radial sections between measurements and computations are conducted.
The experimental data demonstrate that differential diffusion is important only at the first location (x/Lvis = 1/8) [5,7].
The hydrogen-air flame studied here is not influenced by buoyancy because of the high values of the Froude number [28], which are about 17,000.
In this section, both the (IEM, EMST) mixing simulation results for mixture fraction, temperature and major chemical species are presented and compared with the Barlow et al. [29] experimental data.
RANS/TPDF predictions of axial and radial profiles of mean mixture fraction for the IEM and EMST mixing models are shown in Figs. 2, 3 and 4, respectively. The simulation results of the axial profiles are found to be in good agreement with experimental data. The EMST model predicted better the experimental data than the IEM model. The EMST model treats mixing locally in composition space in contrast with the IEM model.  The radial profiles at x/Lvis = 1/4 and 3/4 shows that both models tend to result in a faster rate of mixing as compared to the experiments, with the EMST model predictions in better agreement, compared to the IEM model. Both mixing models underpredict the experimental data at the center axis. This underprediction begins to decrease at further downstream positions. Figure 5 shows the axial evolution of mean temperature. This evolution presents the same behavior as that of experience. The EMST model, which better describes the physical processes of mixing, predicts best the experimental results compared to the IEM model. The maximum value, which is slightly overpredicted in amplitude, is in axial position the same. This axial position (x/Lvis ≈ 0.65) corresponds to the stoichiometric value of the mixture fraction Zst = 0.028.
The radial profiles of mean temperature are shown in Figs. 7 and 8. This evolution shows that the turbulent mixing models, i.e., the IEM and the EMST models, tend to overpredict the experimental data in the region close to the nozzle exit (at x/Lvis = 1/4), however, the temperature peak of the EMST model is in good agreement with experimental data compared to the IEM model. The location of the peak is well predicted by the IEM model. In the far region (at x/Lvis = 3/4), the EMST model slightly overpredicts experimental data while the IEM model tends to underpredict them. The EMST mixing model yields much better results than the IEM model. The difference between numerical simulation and experimental results may be due to the upstream conditions affecting very much the initial zone of the jet and/or to the effect of preferential diffusion [30,31] that characterize the hydrogen turbulent flames. The non-equilibrium chemistry [30] and/or model of turbulence used may also play an important role in this region. It can be also seen on Fig. 6 that the IEM model predicts an extinction of the flame more quickly than the EMST model.
The EMST and IEM predictions of mean profiles of mass fraction of hydrogen, water and oxygen are shown in Figs. 9, 10 and 11, respectively. The numerical results of the EMST model are in good agreement with experimental data compared with those obtained by the IEM model which underpredict the experimental data at large axial distances (x/Lvis ≥ 0.6). Table 3 summarizes the comparison of averaged CPU Time Usage (s), RAM requirements (in Mbytes) and precisions (%) with respect to the experiments for (IEM, EMST) models with (100, 1000) particle numbers in reduced chemistry. An increase in the particle number in the Monte Carlo method improves largely the precision of the results to the detriment of the RAM and the CPU performances. However, the EMST model cures these two disadvantages. Indeed, this approach allows an improvement in the precision and iteration count, definitely better than those carried out by the IEM model. Nevertheless, a light increase in the time computing (CPU) is noted compared to the IEM model. The detailed mechanism implementation with the T-PDF calculation can be very difficult because the RAM becomes very significant exceeding our current data-processing possibilities.

Conclusion and perspectives
A consistent hybrid RANS-TPDF method has been proposed for the simulation of a hydrogen turbulent jet diffusion flame. The flow field has been solved on the basis of the RSM model. A modeled scalar PDF equation has been solved by a Lagrangian Monte Carlo method developed by Pope [17]. This provides a detailed modeling of the turbulence chemistry interaction. The molecular mixing term is modeled by two models IEM and EMST. A twenty tree-step reaction, eleven species kinetic mechanism for hydrogen-Air combustion is employed here. On the basis of the obtained result, we can draw the following conclusions: (i) The EMST model, which treats mixing locally in composition space in contrast with the IEM model, predicts better the experimental data than the IEM model. (ii) Compared to the near field regions of the nozzle exit, the far field regions are characterized by small differences between numerical simulation and experiment data. The effect of preferential diffusion, which characterizes the turbulent hydrogen flames, is present in the regions near the jet nozzle exit as indicated by Meier [30] and Sanders [31]. The effect of non-equilibrium chemistry also plays an important role in these regions, Obieglo [32]. (iii) Finally, we can say that the difference between numerical results and experimental data is due first to the effect of preferential diffusion and non-equilibrium chemistry characterizing the hydrogen flames. The effect of the reaction mechanism (11 species and 23 reactions) and model (RSM-TPDF) used in this study can also favorite this difference. Heat transfer which is very present in hydrogen turbulent diffusion flames is also a factor to consider. The choice of upstream conditions, which have an important influence on the regions near the jet nozzle exit, is another factor not to be neglected. (iv) To better understand and model the physics of reactive flows, the objective and more elaborate formulations are needed for modeling the most important terms (for exampleφ ij andψ(ε)), present in the exact transport equations of Reynolds stresses, dynamic dissipation and transported PDF, for example, the turbulent convection term and the micro-mixing term. The upstream conditions well-defined may also bring improvements to the obtained results.