Multiobjective Genetic Algorithm Optimization for Calculating the

University of Leeds, Leeds LS2 9JT, U.K., and Centre for Aerospace Technology, QinetiQ Group PLC,. Cody Technology Park, Ively Road, Pyestock, ...
0 downloads 0 Views 177KB Size
Download PDF
Ind. Eng. Chem. Res. 2003, 42, 1215-1224

1215

Multiobjective Genetic Algorithm Optimization for Calculating the Reaction Rate Coefficients for Hydrogen Combustion Lionel Elliott,† Derek B. Ingham,† Adrian G. Kyne,‡ Nicolae S. Mera,*,‡ Mohamed Pourkashanian,§ and Christopher W. Wilson| Department of Applied Mathematics, Centre for Computational Fluid Dynamics, and Department of Fuel and Energy, School of Process, Environmental and Materials Engineering, University of Leeds, Leeds LS2 9JT, U.K., and Centre for Aerospace Technology, QinetiQ Group PLC, Cody Technology Park, Ively Road, Pyestock, Farnborough, Hants GU14 0LS, U.K.

It is well recognized that many important combustion phenomena are kinetically controlled. Whether it be the burning velocity of a premixed flame, the formation of pollutants in different stages of a combustion process, or the conversion of NO to NO2 in a gas turbine engine exhaust, it is important that a detailed chemical kinetic approach be undertaken in order to fully understand the chemical processes taking place. In this study, a multiobjective genetic algorithm approach is developed for determining new reaction rate parameters (A’s, β’s, and Ea’s in the three-parameter functional form of the Arrhenius expressions) for the combustion of hydrogen/ air mixtures. The multiobjective structure of the genetic algorithm employed allows for the incorporation of both perfectly stirred reactor and laminar premixed flame data in the inversion process, thus producing more efficient reaction mechanisms. Various inversion procedures based on reduced sets of data are investigated and tested on hydrogen/air combustion in order to generate efficient inversion schemes for future investigations concerning complex hydrocarbon fuels. 1. Introduction Combustion may be modeled using a system of chemical reactions, for which the rates of each reaction are known. Databases holding measurements of the reaction rate parameters (A’s, β’s, and Ea’s) are available. These provide the logarithm of the reaction rate as a function of 1/T, where T is the temperature, but with a large variation observed for some reaction temperatures. However, there is a considerable amount of uncertainty in the rate data, which causes problems to occur when calculating the concentration of species for the combustion products away from their equilibrium values. The Gas Research Institute (GRI) initiated an extensive research program for the optimization of the chemical scheme with regard to natural gas. Precise details of the optimization scheme in that program are discussed in detail by Bowman et al.1 At temperatures above about 1000 K, combustion of simple fuels, such as hydrogen, can be calculated using an eight-step reaction mechanism, for which the reaction rate data are known with confidence.2 However, complex hydrocarbon fuels, such as kerosene, require more than 1000 reaction steps with over 200 species.3 Hence, because all of the reaction rate data are not well-known, there is a high degree of uncertainty in the results obtained using these large detailed reaction mechanisms. We note that various methods have been proposed to find a set of reaction rate parameters that give the best fit to a given set of experimental data. Indeed, once the * Corresponding author. E-mail: [email protected]. Tel: 0044-1133432481. Fax: 0044-1132467310. † Department of Applied Mathematics, University of Leeds. ‡ Centre for Computational Fluid Dynamics, University of Leeds. § Department of Fuel and Energy, School of Process, Environmental and Materials Engineering, University of Leeds. | Centre for Aerospace Technology, QinetiQ Group PLC.

problem is reformulated as an optimization problem, the object function may be minimized/maximized by various methods. However, in the case of complex hydrocarbon fuels, the object function is usually highly structured, having multiple ridges and valleys and exhibiting multiple local optima. For object functions with such a complex structure, traditional gradient-based algorithms4 are likely to fail. Optimization methods based upon linearization of the objective function5,6 fall into the same category. A different method, the solution mapping method, was proposed by Frenklach et al.7 and is based on approximating the surface of the object function by polynomial response surfaces generated by means of a relatively small number of computer experiments. Although efficient from the point of view of computational time, the method of solution mapping has a major drawback because the optimization is based not on the real values of the object function but on a polynomial response surface approximating the real landscape of the object function. For a highly structured surface with multiple valleys and hills, the errors in approximating the object function are large, and these errors are likely to affect the results of the optimization procedure. On the other hand, genetic algorithms (GAs) are particularly suitable for optimizing object functions with complex, highly structured landscapes.8 A powerful inversion technique, based on a singleobjective binary-encoded GA, was developed9 to determine the reaction rate parameters for the combustion of a hydrogen/air mixture in a perfectly stirred reactor (PSR). This technique was further developed by Elliott et al.10 by incorporating constraints associated with the experimental data listed in the National Institute of Standards and Technology (NIST) database. The original binary encoding was also replaced by floating-point encoding. An extension to methane/air combustion was also considered.11 In the above-mentioned studies, GAs were used to determine the reaction rate parameters using measure-

10.1021/ie020501o CCC: $25.00 © 2003 American Chemical Society Published on Web 02/25/2003

1216

Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003

