Comparison of different DRG-based methods for the skeletal reduction of JP-8 surrogate mechanisms

Combustion and Flame xxx (2013) xxx–xxx

Contents lists available at SciVerse ScienceDirect

Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e

Comparison of different DRG-based methods for the skeletal reduction of JP-8 surrogate mechanisms Luca Tosatto ⇑, Beth Anne V. Bennett, Mitchell D. Smooke Department of Mechanical Engineering, Yale University, 15 Prospect Street, New Haven, CT 06511, USA

a r t i c l e

i n f o

Article history: Received 29 December 2012 Received in revised form 27 February 2013 Accepted 19 March 2013 Available online xxxx Keywords: DRG Directed relation graph Mechanism reduction JP-8 Skeletal mechanism

a b s t r a c t Directed relation graph (DRG) techniques are used to generate small skeletal mechanisms capable of accurately simulating the combustion of a two-component surrogate for JP-8 jet fuel. Within the DRG framework, six different reduction techniques are considered, and the effectiveness of different definitions for the connection weights and error propagation is evaluated. The use of DRG reduction techniques for aided sensitivity analysis (DRGASA) and for on-the-fly reduction (flux-based DRG) is studied in detail. An optimal reduction approach based on the sequential use of DRG, DRGASA, and flux-based DRG is proposed. When global reduction is applied to a detailed mechanism of 234 species and 6997 reactions, the six reduction techniques result in very different skeletal mechanisms, but all of them are essentially equivalent in terms of accuracy and number of retained species (82–92 species). Finally, for two-dimensional coflow flame test problems, on-the-fly DRG techniques are investigated. Error-propagation-based methods are found to extinguish the flame artificially and cannot be used in an on-the-fly implementation. Conversely, normal DRG methods greatly improve the mechanism reduction, and accurate solutions are obtained using about 15% of the detailed mechanism. Ó 2013 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction In the last five years, the sizes of chemical mechanisms used in combustion modeling have grown by orders of magnitude [1], increasing from 20–30 species for methane combustion to more than 1000 in some recent detailed mechanisms [2]. This growth is driven by interest in the combustion of ever more complex fuel mixtures and in the formation of pollutants. However, numerical solutions are still mostly limited to a few simplified problems such as adiabatic ignition, freely propagating flames, and counterflow flames, in which the fluid-mechanical portion of the problem can be simplified to 1D or even 0D models. The simulation of truly multidimensional flames is still a great challenge, and the solution of even moderately turbulent flames with simple kinetic models approaches the limit of the world’s largest computational facilities [3]. A way to bridge the gap between computational fluid dynamics and complex chemistry is to replace the full set of chemical species and reactions by a simpler one. The new smaller model is optimized for a given problem and can thus generate results of the same accuracy at a much smaller computational cost. Some techniques try to eliminate unimportant species, obtaining a so-called skeletal mechanism (e.g., sensitivity analysis [4], directed relation ⇑ Corresponding author. Current address: Department of Aeronautics and Astronautics, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139, USA. E-mail address: [email protected] (L. Tosatto).

graph reduction [5], and path flux analysis [6]). Others operate a mechanism reduction in which different timescales present in the chemical reactions are separated and used to eliminate unnecessary degrees of freedom (e.g., quasi-steady-state assumption [7], computational singular perturbation [8], and intrinsic low dimensional manifold [9]). Okino and co-workers survey different reduction techniques in [10], and a recent review by Lu and Law focuses more explicitly on combustion applications [1]. After the pioneering work of Bendtsen et al. [11] who first analyzed the oxidation pathways in combustion using graph structures (see also the earlier related work by Turányi [12]), the idea of using connectivity structures to produce reduction schemes was explored by many authors. In 2005, the original directed relation graph (DRG) method of Lu and Law [5] was proposed, in which a connectivity structure is built to quantify the direct coupling between species in the mechanism. DRG-based methods are known for not providing the optimal (i.e., smallest) reduced mechanism [13], but their conceptual simplicity and low computational cost (compared, for example, with sensitivity analysis approaches) have made this reduction approach widely popular. Many modifications and improvements have followed the initial work of Lu and Law [5].  Pepiot-Desjardins and Pitsch [14] proposed a different definition for the coupling norm that accounts separately for rates of formation and destruction.

0010-2180/$ - see front matter Ó 2013 The Combustion Institute. Published by Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.combustflame.2013.03.024

Please cite this article in press as: L. Tosatto et al., Combust. Flame (2013), http://dx.doi.org/10.1016/j.combustflame.2013.03.024

2

L. Tosatto et al. / Combustion and Flame xxx (2013) xxx–xxx

 In the same publication [14], Pepiot-Desjardins and Pitsch introduced the concept of error propagation in the graph search (DRGEP).  More recently, Luo et al. [15] used a modified definition of the connection edges to increase the robustness of the method when dealing with mechanisms having many isomer molecules. The overall implementation of DRG has also been investigated in the literature. Lu and Law [16] found that restarting the reduction procedure can be beneficial, and Zheng et al. [17] designed the DRGASA algorithm, which uses DRG to aid sensitivity analysis reduction. DRG-based methods are automatic (i.e., they require minimal user input) and very computationally cheap; for this reason they are well suited for on-the-fly implementation. Liang et al. [18] used DRGEP reduction to obtain time-dependent chemical mechanisms adaptively. Recently, their approach has been extended by our research group using a flux-based DRG approach that explicitly considers the contribution of transport to the directed graph and allows both time- and space-dependent on-the-fly reduction [19]. An approach that is somewhat similar to DRG methods is the so-called path flux analysis (PFA) method [6], in which the mass transfer between species is used to assemble a graph structure for the skeletal reduction. Alternatively, the error minimization method of Nagy and Turányi [13] uses the off-diagonal entries in the Jacobian matrix to cluster the species to eliminate. Both PFA and error minimization have been demonstrated to produce slightly smaller mechanisms but require more computationally intensive calculations. In any case, these methods are not considered in the present study because they do not belong to the DRG family of methods in a strictly mathematical sense, i.e., they are not based on the definition of a local coupling norm and the study of its propagation across a directed graph (as explained in Section 2). In this paper, we focus our attention on a particular set of methods for skeletal reduction that are based on the DRG representation of the coupling between chemical species. We provide a unified mathematical formulation for different DRG methodologies and present a careful comparison among the different methods. More specifically, we consider three different definitions of the graph structure, in conjunction with both standard DRG and error propagation (DRGEP) reduction. We also consider aided sensitivity analysis (DRGASA) as an effective way to shrink further the mechanism size for any method of the DRG family. This paper is divided into five sections. In Section 2, different DRG and DRGEP methods are defined, and Section 3 discusses on-the-fly implementations. In Section 4, different reduced mechanisms for the combustion of JP-8 are generated via a sequential reduction approach and the results are examined, as follows.  In Section 4.1, the original DRG reduction method is considered.  In Section 4.2, further reduction is obtained using DRG to aid sensitivity analysis (the so-called DRGASA method).  In Section 4.3, the skeletal mechanisms are verified on the solution of two-dimensional (2D) axisymmetric laminar coflow JP-8 flames.  In Section 4.4, further numerical speedups obtained by on-thefly reduction techniques are analyzed. Finally, Section 5 presents conclusions.

