## Abstract

In the present work, we illustrate the process of constructing a simplified model for complex multi-scale combustion systems. To this end, reduced models of homogeneous ideal gas mixtures of methane and air are first obtained by the novel relaxation redistribution method, and thereafter used for the extraction of all the missing variables in a reactive flow simulation with a global reaction model.

## 1. Introduction and motivation

Solution of the full set of equations, as required in numerical simulations of reactive flows with detailed chemical kinetics, represents a quite challenging task even for modern super-computers. On the one hand, the reason is the large number of kinetic equations needed for tracking each chemical species. On the other hand, detailed combustion mechanisms are typical *multi-scale problems* where different chemical processes, characterized by disparate time scales ranging from seconds down to nanoseconds, are present. As a result, modelling detailed combustion fields comes with a tremendous cost where intensive long simulations are needed to resolve the fastest processes, though one is often interested in the slow dynamics. Thus, simplification methodologies become, other than highly desirable, mandatory in combustion problems where detailed mechanisms for heavy hydrocarbons (with hundreds of chemical species) are used in two- and three-dimensional simulations.

Modern simplification techniques are based on a systematic decoupling of the fast chemical processes from the rest of the dynamics, and are typically implemented by seeking a *low-dimensional manifold*.

Much effort has been devoted to setting up such automated model reduction procedures. The method of invariant grids (MIG) [1], the intrinsic low-dimensional manifold [2], the invariant constrained equilibrium edge preimage curve method [3] and the method of minimal entropy production trajectories [4] are a few popular techniques.

In this work, we introduce an approximated procedure for the fast computation of detailed combustion fields using the novel relaxation redistribution method (RRM) [5] for constructing reduced combustion models, and an application to the mechanism, GriMech 3.0 [6], is discussed.

The paper is organized in sections as follows. The notion of slow invariant manifold (SIM) and the chemical kinetics equations are briefly reviewed in §§2 and 3, respectively. The RRM is discussed in §4. A reduced model for air and methane is used in a planar counter-flow flame simulation in §5, and conclusions are drawn in §6.

## 2. Slow invariant manifold

In this section, we briefly discuss the notion of SIM for a system of autonomous ordinary differential equations in a domain in *R*^{n},
2.1
For more details, the interested reader is referred to the dedicated literature [1]. A manifold is *invariant* with respect to the system (2.1) if the inclusion *y*(*t*_{0})∈*Ω* implies that *y*(*t*)∈*Ω* for all future time *t*>*t*_{0}.

Equivalently, if the tangent space *T*_{y} to *Ω* is defined at *y*, invariance requires *f*(*y*)∈*T*_{y}. In order to transform the latter condition into an equation, it proves convenient to introduce projector operators. For any subspace *T*_{y}, let a projector ** P** onto

*T*

_{y}be defined with image im

**=**

*P**T*

_{y}. Then the invariance condition can be expressed by 2.2 where the left-hand side is often called the

*defect of invariance*

*Δ*. It is worth stressing that, although the notion of invariance discussed above is relatively straightforward,

*slowness*instead is much more delicate. We just notice that invariant manifolds are not necessarily suitable for model reduction (all semi-trajectories are, by definition, one-dimensional invariant manifolds). In this respect, we should also point out that

*slow manifolds*are not uniquely defined in the literature and, in general, different methods deliver different objects. Here, we follow the approach of the MIG [1,7], where slowness is intended as

*stability*, thus an SIM is the

*stable stationary*solution of the relaxation process (

*film equation*) 2.3

We notice that the unknown in both equations (2.2) and (2.3) is the manifold *Ω*, which can be conveniently represented in parametric form, as mapping of a domain (reduced space or parameter space in the following) into the phase space , *Ω* being the image of this mapping: .

## 3. Reaction kinetics equations