ments of the species concentrations in a PSR. However, many practical combustors, such as internal combustion engines, rely on premixed flame propagation. Moreover, burner-stabilized laminar premixed flames are very often used to study chemical kinetics in a combustion environment. Such flames are effectively one-dimensional and can be made very steady, thus facilitating detailed experimental measurements of temperature and species profiles. Therefore, for laminar premixed flames, it is much easier to obtain accurate experimental data to be used in the GA inversion procedure. It is the purpose of this paper to construct such a procedure that takes into account measurements from both laminar premixed flames and PSRs. It was found9 that when using only PSR measurements, the problem of retrieving the reaction rate coefficients has multiple solutions that give equally accurate predictions for species concentrations in PSR simulations. However, further investigations11 have shown that when validated against different experimental data, not all of these solutions are accurate in predicting, for example, the species concentration profiles in laminar premixed flames. Therefore, to obtain a reaction mechanism which can be used for a wide range of practical problems, it is required to develop an inversion procedure which incorporates multiple objectives. Another reason to use this approach is the fact that sometimes in practical situations the amount of experimental data available is limited. Therefore, it is important to use all of the available experimental data, which may come from different measurements for PSRs or laminar premixed flames. Inversion procedures based on a reduced set of available measurements are also considered. In many experimental simulations, some of the species concentrations are subject to large experimental errors both in PSRs and laminar premixed flames. Moreover, in laminar premixed flames, the experimental simulations may provide only the maximum concentration of a particular species rather than the whole profile along the burner. To facilitate the incorporation of various types of experimental data in the GA-based inversion procedures, this study also investigates such GA calculations based on such incomplete or reduced sets of data. 2. Description of the Problem 2.1. Reaction Rate Parameters. In a combustion process, the net chemical production or destruction rate of each species results from a competition between all of the chemical reactions involving that species. In this study, it is assumed that each reaction proceeds according to the law of mass action and the forward rate coefficients are of three-parameter functional Arrhenius form, namely

Rate:

( )

kfi ) AiTβi exp -

Eai

RT

(1)

for i ) 1, ..., NR, where R is the universal gas constant and there are NR competing reactions occurring simultaneously. The rate equation (1) contains the three parameters Ai, βi, and Eai for the ith reaction. These parameters are of paramount importance in modeling the combustion processes because small changes in the reaction rates may produce large deviations in the output species concentrations. In both PSRs and in laminar premixed flames, if the physical parameters required are specified and the reaction rates (1) are

given for every reaction involved in the combustion process, then the species concentrations of the combustion products may be calculated. It is the possibility of the determination of these parameters for each reaction, based upon outlet species concentrations, that is investigated in this paper. Various software packages may be used for the direct calculations to determine the output species concentrations if the reaction rates are known. For these direct calculations, we have used the PSR code12 and the PREMIX code,13 which are described in the next sections. 2.2. PSR and Laminar Premixed Flame (PREMIX) Calculations. The PSR calculations are performed using the PSR FORTRAN computer program that predicts the steady-state temperature and composition of the species in a PSR12 while the laminar premixed flame structure calculations were performed using the PREMIX code.13 We use Xk to denote the mole fraction of the kth species, where k ) 1, ..., K and K represents the total number of species. The temperature and composition that exit the reactor are assumed to be the same as those in the reactor because the mixing in the reactor chamber is intense. differIf PSR calculations are undertaken for NPSR S ent sets of reactor conditions, then the output data that are used in the GA search procedure will consist of a mole concentrations (Xjk, j ) 1, ..., NPSR set of KNPSR S S , k ) 1, ..., K), where Xjk represents the mole concentration of the kth species in the jth set of reactor conditions. The laminar premixed flame structure calculations were performed using the PREMIX code for burnerstabilized flames with a known mass flow rate. The code is capable of solving the energy equation in determining a temperature profile; however, throughout this investigation a fixed temperature profile was preferred. The reason for fixing the temperature was twofold. First, in many flames there can be significant heat losses to the external environment, which are of unknown and questionable origin and thus troublesome to model. Second, the CPU time involved in arriving at the final solution was decreased significantly. If all of the physical operating conditions are specified, then for a given set of reaction rates (1) the profiles of the species concentration and the burning velocity as functions of the distance from the burner’s surface are calculated. If PREMIX calculations are performed for different sets of operating conditions, then the NPSR S output data of the code that is used by the GA in the matching process consist of a set of species concentration profiles [Yjk(x) , j ) 1, ..., NPREMIX , k ) 1, ..., K] and a set S of burning velocity profiles [Vj(x) , j ) 1, ..., NPREMIX ], S where Yjk is the profile of the mole concentration along the burner for the kth species in the jth set of operating conditions and Vj is the burning velocity profile for the jth set of operating conditions. Both PSR and PREMIX programs are not stand-alone programs; they are designed to be run in conjunction with preprocessors from the CHEMKIN library,14 which handles the chemical reaction mechanism and the thermodynamic and transport properties. It is the purpose of this paper to determine the reaction rate coefficients A’s, β’s, and Ea’s in eq 1 given species concentrations Xjk in a PSR and/or profiles of species concentration Yjk and the burning velocity profiles Vj in laminar premixed flames.

Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003 1217

3. Reformulation of the Problem as an Optimization Problem An inverse solution procedure attempts, by calculating new reaction rate parameters that lie between predefined boundaries, to recover the profiles of the species (to within any preassigned experimental uncertainty) resulting from numerous sets of operating conditions. The inversion process aims to determine the unknown reaction rate parameters [(Ai, βi, Eai), i ) 1, ..., NR] that provide the best fit to a set of given data. Thus, first a set of output concentration measurements of species is simulated (or measured). If numerical simulations are performed for NPSR different sets of S different sets PSR operating conditions and/or NPREMIX S of PREMIX operating conditions, then the data will concentration measurements consist of a set of KNPSR S of species for PSR calculations and (K + 1)NPREMIX S profiles of species concentrations and burning velocity for PREMIX calculations. GA-based inversion procedures seek the set of reaction rate parameters that gives the best fit to these measurements. To do so, we consider two objective functions which compare predicted and measured species concentrations for PSR and PREMIX simulations:

