Implementation of the Scalar Dissipation Rate in the REDIM Concept and its Validation for the Piloted Non-Premixed Turbulent Jet Flames

In order to address the impact of the concentration gradients on the chemistry – turbulence interaction in turbulent flames, the REDIM reduced chemistry is constructed incorporating the scalar dissipation rate, which is a key quantity describing the turbulent mixing process. This is achieved by providing a variable gradient estimate in the REDIM evolution equation. In such case, the REDIM reduced chemistry is tabulated as a function of the reduced coordinates and the scalar dissipation rate as an additional progress variable. The constructed REDIM is based on a detailed transport model including the differential diffusion, and is validated for a piloted non-premixed turbulent jet flames (Sandia Flame D and E). The results show that the newly generated REDIM can reproduce the thermo-kinetic quantities very well, and the differential molecular diffusion effect can also be well captured.


Introduction
An accurate modeling of the interaction between mixing processes and the chemical kinetics is an important issue in numerical simulations of turbulent non-premixed flames [1,2,3]. One of the challenges is the modeling of the chemical source terms, since they are strongly introduced by the turbulent mixing process and molecular transport. The Probabilty Density Function (PDF) method [4,5] is a widely used method, because the chemical source term can be solved in a closed form, without any modeling [4,6].
However, since the governing conservation equations for species including detailed chemical kinetics have a high dimensionality due to the large number of species and the high stiffness due to strong non-linearity of elementary reaction rates [7,8], which leads to an extreme complexity and high computational load with respect to the computational time and memory storage, there is an urgent for reduced models for the chemical kinetics. Manifold-based simplified chemistry is one of the major methods for the simplification of chemistry. Examples include The Intrinsic Low-Dimensional Manifolds (ILDM) method [9], the Flame Prolongation of ILDM (FPI) method [10], the Steady Laminar Flamelet Model (SLFM) [11], the Reaction-Diffusion Manifolds (REDIM) method [12], the Flamelet Generated Manifolds (FGM) [13] and many others.
Although the application of the manifold-based simplified chemistry can largely reduce the computational cost, one additional problem arises when they are applied to model turbulent flames: It is well known that the turbulent mixing can significantly affect the chemical reactions [14,15]. It is therefore important to know whether the simplified chemistry can reproduce the effect of the turbulent mixing on the chemical kinetics correctly.
Throughout the whole work, we focus on the REDIM methodology, which takes into account the molecular transport process in the generation of reduced chemistry. Recently in [16] the effect of turbulent mixing process on the chemistry has been successfully modeled in the framework of RE-DIM concept. There the gradients of reduced coordinates are considered as additional independent progress variables in the parametrization of the REDIM reduced chemistry, and these gradients, depending on the turbulence intensity, can be determined simultaneously. However, this algorithm has a drawback that the REDIM can have a high dimension, because for each reduced variable its gradient must be considered as additional progress variable in the tabulation of the REDIM reduced chemistry.
For the non-premixed laminar and turbulent flames, the scalar dissipation rate is a key quantitiy to represent the turbulent mixing process [17], and depends on the scalar gradients. In other words, if we provide the scalar gradient into the REDIM from flame scenarios with different scalar dissipation rates, then the REDIM reduced chemistry can incorporate the influence of turbulent mixing process on the chemical kinetics in terms of the scalar dissipation rates, and the impact of the concentrations gradients on the chemistry -turbulence interaction can be represented.
This work is organized as follows: Section 2 provides a short outline about the models used for the turbulence and the REDIM reduced chemistry taken into account the effect of turbulent mixing in terms of scalar dissipation rate. Section 3 explains in detail the numerical procedure for the generation of the REDIM reduced chemistry and also the coupling between the turbulence model and the REDIM model. Finally, the validation of the proposed algorithm is shown in Section 4. Important issues will then be summarized in the conclusion section.

