658
Ind. Eng. Chem. Res. 2005, 44, 658-667
Reaction Mechanism Reduction and Optimization Using Genetic Algorithms Lionel Elliott,† Derek B. Ingham,† Adrian G. Kyne,*,‡ Nicolae S. Mera,‡ Mohamed Pourkashanian,‡ and Christopher W. Wilson§ Department of Applied Mathematics and Energy and Resources Research Institute, University of Leeds, Leeds LS2 9JT, U.K., and Department of Mechanical Engineering, Mappin Street, University of Sheffield, Sheffield S1 3JD, U.K.
This study describes the development of a new binary encoded genetic algorithm to determine a subset of species and their associated reactions that best represent the full starting point reaction mechanism in modeling a low-pressure stoichiometric 12.5% CH4/25% O2/62.5% Ar, burner-stabilized premixed flame. The number of species in the subset chosen is kept fixed and is specified at the start of the procedure. The genetic algorithm chooses better and better mechanisms on the basis of an objective function which measures how well the new reduced mechanisms predict a set of species’ profiles simulated by the full mechanism. To verify the validity of our approach, a full enumeration was performed on a reduced problem, and it was reassuring to find that the genetic algorithm was able to find the optimum solution to this reduced problem in very few generations. A second step was to take the reduced reaction mechanism and to use a second real encoded genetic algorithm to determine the optimal reaction rate parameters that best model an experimental set of premixed flame species’ profiles. This second step not only improved the reduced mechanism’s ability to model the experimental profiles but also provided a remarkable improvement over the mechanism developed from step 1 in modeling combustion processes outside those used during the mechanism’s development. 1. Introduction The numerical simulation of turbulent reacting flow has made significant progress in recent years, due to advances in physical modeling, development of detailed reaction mechanisms, computer hardware, and numerical algorithms. Nevertheless, the detailed modeling of hydrocarbon combustion within complex 3D geometries is presently out of reach given the tools available today, which accounts for the continued interest in methods for describing the chemical kinetics of hydrocarbon combustion more efficiently. Much effort has been invested in developing simplified two- or three-step reaction mechanisms with rate expressions in generalized Arrhenius form. The reduction of chemical reaction mechanisms constitutes an integral part of modeling, and as such, the reduction strategy and specific techniques depend on the objectives of the modeling. The systematic reduction technique pioneered by Peters and Williams,1 the computational singular perturbation method detailed by Lam,2 and the ILDM approach of Maas and Pope3 are among the most prominent efforts in this area of research. While replacing detailed reaction mechanisms with reduced models alleviates CPU time and memory overheads, the accuracy can often be compromised when applying the reduced schemes to a reaction space outside their range of validity. Schwer et al.4 have addressed this problem by replacing the full chemistry model with an entire library of locally ac* To whom correspondence should be addressed. E-mail:
[email protected]. † Department of Applied Mathematics, University of Leeds. ‡ Energy and Resources Research Institute, University of Leeds. § University of Sheffield.
curate, reduced kinetic models, rather than a single skeletal model. Thus, one reduced model can be used over the large regions of the domain where the flow field is constant and the local chemistry is relatively simple, while another scheme can be used to approximate the full chemistry model in regions where the flow is changing rapidly such that the reactant mixing time scales are of the same order as the chemical reaction time scales. There have been a number of approaches put forward in the literature regarding taking a detailed reaction scheme and finding the minimum subset of species/ reactions that best represent a specific combustion regime. Frenklach et al.5 have developed a detailed reduction method whereby reactions from the full mechanism that correspond to low reaction rates and low enthalpy-weighted reaction rates are eliminated. Vajda et al.6,7 and Brown et al.8 developed a similar method (principal component analysis). This time, rather than reactions being eliminated if their reaction rates fall below a user defined maximum threshold, they are eliminated on the basis of a maximum threshold associated with a sensitivity analysis. The drawback with these reduction strategies5-8 is that they cannot guarantee that the reaction mechanism generated is the most compact mechanism that can perform to the required accuracy. To address this drawback, Bhattacharjee et al.9 have developed a linear integer program which can be solved to guarantee global optimality. Here the smallest reaction subset taken from a detailed mechanism is selected such that the new reduced model satisfies the user’s specified accuracy/validity criteria. A further mechanism reduction methodology that has received considerable attention is the computer assisted reduction method (CARM) developed by Chen.10,11 Here
10.1021/ie049409d CCC: $30.25 © 2005 American Chemical Society Published on Web 01/19/2005
Ind. Eng. Chem. Res., Vol. 44, No. 4, 2005 659
a skeletal mechanism is developed by identifying and eliminating reaction steps based upon a combined sensitivity and rate of production analysis at the conditions being investigated. The next step is to then identify quasi-steady-state (QSS) species and remove those species that fall below a specified cutoff level. Finally, the solution of the coupled and nonlinear set of algebraic equations obtained in the previous steps are used to find the QSS species concentrations together with the reaction rates of the non-QSS species. In this investigation the reduction of chemical reaction mechanisms is treated in general in two steps. First, a set of simplified reactions is selected that describes the combustion process over a given range of operating conditions. Once a reduced mechanism has been constructed, the next step is the selection of suitable reaction rate coefficients. It is necessary not only to establish the best-fit values of the unknown parameters but also to estimate their uncertainties in rigorous statistical terms: i.e., to determine the joint confidence region and to examine the correlations between the parameter values. Both steps represent nonconvex inverse problems, the solution of which is a nontrivial task. Evolutionary algorithms are particularly suitable for optimizing object functions with complex, highly structured landscapes. They are powerful optimization techniques that can be applied to solving optimization problems for which the objective functions have a complex structure and gradient information is not available. Computational intelligence methods, such as neural networks, evolution strategies, and genetic algorithms (GA), are highly adaptive methods originating from the laws of nature and biology. Unlike the traditional gradient-based methods, one of the most important characteristics of computational intelligence techniques is the effectiveness and robustness in coping with uncertainty, insufficient information, and noise. In the first step of our reduction process we used a binary coded genetic algorithm to choose “the best” subset of species and all their associated reactions from a detailed reaction mechanism. The number of species in the subset chosen is kept fixed and is specified at the start of the procedure. In determining “the best” subset of species and their associated reactions the performance is based on an “objective function”, which measures how well the new reduced mechanisms predict a set of species’ profiles simulated by the full mechanism. The second step of the reduction process involved taking the optimal skeletal mechanism generated in step 1 and applying a real encoded genetic algorithm12-16 in order to find “the best” set of the Arrhenius reaction rate coefficients that provide the best fit to an experimental set of species’ profiles. Here the “best set” of reaction rate coefficients is judged using an objective function that examines how effectively the reaction mechanisms predict an experimentally observed set of species’ profiles. The search space available to the GA is defined according to the reaction rate uncertainty associated with the National Institute of Standards and Technology (NIST) chemical kinetic database. Defining the search space in this way leads to reaction mechanisms that can be used outside those conditions used in the optimization process.14 2. Combustion Models Laminar premixed flame structure calculations were performed using the PREMIX code for burner-stabilized
flames with a known mass flow rate (see Kee et al.18) The code is capable of solving the energy equation in determining a temperature profile; however, throughout this investigation an experimental temperature profile was imposed in order to rule out the possibility that the predicted temperature profile leads to an inaccurate species profile prediction. Not solving the energy equation had the additional advantage of a significant CPU time savings. For a given reaction mechanism, having specified the mixture composition, flow rate, temperature profile, operating pressure, and solution controls, the PREMIX code outputs the mixture density, velocity, and product species concentration profiles as a function of distance above the burner’s surface. The PREMIX code was also capable of calculating burning velocities by considering a freely propagating adiabatic flame. This feature was used to test our reaction mechanisms against additional experimental data that lay outside the optimization process. To take advantage of the extensive experimental data that exist for measuring ignition delay times, the SENKIN code of Sandia National Laboratories (see Lutz et al.20) was used to predict the time-dependent chemical kinetic behavior of a homogeneous gas mixture in a closed system. In addition to predicting the species and temperature histories, the program can also compute the first-order sensitivities with respect to the elementary reaction rate parameters. Both PREMIX and SENKIN programs are not standalone programs and are designed to be run in conjunction with preprocessors from the CHEMKIN library (see Kee et al.),19 which handle the chemical reaction mechanism and the thermodynamic and transport properties. For every reaction within any given reaction mechanism the forward rate coefficients are in the three-parameter functional Arrhenius form, namely
( )
Eai kfi ) Ai Tβi exp RT
(1)
for i ) 1, ..., NR, where R is the universal gas constant and there are NR competing reactions occurring simultaneously. The reaction rates (1) contain the three parameters Ai, βi, and Eai for the ith reaction. These parameters are of paramount importance in modeling the combustion processes, since small changes in the reaction rates may produce large deviations when modeling the output species’ profiles in premixed flames or ignition delay times. 3. Optimization Models 3.1. Step 1: Mechanism Reduction. While detailed comprehensive kinetic mechanisms are necessary to fully understand the fundamental chemistry of combustion, they are often impractical for use in industrial codes because of their size. Using GRImech 3.017 as a starting point, the first step of this investigation was to use a binary encoded genetic algorithm (GA) to establish the best reduced subset of species and their associated reactions that best represent a set of species profiles simulated from the full starting point mechanism. In a simple GA, the process of evaluation assigns each individual solution in a population with an associated fitness value indicating the individual’s survivability. Subsequently, the GA must evolve this population of individuals into a new population using the three operations of selection, crossover, and mutation. Selec-
660
Ind. Eng. Chem. Res., Vol. 44, No. 4, 2005
tion probabilistically chooses better individuals for a new generation, while crossover and mutation manipulate candidate solutions to generate new individuals for the selection procedure to process again. After a certain number of generations, when no further improvement is observed, the best chromosome represents an optimal solution. With the full mechanism as the starting point, GRImech 3.0, which contains nrfull ) 325 reactions involving nsfull ) 53 species, an initial starting population of individuals is formed where each individual represents a possible reduced species model. Each individual is represented by a binary vector (chromosome), z, of length nsfull, where zi ) 1 indicates the ith species is included in the reduced model while zi ) 0 indicates that the ith species is not included (i ) 1, 2, ..., nsfull). The number of species, nsreduced, selected in each reduced species model is kept fixed and specified at the start; full thus, Σn1 s zi ) nsreduced. A further constraint to the model is that a number of “definite species”, ndef, are defined such that there are ndef positions in the binary vector, z, where zi must always equal 1. This constraint is used to ensure that fuel and oxidizer species investigated in the PREMIX combustion model are always present in the reduced mechanisms generated by the GA. In addition, major product species or species of interest can be forced to be ever present, thus reducing the search space and subsequently increasing the rate of convergence. The initial starting population consists of nindiv individuals each containing the same definite species with the remaining (nsreduced - ndef) species randomly chosen from the remaining (nsfull - ndef) species of the full mechanism. For example, for an initial population of five individuals where nsfull ) 20, nsreduced ) 10, and ndef ) 3, we could have
individual 1
11100011010111000001
individual 2
10101101010011001100
individual 3
11110000010111100100
individual 4
10101100111001010010
individual 5
10101010100001100111
where the boldface 1’s represent the definite species and the remaining species of the reduced model are randomly distributed along the chromosome. Each reduced species set has, associated with it, a set of possible reactions. Here a reaction is chosen only if all reactant and product species exist within the reduced species set. Also, any third-body reactions will only be included if both the reactant/product and all effective third-body species appear in the reduced species set. The reaction rate coefficients remain the same as those of the full starting point mechanism. Having generated an initial population of individuals, the next step is to devise a criterion for judging how “fit” an individual is. For this process, Chemkin’s PREMIX code is used to generate a set of species profiles for each of the individuals considered. This set of profiles is then compared to a corresponding set of profiles generated with the full starting mechanism using the objective function
{
fPREMIX ) +
(
Nconditions nsreduced
∑ j)1
∑
k)1
)}
|Yjkcalcd - Yjkorig|L2
Wk
|Yjkorig|L2
-1
(2)
where Yjkcalcd represents the mole concentrations of the kth species in the jth set of operating conditions. Yjkorig denotes the corresponding original values that represent the simulated species profiles generated using the full GRI mech 3.0. | |L2 represents the L2 norm of a function, which is numerically calculated using a trapezoidal rule. Wk denotes the normalized weights (lying in the range between 0 and 1) that can be applied to each species, depending on the importance of the species. For those species whose simulated maximum concentration was greater than 10-7, a value of 1.0 was assigned, while the remaining species were considered unimportant and were set to zero. Finally, is a constant added to the fitness function in order to avoid numerical overflow. We note that the maximum value of these object functions is 1/, which indicates a perfect fit to the data. A value of ) 10-8 was used. The validity of the method lies in how it decides what makes one reduced mechanism better than another. This is subjective, depending on what the user wants out of their reaction model, and there are a number of alternative methods of defining the model’s accuracy that could have been used. The accuracy of the reduced mechanisms chosen in this study is given by eq 2. Here an individual is defined to be better than another individual if the sum of the percentage errors divided by 100 in predicting the nsreduced species profiles simulated by the full mechanism is smaller than that of another reduced reaction model. To avoid comparing “less important” species, a condition was added to eq 2 to only consider those species whose peak simulated mole fraction concentration was greater than 10-7. The species whose peak simulated mole fraction concentration is less than 10-7 do not contribute to the fitness but are still important, as they can allow important reactions to be used in the reduced model which may significantly impact the “more important” species profiles. Cases can arise where one reduced model is ranked higher than another but it does not contain a species that the user is interested in or this species is poorly predicted relative to the other species. For such cases, the fitness function can easily be modified to ensure certain species are always chosen by adding these species of interest to the definite species list and/or to bias the reduced mechanisms’ ability to predict these species by modifying the weighting factors, Wk. Having calculated the fitness of each of the individuals, k-tournament selection21 is used to select the fittest individuals to be parents for the next generation. Genetic operators such as crossover and mutation21 were then used to exchange information between the parents in order to form a population of children. To construct the new population, the Elitist model was then used, in which the ne (the elitism parameter) best individuals from the old population were retained in the next population along with the fittest (npop - ne) child chromosomes. These child chromosomes were chosen without repetition, such that no individual would appear twice in the new population. The constraints associated with our problem led us to disregard the more traditional crossover and mutation operators21 in favor of a modified approach. The approach we adopted in exchanging information between parents involved selecting two parents and identifying those species common to both parents. A child is then formed by keeping these common species and randomly choosing the remaining species from the remaining unused positions of the
Ind. Eng. Chem. Res., Vol. 44, No. 4, 2005 661
Figure 1. Crossover operator adopted in this investigation.
chromosome such that the total number of species is equal to nsreduced (see Figure 1). To avoid convergence to local optima, a mutation mechanism has been incorporated. Here a child will be selected for mutation with the probability pmut. If a child is selected, then those positions not containing a definite species have an equal probability of mutating (i.e. changing from a 0 to a 1 or vice versa). The probability of a definite species location being mutated is set equal to zero. The number of mutations carried out along the child chromosome, nmut, is specified at the start of the algorithm. After mutation the number of species in the new child chromosome must equal nsreduced. Thus, for nmut ) 1, mutation involves randomly choosing a nondefinite species location and flipping the value in that position. If the value is a 1, another location containing a 0 is randomly chosen and flipped so that the number of species remains equal to nsreduced. Once the children population has been constructed using crossover and mutation, eq 2 is again used to calculate each child’s fitness. A standard filter function21 is then used to generate a new population by merging the fittest parents and children. From this new population the whole process is repeated with selection, crossover, and mutation used to generate a new children population and so on until convergence is reached where no improvement is observed to the fittest individual. 3.2. Step 2: Reduced Mechanism Optimization. Having found a reduced set of species and their associated reactions that best match a simulated set of species profiles generated using the starting point GRI mech 3.0 scheme, the next step was to optimize this reduced set of reaction rate coefficients. This represents a nontrivial task, as the large number of parameters in the problem creates a solution landscape that is uneven, making it unsuitable for gradient-based methods. Thus, a real coded GA was adopted. A detailed description of the real coded GA approach to reaction rate optimization can be found in the studies of Harris et al.21 and Elliott et al.12-16 The overall methodology is similar to that described in the previous section. The differences lie in how the parent and children chromosomes are defined. This time an individual consists of a set of reaction rate coefficients (real numbers): i.e., Ai, βi, and Eai (eq 1) for each of the reactions included in the reduced model. The initial population consists of a number of sets of these reaction rate coefficients where each coefficient lies between a predetermined bound. In our earlier work12,21 an artificial bound of (25% was placed on each of the reaction rate coefficients. However, it was found that this often led to unrealistic rates that limited the applicability of the reaction mechanisms generated.
Instead, the bounds that were used were based on the uncertainty associated with the experimental findings listed in the National Institute of Standards and Technology (NIST) database. For a given reaction and over a given temperature range, the NIST database contains a number of investigators’ experimental findings for the reaction rate. For certain reactions these findings can be quite conflicting, and Elliott et al.14 have defined a reaction rate search space based on the enclosed area bounded by the maximum and minimum values of these reaction rate findings for a given temperature range. Bounds are then placed on Ai, βi, and Eai such that all calculated reaction rates occupy a region just large enough to contain that of the previously defined search space. A check is then made every time the GA selects Ai, βi, and Eai to ensure that the subsequent reaction rate lies within the NIST search space. If the chosen Ai, βi, and Eai result in a rate that lies outside the NIST search area, then the individual is penalized with a very low fitness value in order to prevent it from surviving the selection process and migrating to the next generation or creating offspring. In adopting this method for defining a physically realistic search space, the reaction mechanisms generated are better equipped for modeling combustion systems that lie outside those used in the optimization process (see Elliott et al.14). Having constructed the initial population, k-tournament selection is used to select parents which undergo the traditional uniform arithmetic crossover and nonuniform mutation (see Harris et al.21). The nonuniform mutation operator (see Michalewicz23) 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 fitness of the individuals is calculated using the same equation (eq 2) as in the previous section; however, this time Yjkorig represents an experimental measurement taken from Bernstein et al.22 of the kth species and jth premix condition. Also, Wk is set to 1 for each of the species where data exists and 0 for those species not measured. 4. Results and Discussion 4.1. Mechanism Reduction. The test problem considered is the combustion of a stoichiometric 12.5% CH4/ 25% O2/62.5% Ar burner-stabilized flame with a mixture flow rate of 0.004 417 (g/(cm2 s)) at a constant reduced pressure of 20 Torr. The reason for these operating conditions will become apparent in the next section, which compares the experimental findings of Bernstein et al.22 obtained under the same conditions. It is these experiments that provide the input temperature profile for our simulations. A genetic algorithm was used to choose the best set of nsreduced species and their associated reactions, as described in section 3.1. The performance of the reduced reaction mechanisms generated is judged on how well they match those output species profiles that have a maximum species concentration above 10-7, as simulated using the full GRI 3.0 mech
662
Ind. Eng. Chem. Res., Vol. 44, No. 4, 2005
Table 1. Sensitivity of pmut and nmut to the Average Number of GA Generations in Choosing the Optimal 4 Species from a Pool of 41, Given 12 Definite Species, 6 Different Starting Seeds, and an Average Taken pmu
nmut
av no. of generations
av no. of function calls
approx CPU time (s)
0.1 0.1 0.1 0.01 0.05 0.2
1 2 3 2 2 2
35 21.7 37.7 44.2 33.2 26.3
1074 675 1035 1350 1020 813
1832 1152 1766 2304 1740 1387
(see eq 2) for the operating conditions described above. To gain confidence in our reduction technique, a full enumeration was conducted to establish the best possible reduced mechanism that exists based on eq 2. The difficulty in a full enumeration is the vast size of the problem. For instance, in finding the optimum 16 species reduced scheme from the full 53 present in GRI 3.0 mech, where 5 definite species are ever present, involves 53-5C16-5 ) 2.26 × 1010 function evaluations, which is computationally unfeasible when one considers that each function evaluation involves running PREMIX. Thus, for the purpose of full enumeration, a much smaller search domain must be tested. To achieve this, the same values for nsreduced ) 16 and nsfull ) 53 were used and the number of definite species being ever present was increased (ndef ) 12), so that the total number of all possible combinations was 53-12C16-12 ) 101 270. While 101 270 calls to PREMIX may still seem excessive, it should be pointed out that many of the reaction mechanisms were small and caused PREMIX to terminate after less than 1 s. The total CPU time taken for the full enumeration on a 3.2 GHz Pentium 4 was 2 days. The 12 ever-present species in the reduced mechanism were CH4, O2, N2, Ar, CO2, H2O, CO, H2, O, OH, H, and CH3. The fitness of all possible combinations of choosing the remaining 4 species from the remaining 41 species of the full mechanism was calculated, and the best fitness and corresponding reduced reaction mechanism was stored. For the same conditions the GA was set the task of finding this optimum reduced reaction mechanism. Table 1 demonstrates the success of the GA in arriving at this optimum solution for a number of cases of interest. In addition to increasing confidence in the ability of the GA to find the optimum solution, the cases considered investigate the sensitivity of the mutation parameters pmut and nmut on the rate of convergence. The departure of our mutation operator from more traditional mutation operators has motivated the interest in these particular parameters. In all cases the GA was able to arrive at the optimum reduced reaction scheme in less than 60 generations. The average number of generations quoted is calculated from a total of 6 runs, where each run used a different starting seed for the random number generator. While the problem represents a relatively simple one for a GA, nevertheless the speed at which the GA finds the solution is reassuring and breeds confidence that the method can be extended to successfully searching a larger parameter space. Following these calculations,
Table 2. Description of GA Generated Reduced Reaction Schemes mechanism descripn Reduced 12sp Reduced 16sp Reduced 20sp
species H, O, O2, OH, H2O, CH2(s), CH3, CH4, CO, CO2, HCO, Ar H2, H, O, O2, OH, H2O, CH2(s), CH3, CH4, CO, CO2, HCO, CH2O, C2H5, C2H6, Ar H2, H, O, O2, OH, H2O, HO2, CH, CH2, CH2(s), CH3, CH4, CO, CO2, HCO, CH2O, CH3O, C2H5, C2H6, Ar
no. of corresponding reacns 16 54 99
the genetic operators and GA parameters were taken to be
population size, npop ) 24 number of offspring, nchild ) 30 tournament selection, tournament size, k ) 2; tournament probability, pt ) 0.8 probability a child is selected for mutation, pmut ) 0.1 number of mutations carried out along the child chromosome, nmut ) 2 elitism, elitism parameter, ne ) 2 Having gained confidence in our methodology, the next task was to apply it to generating the optimal 12, 16, and 20 species reaction mechanisms that best represent the full GRI mech 3.0 under the same premix flame operating conditions. For each case the number of definite species remained the same and the everpresent species were defined to be CH4, O2,Ar, CO2, and H2O. The deficit species were then made up from the remaining 48 species of the full mechanism. In reaching a solution, the GA was allowed to run for 1000 generations, even though the fitness evolution had reached a plateau for all cases after less than 500 generations. Table 2 describes which species were chosen in each of the 3 reaction mechanisms generated and the number of reactions associated with these species. For nsreduced ) 16 and 20 a random initial population was used to start the GA. The difficulty when nsreduced ) 12 was that, because there were relatively few species, there existed a vast number of combinations of species which had associated with them very few reactions and a mechanism consisting of only a handful of reactions would not give a solution. When this happened, a very low fitness was assigned to these poor individuals so that their genes were not passed on to the children population. However, because the vast majority of possible reduced reaction mechanisms contained only a few reactions, when an initial population was established, all the individuals were assigned the same low fitness value and the GA broke down. In order for the GA to proceed, it needed at least one of the individuals in the starting population to have a reasonable fitness value. Thus, to ensure this, a starting population was manually constructed such that there was at least one individual that consisted of the more popular species from the full reaction mechanism. A further test was made on our method in order to add credence to its validity. Taking the slowest case to converge, namely Reduced 20sp (largest search space),
Ind. Eng. Chem. Res., Vol. 44, No. 4, 2005 663
that less than 1% of the CPU time is associated with the GA. In fact, one generation requires running the PREMIX application 30 (number of individuals) × 1000 (number of generations). Although GAs are relatively CPU intensive compared to other gradient-based solution methods, they are very robust, can be applied to many different scenarios (it is thought they could equally well have been applied to a more detailed starting reaction mechanism), and offer a much better guarantee that the solution represents a global and not a local one. 4.2. Reduced Mechanism Optimization. The second step of the investigation uses a GA developed from previous studies12-16 in order to tackle the problem of keeping the same set of species/reactions but finding the best set of reaction rate constants that best fit an experimental set of premix species’ profiles.22 The same premix flame operating conditions of the previous section were used, and the best set of reaction rate coefficients was determined on the basis of the fitness function described in section 3.2. This time the genetic operators and GA parameters were taken to be:
population size, npop ) 24 number of offspring, nchild ) 30 uniform arithmetic crossover, crossover probability, pc ) 0.65 tournament selection, tournament size k ) 2; tournament probability, pt ) 0.8 Figure 2. Performance of three reduced reaction schemes (see Table 2) in predicting GRI 3.0’s simulated profiles of (a) O2, (b) CO2, (c) H2O, (d) CO, (e) OH, (f) O, (g) H, and (h) CH3 in a stoichiometric 12.5% CH4/25% O2/62.5% Ar burner-stabilized flame at a constant reduced pressure of 20 Torr.
a full enumeration was performed to check if the GA had indeed found the best 20 species. It was thought that it would be reasonable to assume that the 16 species that were common with Reduced 16sp would be part of the optimum 20 species, and so these were made definite species. Next, the fitness values of all combinations of choosing the remaining 4 species from the 37 nondefinite species of the full mechanism were calculated and the maximum fitness was verified to be that calculated using the GA. Figure 2 compares the performance of all three reaction mechanisms in predicting the numerous simulated species’ profiles of the full mechanism. As expected, the more detailed the mechanism, the closer the correlation with the simulated data. It can be seen in Figure 2a that the 12-species mechanism significantly overpredicts the oxidizer consumption, which in turn has a knock-on effect with its ability to model the major and intermediate species;this can lead to serious inaccuracies when applied to predicting important combustion phenomena such as ignition delay and burning velocity. From Figure 2 we come to the conclusion that the 16-species mechanism represents the best compromise between accuracy and efficiency and it is for this reason that we adopt this scheme in the next step of the optimization process. In arriving at the optimum 16-species reaction mechanism, the process took approximately 12 h to run on a 3.2 GHz Pentium 4 computer. It is interesting to note
nonuniform mutation, mutation probability, pm ) 0.5 elitism, elitism parameter, ne ) 2 To enhance convergence, the GA inversion procedure uses Reduced 16sp (see Table 2) of the previous section as an initial guess in the starting population. The procedure is unbiased, with each of the six species investigated22 contributing equally to the fitness function of eq 2. A full description of Reduced 16sp GA can be found in Table 3. Figure 3 compares various species’ profiles for the full mechanism (GRI3.0), reduced mechanism (Reduced 16sp), and the reduced mechanism with reaction rate parameters optimized by the GA (Reduced 16sp GA) with the experimental measurements of Bernstein et al.22 (expt). Figure 3a compares the measured and predicted profiles of the CH3 radical. It can be seen that, while both the full and reduced models overpredict the measurements, the GA optimized reduced scheme shows quite a remarkable correlation. Figure 3b displays the concentration profiles of CO, and it can be seen that all predictions lie significantly below those found experimentally. Having improved upon the Reduced 16sp and GRI 3.0 schemes using the GA, an attempt was made to further improve the CO profile by biasing the GA such that the remaining 5 species concentrations are independent of changes in reaction rate. This was achieved by assigning WCO ) 1 and Wk*CO ) 0 in eq 2. Although it has not been shown here, it should be noted that the improved CO prediction was obtained at the significant expense of the remaining species’ predictions. On the basis of the success of the GA in recovering simulated species’ profiles,16 it can be
664
Ind. Eng. Chem. Res., Vol. 44, No. 4, 2005
Table 3 Species Considered temp
1. H2 2. H 3. O 4. O2 5. OH 6. H2O 7. CH2(s) 8. CH3 9. CH4 10. CO 11. CO2 12. HCO 13. CH2O 14. C2H5 15. C2H6 16. Ar
element count
phase
charge
mol wt
low
high
O
H
C
N
Ar
G G G G G G G G G G G G G G G G
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
2.015 94 1.007 97 15.999 40 31.998 80 17.007 37 18.015 34 14.027 09 15.035 06 16.043 03 28.010 55 44.009 95 29.018 52 30.026 49 29.062 15 30.070 12 39.948 00
200.0 200.0 200.0 00.0 200.0 200.0 200.0 200.0 200.0 200.0 200.0 200.0 200.0 200.0 200.0 300.0
3500.0 3500.0 3500.0 3500.0 3500.0 3500.0 3500.0 3500.0 3500.0 3500.0 3500.0 3500.0 3500.0 3500.0 3500.0 5000.00
0 01 1 2 1 1 0 0 0 1 2 1 1 0 0 0
2 1 0 0 0 2 2 3 4 0 0 1 2 5 6 0
0 0 0 0 0 0 1 1 1 1 1 1 1 2 2 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
Reactions Considereda enhanced by 1. 2O + M a O2 + M H2 H2O CH4 CO CO2 C2H6 Ar 2. O + H + M a OH + M H2 H2O CH4
A
b
E
enhanced by
A
b
E
9.07E+16 -0.6 0.0 17. 2H + H2 a 2H2 2.400E+00 18. 2H + H2O a H2 + H2O 5.97E+19 -1.3 0.0 5.52E+20 -2.0 0.0 1.540E+01 19. 2H + CO2 a H2 + CO2 2.000E+00 20. H + OH + M a H2O + M 2.27E+22 -1.6 0.0 7.300E-01 1.750E+00 H2 3.600E+00 H 2O 3.650E+00 3.000E+00 CH4 2.000E+00 8.300E-01 C2H6 3.000E+00 7.67E+17 -0.8 0.0 Ar 3.800E-01 2.000E+00 21. H + CH3 (+M) a CH4 (+M) 1.09E+16 -0.6 641.6 6.000E+00 low-pressure limit: 0.26200E+34, -0.47600E+01, 0.24400E+04 2.000E+00 TROE centering: 0.78300E+00, 0.74000E+02, 0.29410E+04, 0.69640E+04 CO 1.500E+00 H2 2.000E+00 CO2 2.000E+00 H 2O 6.000E+00 3.000E+00 CH4 3.000E+00 C2H6 Ar 7.000E-01 CO 1.500E+00 3. O + H2 a H + OH 2.000E+00 5.02E+04 2.8 7043.9 CO2 4. O + CH2(s) a H2 + CO 3.000E+00 2.26E+13 0.0 0.0 C 2H 6 5. O + CH2(s) a H + HCO 1.49E+13 0.0 0.0 Ar 7.000E-01 5.71E+08 1.5 12334.3 6. O + CH3 a H + CH2O 6.50E+13 0.0 0.0 22. H + CH4 a CH3 + H2 7. O + CH4 a OH + CH3 7.94E+08 1.5 8839.2 23. H + HCO (+M) a CH2O (+M) 1.46E+12 0.6 -283.4 8. O + CO (+M) a CO2 (+M) 1.80E+10 0.0 2399.5 low-pressure limit: 0.24700E+25, -0.25700E+01, 0.42500E+03 low-pressure limit: 0.60200E+15, 0.00000E+00, 0.30000E+04 TROE centering: 0.78240E+00, 0.27100E+03, 0.27550E+04, 0.65700E+04 2.000E+00 H2 2.000E+00 H2 O2 6.000E+00 H 2O 6.000E+00 2.000E+00 H2O 6.000E+00 CH4 CH4 2.000E+00 CO 1.500E+00 CO 1.500E+00 CO2 2.000E+00 CO2 3.500E+00 C 2H 6 3.000E+00 C2H6 3.000E+00 Ar 7.000E-01 Ar 5.000E-01 24. H + HCO a H2 + CO 7.93E+13 0.0 0.0 5.16E+07 2.0 3172.0 9. O + HCO a OH + CO 2.98E+13 0.0 0.0 25. H + CH2O a HCO + H2 10. O + HCO a H + CO2 2.98E+13 0.0 0.0 26. H + C2H5 (+M) a C2H6 (+M) 2.28E+17 -1.1 1693.1 11. O + CH2O a OH + HCO 3.99E+13 0.0 3957.3 low-pressure limit: 0.19900E+42, -0.70800E+01, 0.66850E+04 12. O + C2H5 a CH3 + CH2O 1.95E+13 0.0 0.0 TROE centering: 0.84220E+00, 0.12500E+03, 0.22190E+04, 0.68820E+04 8.83E+07 2.0 7709.3 H2 2.000E+00 13. O + C2H6 a OH + C2H5 14. O2 + CO a O + CO2 2.17E+12 0.0 44185.2 H 2O 6.000E+00 2.000E+00 15. H + O2 a O + OH 2.48E+16 -0.6 16847.6 CH4 16. 2H + M a H2 + M 6.57E+17 -1.1 0.0 CO 1.500E+00 0.000E+00 CO2 2.000E+00 H2 H2O 3.000E+00 0.000E+00 C 2H 6 CH4 2.000E+00 Ar 7.000E-01 CO2 0.000E+00 27. H + C2H6 a C2H5 + H2 1.08E+08 1.9 9629.0 C2H6 3.000E+00 Ar 6.300E-01 6.72E+16 -1.3 0.0
Ind. Eng. Chem. Res., Vol. 44, No. 4, 2005 665 Table 3. (Continued) Reactions Considereda enhanced by
A
b
E
enhanced by
A
b
E
28. H2 + CO (+M) a CH2O (+M) 4.34E+07 1.5 79600.0 44. CH3 + O2 a OH + CH2O 2.57E+12 0.0 16918.7 low-pressure limit: 0.50700E+28, -0.34200E+01, 0.84350E+05 45. 2CH3 (+M) a C2H6 (+M) 6.65E+16 -1.3 653.6 TROE centering: 0.93200E+00, 0.19700E+03, low-pressure limit: 0.34000E+42, -0.70300E+01, 0.27620E+04 0.15400E+04, 0.10300E+05 H2 2.000E+00 TROE centering: 0.61900E+00, 0.73200E+02, 0.11800E+04, 0.99990E+04 H2O 6.000E+00 H2 2.000E+00 CH4 2.000E+00 H 2O 6.000E+00 CO 1.500E+00 CH4 2.000E+00 CO2 2.000E+00 CO 1.500E+00 C2H6 3.000E+00 CO2 2.000E+00 Ar 7.000E-01 C2H6 3.000E+00 29. OH + H2 a H + H2O 1.95E+08 1.4 3418.0 Ar 7.000E-01 30. 2OH a O + H2O 4.31E+04 2.4 -2307.7 46. 2CH3 a H + C2H5 6.40E+12 0.1 11454.5 31. OH + CH2(s) a H + CH2O 2.56E+13 0.0 0.0 47. CH3 + HCO a CH4 + CO 3.34E+13 0.0 0.0 32. OH + CH3 a CH2(s) + H2O 6.93E+17 -1.3 1503.0 48. CH3 + CH2O a HCO + CH4 2.69E+03 3.0 4959.7 33. OH + CH4 a CH3 + H2O 1.01E+08 1.7 2872.1 49. CH3 + C2H6 a C2H5 + CH4 7.18E+06 1.8 10259.1 34. OH + CO a H + CO2 5.48E+07 1.1 67.2 50. HCO + H2O a H + CO + H2O 9.66E+17 -1.1 19809.6 35. OH + HCO a H2O + CO 6.13E+13 0.0 0.0 51. HCO + M a H + CO + M 1.83E+17 -0.8 17838.9 36. OH + CH2O a HCO + H2O 4.41E+09 1.2 -465.3 H2 2.000E+00 37. OH + C2H6 a C2H5 + H2O 4.96E+06 2.2 602.2 H 2O 0.000E+00 38. CH2(s) + O2 a H + OH + CO 1.51E+13 0.0 0.0 CH4 2.000E+00 39. CH2(s) + O2 a CO + H2O 1.19E+13 0.0 0.0 CO 1.500E+00 40. CH2(s) + H2 a CH3 + H 3.66E+13 0.0 0.0 CO2 2.000E+00 41. CH2(s) + CH4 a 2CH3 8.81E+12 0.0 -623.7 C 2H 6 3.000E+00 42. CH2(s) + CO2 a CO + CH2O 6.38E+12 0.0 0.0 52. O + CH3 f H + H2 + CO 3.40E+13 0.0 0.0 43. CH2(s) + C2H6 a CH3 + C2H5 2.31E+13 0.0 -522.3 53. OH + CH3 f H2 + CH2O 1.16E+10 0.7 -2282.3 54. CH2(s) + H2O f H2 + CH2O 6.77E+10 0.3 -927.9 a
k ) ATb exp(-E/RT).
Figure 3. Mole fraction concentration profiles of (a) CH3, (b) CO, (c) O, (d), OH, (e) H, and (f) HCO as function of distance from burner surface for a premixed stoichiometric 12.5% CH4/25% O2/ 62.5% Ar burner-stabilized flame at a constant reduced pressure of 20 Torr (9, experimental data of Bernstein et al.22).
concluded that, within the framework of the mechanism being optimized, it is not possible to recover this CO profile without using nonphysical rate constants that
lie outside the constraints of the NIST database. Parts c and d of Figure 3 compare measured and predicted O and OH profiles, respectively. In both cases the shortcomings in the predictions of the Reduced 16sp and GRI 3.0 mechanisms are eradicated by the GA optimized mechanism, which is in excellent agreement with the measured data. Figure 3e compares the performance of all 3 models in predicting the measured H atom profiles. It can be seen that both the GRI 3.0 and GA generated mechanisms are indistinguishable from each other, both showing an improvement over the H atom overprediction obtained using the Reduced 16sp mechanism. Figure 3f compares measured and predicted profiles of HCO. While the GRI 3.0 mechanism underpredicts the measurements, both the reduced and GA optimized reduced models show a similar level of improvement in capturing the experimental profile. We have demonstrated how a reduced reaction model that contains less than one-fourth of the number of species and less than one-sixth the number of reactions of the full model is still able to outperform it under a certain set of conditions. However, the real test lies in how well it performs in modeling experimental conditions that are situated outside those used in the optimization process. Figure 4 investigates this by comparing all three models against the ignition delay data of Tsuboi and Wagner24 for 0.2% CH4/2% O2/97.8% Ar mixtures for pressures of 21-29 atm with T ) 14002000 K and pressures of 3-4 atm with T ) 1650-2050 K. It is perhaps not surprising that the reduced model developed solely on its ability in recovering premix species’ profiles simulated by a detailed mechanism is not applicable to modeling ignition delay, particularly at pressures in the range 21-29 atm. However, what is remarkable is the improvement in the predicted ignition delay times for both pressure and temperature
666
Ind. Eng. Chem. Res., Vol. 44, No. 4, 2005
5. Conclusions
Figure 4. Ignition delay times in a 0.2% CH4/2% O2/97.8% Ar mixture: (b) experimental data of Tsuboi and Wagner for pressures of 3-4 atm with T ) 1650-2050 K; (9) experimental data of Tsuboi and Wagner24 for pressures of 21-29 atm with T ) 1400-2000 K.
A two-step approach to reaction mechanism reduction has been demonstrated. The first step takes the full GRI mech 3.0 and using a binary encoded GA finds an optimum subset of species and associated reactions that best represents the simulated species’ profiles generated by the full starting mechanism. The number of species present in each reduced subset was defined at the start of the process. Full enumeration was conducted on a reduced problem in order to validate that the GA was indeed finding the optimum reduced reaction set. Three reduced mechanisms containing the optimum 12, 16, and 20 species were then compared to the simulated species’ profiles generated by the full GRI mech 3.0, and on the basis of this comparison the 16-species mechanism was chosen for the second step of the optimization. The second step used a real encoded GA to find the optimum set of reaction rate coefficients of the reduced mechanism generated in step 1 that best represented an experimental set of premix species’ profiles. An excellent agreement was found between the GA optimized reduced scheme and the experimental data. Comparisons with ignition delay and burning velocity data provided evidence that this second step allowed the mechanism to be applied to combustion regimes outside those used in the optimization process. A full description of the final reduced reaction mechanism with optimized reaction rate coefficients can be found in Table 3. Acknowledgment We thank the EPSRC for the funding for this work. Literature Cited
Figure 5. Flame velocity as a function of equivalence ratio for a methane/air mixture at 1 atm and Tu ) 298 K: (b) experimental data of Taylor;25 (9) experimental data of Vagelopoulos et al.26
ranges that occurs when this reduced mechanism is taken and its reaction rates optimized to give the best fit to an experimental set of premix species’ profiles. Figure 5 gives further evidence on how well the three models perform at modeling a different combustion system. Here the mechanisms are applied to predicting the laminar flame velocity measurements of Taylor25 and Vagelopoulos et al.26 Again, a noticeable improvement is apparent in applying step two of our optimization process; however, this time Reduce 16sp GA is unable to outperform the full mechanism in capturing the peak burning velocity profile. It should be pointed out that for modeling this particular problem it was necessary to add N2 to the reduced mechanisms, as air was present in the oxidizer. No associated reactions were added.
(1) Peters, N.; Williams, F. A. The asymptotic structure of stoichiometric Methane-Air Flames. Combust. Flame 1987, 86, 185. (2) Lam, S. L. Using CSP to understand complex chemical kinetics. Combust. Sci. Technol. 1993, 89, 375. (3) Mass, U.; Pope, S. B. Simplifying chemical kinetics: intrinsic low-dimensional manifolds in composition space. Combust. Flame 1992, 88, 239. (4) Schwer, D. A.; Lu, P.; Green, W. H. Adaptive chemistry approach to modelling complex kinetics in reacting flows. Combust. Flame 2003, 33, 451. (5) Frenklach, M.; Wang, H. Detailed reduction of reaction mechanisms for flame modelling. Combust. Flame 1991, 87, 365. (6) Vajda, S.; Valko, P.; Turanyi, T. Principal component analysis of kinetic models. Int. J. Chem. Kinet. 1985, 55, 81. (7) Turanyi, T.; Berces, T.; Vajda, S. Reaction rate analysis of complex kinetics systems. Int. J. Chem. Kinet. 1989, 21, 83. (8) Brown, N. J.; Li, M. L.; Koszykowski, Mechanism reduction via principal component analysis. Int. J. Chem. Kinet. 1997, 393. (9) Bhattacharjee, B.; Schwer, D. A.; Barton, P. I.; Green, W. H. Optimally reduced kinetic models: reaction elimination in large-scale kinetic mechanisms. Combust. Flame 2003, 135, 191. (10) Chen, J.-Y. A general procedure for constructing reduced reaction mechanisms with given independent reactions. Combust. Sci. Technol. 1988, 57, 89. (11) Chen J.-Y. Development of reduced mechanisms for numerical aspects of reduction in chemical kinetics. CERMICS-ENPC Cite Descartes, Champus sure Marne, France, 1997. (12) Elliott, L.; Ingham, D. B.; Kyne, A. G.; Mera, N. S.; Pourkashanian, M.; Wilson, C. W. A real coded genetic algorithm for the optimization of reaction rate parameters or chemical kinetic modelling in a perfectly stirred reactor. Genetic and Evolutionary Computation Conference; Late Breaking Papers; GECCO: New York, 2002; pp 138-144. (13) 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.
Ind. Eng. Chem. Res., Vol. 44, No. 4, 2005 667 (14) 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, 175, 619. (15) Elliott, L.; Ingham, D. B.; Kyne, A. G.; Mera, N. S.; Pourkashanian, M.; Wilson, C. W. Multi-objective genetic algorithms optimisation for calculating the reaction rate coefficients for hydrogen combustion. Ind. Eng. Chem. Res. 2003, 42, 1215. (16) Elliott, L.; Ingham, D. B.; Kyne, A. G.; Mera, N. S.; Pourkashanian, M.; Wilson, C. W. Genetic algorithms for optimisation of chemical kinetics reaction mechanisms. Prog. Energy Combust. Sci. 2004, 30, 297. (17) Smith, G. P.; Golden, D. M.; Frenklach, M.; Moriarty, N. W.; Eiteneer, B.; Goldenberg, M.; Bowman, C. T.; Hanson, R. K.; Song, S.; Gardiner, W. C.; Lissianski, V. V.; Qin, Z. http:// www.me.berkeley.edu/gri-mech. (18) 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, Sandia, NM, 1985. (19) Kee, R. J.; Miller, J. A.; Jefferson, T. H. CHEMKIN: A general-purpose, problem-independent, transportable, FORTRAN chemical kinetics code package. Sandia Report SAND80-8003; Sandia National Laboratories, Sandia, NM, 1980. (20) Lutz, A. E.; Kee, R. J.; Miller, J. A. SENKIN: A FORTRAN program for predicting homogeneous gas-phase chemical kinetics with sensitivity analysis. Sandia Report SAND87-8248; Sandia National Laboratories, Sandia, NM, 1987.
(21) 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. (22) Bernstein, J. S.; Fein, A.; Choi, J. B.; Cool, T. A.; Sausa, R. C.; Howard, S. L.; Locke, R. j.; Miziolek, A. W. Laser-based flame species’ profile measurements-a comparison with flame model predictions. Combust. Flame 1993, 92, 85. (23) Michalewicz, Z. Genetic Algorithms/Data Structures/ Evolution Programs, 3rd ed.; Springer, Berlin, 1996. (24) Tsuboi, T.; Wagner, H. G. Homogeneous thermal oxidation of methane in reflected shock waves. Fifteenth Symposium (International) on Combustion; The Combustion Institute, Pittsburgh, PA, 1974; p 883. (25) Taylor, S. C. Burning velocity and the influence of flame stretch. Ph.D. Thesis, University of Leeds, 1991. (26) Vagelopoulos, C. M.; Egolfopoulos, F. N.; Law, C. K. Further considerations on the determination of laminar flame speeds with the counterflow twin-flame technique. Twenty-fifth Symposium (International) on Combustion; The Combustion Institute, Pittsburgh, PA, 1994; p 1341.
Received for review July 6, 2004 Revised manuscript received October 7, 2004 Accepted October 26, 2004 IE049409D