(1) Define a norm (or connection weight) to quantify direct coupling between pairs of unknowns. (2) Use the norm to define a directed relation graph structure that connects the unknowns. In the graph, every chemical species maps to a vertex, and an edge is present between two species if and only if direct coupling exists (i.e., if the species have a reaction in common). (3) Define a set of target species that are necessary for the accurate solution of the problem. A typical choice for combustion problems is the H radical; alternatively, a combination of fuel molecules, oxygen, and combustion products can be used. Pollutants and secondary products should be added to the list of target species if their accurate prediction is of interest. (4) Starting from the target species, search the directed graph and attribute to each of the vertices an importance coefficient. This coefficient quantifies how strongly a given species is connected to the target species. (5) Eliminate from the mechanism any species whose importance coefficient is below a user-defined threshold, because such species are only loosely connected to the main combustion pathway. DRG-based methods differ in their definitions of connection weights and importance coefficients. 2.1. Definition of connection weights The original DRG method of Lu and Law [5] defines the connection weight from species i to species j as a ratio between the chemical activity of the couple i,j compared to the total chemical activity of species j, namely ðLuÞ Ri!j

P

a2Cði;jÞ jmia r a j ¼P ; a2RðiÞ jmia r a j

where RðiÞ is the set of reactions that pertain to species i; Cði; jÞ ¼ RðiÞ \ RðjÞ is the set of reactions in which both species i and j participate, mia is the stoichiometric coefficient of species i in reaction a, and ra is the net reaction rate (forward minus reverse). Pepiot-Desjardins and Pitsch [14] noticed that Eq. (1) does not distinguish between reactions that create or destroy species i, which they argued could result in a low extent of reduction. They suggest a possible alternative definition: ðPep:Þ

Ri!j

P     a2Cði;jÞ mia r a  P ; ¼ P þ  max a2RðiÞ ðmia r a Þ ; a2RðiÞ ðmia r a Þ

ðLuoÞ

All DRG-based reduction methods can be summarized in five main steps.

ð2Þ

where the operator ()+ selects only the positive terms in the summation and the operator () selects only the negative terms and makes them positive. Eq. (2) quantifies the contribution of the species couple i, j to the total rate of formation or destruction of species i. Note that all forward and backward rates must be considered as a single reaction when using the connection weights (2), or else partial equilibrium reactions could result in artificially low connection weights. Finally, in recent work, Luo et al. [15] noticed that the original formulation (1) can be improved if the detailed mechanism contains many isomers. For this reason, a variation of the method is introduced, in which the maximum norm is used instead:

Ri!j ¼ 2. DRG methods

ð1Þ

maxa2Cði;jÞ jmia r a j : maxa2RðiÞ jmia r a j

ð3Þ

Eqs. (1)–(3) enable the definition of a directed graph structure in which every species maps to a vertex of the graph and a directed edge (connection) is present only if its R coefficient is not null. A

Please cite this article in press as: L. Tosatto et al., Combust. Flame (2013), http://dx.doi.org/10.1016/j.combustflame.2013.03.024

L. Tosatto et al. / Combustion and Flame xxx (2013) xxx–xxx

3

This iterative procedure is also known as Dijkstra’s graphsearch algorithm, which is attractive for its simplicity. Other equivalent graph-search algorithms are possible in the context of DRG reduction [5,20]. As an illustrative example, the DRG importance coefficients are reported in black in Fig. 1.

Fig. 1. Schematic example of a directed relation graph. The R connection weights are shown in blue, the DRG importance coefficients relative to the target species A are shown in black, and below them the DRGEP coefficients are shown in red. Note that the right portion of the graph is assigned a zero value because it is disconnected from A (E connects to A but not vice versa). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

2.2.2. DRGEP importance coefficients The definition (4) completely neglects the ‘‘distance’’ of a species from a target. For example, in Fig. 1, species C is assigned the same importance as species B despite having a longer chain of connections necessary to reach target species A. In an attempt to include a penalty for species with a long connecting chain, Pepiot-Desjardins and Pitsch [14] formulated the DRGEP method (directed relation graph with error propagation), in which I(DRGEP) is obtained by multiplying all the connections in a path: ðDRGEPÞ Ii

simple sketch of a DRG structure is shown in Fig. 1. Each of the connection weights quantifies the direct error that would be introduced in the solution of species i if species j were to be eliminated from the mechanism. 2.2. Definition of importance coefficients The pair-wise connection weights are not sufficient to quantify the importance of each individual species; in fact, importance coefficients must be considered to account for the error induced by chemical species that are not directly coupled but are connected by a chain of dependence on the graph. For example, referring to Fig. 1, if species A is a target species that needs to be solved accurately, the elimination of species D could introduce a significant error on species B, which would in turn affect A. The importance coefficients take into account connection weights and position in the graph; thus, they estimate the error induced on a target when a chemical compound is removed from the mechanism. 2.2.1. DRG importance coefficients The original DRG method of Lu and Law [5] defines the importance coefficient of species i as1 ðDRGÞ Ii

8 < 1 if species i is a target species    ¼ ðDRGÞ otherwise; : max min Rj!i ; Ij

ð4Þ

j2S

8 < 1 if species i is a target species   ¼ ðDRGEPÞ otherwise: : max Rj!i  Ij

Figure 1 shows that the I(DRGEP) coefficients tend to become smaller when moving away from the target species, contrary to the I(DRG) coefficients. Thus DRGEP is a more ‘‘aggressive’’ technique than DRG, meaning that it tends to assign a low value of I to a larger number of chemical compounds, thus making them more prone to elimination. Each definition of the connection weights (1)–(3) can be used in combination with the standard DRG (4) or with DRGEP (5). Therefore, six different possible reduction methods will be investigated. 2.3. Elimination procedure Any DRG-based reduction approach is based on the solution of test problems that are relevant to combustion applications but that can be easily solved using a detailed mechanism. For example, later we will consider a reduced mechanism generated from stirred reactor and adiabatic ignition calculations for different values of pressure and equivalence ratio. Different directed graph structures are assembled by sampling the test problems for different values of the parameters. For example, for an ignition problem, different DRG structures must be built for different values of time, pressure, and equivalence ratio. A global value for the importance coefficient can be obtained taking the maximum across all the samples

Iglobal;i ¼ max ðIi Þ; samples