fPSR[(Ai,βi,Ei)i)1,NR] )

{

10

-8

NS K

+

∑ ∑ Wk j)1k)1

}

orig |Xcalc jk - Xjk |

Xorig jk

-1

(2)

fPREMIX[(Ai,βi,Ei)i)1,NR] )

{

10

NS

-8

+

∑ j)1

(

|Vcalc - Vorig j j |L2

+

|Vorig j |L2

)}

K

orig |Ycalc jk - Yjk |L2

k)1

|Yorig jk |L2

∑ Wk

-1

(3)

calc calc represent the mole conwhere (i) Xcalc jk , Yjk , and Vj centrations of the kth species in the jth set of operating conditions and the burning velocity in the jth set of operating conditions using the set of reaction rate orig parameters [(Ai, βi, Eai), i ) 1, ..., NR], (ii) Xorig jk , Yjk , orig and Vj are the corresponding original values which were measured or simulated using the exact values of the reaction rate parameters, (iii) | |L2 represents the L2 norm of a function which is numerically calculated using a trapezoidal rule, and (iv) Wk are different weights that can be applied to each species depending on the importance of the species. Initially all of these weights are set to be equal to 1.0. It should be noted that the fitness function (2) is a measure of the accuracy of species concentration predictions obtained by a given reaction mechanism for PSR simulations. The second fitness function is a measure of the accuracy in predicting species concentrations and burning velocity profiles in the case of laminar premixed flames. We note that the maximum value of these object functions is the constant 108, which indicates a perfect fit to the data. The constant has been added to the fitness function to avoid numerical overflow. A high value was used to increase the selective pressure applied by the GA. However, similar results are obtained with various other values.

It is worth noting that when working with accurate experimental data or numerically simulated data the two fitness functions (2) and (3) have a common global optimum given by the real values of the reaction rate coefficients or the values that were used to generate the data, respectively. Therefore, for the problem considered the two objective functions are not conflicting, and the method we have proposed can be used although it searches for a single common optima rather than a set of nondominated solutions. It is worth noting that if the experimental data are corrupted with a large level of noise, then the two objective functions (2) and (3) may be conflicting and then Pareto-oriented techniques15 have to be used in order to construct the front of nondominated solutions, i.e., a set of possible optima. In this situation other considerations have to be used to make a decision on what is the best location of the optimum. 4. Multiobjective GA Inversion Technique There is a vast variety of GA approaches to multiobjective optimization problems. GAs are, in general, a very efficient optimization technique for problems with multiple, often conflicting, object functions.15 However, most GA approaches require that all the object functions are evaluated for every individual generated during the optimization process. For the problem considered in this paper, the object function fPREMIX requires NPREMIX S executions of the PREMIX code, which has a long run time. For example, in the hydrogen case a PREMIX run was found to take an average of 7.2 s CPU time on a Pentium IV processor running at 1.7 GHz. This leads to a running time of more than 5 days for evaluating say 1000 GA generations with 60 individuals generated at every iteration. For complex hydrocarbon fuels, the CPU time required for a single individual is substantially increased and the overall run time required to achieve convergence becomes prohibitive. On the other hand, PSR calculations are much faster. Therefore, in this paper we develop a multiobjective GA that avoids PREMIX evaluations at every generation. The algorithm is based on alternating the fitness functions that are used for selecting the survivors of the populations generated during the GA procedure. Under the assumptions that the optimization problem considered has n objective functions, the multiobjective GA employed in this paper consists of the following steps: 1. Construct n fitness functions f1, f2, ..., fn, and assign a probability for each of these functions p1, p2, ..., pn. 2. Construct an initial population. 3. Select a fitness function fi with the probability pi using a roulette wheel selection. 4. Using the fitness function selected at step 3, construct a new population. (a) Use the fitness function fi to select a mating pool. (b) Apply the genetic operators (crossover and mutation) to create a population of offsprings. (c) Use the fitness function fi to construct a new population by merging the old population with the offspring population. 5. Repeat steps 3 and 4 until a stopping condition is satisfied, i.e., until convergence is achieved. For the problem considered in this chapter, we have n ) 2, f1 ) fPSR, and f2 ) fPREMIX, and three different sets of selection probabilities for the object functions were investigated in the following GA formulations:

1218

Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003