Mathematical models
The method for the numerical simulation of the turbulent reacting flows used in this is a hybrid CFD/transported-PDF model. This model consists of two parts which are coupled with each other: one is the CFD program solving the Reynolds-averaged Navier-Stokes (RANS) equation [2,18,19] to obtain the hydrodynamic quantities; the other part is the PDF program solving the transported-PDF equation [4] to obtain the thermo-kinetic quantities of the flow. Furthermore, in order to reduce the computational effort caused by the detailed chemical kinetics, the Reaction-Diffusion Manifolds (REDIM) method [12], one of the model reduction for the chemical kinetics, is used. In the following, all three methods will be briefly outlined. The coupling between these methods are presented in the next section.

Reynolds-averaged Navier-Stokes (RANS) method
The basic idea to derive a Reynolds-averaged Navier-Stokes (RANS) equation is the Favreaveraging, in which a quantity q(x, t) can be split into a density-weighted Favre-averaged quantity and the corresponding fluctuation q′′(x,t) [2,18,20]: By applying the Favre-averaging, the Favreaveraged continuity equation and momentum equation for RANS can be written as [2,18,20]: where i,j = 1,2,3, and the Einstein summation convention is used. In both equations, t is the time, ρ the density of mixture and p the pressure. The u = (u 1 , u 2 , u 3 ) T is the vector of the components of the flow velocity, and τ ij is the stress tensor [19]. δ ij denotes the Kronecker delta.
The Reynolds stresses appear in an unclosed form and need to be modeled. In the present work, the Reynolds stresses are obtained from the PDF part. Moreover, we consider the system as adiabatic. The Favre-averaged temperature is not determined by solving an Favre-averaged conservation equation for the energy, but from the PDF part (see below for details).

Transported Probability-Density-Function (TPDF) method
In the present work, a one-point one time joint Probability Density Function (PDF) [4] of velocity, composition and turbulent frequency f ωuΨ (Θ, v, Φ; x, t) [3] is employed. The motivation to obtain the f ωuΨ (Θ, v, Φ; x, t) is that important information such as average, variance and co-variance can be obtained [4,6]. The PDF of a flow is a-priori unknown, but can be obtained by solving a transported-PDF equation to obtain the instantaneous PDF at every position in the flow field [4,6].
Eq. 3 In the derived transported-PDF equation, the thermo-kinetic source term appears in a closed form [4].

� � ′′ ′′
To solve the integral-differential transported-PDF equation efficiently, the Monte-Carlo particle method was suggested [4], and the PDF is represented by an ensemble of notional particles evolving according to the stochastic processes described by either ordinary differential equations (ODEs) or stochastic differential equations (SDEs): • The position of each notional particle evolves by its own velocity in form of ODE. • The evolution of particle velocity is modeled by the simplified Langevin model (SLM) in form of SDE [4,19]. • The turbulent frequency of each notional particle is modeled based on the gamma-distribution model [19] in form of SDE. • The evolution of the thermo-kinetic state is calculated in form of ODE as: where i = 1,2, ···, N p (N p : number of notional particles per CFD cell). S(Ψ *,i ) is the thermo-kinetic source term, and M describes the molecular mixing process. It should be mentioned here that the model for the molecular mixing process M affects the accuracy of the overall numerical simulation significantly. Various mixing models have been proposed including, but not limited to, the Interaction by Exchange with the Mean (IEM) [21], the Curl's coalescence and dispersion Model (CD) [22], the Euclidean Minimum Spanning Trees (EMST) [23], the Conditional Moment Closure (CMC) [24] and the Multiple Mapping Conditioning (MMC) [25]. A review on different mixing models and the discussion on their advantages and disadvantages can be found in e.g. [6,26]. Although the PDF method deals with the chemical kinetic in an exact way without any modeling, the needed system of the ODEs (Eq. 4) using detailed chemical kinetics can become very complex in terms of dimensionality and non-linearity [7,8], leading to a high computational cost in the numerical integration for the evolution of thermo-kinetic state due to chemical reaction. Therefore, in the next section, method for the simplification of the chemical kinetics are introduced, so that the computational effort due to chemical kinetics can be significantly reduced without losing accuracy.