where S is the full set of chemical species, Rj?i can be any of the connection weights defined in Eqs. (1)–(3), and it is implicitly assumed that if two species are not connected then Rj?i = 0. Eq. (4) defines the importance coefficient for species i as the smallest connection on any path towards a target species. This problem is well known in computer science and can be solved efficiently using a minimum-cost graph search algorithm [20]. More specifically, I(DRG) is calculated iteratively, as follows.  Starting from the target species, the nearest neighbors in the graph can be assigned a temporary value for I(DRG) = Rj?i.  The nearest neighbors of the nearest neighbors are considered, and a new group of I(DRG) values can be calculated (eventually overwriting the previously calculated values).  This process of updating of the values of I is propagated through the graph until all the connections are explored. 1 In the original work of Lu and Law [5], there is no explicit notion of importance coefficients. However, the definition formulated here is strictly equivalent to their method and enables the DRG and DRGEP methods to be presented in a unified mathematical formulation.

ð5Þ

j2S

ð6Þ

where I signifies either I(DRG) or I(DRGEP), depending on whether Eq. (4) or (5) is used. It must be noticed that Eq. (6) is a very conservative approach to the problem in which a species can be assigned a large importance even though it is effectively important for only a single test case. Thus, this approach can result in mechanisms that are not smallest in size if only some (but not all) major system parameters are involved in the test problems. However, this approach is very common in mechanism reduction [13,15,21], because it can ensure that the accuracy requirements are satisfied for all the test problems. Once the global importance coefficients are defined, a simple reduction approach is to eliminate all the species whose importance is smaller than a user-defined threshold e. This elimination procedure was suggested in the first graph reduction paper of Lu and Law [5] and will be indicated here with the generic term of DRG(EP). The shorthand DRG(EP) refers to DRG when Eq. (4) is used, and it refers to DRGEP when Eq. (5) is used. It should be noted that the simple DRG(EP) approach estimates the error associated with the elimination of a single species but not

Please cite this article in press as: L. Tosatto et al., Combust. Flame (2013), http://dx.doi.org/10.1016/j.combustflame.2013.03.024

4

L. Tosatto et al. / Combustion and Flame xxx (2013) xxx–xxx

of clusters of species; thus, if many species are eliminated, the importance coefficients I can become rough approximations of the elimination error. For example, it can happen that the elimination of a species completely blocks an unimportant chemical path and, as a consequence, it enhances the activity of other chemical reactions, thus increasing the importance of some other species in the graph. In principle, these effects are small, because the importance coefficients I are explicitly designed to identify unnecessary species whose elimination does not affect the overall combustion chemistry. However, if a large enough number of species is eliminated, these secondary effects can accumulate and introduce significant inefficiencies in the elimination procedure. This phenomenon is well explained by Lu et al. [16], who also introduced a very effective correction through the so-called multi-stage DRG(EP) approach, in which after a first elimination step another DRG(EP) procedure is executed. The rationale behind the restart approach is that, after the first elimination, the graph structure is modified and some moderately important species may disconnect from the main oxidation pathways; the second DRG(EP) reduction can identify these species and safely eliminate them from the mechanism. An approach to reduce further the skeletal mechanisms obtained by DRG(EP) methods is to use sensitivity analysis, as implemented in DRGASA [17] and DRGEPSA [22]; we will refer to these two methods collectively as DRG(EP) ASA. In these methods, all species whose importance coefficients are less than the user-defined threshold e are tested one by one, starting from the smallest value of Ii. If the species elimination yields a small error in the test problem simulations, then the chemical compound is removed from the mechanism; otherwise, it is kept regardless of its importance value. This approach has two major advantages: it relates the elimination of a species to the actual simulation error, and it allows larger values of e, which generate smaller skeletal mechanisms. However, improved performance comes at a much higher computational cost, since many numerical simulations are necessary to control the error. It makes sense to combine the different techniques to obtain a small skeletal mechanism in a short computational time. For example, Luo et al. [15] suggested the sequential application of a DRG(EP) reduction and a DRG(EP) ASA reduction. We employ a similar strategy in the present work.

3. On-the-fly flux-based reduction DRG(EP) methods are well suited for on-the-fly reduction, in which the skeletal model is generated automatically during the solution of the problem without the need to resort to test problems. These methods do not allow direct control over simulation error (like in the DRGASA method) but can generate skeletal mechanisms that are tailored to a specific problem, thus offering very high reduction efficiencies. We have applied the six combinations of connection weights (1)–(3) and importance coefficients (4) and (5) using an on-thefly flux-based DRG implementation. This method is particularly useful when solving multidimensional problems, because it operates in a spatially inhomogeneous fashion and generates different chemical mechanisms at different spatial locations. A complete description of this method appears in Tosatto et al. [19]; we briefly recall some basic concepts here. The on-the-fly flux-based DRG method maps each chemical species at each point of a computational grid to a different vertex in the directed graph structure. Within each grid point, different chemical species are connected by the chemical reactions; moreover, graph edges that connect the same chemical species across different grid points are also present, to account for mass transport

effects. For use within the on-the-fly flux-based DRG framework, the graph structure presented earlier must be extended by introducing two different connection weights: RI, i?j, which represents the connection from species i to species j at grid point I, and RI\J,i, which quantifies the connection due to the transport of species i across the interface between points I and J. The R(Lu) connections can be extended to account for transport in the following way [19]: ðLuÞ

RI;i!j ¼

X

jmia r a j=DI;i

a2Cði;jÞ

    ðLuÞ RI\J;i ¼ ðW i V I Þ1 UI\J i =DI;i

!  X  UI\J  DI;i ¼ max jmia ra j;  i  ; W i V I  a2RðiÞ J2N ðIÞ

ð7Þ

X

where Wi is the molecular weight of species i, VI is the volume associated with grid point I; UI\J is the mass flux at the interface bei tween points I and J (including both convection and diffusion), and N ðIÞ is the set of grid points that are close to I. Despite the more complex notation, Eq. (7) is analogous to Eq. (1). The first line quantifies the connection between two species as the ratio between the chemical activity of the couple i, j with respect to the total transport or reaction of species i—whichever is bigger. In the second line, the connection of species i between two neighboring cells in the computational domain is defined as the total mass flux of i divided by the chemical or transport activity. In other words, the on-the-fly flux-based approach treats the mass flux across the cells of a finite volume mesh as if they were ‘‘reactions’’ that, instead of transferring mass between two chemical species, exchange mass between computational cells. This approach is qualitatively analogous to what was done by Valorani et al. [23] and Goussis et al. [24] when studying time scale separation in reactive systems (see also [25]). Following this line of reasoning, the R(Pep.) connections can be generalized according to the following rules:

,   X  ðPep:Þ RI;i!j ¼  mia ra  DI;i  a2Cði;jÞ  . ðPep:Þ 1  I\J  RI\J;i ¼ ðW i V I Þ Ui  DI;i

ð8Þ

 !  X X  UI\J   i DI;i ¼ max jmia ra j ;  ;  W i V I  a2RðiÞ J2N ðIÞ P  P P where j  j ¼ maxð j  jþ ; j  j Þ is an operator that selects the positive or negative terms, similar to what is done in Eq. (2). Finally, the on-the-fly flux-based generalization of the R(Luo) coefficients of Eq. (3) is ðLuoÞ