(1) GA_PSR (pPSR ) 1.0, pPREMIX ) 0.0). This formulation is equivalent with the single objective optimization procedure developed in ref 10 and is based only on PSR calculations. (2) GA_PREMIX (pPSR ) 0.0, pPREMIX ) 1.0). This formulation is based only on PREMIX calculations and is computationally expensive. (3) GA_MIXED (pPSR ) 0.9, pPREMIX ) 0.1). This formulation uses both PSR and PREMIX calculations. On average, out of 10 generations, 9 are evaluated using fPSR and 1 is evaluated using fPREMIX. Different values may be set for pPSR and pPREMIX. In terms of computational time, GA_PREMIX is the most time-consuming procedure, requiring an average of 7.2 min/generation, while GA_PSR only requires 0.66 min/generation. The CPU time required by the GA_ MIXED depends on the values set for the probabilities pPSR and pPREMIX. For the values considered in this paper, the GA_MIXED was found to require an average of 1.32 min/generation. This reduction in computational time in comparison with GA_PREMIX will become very important when optimization procedures need to be run for a large number of generations for complex hydrocarbon fuels. It is worth noting that the fitness selection probabilities pPSR and pPREMIX can be used to bias the convergence of the algorithm toward accurate predictions in PSR or PREMIX, respectively. However, they do not act as weights because at every generation only one of the objective functions (2) and (3) is evaluated and not a linear combination of both of them. Instead, these objective functions are alternated during the evolution process. This novel method maintains the multiobjective character of the optimization process while avoiding the evaluation of the expensive fitness function (3) during every generation. We note that the GA employed is an extension for multiobjective optimization of the floating point encoded GA developed by Michalewicz8 and used for reaction rate retrieval.10 The genetic operators and the GA parameters were taken to be as follows: (i) population size npop ) 50; (ii) number of offspring nchild ) 60; (iii) uniform arithmetic crossover, crossover probability pc ) 0.65; (iv) tournament selection, tournament size k ) 2, tournament probability pt ) 0.8; (v) nonuniform mutation, mutation probability pm ) 0.5; (vi) elitism, elitism parameter ne ) 2. We note that the nonuniform mutation operator8 is an operator that enables fine local tuning by searching the space uniformly during the first generations and locally at later stages. This enables a high level of exploration during the first few generations when the GA is searching for promising regions in the parameter space, while at later stages the solutions found are refined by analyzing the effects on the fitness function of small perturbations applied to the reaction rate parameters. The new sets of rate constants generated every iteration are constrained to lie between predefined boundaries that represent the uncertainty associated with the experimental findings listed in the NIST database.10 5. Numerical Results To test the multiobjective GA inversion procedure and to compare various reduced versions of it, the test problem considered is the recovery of the species concentration profiles predicted by a previously vali-

Figure 1. Fitness function fPSR as a function of the number of generations performed for i ) 1000 generations.

dated hydrogen/air reaction mechanism for a wide range of operating conditions. The reaction scheme used to generate the data involves 9 species across 19 reactions.9 Ideally, experimental data would have been preferred for the optimization process; however, because this is difficult to find for all of the conditions at which the combustion scheme should be tested, it was decided to undertake numerical computations on a set of numerically simulated data. Because the purpose of this study is to develop and to assess the technical feasibility of the multiobjective GA techniques for the development of new reaction mechanisms and to compare different formulations in order to employ the most efficient schemes for future investigations concerning complex hydrocarbon fuels, for this study numerically simulated data are as efficient as experimental data. The test problem considered is the combustion of a hydrogen/nitrogen/oxygen mixture at a constant atmo) 20 spheric pressure, p ) 1 atm. A set of NPSR S different operating conditions have been considered for PSR runs corresponding to various changes to the inlet temperature T ) 1000, 1250, 1500, 1750, and 2000 K and equivalence ratios F ) 0.5, 1.0, 1.5, and 2.0. The reactor has a volume V ) 67.4 cm3, and the nominal residence time is τ ) 3.0 × 10-5 s. The inlet composition is defined by a specified fuel/air equivalence ratio, Φ ) 1.0. The fuel is 80% H2 and 20% N2, and the oxidizer is air, i.e., 79% N2 and 21% O2. We assume that there is no heat loss from the reactor, i.e., Q ) 0, and the temperature profile is specified rather than the energy equation being solved. ) 5 operFor PREMIX calculations, a set of NPREMIX S ating conditions, corresponding to various values of the equivalence ratios, namely, Φ ) 0.5, 0.75, 1.0, 1.25, and 1.5 for a steady-state burner-stabilized premixed laminar flames with specified temperature, pressure p ) 1 atm, and mass flow rate ) 4.63 × 10-3 (g/cm2/s). 5.1. Inversion Procedures for Multiobjective Optimization and Comparison with Single-Objective GA Procedures. It is the purpose of this section to investigate the convergence and efficiency of the multiobjective GA proposed, GA_MIXED, and to compare it with the single-objective GAs GA_PSR and GA_PREMIX. Figure 1 presents the values of the fPSR

Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003 1219

Figure 2. Reaction rate k for the reaction H + O2 ) O + OH given by the inversion procedures employed in various GA formulations, namely, in comparison with the original reaction rate and the NIST constraints imposed on the reaction rate.

fitness function for the fittest individual in the population as a function of the number of generations. We note that, although an elitist model is employed, the fitness function value is not monotonic with respect to the number of generations performed. This is due to the fact that at about every 10 generations the fitness function used to select the new population is changed and the selection is done using the fPREMIX fitness function. If this is the case, then the value of the fPSR value of the fittest individual in the population may decrease if that individual does not give a high fitness value for the fPREMIX fitness function and, hence, it does not necessarily propagate to the next generation. This technique eliminates those individuals that artificially give a high value for one of the fitness functions but fail to give a high value for the other. It was found that such individuals are numerically constructed pseudosolutions that give reaction rates different from the original one yet still provide good fitness to the finite set of data used in the optimization process. However, they fail to produce accurate results for a wider range of problems, under different operating conditions. The proposed multiobjective GA constructs and evolves populations based on an object function, and then the solutions generated are validated against another object function. This filters out individuals that generate reaction mechanisms with restricted ranges of operation and evolves populations of solutions that will give accurate predictions for various types of problems under various operating conditions. We also note in Figure 1 that after about 800-900 generations the fitness function does not increase substantially, suggesting that convergence has been achieved. Although not presented here, it is reported that the fitness function fPREMIX has a similar behavior. Therefore, each of the three formulations considered was run for 1000 generations to produce the reaction mechanisms GA_PSR, GA_MIXED, and GA_PREMIX, which are compared next. Figure 2 presents the reaction rate retrieved for the first reaction, namely, H + O2 ) O + OH, by the three