Reduction of Chemical Kinetics: Reaction-Diffusion Manifolds (REDIM)
We begin with the general form to describe the evolution of the thermo-kinetic states (enthalpy h, pressure p, species mole number ϕ i ) in a reacting flow in terms of state vector Ψ = (h, p, ϕ 1 , ϕ 2 , ..., ϕn sp ) T [7,8]: where ρ is the mixture density, ϕ i = w i /M i the ratio between mass fraction and molar mass of i-th species, and n sp the number of species. D is the n × n -dimensional transport matrix including molecular diffusion and heat transport F(Ψ) the n -dimensional source term accounting for the chemical reactions, and u the vector of velocity. The high dimensionality and stiffness of the governing conservation equations for species including detailed chemical kinetics increase significantly the computational effort [8]. However, as shown in [27], analysis based on Direct Numerical Simulation (DNS) reveals two major observations. Firstly, not all composition space is accessed. Secondly, the whole dynamic system is constrained around some low-dimensional manifolds due to the existence of the characteristic time-scales within a combustion system differing by orders of magnitude typically from 1 s to 10 -10 s [9]. Following both observations, one efficient way for the reduction of chemical kinetics is to find out such low-dimensional manifolds that can be used to approximate the full system dynamics: In this work, the Reaction-Diffusion Manifolds (REDIM) method [12] is used for the simplification of chemical kinetics, which takes into account transport processes (e.g. convection, molecular diffusion, heat transport). In the following, we will first briefly review the concept of the REDIM method. More details about the theory can be found in e.g. [12]. According to the [12], by solving the following evolution equation: its steady solution yields the REDIM. In this equation, Ψ θ is the matrix of partial derivatives of the thermo-kinetic state vector Ψ with respect to the reduced coordinates θ. is the pseudo-inverse of Ψ θ fulfilling the condition: . The most important feature in this REDIM evolution equation (Eq. 7) is that besides the initial and boundary conditions, the integration of this evolution equation requires additional information about the gradient estimate γ(θ) = grad(θ). As shown in [28,29], for both premixed and non-premixed combustion system, the gradient estimate can be important for the REDIM with low dimensions (usually one-dimensional and two-dimensional), but the REDIM becomes less sensitive to the gradient estimate if the REDIM is three-dimensional or higher [28,29]. Moreover, it has been shown in [30] that the REDIM reduced chemistry does not depend on the choice of the reduced coordinates θ, since the REDIM is an invariant manifold.
Different approaches have been proposed in previous works, including 1) estimating the gradient by constant values [30]; 2) taking the gradient from typical flame scenarios such as burning flames and mixing process [31]; 3) taking the gradient from the DNS data if available [32].
In this work, we use the second approach, in other words, the REDIM is built by taking the gradient from detailed chemistry calculations of steady laminar counterflow flame scenarios. This has the advantage that these flame scenarios contain information about their scalar dissipation rate defined as [33]: Eq. 8 at the location of stoichiometric mixture fraction ξ st , and D ξ is the diffusion coefficient of mixture fraction ξ, and in the framework of turbulent flow it can be modeled according to [18,11] based on the mixture fraction variance ξ′′ and the turbulent frequency ω as: where the model constant C χ is usually set to be 2 [11].
The scalar dissipation rate is an important variable in turbulent combustion systems [1,34] which governs the interaction between the mixing process and the chemical reaction process. If one takes the gradient from these flame scenarios, the REDIM incorporates the scalar dissipation rate from turbulent flows automatically, and the effect of the scalar dissipation rate χ on the chemical kinetics is included in the REDIM reduced chemistry kinetics.
In the tabulation of the REDIM chemistry, the scalar dissipation rate is then an additional variable besides reduced coordinates θ. In this case, the REDIM reduced chemistry is tabulated as a function of θ and χ: Eq. 10 and the effect of the scalar gradient on the turbulence-chemistry interaction can be represented in terms of the variable χ. In the next section, the generation of REDIM integrated with the scalar dissipation rate will be discussed in details.