RI;i!j ¼ max jmia ra j=DI;i a2Cði;jÞ  .   ðLuoÞ RI\J;i ¼ ðW i V I Þ1 UI\J i  DI;i

!   UI\J   i  DI;i ¼ max maxjmia ra j; max  : a2RðiÞ J2N ðIÞW i V I 

ð9Þ

In the present research, the flux-based graph definitions (7)–(9) have been employed in on-the-fly reduction using a standard DRG(EP) procedure nested inside a fully implicit, reactive flow solver [26]. At the start of each numerical iteration, a flux-based DRG(EP) reduction is performed to reduce the number of unknowns to be calculated. Overall, we implement a three-step reduction technique; first the two global reduction techniques, DRG(EP) and DRG(EP) ASA are used, and then the reduced mechanism is used

Please cite this article in press as: L. Tosatto et al., Combust. Flame (2013), http://dx.doi.org/10.1016/j.combustflame.2013.03.024

5

L. Tosatto et al. / Combustion and Flame xxx (2013) xxx–xxx

as an input for an on-the-fly reduced simulation using flux-based DRG. More detail appears in the next section. 4. Results In order to compare the reduction methods defined in the previous sections, we now apply each to a detailed mechanism for a two-component surrogate for JP-8 jet fuel (Aachen surrogate [27], whose molar composition is 77% decane, 23% trimethylbenzene), with the goal of generating skeletal mechanisms for the simulation of coflow diffusion flames of JP-8 surrogate. The detailed mechanism comes from Ranzi et al. [28] and is composed of 234 species and 6997 reactions. This kinetic scheme [28] was built hierarchically, starting from detailed sub-mechanisms for H2 and C1–C4 species. Reaction subsystems that describe the oxidation of higher hydrocarbons were added using a modular structure in which the chemical reactions are divided into homologous classes, so that only a few main classes of fundamental kinetic parameters were needed to extend the scheme to heavier species. This modular approach provided an organized mechanism in which lumping procedures could be utilized to reduce the dimension of the scheme. The overall resulting scheme [28] is quite compact (when compared to detailed mechanisms that deal with similar classes of molecules [29,30]) and has few isomer molecules that may require a specific treatment during the reduction phase [21]. For these reasons, it provides an ideal benchmark on which to run parametric studies on reduction methods, because it makes it possible to approach complex oxidation problems in an affordable computational time. We define two test problems: adiabatic ignition with an initial temperature of 1000 K, and a perfectly stirred reactor (PSR) with an inlet temperature of 500 K. The PSR tests used a continuation method that can resolve the unstable branch below the extinction point [31]. A finely spaced mesh for equivalence ratios and pressures was used to cover all the chemical states normally found in a flame; in particular, we considered pressures in the set P = {1, 5, 10, 15, 20, 25, 30, 35, 40} and ten equally spaced values of equivalence ratio in the range / 2 [0.5, 2]. A relative tolerance of 5% was imposed on ignition time, ignition temperature, extinction residence time, and extinction temperature. Any skeletal mechanism that violated this threshold was considered inaccurate and was rejected. We observe here that this choice of test problems is not the only possible one; for example, freely propagating premixed flames or even counterflow flames could be considered [25]. In this study, we followed closely what was proposed by

20

4.1. DRG and DRGEP reductions First, a parametric study of DRG and DRGEP reductions was conducted using the different connection weights (1)–(3) for various values of e. Results are shown in Fig. 2, which reports the maximum relative error in the test cases as a function of the number of species in the mechanism. All methods give similar results, and they all produce skeletal mechanisms of similar size, as shown in the ‘‘DRG(EP)’’ column of Table 1. (The rest of the table is discussed in later sections.) As the number of species decreases (moving from right to left on the plots in Fig. 2), none of the methods generates strictly growing error curves; for example, the DRG method with Pepiot-Desjardins weights has a region (between 122 and 105 species) in which the simulation error decreases significantly while removing species from the mechanism. This result shows that for all methods, the importance coefficient Ii is only a coarse representation of the actual elimination error; in fact, ideally, an optimal reduction method would always produce increasing simulation errors for decreasing mechanism size. The nonmonotonicity of the error curve is a natural feature not only of DRG but of most reduction methods that are based on the inspection of reactive states. The nonlinearity of the combustion problem makes it impossible to estimate the error generated by species elimination within acceptable bounds. Formally, it is not impossible to impose strictly monotonic error curves using adaptive elimination techniques, such as the methods proposed by Nagy and Turányi [13]. However, in such cases, strictly growing error curves are obtained by a ‘‘trial and error’’ approach, in which different skeletal mechanisms are rejected or retained depending on their simulation error; thus, a very significant computational overhead is introduced for large mechanisms. There are several reasons that justify nonmonotonic curves in Fig. 2.

20

DRG Lu DRG Pep. DRG Luo

18 16

16

14

14

12

12

10 8

10 8

6

6

4

4

2

2

0

100

150

200

Number of species

DRGEP Lu DRGEP Pep. DRGEP Luo

18

Error [%]

Error [%]

Luo et al. Luo et al. [15] and thus concentrated on extinction and ignition problems. Furthermore, we decided to focus on 0D problems, because their low computational cost allowed us to perform extensive parametric tests. Preliminary computational experiments for methane combustion were also conducted on a 171-species 1001-reaction mechanism that models C1–C4 combustion (the so-called SERDP mechanism [32]). The results obtained are in qualitative agreement with the JP-8 surrogate case. For this reason, our discussion will focus only on the jet-fuel case.

0

100

150

200

Number of species

Fig. 2. Maximum relative error in the ignition time, ignition temperature, extinction residence time, and extinction temperature as a function of the number of species in the mechanism. DRG methods are shown on the left, and DRGEP methods are on the right. The thick horizontal line indicates an error threshold of 5%.

Please cite this article in press as: L. Tosatto et al., Combust. Flame (2013), http://dx.doi.org/10.1016/j.combustflame.2013.03.024

6

L. Tosatto et al. / Combustion and Flame xxx (2013) xxx–xxx

Table 1 Number of species retained in the reduced mechanism when applying different reduction methods. The leftmost column gives the names of the methods along with the relevant equation numbers; ‘‘Lu’’, ‘‘Pep.’’, and ‘‘Luo’’ indicate the type of connection weights used. The remaining three columns give mechanism size under various circumstances. Column ‘‘DRG(EP)’’ applies DRG and DRGEP methods to stirred reactor and ignition problems; see Section 4.1. Column ‘‘DRG(EP) ASA’’ applies DRGEP and DRGEPSA methods to stirred reactor and ignition problems; see Section 4.2. Column ‘‘On-the-fly DRG(EP)’’ applies DRG and DRGEP methods in an on-the-fly flux-based reduction of the solution of a axisymmetric laminar coflow flame; see Section 4.4. For on-the-fly reduction with DRGEP, the coflow flame was artificially extinguished. Method name

DRG(EP)

DRG(EP) ASA

On-the-fly DRG(EP)