Figure 3. Mole fraction of O (a) and H2O2 (b) for various temperatures, as calculated by PSR using the reaction mechanisms generated by the inverse GA formulations considered.

inversion formulations considered. It can be seen that all three formulations are efficient in retrieving the reaction rate because the graphs are almost indistinguishable from the original reaction rate. However, a more efficient way to compare the reaction mechanisms generated is by analyzing the predictions that they provide for the various chemical processes because, for example, very small variations in the reaction rates may produce large deviations in the predicted species concentrations. Figure 3 presents the species concentration calculated with the PSR program, using the reaction mechanisms generated in the GA inversion procedure and the original mechanism for various temperatures and the air/fuel ratio given by Φ ) 1.0. It can be seen that GA_PSR outperforms GA_PREMIX in predicting the PSR output species. This can obviously be explained by the fact that GA_PSR is optimized using PSR species measurements, while GA_PREMIX is optimized for laminar premixed flame predictions. The errors of

1220

Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003

Figure 4. Profiles of the mole fractions along the burner, for various species, as calculated by PREMIX using the reaction mechanisms generated within various inversion formulations, namely, GA_PSR (9), GA_MIXED (×), and GA_PREMIX (2), in comparison with the profiles obtained with the original mechanism (s).

GA_PREMIX predictions of PSR output species are large, particularly for the unstable species; see, for example, H2O2 in Figure 3. However, GA_MIXED, the mechanism generated on both PSR and PREMIX simulations, is very efficient in predicting all of the PSR output species. If the profiles of the mole fractions in the PREMIX calculations are investigated, then as expected GA_ PREMIX outperforms GA_PSR; see Figure 4. However, again GA_MIXED gives accurate predictions for both stable and unstable species, while GA_PSR only gives accurate results for the stable species. We note that GA_MIXED equals GA_PREMIX in species profile predictions, although it is obtained by a much faster inversion procedure, as mentioned in section 4. Figures 3 and 4 presented the predictions given by the constructed reaction mechanism for operating conditions which were included in the optimization process.

Figure 5. Burning velocity and main species mole fraction profiles for PREMXI calculations at an air/fuel ratio Φ ) 2.0 for the reaction mechanism generated by various GA inversion procedures, namely, GA_PSR (9), GA_MIXED (×), and GA_PREMIX (2), in comparison with the profiles obtained with the original mechanism (s).

However, the aim of the inversion procedure is to develop reaction mechanisms that are valid over wide ranges of operating conditions. Therefore, next the reaction mechanisms are tested at operating conditions outside the range included in the optimization procedure. Table 1 presents the species concentrations obtained by PSR calculations for various reaction mechanisms, in comparison with those generated by the original mechanism at a pressure given by p ) 2 atm. Again it can be seen that GA_PSR and GA_MIXED give accurate predictions, while GA_PREMIX fails to predict accurately some of the species. GA_MIXED was found to give good predictions also for PREMIX operating conditions outside the ranges used for the optimization process. Figure 5 presents the predicted burning velocity profile and the main species concentration profiles obtained at an air/fuel ratio Φ ) 2.0, which is outside the range of operating conditions used for optimization. We note (see Figure 4) that all of the mechanisms are accurate in predicting the burning velocity and the main species concentrations. However, for unstable species

Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003 1221 Table 1. Mole Fractions of the Output Species at Operating Conditions That Are Not within the Range Covered by the Data Used in the Inversion Process As Calculated by PSR Using the Reaction Mechanisms Generated by the Inverse GA Formulations Considered GA_PSR original O O2 H H2 OH HO2 H2O H2O2 N2

10-3

4.80 × 2.77 × 10-2 2.86 × 10-2 4.89 × 10-2 5.86 × 10-3 2.40 × 10-5 2.37 × 10-1 3.98 × 10-5 6.47 × 10-1

calcd

GA_PREMIX % error

10-3

4.82 × 2.81 × 10-2 2.83 × 10-2 4.98 × 10-2 5.73 × 10-3 2.41 × 10-5 2.37 × 10-1 4.18 × 10-5 6.47 × 10-1

0.36 1.39 1.11 1.81 2.18 0.36 0.3 4.98 0.02

calcd

GA_MIXED

% error

10-3

4.34 × 2.83 × 10-2 2.77 × 10-2 5.00 × 10-2 5.52 × 10-3 2.25 × 10-5 2.37 × 10-1 6.06 × 10-5 6.47 × 10-1

9.69 2.22 3.15 2.17 5.77 6.41 0.13 52.14 0.05

calcd 10-3

4.72 × 2.80 × 10-2 2.86 × 10-2 4.94 × 10-2 5.75 × 10-3 2.37 × 10-5 2.37 × 10-1 4.04 × 10-5 6.47 × 10-1

% error 1.78 1.11 0.15 0.93 1.77 1.19 0.21 1.31 0.02