Here, we consider closed reactive systems with fixed mixture-averaged enthalpy , and total pressure *p*, where *n* chemical species and *d* elements participate in a complex network of elementary reactions. Species compositions are represented in terms of mass fractions *Y*_{i}=*m*_{i}/*m*_{tot}, with *m*_{i} and *m*_{tot} denoting the mass of species *i* and the total mass, respectively. The mixture enthalpy, at a temperature *T*, can be expressed as , while the governing equations in a closed reactor take the form [8]
3.1
where is the mixture density, while *W*_{i}, and *h*_{i}(*T*) denote the molecular weight, molar concentration rate and specific enthalpy of species *i*, respectively. Following ChemKin [9], specific enthalpy can be approximated by a polynomial fit as follows:
3.2
where *a*_{ji} are the tabulated Nasa coefficients and *R* is the universal gas constant. Molar concentration rates take the explicit form
3.3
where *α*_{ij} and *β*_{ij} are the stoichiometric coefficients of the *j*th elementary reaction . The *j*th reaction rate *r*_{j} is expressed using the popular *mass action law*:
3.4
with [*X*_{i}] denoting the molar concentration of the *i*th species and the *j*th reaction rate constants typically expressed in the Arrhenius form [8].

## 4. Model reduction technique

In the following, we use a discrete representation for manifolds referred to as *grids* [1], consisting of a set of *interconnected* nodes (i.e. it is assumed that the nearest neighbours of an arbitrary node *y* can be identified). A grid is defined by the restriction of mapping *F* on the discrete subset of the parameter space into the phase space , whereas an *invariant grid* satisfies the grid version of the invariance equation: *f*(*F*(*ξ*))−*P**f*(*F*(*ξ*))=0, for all [1]. Notice that, owing to the node *connectivity*, it is possible to compute the local tangent space, hence the projector ** P** (e.g. using approximated differential operators).

### (a) Relaxation methods

Construction of one-dimensional invariant grids is accomplished by the RRM, which can be considered an efficient approach for solving the film equation (2.3) starting from an initial grid (for details, see [5]). Referring to figure 1, for simplicity, in this work is chosen regular in terms of the parameter *ξ*.

Let a numerical scheme (Euler, Runge–Kutta, etc.) be chosen for solving the system of kinetic equations (3.1), and let all the grid nodes relax towards the SIM under the detailed dynamics *f* during one time step. The fast component of *f* brings a grid node closer to the SIM, while, at the same time, the slow component causes a contraction towards the steady state of equation (3.1). As a result, the grid becomes dense in a neighbourhood of the steady state and coarse far from it, when keeping relaxing. The slow motion can be neutralized by a node redistribution after the grid relaxation (thus mimicking the term −*P**f* in equation (2.3)). In other words, as illustrated in figure 1, the relaxed states are redistributed on a regular grid in terms of the parameter *ξ* via linear interpolation.

Notice that all intermediate grids are by construction regular in terms of *ξ* and, for invariant grids, the overall effect due to relaxation and redistribution is null (and the invariance condition (2.2) is satisfied).

In our study, we consider a detailed combustion mechanism for air and methane (GriMech 3.0), where *n*=53 chemical species and four elements participate in 325 elementary reactions [6].

Here, at a fixed mixture-averaged enthalpy and pressure *p*, is the *mixing line* between the two states *y*^{fresh} (stoichiometric fresh mixture) and *y*^{eq} (stoichiometric chemical equilibrium state) discretized by *N*=200 nodes (figure 2). Iterations are carried out until at every grid node, where *δy* is the overall movement due to the relaxation and redistribution while is the movement due to relaxation alone of an arbitrary node *y* with a tolerance *ε*=0.001 [5]. To the end of constructing a two-dimensional invariant grid parametrized with respect to *ξ*=*Y*_{CO2} and , the above construction is performed over a range of enthalpies (kJ kg^{−1}) −415 with a step (kJ kg^{−1}) . In figure 3, a projection of the above invariant grid in the three-dimensional subspace –*Y*_{CO2}–*Y*_{CO} is reported.

## 5. Reactive flow simulation

Let us consider the planar stagnation point flow, where a premixed stoichiometric mixture of air and fuel, initially at room conditions (*T*=300 K, *p*=1 bar), impinges against a stream of hot products. Owing to symmetry, a flat flame can be established in this flow at 0<*y*<*L* as schematically depicted in figure 4. Although the above is effectively a two-dimensional problem, under the assumptions of symmetry, boundary layer approximation and low Mach number regime, it is possible to consider the following one-dimensional system of governing equations imposing conservation of mass, momentum, energy and chemical species, respectively, along the symmetry line (*x*=0).
5.1
5.2
5.3
5.4
where the ideal gas law, , can be used for the closure. Moreover, *μ*, *λ*, *D*_{i} and are the dynamic viscosity, thermal conductivity, diffusion coefficient of species *i* and density of the fresh mixture at the inlet, respectively. The mean specific heat (under constant pressure) and the mean molecular weight take the explicit forms:

with *c*_{pi} being the specific heat of species *i* (mass unit). Let *u*, *v* and *ϵ* be the two velocity components of two-dimensional flow field along the *x*- and *y*-axes and the flame strain rate, respectively. We note that the above set of equations (5.1)–(5.4) are conveniently expressed in terms of the quantities and , with . The latter problem is solved by imposing the fresh mixture condition at the inlet,
5.5
while chemical equilibrium and zero-flux condition can be chosen at the outlet,
5.6
where is the density of the fully burned mixture. The detailed derivation of the set of equations (5.1)–(5.4) along with the boundary conditions (5.5) and (5.6) can be found in Marzouk [10]. In this study, spatial derivatives in equations (5.1)–(5.4) have been approximated by finite differences (upwind for convective terms while central differences for diffusive terms), and the corresponding ordinary differential equations system has been solved by a numerical stiff solver [11] readily available in Matlab (*ode15s*). Moreover, we first consider *n*=4 reactive chemical species (CH_{4}, CO_{2}, O_{2}, H_{2}O) with an abundant inert (N_{2}) participating in the one-step global oxidation (*s*=1):
5.7
whose rate *r*_{1}, according to Turns [8], can be evaluated as
5.8

### (a) Retrieval of the missing fields

Upon solving the equations (5.1)–(5.4) in combination with the global reaction mechanism (5.7), five chemical fields, the temperature, mixture density and flow velocity are available along the symmetry line (*x*=0) in figure 4. On the one hand, the computational cost of the latter simulation is drastically reduced compared with a reactive simulation with a detailed combustion mechanism, where a much stiffer and larger set of species equations (5.4) must be taken into account. However, on the other hand, global approaches inevitably come with a significant lack of information concerning all intermediate and minor chemical species, which underlie a complex phenomenon such as the one represented by equation (5.7).

In order to fill this gap, here we perform an *a posteriori* retrieval of all the missing variables in the above computation (e.g. minor species such as radicals) via linear interpolation from the table describing the invariant grid evaluated (once and forever) in §4. To this end, the two variables *Y*_{CO2} and can be extracted from the above simulation and used to access any of the coordinates of the invariant grid (see also figure 3). Some of the interpolated variables are reported in figure 5. Here, a comparison has been carried out between retrieved data and detailed solution by Goodwin using Cantera [12]. Major species (e.g. H_{2}O, O_{2}) and the temperature are in good agreement whereas for most of the radicals (e.g. OH, O) deviations up to 35 per cent are observed. However, prediction of minor species with *Y*_{i}<10^{−6} (e.g. CN, HCN) is out of the reach (and beyond the scope) of the present approximated methodology. We note that the above deviations are expected and mainly due to the effect of diffusion, which, in the reduced model, is taken into account only through the diffusion of manifold parameters. Finally, the implementation of the method does not require significant changes of the flow solver (only fewer species equations (5.4) are to be solved), hence the saving in terms of computational time can be directly linked to the ratio between the number of species in detailed and global models. The latter estimate is quite conservative as the computation of rates in a detailed model is more time-consuming.

## 6. Conclusion

In this work, we first demonstrate that a recently introduced model reduction technique (RRM) is suitable for handling a complex multi-scale combustion mechanism for hydrocarbons. Moreover, on the basis of the latter simplified model, we introduce and test an embarrassingly simple method for the computation of minor chemical species upon a reactive flow simulation efficiently performed with a global (not necessarily one-step) reaction. The present study is the first step towards fast computation of detailed combustion fields for heavy hydrocarbons in two- and three-dimensional problems. The adopted technique is only one possible use of the RRM for reactive flow simulations. Other approaches are discussed, for example, in Chiavazzo *et al*. [13].

## Footnotes

One contribution of 25 to a Theme Issue ‘Discrete simulation of fluid dynamics: applications’.

- This journal is © 2011 The Royal Society