DRG Lu DRG Pep. DRG Luo

Eqs. (1) and (4) Eqs. (2) and (4) Eqs. (3) and (4)

127 137 150

82 84 89

avg = 30.56, max = 80 avg = 30.78, max = 83 avg = 31.39, max = 87

DRGEP Lu DRGEP Pep. DRGEP Luo

Eqs. (1) and (5) Eqs. (2) and (5) Eqs. (3) and (5)

136 133 124

85 85 92

Extinguished Extinguished Extinguished

 The importance coefficients Ii are representative of the magnitude of the elimination error but not of the sign. For example, the elimination of a species could generate an error of +3% in the ignition time and a subsequent elimination could generate a 3% error, thus generating an overall negligible difference with respect to the full mechanism. The DRG(EP) methods cannot compensate for these error cancellation effects.  The DRG graph structure is affected by the elimination of the species. After a few species are eliminated, it is possible that the overall oxidative structure of the mechanism is changed, so elimination errors can occur because the directed graph is not re-evaluated. This fact is consistent with observations; the error curves in Fig. 2 are monotonic if a small number of species is eliminated (mechanism size of 180–230 species), while nonmonotonic trends become more probable if small mechanisms are sought (120–160 species).  Finally, DRG methods are captivating because of their physically sound but simple construction, but they are partly founded on intuition. For example, the definitions of the R connection weights (7)–(9) are all well posed and phenomenologically sound, but there is no theorem that effectively relates any of these formulations to extinction or ignition temperatures.

Table 2 The algorithm used for aided sensitivity analysis reduction. A DRG(EP) reduction step is considered at the beginning, and then a new reduction graph is evaluated at the start of the DRG(EP) ASA reduction. The chemical species are sorted according to their importance in the graph structure and are then tested individually to determine whether they can be eliminated.

When comparing the DRG and DRGEP methods, it is evident that the error propagation is advantageous if a small reduction extent is sought; Fig. 2 shows that DRGEP can provide more accurate mechanisms than DRG if only a handful of species needs to be removed (mechanism size of 150–200 species). However, no evident difference is present among the different methods in the region where the error crosses the 5% line; the error propagation (5) brings significant improvements for the R(Luo) graph weights, but it is actually slightly detrimental to the other methods. Finally, it can be seen that the connection weights defined in Eq. (2) provide better performance in the portion of the graph where the error is very high. This result signifies that R(Pep.) better quantifies the connection between the important species in the mechanism.

The full results for the DRGASA reduction are shown in Fig. 3 and are summarized in the ‘‘DRG(EP) ASA’’ column of Table 1. Starting from the mechanisms obtained using DRG(EP), the use of sensitivity analysis significantly reduced the size of each mechanism to 82–92 species without significantly increasing the simulation error in the test problems. Similar to what was observed for the DRG(EP) tests of Fig. 2, the DRG(EP) ASA tests of Fig. 3 also show that all six reduction methods are rather equivalent in terms of performance. All the investigated methods are effective and generate small skeletal mechanisms, and all of them greatly benefit from the aided sensitivity analysis to retain some of the important species in the sorted list. Indeed, the bottom plots in Fig. 3 show that some species are tested, but not eliminated, even at the very early stages of the DRG(EP) ASA process (100–130 species). Thus, the DRG(EP) ASA procedure shown in Table 2 improves the original skeletal mechanism of DRG(EP) because it corrects for the importance coefficients I being only approximations of the elimination error. It must be noticed that the more proficient reduction of DRG(EP) ASA is achieved by increasing the computational cost of the reduction procedure, since a number of test problem simulations (equal to the number of species in the sorted list) must be performed. Comparing the DRGASA and DRGEPSA methods (left and right halves of Fig. 3, respectively), it appears that error-propagationbased methods require fewer trials to perform a species elimination. This result implies that the error propagation approach of Eq. (5) is more proficient at sorting the species; in other words, I(DRGEP) is a better estimator of the elimination error than I(DRG) in

4.2. Aided sensitivity analysis Using the mechanisms reduced via DRG(EP) (see Table 1) as input, a second reduction step was performed using DRG(EP) ASA. In DRGASA (DRG-aided sensitivity analysis) or DRGEPSA reduction, the directed graph structure is used to sort species by importance. The species with moderate and low importance are then progressively tested, and only the species whose simulation error is below a specified threshold (in our case, 5%) are eliminated. The algorithm for DRG(EP) ASA appears in Table 2. Compared to the ‘‘simple’’ DRG(EP) approach discussed in Section 4.1, DRG(EP) ASA is much more computationally expensive, but it can provide significantly smaller skeletal mechanisms and allows direct control of the simulation error.

! First reduction: DRG(EP) Test the full mechanism. Assemble the directed graph. Calculate importance coefficients. Eliminate species if I < e. ! Second reduction: DRG(EP) ASA Test the DRG(EP)-reduced mechanism. Assemble the new directed graph. Calculate importance coefficients. Sort species by importance. For every species in the sorted list Eliminate the species. Test the reduced mechanism. If the error is large Reintroduce the species. Else Permanently eliminate the species. End End

Please cite this article in press as: L. Tosatto et al., Combust. Flame (2013), http://dx.doi.org/10.1016/j.combustflame.2013.03.024

7

5

4

4

3 2 1 0 80

Number of trials

Error [%]

5

90

100

110

120

Number of species

20

130

DRGASA Lu DRGASA Pep. DRGASA Luo

15 10 5 0 80

90

100

110

120

Number of species

130

140

3 2 1 0 80

140

Number of trials

Error [%]

L. Tosatto et al. / Combustion and Flame xxx (2013) xxx–xxx

90

100

110

120

Number of species

20

130

140

DRGEPSA Lu DRGEPSA Pep. DRGEPSA Luo

15 10 5 0 80

90

100

110

120

Number of species

130

140

Fig. 3. Maximum relative error in the ignition time, ignition temperature, extinction residence time, and extinction temperature as a function of the number of species in the mechanism for the DRGASA reduction (top left plot) and the DRGEPSA reduction (top right plot). The bottom plots show the number of species that had to be tested (i.e., the number of trials made) before finding a chemical compound suitable for elimination.