Table 2. Mole Fractions of the Output Species As Calculated by PSR for Φ ) 1.0, T ) 1500 K, and p ) 1 atm Using the Reaction Mechanisms Generated by the Inverse GA Formulations Based on Incomplete Sets of Data Considereda REDUCED1 O O2 H H2 OH HO2 H2O H2O2 N2 a

REDUCED2

GA_MIXED

original

calcd

% error

calcd

% error

calcd

% error

2.23 × 10-2 3.33 × 10-2 7.66 × 10-2 6.20 × 10-2 2.29 × 10-2 5.35 × 10-6 1.74 × 10-1 9.30 × 10-6 6.09 × 10-1

2.22 × 10-2 3.34 × 10-2 7.68 × 10-2 6.21 × 10-2 2.30 × 10-2 5.38 × 10-6 1.74 × 10-1 1.01 × 10-5 6.09 × 10-1

0.18 0.28 0.21 0.17 0.33 0.55 0.16 8.14 0.02

2.23 × 10-2 3.35 × 10-2 7.71 × 10-2 6.23 × 10-2 2.29 × 10-2 5.87 × 10-6 1.73 × 10-1 1.02 × 10-5 6.09 × 10-1

0.04 0.73 0.58 0.41 0.15 9.82 0.35 9.31 0.05

2.22 × 10-2 3.36 × 10-2 7.69 × 10-2 6.24 × 10-2 2.28 × 10-2 5.19 × 10-6 1.73 × 10-1 9.23 × 10-6 6.09 × 10-1

0.46 0.97 0.35 0.55 0.6 2.87 0.29 0.73 0.03

For every mechanism the italicized cells mark the species that were not measured in the corresponding inversion procedure.

Table 3. Mole Fractions of the Output Species As Calculated by PSR for Φ ) 1.0, T ) 1500 K, and p ) 2 atm Using the Reaction Mechanisms Generated by the Inverse GA Formulations Based on Incomplete Sets of Data Considereda REDUCED1 O O2 H H2 OH HO2 H2O H2O2 N2 a

REDUCED2

GA_MIXED

original

calcd

% error

calcd

% error

calcd

% error

4.80 × 10-3 2.77 × 10-2 2.86 × 10-2 4.89 × 10-2 5.86 × 10-3 2.40 × 10-5 2.37 × 10-1 3.98 × 10-5 6.47 × 10-1

4.65 × 10-3 2.78 × 10-2 2.83 × 10-2 4.91 × 10-2 5.94 × 10-3 2.45 × 10-5 2.37 × 10-1 6.81 × 10-5 6.47 × 10-1

3.19 0.23 1.05 0.4 1.41 2.23 0.02 70.9 0.02

4.73 × 10-3 2.80 × 10-2 2.86 × 10-2 4.94 × 10-2 5.80 × 10-3 2.72 × 10-5 2.37 × 10-1 7.35 × 10-5 6.47 × 10-1

1.5 0.86 0.06 0.87 0.93 13.34 0.2 84.41 0.01

4.72 × 10-3 2.80 × 10-2 2.86 × 10-2 4.94 × 10-2 5.75 × 10-3 2.37 × 10-5 2.37 × 10-1 4.04 × 10-5 6.47 × 10-1

1.78 1.11 0.15 0.93 1.77 1.19 0.21 1.31 0.02

For every mechanism the italicized cells mark the species that were not measured in the corresponding inversion procedure.

GA_MIXED was again found to outperform GA_PSR and to be approximately equal to GA_PREMIX. Similar conclusions have been obtained for a wide range of operating conditions. Overall, the GA_MIXED reaction mechanism generated by the multiobjective GA proposed in this study was found to be very efficient in predicting the output species concentrations and burning velocity for both PSR and PREMIX simulations. Thus, we may conclude that the inversion procedure proposed efficiently generates in a reduced CPU time a general reaction mechanism that can be used to predict the combustion products for various combustion problems at various operating conditions. 5.2. Multiobjective Inversion Procedures Based on Incomplete Sets of Data. The reaction mechanism investigated in the GA formulations in the previous section were generated using complete sets of data; i.e., the concentrations of all of the output species were assumed to be known for a set of operating conditions. However, in many practical situations, measurements are available only for the most important species. Therefore, it is the purpose of this section to test the GA inversion procedure proposed for the case of incomplete sets of data. Two reaction mechanisms are gener-

ated using two reduced sets of data simulating incomplete sets of measurements as follows: (1) REDUCED1 is generated using multiobjective optimization based on PSR and PREMIX data only for the major species, namely, O2, H2, H2O, and N2. (2) REDUCED2 is generated using multiobjective optimization based on PSR and PREMIX data for the major species O2, H2, H2O, and N2 and for the free radicals O, H, and OH. These two reaction mechanisms are compared with the GA_MIXED reaction mechanism that was developed in the previous section, which assumed that measurements were available for all of the species, including the minor species H2O2 and HO2. Table 2 presents the species concentration predictions obtained for PSR calculations with the above-mentioned reaction mechanisms at T ) 1500 K, Φ ) 1.0, and p ) 1 atm. For every mechanism the cells marked in gray indicate the species for which measurements were not used in generating that reaction mechanism. By investigation of the results obtained for the mechanism REDUCED1, it can be seen that the free radicals O, H, and OH are accurately predicted even if measurements for these species are not used at all. However, the minor

1222

Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003

