A Unified Reduced Model for Auto-Ignition and Combustion in Premixed Systems

In this paper, two complementary chemistry model reduction methods for combustion simulations are further developed and combined. A progress variable model (PVM), which follows the idea of trajectory generated manifolds (TGLDM), is tailored for describing auto-ignition in situations where the influence of molecular transport on chemical reaction is weak, like auto-ignition in media with weak scalar gradients. The other model using the reaction diffusion manifold approach (REDIM) is designed for situations where the interaction of chemistry with molecular transport is essential. The formulation of both models is discussed and implementational issues of each single model are given. Also, each model is tested in its respective range of applicability (quasi-homogeneous combustion under steady/unsteady physical boundary conditions for the PVM, combustion in fields with essential scalar gradients for REDIM). The coupling of the two models into a unified model, which covers combustion in both regimes and during the transitions between regimes, is discussed, based on the global quasi-linearization concept (GQL).


Introduction
A common way to deal with combustion in CFD/ LES simulations is to use simplified, reduced models for chemical reaction instead of the computationally prohibitive detailed chemical kinetics [1]. The reduced models typically introduce two simplifications: On the one hand, the number of variables to be solved for is strongly reduced, mainly because of the abandonment of the multitude of chemical species in the detailed description; typically, some ten to thousands of species are replaced by a few variables that describe chemistry. On the other hand, further simplification consists in the removal of strongly disparate timescales governing the system's development, which also aggravate the detailed description [2].
Formally, the simplification can be described as follows: Instead of solving the full n-dimensional system of PDEs (where n denotes the dimension of the detailed thermo-chemical state vector Ψ, e.g. comprising temperature, pressure and species), the system of PDEs is reduced to m < n dimensions (where m denotes the dimension of the reduced state vector θ). Hence, the m-dimensional system of conservation equations for the reduced state vector θ is given by [3,4]: with v denoting the flow velocity, ρ the density and D * the generalized diffusion matrix [5] and P a projection operator onto the reduced space. The appropriate chemical source terms S as well as the projection matrix and the diffusion matrix for the conservation equations of eq. (1) must be pre-calculated and provided by the model as a function of the reduced state θ. To solve this equation, initial and boundary conditions have to be known. Many reduced approaches have been developed in the past to calculate and store the chemical source terms (e.g. ILDM [2,6], REDIM [3,7], GQL [8], FGM [9,10], flamelet concept [11], flamelet-prolongation of ILDM [12], CSP [13], QSSA [14] or progress variable models [15,16,17,18], and all of them have their advantages and drawbacks, which are e.g. discussed in [19] where also further references to reduction methods can be found. Models normally are tailored for a certain range of applicability (i.e. combustion scenarios). For instance, the flamelet-like models implicitly contain a coupling of chemistry with molecular transport, and their natural field of application is therefore combustion scenarios where chemical reaction is accompanied by transport, e.g., flame-like combustion or auto-ignition at mixing layers of hot gas with unburned gas [20,21]. For describing situations where, e.g. chemical reaction occurs without essential influence of transport, many concepts are available [17,19].
But in practice such scenarios appear either in parallel or sequentially. A typical case is the onset of combustion in compression-ignition engines, which often starts as auto-ignition at some site where the local conditions favor auto-ignition, e.g., by slightly elevated temperatures (hot-spots) [22]. The transport processes at that site are often weak, so that the process before and during auto-ignition is dominated almost exclusively by the reaction source terms. After ignition, fast chemical reactions produce combustion products and high temperatures, thus creating strong gradients with respect to the surroundings, which, in turn, cause strong transport processes at the hot-spot's boundary, giving rise to a flamelike propagating reaction front. Therefore, while the initiation of combustion would be described by a model tailored for homogeneous auto-ignition, the subsequent flame with its transport driven propagation mode would call for a flamelet-like model.
There is a lack of models for such situations with «overlapping effects», so that in practical simulations, often some manual selection and adjustment of models has to be performed.
In an attempt to improve the situation, we combine two reduced simplified models for different combustion scenarios to a unified model: A model tailored for auto-ignition without transport (progress variable model, PVM), and the REDIM (Reaction Diffusion Manifold) model for combustion with diffusive processes.
We demonstrate the performance of these models within their respective domain of validity, by comparing predictions of the reduced models with the outcome of detailed numerical simulations for the same conditions. In particular, for the progress variable model, homogeneous auto-ignition for steady and unsteady conditions with respect to pressure and enthalpy is studied. For the REDIM model, the propagation of a premixed laminar flame is computed.
The PVM and REDIM both describe low-dimensional subsets (manifolds) in state space to which the system's dynamics is confined. A coupling of the two models corresponds to a connection of the two manifolds. A concept for accomplishing such a connection based on the concept of global quasi-linearization (GQL) is developed.
The final outcome is a unified model, which allows an efficient model reduction with a large range of application.

Model Reduction Strategies
As noted above, in this work we combine two model reduction concepts, which are favorable for different combustion regimes: a PVM and a REDIM concept. These will be described in the following, where we focus only on the aspects relevant for the unified model.

Progress Variable Model (PVM)
To treat auto-ignition processes with simplified chemistry, a progress variable model has been developed [15]. It is based on mapping information about the dynamics of chemical reactions from detailed calculations of homogeneous, adiabatic, isobaric reactors onto a single variable, the progress variable. The model follows the concept of the trajectory generated low dimensional manifolds (TGLDM) [23]. A large range of physical conditions under which reaction occurs can be incorporated into the model by computing homogeneous reaction systems for different initial temperatures, pressures and mixture compositions. The calculations for the detailed trajectories were performed with the in-house code HOMREA [24].

Mathematical Model
The thermochemical state of the homogeneous reaction system is described by the state vector Ψ = [h,p, Y 1 ,...,Yn s ] T as a point in the (n s + 2)-dimensional state space. Here h is the specific enthalpy of the mixture, p the pressure and Y i are the mass fractions of the chemical species i (with n s as the number of the chemical species). The simplifications for an isobaric, adiabatic and homogeneous system lead in vector form to with the chemical source term F (Ψ) = [0,0,M i ω i / ρ,...,Mn s ωn s / ρ] T and the density ρ, the molar mass M i of species i and the molar production rate ω i of species i.
Solving this system of ordinary differential equations for a given initial state Ψ 0 = Ψ (t = 0) the temporal evolution of the system is obtained. At any given point in time, the detailed state vector, and therefore also the chemical source, is then known.
A set of different initial values defines a manifold generator [23], and if the manifold defined by the initial conditions is m-dimensional, then the resulting low-dimensional manifold is (m + 1)-dimensional with the reaction progress being the additional coordinate.
Parametrizing a trajectory by some suitable progress variable χ (see below), a detailed state Ψ on a trajectory is completely specified by initial conditions and a value of χ.
The initial conditions for a trajectory can be specified by the initial enthalpy h of the mixture, the pressure p and the unreacted mixture composition Y i,0 , which, in turn, can be expressed by the mixture fraction Z [14], if all admissible initial mixtures result from mixing of fuel and air in different ratios. This can be generalized to incorporate also other quantities like exhaust gas content, exhaust gas composition, statistical moments of these quantities if necessary. When a set of trajectories for different initial conditions is present, a state on the subspace formed by all trajectories can be specified by a reduced state vector θ: Precalculating trajectories, one can compute and tabulate the chemical source term of the progress variable for use in eq. (1) as a function of the reduced state vector.

Definition of the Progress Variable
Various models rely on the notion of some measure for chemical progress, and several approaches have been taken to define a chemical progress variable [15,17,18,25,26]. Some are based on temperature or the thermal energy released by chemical reaction, others on chemical species or combinations thereof. The importance of temporal monotonicity of that progress variable is stressed in many of these works. In general, one may impose this and some additional requirements on a progress variable [15,27], namely monotonicity in time during reaction ( ), positivity and boundedness, such that a normalization is possible (χ = 0 for unreacted fuel/air mixture, χ = 1 for chemical equilibrium), sensitivity to important stages of chemical reaction (ignition process, combustion) and universality (applicability to all reaction systems).
Finding such a variable is nontrivial [15,27]. The entropy of the mixture is guaranteed to be monotonically increasing in a closed, adiabatic system by virtue of the second law of thermodynamics and the irreversibility of the ignition and combustion process. Furthermore, entropy also displays good sen-sitivity to chemical activity, including early stages of ignition, which results in a good resolution of all phases of the ignition process.
Based on the rate of production of the specific entropy ṡ of the mixture the progress variable χ and its time derivative are defined as: with (assuming an ideal gas mixture) and s i denoting the specific entropy, c p,i the specific heat capacity at constant pressure and p i the partial pressure of species i. Furthermore, T denotes the temperature and the universal gas constant.
This definition fulfills all requirements specified for the progress variable above, and is therefore used in this work.

Reaction Diffusion Manifold (REDIM)
The REDIM method has been described in detail in [3,7]. In the following, only a short account of the REDIM method will therefore be given.
Using the thermochemical state vector Ψ = [h,p, Y 1 ,...,Yn s ] T as introduced above, the conservation equations for a reacting system are given as: with v denoting the flow velocity and D the generalized diffusion matrix [5].
Applying the framework of invariant system manifolds, it is assumed that the evolution of the system is part of the m-dimensional subspace, namely the manifold, which is parameterized by the reduced coordinates θ. This means that the vector field, described by the right-hand side of eq. (6) belongs to the tangent space of the manifold and the projection into the normal space of the manifold has to be zero: using the unit matrix I, the derivative Ψ θ of Ψ with respect to θ and its Moore-Penrose pseudo-inverse [28].
of eq. (7) to a system of PDEs [3] leads to the evolution equation of the manifold: For the extension of the REDIM method to systems with detailed transport models with a generalized diffusion matrix D see [29].
Starting from an initial guess and solving the system of PDEs of eq. (8) by using Dirichlet-type boundary conditions, one obtains the manifold as the stationary solution. For detailed procedure, see [3,30]. Finally, the chemical source terms to use in eq. (1) can be tabulated as a function of the reduced state vector θ.
Although the concept is quite simple, the solution of the evolution equation is challenging. To our knowledge this concept has only been applied to H 2 , CH 4 and syngas combustion. For this work we have improved the numerical solution method such that it allows the treatment of large reaction mechanisms, and the example shown below is a first extension to a combustion system relevant for engine combustion.

Tabulation Strategies
To use the models described above for calculations with simplified chemical kinetics, one has to tabulate and provide means for efficient evaluation of the chemical source terms of the reduced variables as a function of the reduced state vector θ. The strong non-linearity of the chemical source terms is a particular challenge for both the tabulation strategy and the lookup (function evaluation) routine.
Several tabulation strategies exist; some of them work in situ during a calculation [31]. Here we focus on pre-calculating tabulations.
A straightforward approach for multivariate function tabulation is to use a rectilinear grid point system with equidistant spacing along each dimension. The resulting grid structure is independent of the function to be represented, but has to be specified a-priori by the user. The table-lookup (by interpolation) is simple and computationally fast.
However, the strongly non-linear nature of the chemical source terms often requires a locally very fine grid, while in other regions a coarse grid would be sufficient. The uniform grid spacing has to be adapted globally to the finest required resolution.
For a large number of tabulation variables (high-dimension of the table), the equidistant approach can become infeasible due to the large memory requirements, which increase exponentially with dimension (curse of dimensionality [32]).
An interesting approach to overcome the curse of dimensionality is the application of adaptive grids. Based on the framework of sparse grids [33,34,35] an adaptive sparse grid algorithm [36] has been implemented. Starting with the boundary grid points of a user-defined rectangular tabulation domain, the step size is successively halved per dimension. For each originating grid point the tabulation accuracy and the interpolation error, respectively, is controlled by checking a defined criterion. If necessary, the grid point is added to the dataset of the tabulation. The refinement continues until the required accuracy or the maximum refinement level is reached.
A requirement of the adaptive grid algorithm is the representation of the function (i.e. the chemical source terms) based on discrete supporting points in hierarchical basis [37,38]. In contrast to the representation in classical nodal basis a hierarchy on the supporting points is defined by successively refining the grid and halving the step size, respectively.
For simplicity, the following description is limited to the one-dimensional case. Let f (x) be a univariate function on the interval [0,1]. Then the interpolant based on n supporting points can be represented by a linear combination of the basis function β with the coefficients u i : using and Δx denoting the spacing of the supporting points.
In hierarchical basis the spacing is given to Δx = 2 -l with the hierarchical level l. In this case, the basis function depends on the hierarchical level l and the interpolant is given by the sum over all levels:  Figure 1 shows the approximation of a one-dimensional example function for hierarchical basis. It can be seen that the mesh size of the hierarchical basis function depends on the level l.
The transition to the multi-dimensional case is via a tensor product approach [34], applied separately to each basis function. This results in the multivariate basis function: There are many possibilities to choose the onedimensional basis functions β (like polynomials and wavelets [39]). In this work we will use the hatfunction as basis function for linear interpolation: β because this allows a simple, fast and efficient implementation.

Coupling Strategy
Both methods used in this work rely on lowdimensional manifolds in composition space. However, they inevitably are based on different model reduction assumptions, and therefore the low-di-mensional manifolds do not coincide. Therefore it is necessary to have a smooth transition between the two manifolds. To combine the PVM and REDIM approach, a transition between the two models is implemented, based on the global quasi-linearization approach (GQL) [8].
A simple connection can be established by identifying points on the two manifolds that feature a correspondence with respect to the fast chemical subspace. For instance, given a reduced vector θ PVM which describes a state on the manifold described by the progress variable model, a corresponding reduced state θ REDIM on the REDIM manifold can be found by the condition that both points are on the same trajectory of the fast subspace, or in other word, that the coordinates in the slow subspace are equal: with the matrix M s representing the global slow chemical directions identified by the GQL. By this, the reduced state θ PVM on the PVM manifold is projected to the reduced state θ REDIM on the REDIM manifold in direction of the fast chemical processes. The slow subspace matrix M s is obtained via the standard GQL technique and the mapping between the two manifolds is stored together with the two tables. Based on an identification of the respective combustion regime, the transition from one to the other manifold is performed during the application of the reduced schemes in the CFD calculation.

Implementation Details
For comparison and validation, we computed different combustion systems with detailed and simplified chemistry.
For all detailed calculations, a reaction mechanism for toluene reference fuel (TRF) [40]  For domains where chemistry governs the process, a PVM using an adaptive sparse grid was built. The parameter range in progress variable was set to 0.0 < χ < 1.0, in pressure to 5 bar < p < 55 bar and mixture fraction to 0.25 < Z < 0.013. Instead of the specific enthalpy, the temperature T * of the corresponding unreacted mixture was used as a tabulation parameter in the range of 700 K < T * < 1250 K.
To control the refinement, the local interpolation error of the source term of the progress variable was checked:  This example shows clearly that the tabulation strategy based on adaptive sparse grids allows a very efficient and accurate tabulation.
The accuracy could be increased further by using a finer grid for the tabulation; however, the observed deviation is already much smaller than, for example, the typical deviation between predictions of two different detailed reaction mechanisms, so that further refinement would be hardly worthwhile in practice.

Assessment of the Generality of the PVM
The test case shown above validates the tabulation approach. It does however treat only systems with constant T * , p and Z, which had been used for ( ) using a relative tolerance of ε rel = 0.1, an absolute tolerance of ε abs = 10 -4 and a maximum refinement level of 20+d. Furthermore, the refinement was limited to Δχ > 10 -4 , ΔT * > 2 K, Δp > 2 bar and ΔZ > 10 -2 .
For reduced calculations of domains governed by reaction and diffusion, the REDIM method was used with four parameters, namely the mixture fraction Z, pressure p, temperature of the unreacted mixture T * and a chemical variable, on an equidistant grid for the reduced variable with 51 grid points for constant T * , p and Z. The gradient estimate [3] was taken from detailed calculations of laminar, adiabatic, isobaric, pre-mixed flames.
Besides the chemical source terms, other properties of the mixture (e.g. mean molar mass , temperature T, density ρ, specific heat capacity at constant pressure c p and rate of temperature increase by chemical reaction Ṫ chem ) are included in the tabulations as well. For analysis and comparison with detailed calculations, mass fractions of chemical species are also part of the tabulations.

Assessment of the Tabulation Quality of the PVM
One application of the PVM with stationary conditions is the reduced calculation of ignition processes in a homogeneous, adiabatic, isobaric reactor.
Given the initial conditions the system of PDEs of eq. (1) yield in this case: This test case assesses the quality of the PVM exactly for the conditions it was built for and therefore validates the tabulation strategy (interpolation error). Figure 2 compares the development of temperature and mass fraction of an intermediate species (CH 2 O) for detailed and reduced calculation, for χ 0 = 0, * 0 T = 800 K, p 0 = 10 bar and mixture fraction according to equivalent ratio ϕ 0 = 2. One can see that multistage ignition, which is an important phenomenon in ignition processes, is represented realistically by the reduced model. For rich mixtures, the drop of temperature in the later stages of ignition is also reproduced. Besides the temperature also the intermediate species is close to the detailed simulations.
ign the manifold generation. Therefore it is interesting to see whether the model can handle conditions with varying T * , p and Z, which had not been used for the manifold calculation.
To validate the PVM for instationary conditions with respect to enthalpy and pressure a second test case was chosen, namely model calculations of homogeneous engine cycles. For this, the temporal evolution of the system's volume V was prescribed as: where b is the bore (100 mm), s the stroke (83 mm) and r the length of the conrod (149 mm). The compression ratio used was 14:1.
The crank angle α changed with time according to: yielding its temporal derivative: Here, n denotes the engine speed. Eq. (19) imposes a temporal boundary condition onto the system, rendering it non-autonomous.
For homogeneous engine cycle calculations the system of PDEs of eq. (1) yield the following system of ODEs: (25) with m denoting the total mass, the specific heat capacity at constant pressure for the unreacted mixture, the mean molar mass and κ the isentropic exponent. The «source» terms of temperature Ṫ chem and mean molar mass are provided by the model tabulation. The initial conditions were set to χ 0 = 0, * 0 T = 500 K, p 0 = 1 bar and mixture fraction according to equivalent ratio ϕ 0 = 1.  Figure 4 shows the evolution of pressure for both detailed and reduced calculations for a variation of engine speed. Even though the PVM is based on constant-pressure adiabatic trajectories (with pressure appearing as a tabulation parameter), it still reproduces the chemical dynamics in a system with an unsteady pressure evolution very well. Especially when ignition appears before the top dead center (TDC), the deviation in ignition time between detailed and reduced calculations is less than 1 CAD. Even for higher engine speed and later ignition points, respectively, the pressure curve fits quite well to the detailed results (deviation in ignition time less than 2 CAD). For the presented, engine relevant parameters, the onset of auto-ignition is captured well. The peak pressure near TDC differs by maximum 3.3 %.
This model could e.g. also be used in a stochastic reactor model for engine combustion [41]. Fig. 4. Variation of engine speed for calculation of homogeneous engine cycle using the reduced model with an adaptive sparse grid, T = 500 K, p = 1 bar, ϕ = 1, square: n = 1000 min -1 , triangle: n = 3000 min -1 , circle: n = 6000 min -1 , solid: detailed calculation, dashed: reduced calculation.

Assessment of the Model Quality in the Domain of Coupling of Reaction and Diffusion
A third test case is flame propagation. Because the model couples a PVM with a REDIM model, it switches to the REDIM based table if coupling between reaction and diffusion becomes important. A test case for this domain are laminar flames. The system of PDEs yielding for the reduced calculation of flamelet trajectories are described in [3]. Note that the reduced model was designed to handle premixed or slightly stratified systems due to the choice of the gradient estimates for the REDIM construction. In the case of non-premixed flames, the reduced model could be easily adapted by using other gradient estimates [30].

Coupling of the Models
Finally, a case for an auto-ignition at some hotspot followed by flame propagation is presented.
This one-dimensional simulation starts with a stoichiometric fuel/air mixture at 10 bar. The initial temperature profile features a 500 K baseline, with a superimposed hot-spot with a temperature elevation of 7 K and a radius of 5 mm. The system then reacts under adiabatic-isobaric conditions. Since initially there are practically no spatial gradients, this process is modeled by the progress variable model. After a certain time (corresponding to the ignition delay time of a homogenous reactor at the given conditions), the hot-spot auto-ignites. This homogeneous reaction is depicted by red lines in Fig. 7, which show the temporal development of the spatial CO 2 profile. With time, the reaction leads to a steep temporal rise in CO 2 inside the hot-spot, while the surroundings remain in a practically unreacted state. The chemical reaction in the hot-spot finally causes strong gradients with respect to temperature and species at the hot-spot's boundary. To simulate the subsequent development, therefore, the REDIM model is used.
The coupling between PVM and REDIM in this case is simply performed by mapping the CO 2 values of the last PVM time-step to the REDIM model (using a matrix M s in eq. (16) consisting of one line of all zeros, with only the column corresponding to CO 2 being set equal to 1).
This further development is indicated by black lines in the figure. The system continues to react (with the chemical source terms now governed by the REDIM model), leading to a further increase of CO 2 , and the hot-spot expands. The gradients at the boundary continue to increase with time. With further development in time, a flame would develop at the boundary (not visible here) and propagate.
The combined effect of auto-ignition with subsequent flame development is captured by the unified PV/REDIM model. The laminar flame speed is an important characteristic of premixed flame propagation. Figure 6 compares the laminar flame speeds obtained by detailed and reduced chemistry. The results from the reduced calculations show a slightly larger laminar flame speed. The deviation increases from 12% to 24% in the direction of richer mixtures.

Conclusions
In this paper, two reduced combustion models are presented and combined, which aim at describing different complementary combustion scenarios. One of the models is a progress variable model, designed for describing homogeneous reactions. It is based on a generic definition of a progress variable (derived from the specific entropy of the reacting mixture). Even though this model is based on chemical reaction trajectories of adiabatic constant-pressure systems, it has performed well also in situations with an unsteady evolution of pressure, namely in representative simulations of homogeneous engine cycles. On the other hand, for combustion with strong coupling of reaction with molecular transport, the REDIM-model is used. Both models are implemented for a complex fuel mixture (toluene reference fuel).
For situations where both homogeneous reaction and reaction/transport-dominated modes occur simultaneously or in temporal sequence, an unified model which couples PVM and REDIM is introduced and successfully applied to a simple test case (hot-spot auto-ignition and expansion).