the context of aided sensitivity analysis. However, it can also be noticed that DRGEPSA methods generate reduced mechanisms that are slightly larger than those produced by DRGASA. This apparent contradiction can be explained in terms of species pairing effects. In a chemical mechanism that consists of many interacting, but competing, oxidation pathways, the elimination of one species may cause an entire pathway to be eliminated. The DRGASA and DRGEPSA procedures are based not only on connection weights that account for the direct coupling between species but also on importance coefficients that estimate the error due to the elimination of a single species; thus, species pairing effects are never accounted for in the elimination procedures. For this reason, the final size of the mechanism can be highly influenced by species whose elimination early in the elimination procedure caused the removal or amplification of important oxidation paths later in the procedure. Further observations about DRG(EP) ASA reduction can be made by looking at the species that are retained in the skeletal mechanisms. Despite their similarity in size, the six mechanisms have only about 40 species in common. Accounting for the fact that the fundamental oxidation pathways for hydrogen and C1–C2 molecules cannot be eliminated (about 30–35 species), we must conclude that the six mechanisms are actually very different from each other. This finding is quite important. It implies that in the DRG approach, there is a range of equivalent solutions to the reduction problem for the set of test problems and parameters considered here (ignition time, ignition temperature, extinction residence time, and extinction temperature). The existence of these equivalent solutions is probably due to one or both of the following reasons. (1) The chosen test cases are not sufficiently ‘‘strict’’ to identify a skeletal mechanism unequivocally. This hypothesis may signify that DRG could act as an efficient selector for different oxidation pathways, but that the different oxidation pathways have equal sensitivity to the test problems. (2) The DRG methods are all equally effective. This hypothesis would signify that all DRG-based reduction can provide efficient skeletal mechanisms with small computational effort, but the obtained mechanisms preserve different oxidation pathways.

The generation of equivalent solutions to the reduction problem is actually not surprising, because the different methods can operate very different species selection. For example, comparing the DRG and DRGEP formulations in Eqs. (4) and (5), very different importance coefficients are attributed to species that are far away from the target species. DRGEP retains only the important reactions for the species that are close to the starting species and tends to eliminate species that are far from the starting species in the graph, regardless of their connection strength. Therefore, the reduced mechanisms obtained with DRG and DRGEP can retain different oxidation pathways and yet still correctly satisfy the test problems.

4.3. 2D calculations Before considering on-the-fly reduction techniques, we now apply the skeletal mechanisms obtained using DRGASA and DRGEPSA to the solution of 2D coflow flames. These simulations are meant to verify that the ignition and extinction test problems used in the reduction are well suited for more complex combustion problems. Coflow, axisymmetric, laminar diffusion flames are considered in this study. Specifically, we model a burner with two concentric feeds of 0.241 cm diameter for the fuel and 2.54 cm diameter for the oxidizer coflow. The fuel and oxidizer streams flow at average speeds of 79.0 cm/s and 68.7 cm/s, respectively. A parabolic velocity profile is assumed for the inner (fuel) stream, while a top-hat profile is imposed on the outer (oxidizer) stream. The same fuel used in the ignition and PSR problems (Aachen JP-8 surrogate of composition 77% decane, 23% trimethylbenzene [27]) is considered. Both streams are diluted with nitrogen. The mole fraction of fuel in the inner stream is 0.0367, and the outer stream is oxygen-vitiated air with 31% oxygen. This configuration was chosen because experimental measurements are available. The use of experimental data in the context of validation of mechanism reduction is important because it can provide error bounds for the reduction, meaning that we can consider sufficiently accurate any reduced mechanism that is capable of predicting the flame structure of the detailed mechanism within the error bounds of the experimental measurements. The coflow flame is simulated by simultaneous solution of

Please cite this article in press as: L. Tosatto et al., Combust. Flame (2013), http://dx.doi.org/10.1016/j.combustflame.2013.03.024

8

L. Tosatto et al. / Combustion and Flame xxx (2013) xxx–xxx

Fig. 4. Simulations of 2D coflow flames of JP-8 using the skeletal mechanisms for JP-8 surrogate combustion that were output by DRGASA and DRGEPSA. Each of the six subfigures, consisting of two plots, contains the following data: the experimental measurements (left plot’s left half); the numerical solution using the skeletal mechanism (left plot’s right half); the numerical solution using the full 234-species 6997-reaction mechanism (right plot’s left half); and the numerical solution using the skeletal mechanism again (right plot’s right half). Here, ‘‘skeletal mechanism’’ refers to the output of DRGASA (three subfigures in the left column) or DRGEPSA (three subfigures in the right column). The white triangle in the experimental temperature map covers the sooting region of the flame where the measurement is inaccurate.

 the Navier–Stokes equations, which are recast into their vorticity-velocity form;  one equation for the conservation of energy; and  a convective-diffusive-reactive equation for each of the chemical species. Through this formulation, a numerical solution can be found without resorting to staggered grids [33,34]. Dirichlet conditions are imposed at the lower boundary of the domain to fix velocity, temperature, and mixture composition. Symmetry conditions are imposed at the burner centerline, while homogeneous Neumann conditions are used to model the ‘‘outflow’’ conditions at the top of the domain. The computational domain is of size 7  20 cm2, which is about 20 times the size of the flame (0.7  1 cm2), thus ensuring that the

velocity field in the flame region is not affected by the boundary conditions. Finite difference discretizations over a nonuniform structured grid of size 128  128 are used to approximate the differential operators. More specifically, about 70  90 points are used to discretize the actual flame region while the remaining points are used to model the cold flow region outside the flame. Several grid refinement steps were used so that the final numerical solution is independent of the computational domain’s spatial discretization. The discretized governing equations generate a system of nonlinear equations, which is solved using Newton’s method. Pseudo-timestepping is implemented to ease convergence from an arbitrary starting estimate. The Jacobian matrix is assembled at the beginning of each timestep by numerical differentiation of the residual function, and the resulting linear system is solved using GMRES (m) [35].

Please cite this article in press as: L. Tosatto et al., Combust. Flame (2013), http://dx.doi.org/10.1016/j.combustflame.2013.03.024

L. Tosatto et al. / Combustion and Flame xxx (2013) xxx–xxx

9

Fig. 5. Simulations of 2D coflow flames of JP-8 using on-the-fly reduction after the DRGASA reduction. In the leftmost column, Rayleigh scattering temperature measurements (left half) are compared to the on-the-fly numerical simulations (right half). The white triangle in the experimental temperature map covers the sooting region of the flame where the measurement is inaccurate. In the center column, the temperature obtained with the full 234-species 6997-reaction mechanism (left half) is compared with that obtained using on-the-fly flux-based DRG (right half). In the rightmost column, the temperature obtained using on-the-fly flux-based DRG is again displayed (left half), along with the number of active species at each location in the flame (right half).

Due to the large number of species and reactions considered in the detailed chemical mechanism, parallel solution methods are necessary. The domain has been partitioned among different processes, and both the Jacobian assembly and the linear system solution were performed in parallel using the OpenMPI library for network communication. The overall parallel solution algorithm is explained in [26] and draws on previous work of Ern et al. [36]. Experimental 2D temperature maps are available by laser-induced Rayleigh scattering measurements. Details about the measurement technique used for this dataset are explained in [37], and a more general review of the experimental technique is given by Long [38] and Sutton et al. [39]. The experimental temperature maps shown in Fig. 4 can be considered accurate to within about

6.5% because of uncertainties in the feed composition and measurement noise. Furthermore, a minimal uncertainty exists in the definition of the boundary conditions for the flame, since the plug-flow condition for the oxidizer feed may not be a good representation of the effective flow that is produced using a honeycomb mesh. This uncertainty in the boundary condition can slightly affect the liftoff height of the flame. Figure 4 compares the reduced mechanism results both to the detailed mechanism results and to Rayleigh temperature measurements. It can be observed that the reduced mechanisms lead to minimal but observable differences from the solution with the detailed mechanism. More specifically, the smallest skeletal mechanisms—those obtained using the connection weights R(Lu) (1) and