Numerical Procedure
In this section, the numerical procedure for the generation of the REDIM reduced chemistry taking into account the scalar dissipation rate will be discussed. Then the numerical algorithm for the hybrid CFD/transported-PDF model coupled with REDIM will be outlined.

Implementation of scalar dissipation rate in the REDIM concept
To illustrate the implementation of scalar dissipation rate in the REDIM, a non-premixed CH 4 diffusion flame in counterflow configuration is considered, corresponding to the Sandia Flame inlet conditions [35]: are performed using the in-house HOMREA and INSFLA code [36]. A detailed transport model based on Curtiss-Hirschfelder approximation [37] including the Soret effect [38] is used. The GRI 3.0 mechanism [39] is used. A detailed description for the REDIM generation procedure about the numerical issues including the initial conditions, boundary conditions and the grid sensitivity with respect to the REDIM accuracy can be found elsewhere such as [28,40,41]. In the present work, two REDIMs with different dimensions are generated: • One is the two-dimensional (2D) REDIM given as: This 2D REDIM is a bunch of one-dimensional REDIMs parametrized by one reduced variable θ 1 and different scalar dissipation rate χ.
The initial solution for every one-dimensional REDIMs uses simply the Burke-Schumann solution [42], and the boundary conditions are set as fixed boundaries determined by two points: the thermo-kinetic states of fuel side and of oxidizer side, as represented in Fig. 1(a). Then the gradient estimate γ(θ) in Eq. 7 are taken from different detailed flame solutions with the corresponding different scalar dissipation rates, as shown in Fig. 1(b).
The other one is the three-dimensional (3D) given as: This 3D REDIM is a bunch of two-dimensional REDIMs parametrized by two reduced variables θ = (θ 1 , θ 2 ) T different scalar dissipation rate χ.
The boundary conditions for each two-dimensional REDIMs are the same. The very left and the very right points represent the thermo-kinetic states of the fuel and oxidizer sides, which are set as fixed values. The upper boundary (red line in Fig. 2(a)) is determined by the Burke-Schumann solution [42] and the lower boundary by the pure mixing process without reaction. Both upper and lower boundaries are set to be moving boundaries and can be adapted according to the algorithm proposed in [28].
The gradient estimates γ(θ) in Eq. 7 are also taken from different detailed flame solutions with different scalar dissipation rates. And the gradient estimates are extrapolated to the upper and lower boundaries of the whole domain, as shown in Fig. 2(b).   Figure 3 shows representative 2D REDIM (a) and 3D REDIM (b) surface for different values of scalar dissipation rates χ. Note that for the 2D REDIM, the REDIM must coincide with the detailed flame solution with the same scalar dissipation rate, if the exact gradient is supplied. From this figure, the effect of the scalar gradient on the reduced chemical kinetics can be observed clearly. For larger scalar dissipation rate, corresponding to a higher scalar gradient (c.f. Fig. 2(b)), the reactants mix with each other in a less complete degree and the system departures away from the equilibrium state. This phenomena indicates that the influence of the scalar gradient on the state space can be well represented in the framework of REDIM concept with scalar dissipation rate as an additional variable. Figure 4 outlines the whole numerical procedure for the hybrid RANS/transported-PDF coupled with REDIM reduced chemistry. In the RANS part (blue color), the Favre-averaged continuity (Eq. 2) and Navier-Stokes equations (Eq. 3) are