species H2O2 and HO2, which are in very low concentrations, are less accurately predicted, even if all of the other species are measured; see the results for mechanism REDUCED2 in Table 2. The errors in predicting the minor species are even larger if the mechanisms are tested at a set of operating conditions that was not used in the optimization procedure; see Table 3, which presents the same mole fraction predictions but at a different pressure, namely, p ) 2 atm. Again the main species and the major free radicals are accurately predicted, even if they are not included in the set of data used for the inversion procedure. Moreover, we note that also the minor species are accurately predicted by GA_MIXED; see Tables 2 and 3. Thus, if a reaction mechanism is generated using measurements for the minor species at some operating conditions, then this reaction mechanism accurately predicts all of the species concentrations also at other operating conditions outside the range of operating conditions used for optimization. However, if no measurements are available for the minor species, then the reaction mechanism generated will only give accurate predictions for the major species and the major free radicals. Similar conclusions are obtained for PREMIX predictions; see Figure 6, which presents the profiles of the mole concentrations of free radicals predicted by various mechanisms. Again we note that the major radicals are accurately predicted by all of the mechanisms, while the minor species H2O2 and HO2 are only predicted accurately by GA_MIXED, which used measurements of these radicals in the inversion process. 5.3. Multiobjective Inversion Procedures Based on Different Types of Data. In the previous sections, the reaction mechanisms tested were generated by matching the whole profile of species concentrations and burning velocity generated by a potential solution to the corresponding measured or numerically simulated profile. Next we generate a mechanism GA_PEAK by matching only the peaks of the species concentration profiles given by PREMIX calculations. Such an approach is useful if the whole profile of the species concentrations is not available. The mechanism GA_PEAK was generated using multiobjective optimization based on the set of operating conditions used in the previous sections, but the premix fitness function fPREMIX is changed to

{

fPREMIX ) 10

NS

-8

+

∑ j)1

(

orig |max(Vcalc j ) - max(Vj )|

max(Vorig j )

)}

K

orig |max(Ycalc jk ) - max(Yjk )|

k)1

max(Yorig jk )

∑ Wk

+ -1

(4)

to generate a reaction mechanism based only on the information regarding the peaks of the species concentration profiles. Figure 7 presents the profiles of the species concentration for the free radicals at Φ ) 1.0 and p ) 1 atm, calculated with the GA_PEAK mechanism in comparison with those generated by the GA_MIXED mechanism developed in section 5.1 and the original mechanism. We note that accurate predictions are obtained also for the mechanism GA_PEAK all the way along the burner, even if in generation of this mechanism only information

Figure 6. Profiles of the mole fractions along the burner, for various species, as calculated by PREMIX using the reaction mechanisms generated by the inversion procedures based on incomplete sets of data REDUCED1 (2) and REDUCED2 (9), in comparison with the profiles obtained with GA_MIXED (×) and the original mechanism (s).

about the peaks of the profiles was used. Although not presented here, it is reported that the mechanism GA_PEAK produced accurate predictions for both PSR and PREMIX calculations for a wide range of operating conditions. Therefore, we may conclude that if information about the whole profile of the species concentrations along the burner is not available, then efficient reaction mechanisms may be developed by using only the values of the peak concentrations. 5.4. Validation of the Reaction Mechanisms Generated by GAs. To provide additional validation of the results obtained by the inversion procedures presented in this paper, the predicted ignition delay times are computed using all of the reaction mechanisms generated in this study. In tests of the ignition delay time, the SENKIN code of Sandia National Laboratories16 was used to predict the time-dependent chemical kinetic behavior of a homogeneous gas mixture in a closed system. In addition to prediction of the

Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003 1223 Table 4. Comparison of the Ignition Delay Time (s) Predicted by the Reaction Mechanism Generated by Various GA Inversion Procedures Considered in This Paper for Different Initial Temperatures T ) 1250 K

T ) 1500 K % error

original GA_PSR GA_PREMIX GA_MIXED REDUCED1 REDUCED2 PEAK

4.02 × 10-5 3.68 × 10-5 3.88 × 10-5 3.82 × 10-5 3.58 × 10-5 3.74 × 10-5 3.89 × 10-5

8.32 3.50 4.86 10.97 6.97 3.13

species and temperature histories, the program can also compute the first-order sensitivities with respect to the elementary reaction rate parameters. For the purpose of this paper, the system was considered to be adiabatic under conditions of constant volume and atmospheric pressure, with the initial mole fractions such that the air/fuel ratio is Φ ) 1.0. Table 4 compares the ignition delay times predicted by the reaction mechanisms generated in this study at three different initial tem-

T ) 1750 K % error

2.29 × 10-5 2.19 × 10-5 2.20 × 10-5 2.29 × 10-5 2.14 × 10-5 2.25 × 10-5 2.32 × 10-5

4.66 3.89 0.06 6.48 1.97 1.20

% error 2.33 × 10-5 2.24 × 10-5 2.14 × 10-5 2.37 × 10-5 2.28 × 10-5 2.35 × 10-5 2.41 × 10-5

3.73 8.17 1.60 2.27 0.99 3.23