Please cite this article in press as: L. Tosatto et al., Combust. Flame (2013), http://dx.doi.org/10.1016/j.combustflame.2013.03.024

L. Tosatto et al. / Combustion and Flame xxx (2013) xxx–xxx

4.4. On-the-fly reduction in 2D calculations Further reduction was obtained by applying on-the-fly fluxbased DRG reduction to the mechanisms obtained with DRGASA and DRGEPSA. This reduction approach does not formally reduce the size of the chemical mechanism, but it provides computational speedup by selectively deactivating species, reactions, and diffusive transport in appropriate regions of the computational domain. Each row in Fig. 5 shows the results for the on-the-fly reduced simulations using a different definition for the connection weights, chosen from among Eqs. (1)–(3). For the three methods presented in the figure, there is no appreciable degeneration in the accuracy of the solution; indeed, similar to what is shown in Fig. 4, there is a minimal difference in liftoff height between the skeletal solutions and the full mechanism solution. The comparison between numerical and experimental data is very good, both for the solution using the detailed mechanism and for the solutions using the on-the-fly mechanism reduction (see the first two columns of Figs. 4 and 5). The rightmost column in Fig. 5 shows how on-the-fly flux-based DRG works. Starting from the skeletal mechanisms obtained using the DRGASA procedure explained in Section 4.2 (82–89 species), the flux-based DRG approach further reduces the model by solving most of the reactive flow using only 30–60 species, while a more accurate chemical representation is used only in a narrow spatial region. Furthermore, the cold flow that shrouds the flame (not shown in the figure) is modelled using only a handful of chemical compounds. This kind of spatially inhomogeneous reduction enables an accurate solution of the flame using only a small fraction of the DRGASA reduced mechanism. A quantitative comparison is given in Table 1, where the average number of species used in each computational cell is shown in the rightmost column, along with the maximum number of species used in the ignition region at the base of the flame; it can be seen that all the DRG methods perform equivalently well. It should be mentioned that we did indeed implement all six methods considered thus far (see Table 1) in the context of onthe-fly reduced simulations, but DRGEP-based methods (Eq. (5)) artificially extinguished the flame; therefore, they do not appear in Fig. 5. This peculiar failure of DRGEP-based reduction was investigated, and the culprit was found to lie in the use of Eq. (5) in conjunction with transport connections in the directed graph. More specifically, the effect of transport connections was found to decay very rapidly in densely gridded regions of the flame, which caused important species to be eliminated erroneously. In order to quantify better the performance of the on-the-fly flux-based DRG, the relative difference in liftoff height between the reduced calculation and the detailed mechanism is used to measure the error. (Other metrics, such as relative difference in flame length and the L2-norm of the difference between the calculated and measured temperature profiles, were considered and produced similar results.) Furthermore, we will use the average number of species in the simulation to measure the efficiency of the on-the-fly reduction. Figure 6 shows that on-the-fly flux-based DRG reduction can significantly reduce the number of species used in the calculation while keeping the error essentially unchanged. More specifically, starting from a DRGASA reduction that retains

30

Error [%]

R(Pep.) (2) and the DRG–DRGASA importance coefficients (4)—show an appreciable difference in liftoff height with respect to the detailed mechanism. Conversely, the skeletal mechanisms obtained using the connection weights R( Luo) (3) and the DRGEP–DRGEPSA importance coefficients (5) are slightly larger but result in temperature maps that are indistinguishable from the full mechanism. Regardless of these minor differences, the 2D computations using all six skeletal mechanisms are in very good agreement with the experimental results.

20 10 0

0

0.02

0.04

0.06

Epsilon

0.08

50

Error [%]

10

0.1

DRG Lu DRG Pep. DRG Luo

40 30 20 10 0

0

0.1

0.2

0.3

Fraction of active species

0.4

0.5

Fig. 6. Error in liftoff height for the coflow JP-8 flames as a function of the threshold e (top), and as a function of the average fraction of active species with respect to the detailed 234-species mechanism. The on-the-fly flux-based DRG calculations were started using the DRGASA skeletal mechanisms.

about 40% of the full mechanism, the flux-based DRG rapidly reduces the number of species to less than 15% of the total when averaged over all the grid points in the domain. It is interesting to notice that the reduced mechanism generated using the connection weights in Eq. (3) contains slightly more species, but it also provides more accurate results. Worth noting is the fact that a large quantity of species are eliminated even for a low value of the threshold e, and very little error is introduced by the on-the-fly flux-based DRG. However, if large values of e are imposed, a significant error is introduced, even though the number of species considered is mostly unchanged. This behavior is possibly due to the fact that the mechanism was previously reduced using DRGASA. Finally, it is interesting to observe that there are chemical species that the on-the-fly flux-based DRG eliminates from the entire computational domain. More specifically, light alkanes (propane and butane) are included in some of the skeletal mechanisms but are never selected by the on-thefly flux-based DRG as important species. This observation hints once again that the ignition and extinction test problems we considered, although adequate at generating an accurate and small skeletal mechanism, may not be entirely accurate at representing the flame structure at a deeper physical level.

5. Conclusions Different methodologies based on directed relation graph reduction have been applied to the generation of skeletal mechanisms. Two types of reduction have been investigated: global reduction to produce skeletal mechanisms, and on-the-fly reduction subsequently applied to those skeletal mechanisms. All of the global reduction techniques considered generate mechanisms of similar size: 82–92 species (about 35% of the full mechanism). When on-the-fly reduction is considered, a very significant improvement in reduction efficiency is obtained. In that setting, flux-based DRG methods yield accurate solution of flame structures using on average only about 15% of the detailed mechanism, and using about 30% of the full mechanism in a small region of the domain where the flame is located. While all of the methods investigated are very effective at generating skeletal mechanisms, a more in-depth analysis has shown that the reduced mechanisms produced by the DRG and DRGEP methods are rather different in terms of species selection. We have

Please cite this article in press as: L. Tosatto et al., Combust. Flame (2013), http://dx.doi.org/10.1016/j.combustflame.2013.03.024

L. Tosatto et al. / Combustion and Flame xxx (2013) xxx–xxx