Numerical Procedure for the coupling of hybrid RANS/transported-PDF/ REDIM models
The Favre-averaged temperature for one CFD cell (R g : the gas constant) is calculated from all notional particles T * (θ * ,χ) in the same CFD cell in the PDF part. The coupling between the PDF part and the RANS part is achieved by determining the averaged pressure for the CFD cell: solved to get the averaged density and velocity for each CFD cell, and will be fed into the PDF part. In the PDF part (orange color), N p notional particles are distributed in each CFD cell. Then the evolution equation for the spatial position X * , the velocity fluctuation u′′ ,* ) and the turbulent frequency ω * are numerical integrated, consistent with the algorithm used in [43,44]. After that the change of the reduced coordinates for each notional particle is determined by the mixing model.
For the reaction process, the REDIM subroutine will be called. The scalar dissipation rate of each notional particle will be calcuated using Eq. 9. Based on the information of the reduced coordinates θ and the scalar dissipation rate χ, the reaction rates RR(θ * ,χ) and the temperature T * (θ * ,χ) can be retrieved from the REDIM table via a table look-up. Then the evolution of the reduced coordinates due to chemical reaction can be calculated as: Eq. 12 and the Reynolds stresses determined from the SLM model in the PDF part.

Results and discussion
To validate the generated REDIM reduced chemistry and to investigate whether the effect of the turbulent mixing process on the REDIM can be re-produced if one considers the scalar dissipation rate as additional progress variable, the wellknown piloted non-premixed methane/air diffusion flames, Sandia Flames [35], are selected as representative validation flames. For this flame series, species such as H 2 need to be especially well predicted to capture the effect of the turbulent mixing due to its high diffusivity [45]. The Sandia Flame series consist of • a main jet with diameter D j = 7.2 mm and mixture composition of 25% methane and 75% air by volume with temperature T j = 294 K, and • a co-axial pilot with an inner diameter of 7.7 mm and an outer diameter of 18.2 mm and a mixture of C 2 H 2 , air, CO 2 and N 2 under lean condition (Φ = 0.77).
A more detailed description for the experimental configuration can be found in [35]. The Sandia Flame D (Re = 22000) and Sandia Flame E (Re = 36000) [35] are investigated in the present work. The mixture fraction is defined consistently with the definition used in the experimental measurements [35]:  Table. The particle cloning and clustering algorithm [6] is used to ensure a high computational efficiency, avoiding excessive number of notional particles in one CFD cell.
The boundary conditions for the mean velocity, the Reynolds stresses, the turbulent frequency and the turbulent kinetic energy are the same as the numerical setup used in [46].
For the chemical kinetics, two REDIM reduced chemistries are applied: • N 2 -χ -REDIM: this is a 2D REDIM tabulated as: • N 2 -θ 2 -χ -REDIM: this is a 3D REDIM tabulated as: where the second reduced coordinate θ 2 is given as: The choice of this combination is based on the investigation on the coupling of the turbulent mixing processes with manifold-based simplified chemistry (more details in [47]).
For both cases, the ϕN 2 is selected as the first reduced variable because it is chemical neutral and represents the mixing process. Note again that the generation of the REDIM reduced chemistry is invariant with respect to the choice of θ, while the application of the REDIM coupled with the mixing process in the turbulent simulation requires a suitable choice of θ (more details in [47]). Figure 5 shows the predicted burning index of the temperature BI(T) [48] introduced to describe the degree of local extinction by using the N 2 -χ -REDIM chemistry and the N 2 -θ 2 -χ -REDIM chemistry along the centerline. Remember that BI(T) = 1 means complete burning solution, and BI(T) = 0 means complete extinction. It is seen that the BI(T) for Sandia Flame D is almost equal to 1 everywhere, indicating that the Sandia Flame D is a stable flame and low degree of local extinction is observed. The predicted BI(T) using the N 2 -χ -REDIM chemistry agrees very well with the experimental measurements [35]. However, C χ 2.0 Scalar dissipation rate (Eq. 9) [11,18] for Sandia Flame E, the lowest BI(T) is at position x/D j = 7.5 and increases to around 1 at positions x/D j ≥ 30 along the centerline. This shows that the a moderate degree of local extinction can be observed at positions x/D j = 7.5 and 15, and a re-ignited burning solution is obtained at position x/D j = 30. The numerical calculation shows that using the N 2 -θ 2 -χ -REDIM chemistry can well capture this phenomenon, while using the N 2 -χ -REDIM chemistry overpredictes the burning index of temperature, indicating an underprediction of local extinction. Such inaccurate prediction of the local extinction can cause the inaccurate prediction for the thermo-kinetic quantities, as shown in the following text. Figure 6 shows the conditional Favre-averaged mass fraction of CO ( Fig. 6(a)) and H 2 ( Fig. 6(b)) over mixture fraction ξ at three different locations (x/D j = 7.5, 15 and 30) for Sandia Flame D. It is observed that for Sandia Flame D, the numerical calculations using N 2 -χ -REDIM and N 2 -θ 2 -χ -REDIM provide almost the similar results. This is because, as already mentioned above, Sandia Flame D has a low degree of local extinction and is a highly stable flame (c.f. Fig. 5). Therefore, the flame front of the Sandia Flame D can be considered as an ensemble of laminar steady flamelets, and the N 2 -χ -REDIM (a bunch of various 1D REDIM depending on the χ -values) is already accurate enough. The only inaccuracy predicted by using the N 2 -χ -REDIM is the conditional mass fraction of H 2 at position x/D j = 30. At this downstream position, the H 2 value predicted by N 2 -χ -REDIM over-predicts the differential molecular diffusion effect, which has also been confirmed in [49] where the flamelet model for treating differential diffusion effect was applied. However, use of N 2 -θ 2 -χ -REDIM can predict the H 2 value at both upstreams and downstreams very well, indicating that the effect of the differential diffusion is well captured.
In the Fig. 7, we show furthermore the conditional Favre-averaged mass fraction of CO ( Fig.  7(a)) and H 2 ( Fig. 7(b)) over mixture fraction ξ at three different locations (x/D j = 7.5, 15 and 30) for Sandia Flame E, in which a moderate degree of local extinction and re-ignition can be observed. The under-prediction of local extinction at positions x/D j = 7.5 and 15 by using the N 2 -χ -REDIM (c.f. Fig. 5(b)) causes an over-prediction of CO and H 2 mass fractions, while an accurate prediction of local extinction by using the N 2 -θ 2 -χ -REDIM chemistry leads to a good agreement of CO and H 2 mass fractions compared to the experimental measurements [35]. At position x/D j = 30 where one has a burning solution (BI(T) ≈ 1, c.f. Fig. 5(b)), the result of CO mass fraction using the N 2 -χ -REDIM again have a good agreement with the experimental measurements [35], while the overestimate of the differential diffusion is again observed in H 2 mass fraction (same as in Sandia Flame D, see Fig. 6(b)). It should be mentioned that a slight over-estimation of CO mass fraction has also observed in [46] where the GRI 3.0 [39] detailed chemistry model is used.
The remaining but also interesting issue is to check the predicted scalar dissipation rate. Although the experimental measurements for the scalar dissipation rates in the Sandia Flame D and E are not provided, Pitsch and Steiner investigated in [50] the conditional Favre-averaged scalar dissipation rates for Sandia Flame D based on the Large-Eddy-Simulation (LES) method using the Flamelet model. It can be observed from Fig. 8 (a) (b) that the in the present work using N 2 -θ 2χ -REDIM chemistry gives the same qualitative (or even quantitative) behavior as in [50]: 1) the decreases far from the jet outlet, and 2) the high values of at x/D j = 7.5 indicate a higher  degree of local extinction (lower value of burning index, c.f. Fig. 5(a)), while the low values of at x/D j = 30 indicate a burning solution (higher value of burning index, c.f. Fig. 5(a)).