peratures. We note that all of the mechanisms based on complete sets of data reproduce the same ignition delay times as the original mechanism, and this indicates that the GA does not just generate a mechanism whose validity only extends as far as those operating conditions included in the optimization process. However, it can be seen that overall the GA_MIXED mechanism outperforms the GA_PSR and GA_PREMIX mechanism which are single optimization procedures. GA_PEAK also gives a good prediction for the ignition delay time. The accuracy in predicting the ignition delay time is slightly decreased if the inversion procedures are based on incomplete sets of data. However, also the mechanisms REDUCED1 and REDUCED2 produced good estimates for the ignition delay time. 6. Conclusions In this study we have presented a multiobjective GA inversion procedure that can be used to generate new reaction rate coefficients that successfully predict the profiles of the species produced by a previously validated hydrogen/air reaction scheme. Various formulations have been tested and compared to investigate the possibility of adapting the inversion procedure to various types of experimental data that may be available. The new reaction rate coefficients generated by GAs were found to successfully match mole fraction values in PSR simulations and burning velocity and species profiles in PREMIX simulations, as well as the ignition delay time predicted by the original starting mechanism over a wide range of operating conditions, despite their origin being based solely on measurements for a restricted range of operating conditions. Overall, the reaction mechanism developed by multiobjective GA proposed in this study was found to outperform the reaction mechanisms generated by the single-objective GAs previously employed. Acknowledgment The authors thank the EPSRC and the Corporate Research Program from the MOD and DTI/CARAD Program for funding and permission to publish this work. Literature Cited

Figure 7. Profiles of the mole fractions along the burner, for various species, as calculated by PREMIX using the reaction mechanisms generated within the inversion formulation GA_ PEAK (2), in comparison with the profiles obtained with GA_ MIXED (×) and the original mechanism (s).

(1) Bowman, C. T.; Hanson, R. K.; Gardiner, W. C.; Lissianski, V.; Frenklach, M.; Goldenberg, M.; Smith, G. P. GRI/MECH 2.11. An optimised detailed chemical reaction mechanism for methane combustion and NO formation and reburning; GRI Technical Report 97/0020; GRI: 1997. (2) Dixon-Lewis, G.; Goldsworthy, F. A.; Greenberg, J. B. Proc. R. Soc., London 1975, A346, 261. (3) Dagaut, P.; Reuillon, M.; Boetner, J.-C.; Cathonnet, M. Kerosene combustion at pressures up to 40 atm: experimental

1224

Ind. Eng. Chem. Res., Vol. 42, No. 6, 2003

study and detailed chemical kinetic modelling. Proc. Combust. Inst. 1994, 25, 919-926. (4) Rabitz, H.; Kramer, M.; Dacol, D. Annu. Rev. Phys. Chem. 1983, 34, 419. (5) Milstein, J. The inverse problem: estimation of kinetic parameters. In Modelling of Chemical Reaction Systems; Ezbert, K. H., Deuflhard, P., Jager, W., Eds.; Springer: Berlin, 1981. (6) Bock, H. G. Numerical treatment of inverse problems in chemical reaction kinetics. In Modelling of Chemical Reaction Systems; Ezbert, K. H., Deuflhard, P., Jager, W., Eds.; Springer: Berlin, 1981. (7) Frenklach, M.; Wang, H.; Rabinowitz, J. Optimisation and analysis of large chemical kinetic mechanisms using the solution mapping method- combustion of methane. Prog. Energy Combust. Sci. 1992, 18, 47. (8) Michalewicz, Z. Genetic Algorithms/Data Structures/Evolution Programs, 3rd ed.; Springer: Berlin, 1996. (9) Harris, S. D.; Elliott, L.; Ingham, D. B.: Pourkashanian, M.; Wilson, C. W. The optimisation of reaction rate parameters for chemical kinetic modelling of combustion using genetic algorithms. Comput. Methods Appl. Mech. Eng. 2000, 190, 1065. (10) Elliott, L.; Ingham, D. B.; Kyne, A. G.; Mera, N. S.; Pourkashanian, M.; Wilson, C. W. Incorporation of Physical Bounds on Rate Parameters for Reaction Mechanism Optimisation Using Genetic Algorithms. Combust. Sci. Technol. 2003, in press. (11) Elliott, L.; Ingham, D. B.; Kyne, A. G.; Mera, N. S.; Pourkashanian, M.; Wilson, C. W. The optimisation of reaction rate parameters for chemical kinetic modelling using genetic algorithms. GT-2002 30092 ASME TURBO EXPO, 2002.

(12) Glarborg, P.; Kee, R. J.; Grcar, J. F.; Miller, J. A. PSR: A FORTRAN program for modelling well-stirred reactors; Sandia National Laboratories Report SAND86-8209; Sandia National Laboratories: Albuquerque, NM, 1988. (13) Kee, R. J.; Grcar, J. F.; Smooke, M. D.; Miller, J. A. A Fortran Program for Modelling Steady One-Dimensional Premixed Flames; Sandia Report SAND85-8240; Sandia National Laboratories: Albuquerque, NM, 1985. (14) Kee, R. J.; Miller, J. A.; Jefferson, T. H. CHEMKIN: A general-purpose, problem-independent, transport table, FORTRAN chemical kinetics code package; Sandia National Laboratories Report SAND80-8003; Sandia National Laboratories: Albuquerque, NM, 1980. (15) Van Velhuizen, D.; Lamont, G. Multi-objective Evolutionary Algorithms: Analyzing the State-of-the-Art. Evol. Comput. 2000, 8 (2), 125. (16) Lutz, A. E.; Kee, R. J.; Miller, J. A. SENKIN: A FORTRAN program for predicting homogeneous gas-phase chemical kinetics with sensitivity analysis; Sandia National Laboratories Report SAND87-8248; Sandia National Laboratories: Albuquerque, NM, 1987.

Received for review July 3, 2002 Revised manuscript received December 3, 2002 Accepted January 21, 2003 IE020501O