discussed a few reasons for this difference, which can be summarized in the following four points. (1) The ignition and extinction test problems that are normally used are not a perfect replica of the flame conditions. (2) The species pairing effect, in which the elimination of an unnecessary species in the early stages of the reduction can cause the elimination of a necessary oxidation pathway in a later stage, is not accounted for by the DRG and DRGEP approaches. (3) The elimination of one or more species alters the directed graph structure. This undesired effect can be mitigated by graph re-evaluation (‘‘restart’’). (4) The definitions of the connection weights R in Eqs. (7)–(9) and of the importance coefficients I in Eqs. (4) and (5) are based on intuition, and for this reason they select different species for elimination. In the authors’ opinion, these remarkable differences among the skeletal mechanisms raise the (unanswered) questions of whether an optimal mechanism (either in terms of mechanism size or robustness) could be selected and whether there is a possibility to improve the DRG approach in this sense. A final interesting observation relates to the failure of DRGEP methods when used in an on-the-fly flux-based DRG context. This failure is tied to the definition of the importance coefficient as the product of the connections on the path (Eq. (5)). This definition provides more accurate (although slightly larger) skeletal mechanisms during the global reduction phase, but it fails at identifying important species in regions with fine gridding when used in multidimensional calculations. Acknowledgments The authors would like to acknowledge Professor Tianfeng Lu of the University of Connecticut for valuable discussions and suggestions, and Mr. Federico Mella and Professor Marshall Long of Yale University for the Rayleigh temperature measurements. The financial support of the US Air Force Office of Scientific Research under Grant #FA9550-06-1-0018 (Dr. Julian Tishkoff, program manager), the US Department of Energy Office of Basic Energy Sciences under Grant #DE-FG02-88ER13966 (Dr. Wade Sisk, program manager), and the National Science Foundation under Grant #CTS-0328296 is gratefully acknowledged. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.combustflame. 2013.03.024.

11

References [1] T. Lu, C.K. Law, Prog. Energy Combust. Sci. 35 (2009) 192–215. [2] C.K. Westbrook, W.J. Pitz, O. Herbinet, H.J. Curran, E.J. Silke, Combust. Flame 156 (2009) 181–199. [3] D.O. Lignell, J.H. Chen, P.J. Smith, Combust. Flame 155 (1–2) (2008) 316–333. [4] Y. Reuven, M.D. Smooke, H. Rabitz, J. Comput. Phys. 64 (1986) 27–55. [5] T. Lu, C.K. Law, Proc. Combust. Inst. 30 (2005) 1333–1341. [6] W. Sun, Z. Chen, X. Gou, Y. Ju, Combust. Flame 157 (7) (2010) 1298–1307. [7] F.A. Williams, Combustion Theory, Benjamin/Cummings, 1984. pp. 565–575. [8] S.H. Lam, D.A. Goussis, Int. J. Chem. Kinet. 26 (1994) 461–486. [9] U. Maas, S. Pope, Proc. Combust. Inst. 24 (1992) 103–112. [10] M.S. Okino, M.L. Mavrovouniotis, Chem. Rev. 98 (1998) 391–408. [11] A.B. Bendtsen, P. Glarborg, K. Dam-Johansen, Comput. Chem. 25 (2) (2001) 161–170. [12] T. Turányi, New J. Chem. 14 (1990) 795–803. [13] T. Nagy, T. Turányi, Combust. Flame 156 (2009) 417–428. [14] P. Pepiot-Desjardins, H. Pitsch, Combust. Flame 154 (2008) 67–81. [15] Z. Luo, T. Lu, M.J. Maciaszek, S. Som, D.E. Longman, Energy Fuels 24 (2010) 6283–6293. [16] T. Lu, C.K. Law, Combust. Flame 144 (1–2) (2006) 24–36. [17] X. Zheng, T. Lu, C. Law, Proc. Combust. Inst. 31 (1) (2007) 367–375. [18] L. Liang, J.G. Stevens, J.T. Farrell, Proc. Combust. Inst. 32 (1) (2009) 527–534. [19] L. Tosatto, B.A.V. Bennett, M.D. Smooke, Combust. Flame 158 (5) (2011) 820– 835. [20] D.E. Knuth, The Art of Computer Programming, vol. 1, Addison-Wesley, 1968. [21] T. Lu, C.K. Law, Combust. Flame 154 (1–2) (2008) 153–163. [22] K.E. Niemeyer, C.-J. Sung, M.P. Raju, Combust. Flame 157 (9) (2010) 1760– 1770. [23] M. Valorani, H.N. Najm, D.A. Goussis, Combust. Flame 134 (1–2) (2003) 35–53. [24] D.A. Goussis, M. Valorani, F. Creta, H.N. Najm, Prog. Comput. Fluid Dyn. 5 (6) (2005) 316–326. [25] M. Valorani, F. Creta, D.A. Goussis, J.C. Lee, H.N. Najm, Combust. Flame 146 (1– 2) (2006) 29–51. [26] L. Tosatto, B.A.V. Bennett, M.D. Smooke, Combust. Theory Modell. 15 (4) (2011) 455–486. [27] S. Honnet, K. Seshadri, U. Niemann, N. Peters, Proc. Combust. Inst. 32 (1) (2009) 485–492. [28] E. Ranzi, A. Frassoldati, S. Granata, T. Faravelli, Ind. Eng. Chem. Res. 44 (14) (2005) 5170–5183. [29] O. Herbinet, W.J. Pitz, C.K. Westbrook, Detailed chemical kinetic oxidation mechanism for a biodiesel surrogate, in: Fall Technical Meeting of the Western States Section of the Combustion Institute, 2007. [30] C. Westbrook, W. Pitz, J. Boercker, H. Curran, J. Griffiths, C. Mohamed, M. Ribaucour, Proc. Combust. Inst. 29 (1) (2002) 1311–1318. [31] V. Giovangigli, M. Smooke, Appl. Numer. Math. 5 (4) (1989) 305–331. [32] X. You, A.V. Joshi, S.G. Davis, A. Laskin, F. Egolfopoulos, C.K. Law, H. Wang, USC Mech Version II. High-Temperature Combustion Reaction Model of H2/CO/C1– C4 Compounds, ignis.usc.edu/Mechanisms/USC-Mech II/USC_Mech II.htm (2007). [33] T.B. Gatski, C.E. Grosch, M.E. Rose, J. Comput. Phys. 48 (1) (1982) 1–22. [34] A. Ern, M.D. Smooke, J. Comput. Phys. 105 (1) (1993) 58–71. [35] H. van der Vorst, Iterative Krylov Methods for Large Linear Systems, Cambridge University Press, 2003. [36] M.D. Smooke, V. Giovangigli, Int. J. Supercomput. Appl. 5 (4) (1991) 34–49. [37] L. Tosatto, F. Mella, M.B. Long, M.D. Smooke, A study of JP-8 surrogate coflow flame structure by combined use of laser diagnostics and numerical simulation, Combust. Flame. 159 (10) (2012) 3027–3039. [38] M.B. Long, Multidimensional imaging in combusting flows by Lorenz-Mie, Rayleigh, and Raman scattering, in: A.M.K.P. Taylor (Ed.), Instrumentation for Flows with Combustion, Academic Press, London, 1993. [39] G. Sutton, A. Levick, G. Edwards, D. Greenhalgh, Combust. Flame 147 (1–2) (2006) 39–48.

Please cite this article in press as: L. Tosatto et al., Combust. Flame (2013), http://dx.doi.org/10.1016/j.combustflame.2013.03.024