Kinetic Parameters Estimation for Three Way ... - ACS Publications

ARTICLE pubs.acs.org/IECR

Kinetic Parameters Estimation for Three Way Catalyst Modeling Karthik Ramanathan*,† and Chander Shekhar Sharma‡,§ † ‡

India Science Lab, General Motors Global Research and Development Center, Bangalore, Karnataka 560066, India General Motors Powertrain, Bangalore, Karnataka 560066, India ABSTRACT: One of the critical needs of a Three Way Catalyst (TWC) model is to be able to predict light-off. This is crucial for application studies and vehicle architectural studies because most of the emissions from a TWC occur before light-off (called coldstart emissions). Laboratory experiments give detailed insights to the reaction mechanism and analytical forms of the rate expressions as they are well-controlled and well-behaved as compared to vehicle tests. However, to predict emissions on a vehicle test, the laboratory-estimated kinetic parameters are not entirely capable because of the various uncertainties in the vehicle tests. In this work, six different vehicle data sets are used to calibrate and validate the TWC global kinetic model. Our emphasis in this work is restricted to predicting the light-off (cold-start emissions) in TWC. The kinetic model is calibrated using 4 vehicle data sets (which use the FTP drive cycle) using iSIGHT software package. The kinetic parameters of the various reactions occurring in the TWC are estimated to match the experimental data through exploratory and local optimization methods. A systematic approach (with increasing complexity) is used to estimate the kinetic parameters. The estimated parameters are then used to validate the model on two different vehicle data sets (one NEDC drive cycle and one FTP drive cycle) with different catalyst compositions and engine power (and hence different engine out exhaust compositions). The model with estimated kinetic parameters predicts the light-off reasonably well for the new data sets. The parameter estimation approach in this work is kept as generic as possible to exhaust aftertreatment devices, and a set of guidelines for parameter estimation (specifically for use in exhaust aftertreatment devices) is presented (in the Appendix).

’ INTRODUCTION AND LITERATURE REVIEW Catalytic monoliths are one of the most widely used forms of catalytic reactors. One of the major advantages of such reactors over traditional packed bed reactors is the lower pressure drops that they afford for the high flow rates that are encountered in many industrial applications. One of the most common applications of a monolith reactor is in the control of automobile emissions (CO, hydrocarbons, and NOx). The monolith reactor is a unitary structure usually made of a ceramic material and consisting of a large number of thin-walled, narrow, parallel flow channels that makes up the substrate. The catalyst is deposited on the walls of the narrow channels either as a thin layer or as a porous washcoat which is supported on the monolith substrate. The gas that enters the reactor is transported by convection and diffusion in parallel flow streams through these narrow channels and reactions occur on the active catalytic surfaces present in the washcoat. There is extensive literature dealing with experimental, modeling, and simulation results of catalytic monoliths. The various review articles/books113 summarize the recent progress in catalytic monolith reactors. Catalytic monoliths (or catalytic converters) used in automobile emission control operate under highly transient conditions, and their performance during the warm-up period is an important factor in the design of catalytic converters. Because of the time-varying inlet conditions and the complex coupling of the physical and chemical processes in the monolith, improving the design requires a comprehensive understanding of the influence of various design and operating variables on the performance of the monolith. Mathematical modeling and simulation are useful in identifying the optimal designs and reducing the amount of experimentation needed to achieve this goal. Mathematical r 2011 American Chemical Society

models of varying degrees of complexity (in terms of the geometric dimensionality of the model, flow modeling, washcoat diffusion modeling, and chemical kinetics modeling) have been in use1427 over the last 23 decades. In particular, the onedimensional two-phase (gassolid) model,17,28 which differentiates between the gas and solid phase concentrations and temperatures has been found to describe these processes in a monolith reactor quite satisfactorily.29,30 One of the main challenges in modeling catalytic converters is in obtaining the right kinetic mechanisms and expressions for the various heterogeneous reactions (surface reactions on the active metal surface) occurring in the exhaust aftertreatment device. Using the right kinetic model is very crucial to improve the predictability of the model (particularly predicting the cold-start emissions and light-off). The simplest approach to kinetic modeling is the look-up table approach,31 where the model is populated with a database of conversions for various species as a function of temperature and space velocity, from which conversions can be predicted by interpolation. Adapting to another catalyst formulation is relatively quick as it just requires measuring new conversion curves. These models are very fast in terms of computational time. However, since it does not describe the underlying chemical processes, the range of applicability is limited. On the other extreme are microkinetic models where all intermediate reaction steps are accounted for, and rate expressions for these reaction steps are proposed and the various Received: April 6, 2011 Accepted: June 28, 2011 Revised: June 9, 2011 Published: July 21, 2011 9960

dx.doi.org/10.1021/ie200726j | Ind. Eng. Chem. Res. 2011, 50, 9960–9979

Industrial & Engineering Chemistry Research reaction rate parameters are obtained from experiments16,3235 or theoretical studies.36 These models are computationally intensive and hence cannot be easily ported to three-dimensional models and reduced order models for control applications which requires a small set of reactions/species. Additionally, it usually take a lot of time (years) to develop these models. With the limitations in the measurements of various species compositions in a real world engine data (i.e., on-vehicle), obtaining kinetic mechanisms at the microkinetic level of detail is also limited. A middle path is often chosen by modelers where a single global rate expression is written for each reaction. Toward this goal, one of the practical and increasingly used17,23,25,3741 approaches to determine the reaction rates is to use a class of global kinetic models with LangmuirHinshelwood expressions. These models typically have good predictive capability of the model and help in rapid design evaluations and also help in control algorithms. The LangmuirHinshelwood type of kinetics used in automotive catalytic converters has its basis in the pioneering work of Voltz et al.,42 who studied the kinetics of oxidation of CO and C3H6 on Pt/Alumina catalyst and reported the inhibition due to CO, C3H6, and NO. However, they considered NO to be a nonreacting species. This form of the kinetics has been extended to Pd/Rh catalysts and also to include H2 oxidation reactions. This category of models is semiempirical which involve tunable parameters, and all phenomena that are not explicitly accounted for (like adsorption, desorption, etc.) are lumped into the values of the tunable parameters in the kinetic expressions. Laboratory experiments give detailed insights to the reaction mechanism and analytical forms of the rate expressions as they are well-controlled and well-behaved. However, these estimated kinetic parameters do not translate directly to engine/vehicle scale data because of the following reasons.4345 First, kinetic parameters are very sensitive to washcoat and catalyst formulations. Second, catalytic processes suffer from deactivation and aging which depend on the conditions at which it is aged and typically in automotive exhaust catalyst it is very difficult to get to the same level of aging even if we start with identical washcoat formulations. To add to this complexity, the aging/deactivation of the catalyst varies along the converter length as well, which means that the number of active sites would vary along the converter length. Third, the reaction path depends on the possible presence of species in real exhaust gas that are not present in controlled experimental studies (for example, different hydrocarbons). Additionally, water and CO2 play a promoting and inhibiting role for a few of the reactions, and if this is not considered in controlled experimental studies, the estimated kinetic parameters cannot be generalized. Fourthly, the kinetic parameters are very sensitive to the feed composition and temperature, which are very transient in typical engine exhaust and varies with engine type, velocity, outdoor temperatures, engine load, the drive cycle, and the variations due to the type of fuel and its additives. In short, the reaction kinetics in a catalytic converter may be affected by the engine exhaust gas, the washcoat formulation, and the aging history of the catalyst. This necessitates the need to validate/tune the kinetic mechanism and parameters using complete exhaust compositions and engine scale data. The kinetic parameters estimated from laboratory scale experiments or from literature are used as a starting point for such a validation/tuning. Parameter estimation using optimization methods (that minimize the error between the measured/experimental and simulated data) is very useful to estimate the kinetic parameters.

ARTICLE

One of the first reported computer aided tuning/estimation of kinetic parameters in the literature in the context of the Three Way Catalytic converter (TWC) based on the method of conjugate gradients was by Montreuil et al.46 They generated an experimental database of steady state conversion efficiencies of CO, NO, C3H8, C3H6, H2, and O2 for different inlet feed compositions and temperatures over Pt/Rh and Pd/Rh catalysts. Using a model similar to that of Oh et al.,17 they considered 13 different reactions with kinetic expressions of the LangmuirHinshelwood form and estimated 95 parameters using a multidimensional conjugate gradient optimization method. However, this work was based on data obtained under controlled steady-state conditions in an experimental lab scale setup. Dubien and Schweich47 proposed a theoretical methodology to determine the pre-exponential factor and activation energy of simple rate expressions from light-off experiments and compared these values to parameters estimated with an optimization method (Downhill Simplex method). Their theoretical approach, however, is applicable only for noncompeting reactions. And hence could not be used directly for parameter estimation for a realistic catalytic converter where many competing reactions take place. The transient behavior for a real automobile exhaust system has been analyzed by Pontikakis and Stamatelos.48 They used a two-dimensional model for the solid phase energy balance in the channel. They tuned about 12 kinetic parameters using conjugate gradient based optimization method coupled with a Bates and Watts transformation49 to convert a constrained optimization problem into an unconstrained one. Pontikakis et al.44 used a slightly different approach to tune about 12 parameters including adsorption parameters using conjugate gradient based optimization method coupled with a penalty function approach to handle constraints. The conjugate gradient method though accurate and fast has a tendency to converge to a minimum close to the initial guess. However, in systems where there are multiple kinetic parameters (and where multiparameter optimization is performed), the system can possess multiple minima. It is hence necessary to obtain a global minimum to get the lowest objective function (for physically realistic kinetic parameters). Using exploratory/ evolutionary optimization algorithms like the Genetic Algorithm (GA), Adaptive Simulated Annealing, and Downhill Simplex can help to obtain a global minimum. One of the early works in the literature to use GA in the context of this problem was by Glielmo and Santini.50 They developed a fast, two time scale numerical algorithm to solve the model equations. The reactions they considered include oxidation of CO, CH4, C3H6, H2, and reduction of NOx. They tuned a total of 7 parameters simultaneously. Pontikakis et al.51 used a Genetic Algorithm based method to tune around 10 parameters. In summary, the two methods of optimization used in the literature for kinetic parameter estimation are the conjugate gradient method and the genetic algorithm (local and exploratory methods, respectively). In this work, we use global kinetics to represent the reactions in a TWC and estimate the kinetic parameters for these reactions from engine/vehicle drive cycle data using both exploratory and local optimization methods and validate them on different vehicle driving cycles and different catalyst compositions. Our emphasis in this work is restricted to predicting the light-off (cold-start emissions) in TWC. The parameter estimation approach in this work is kept as generic as possible to exhaust 9961

dx.doi.org/10.1021/ie200726j |Ind. Eng. Chem. Res. 2011, 50, 9960–9979

Industrial & Engineering Chemistry Research

ARTICLE

Initial conditions.

aftertreatment devices and a set of guidelines for parameter estimation (specifically for use in exhaust aftertreatment devices) is presented (in the Appendix).

’ MATHEMATICAL MODEL A single-channel (i.e., channel-to-channel variations and interactions are ignored) one-dimensional model is used for describing the monolith converter.17,43 The axially varying gas phase temperature, velocity, and reactant concentrations are to be interpreted as cross-sectional averages within each channel. Axial diffusion of mass and heat is also assumed to be negligible in the gas phase. (Axial Peclet number, defined as the ratio of axial diffusion time to the axial convection (residence) time, is .1000 for both heat and mass transport, indicating dominance of the convective heat and mass transfer). The wall of the individual channel is coated with a porous, catalytically active layer (i.e., washcoat) of uniform thickness, and it is also assumed that the washcoat is sufficiently thin (typically 50 μm or less) so that its curvature is unimportant and diffusion resistance within the washcoat can be neglected.16,52,53 The governing model equations derived under these assumptions are given below: Solid-phase energy balance.

ðfsb Fsb cp, sb þ fwc Fwc cp, wc Þ ¼ fsb

∂ ∂Ts λsb ∂z ∂z

∂Ts ∂t

nrxn

∑ aj ðzÞðΔHÞjRj

j¼1

þ hSðTg Ts Þ

Ts ðz, t ¼ 0Þ ¼ Tinit ðzÞ θk ðt ¼ 0Þ ¼ θk, init ðzÞ k ¼ 1, 3 3 3 , ndfq

Here, Ts is the solid phase temperature, Tg is the gas temperature, and xg and xs are the gas and solid phase mole fractions of species ‘i’. Here, nrxn represents the number of reactions, and nsp represents the number of reacting species. The subscripts ‘wc’ and ‘sb’ represent the washcoat and the substrate, respectively. The stoichiometric coefficients are given by si,j. In the above equations, θk are the coverage parameters, and Fk is a function of the solid temperature, solid phase mole fractions, and the coverage parameters. Depending on the system under study, the coverage parameters are governed by differential equations or algebraic equations. Here, ndfq represents the number of surface coverage parameters determined by differential equations, and nsrf represents the total number of surface coverage parameters. In the above equations, S represents the geometric surface area per converter volume. These, as well as all other variables, are defined in the Nomenclature section. The solid-phase energy balance takes into account the axial heat conduction (term 1 in the right-hand side (r.h.s.)), the heat generated by reactions (term 2), and heat transfer between the gas and the solid (term 3). The boundary condition for the solidphase energy balance is given by eq 7. The heat and mass transfer coefficients in the above model equations are determined by the following equations and correlations

ð1Þ

h¼ km, i

Gas-phase energy balance.

∂Tg w cp, g ¼ hSðTs Tg Þ A ∂z

ð2Þ

nrxn

∑ ajðzÞRj si, j j¼1

ð4Þ

Fk ðTs , xs , θÞ ¼ 0

ð5Þ

Boundary conditions.

xg, i ðz ¼ 0, tÞ ¼ xi, in ðtÞ Tg ðz ¼ 0, tÞ ¼ Tin ðtÞ

λsb

∂Ts ¼0 ∂z

@z ¼ 0, L

ð6Þ

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 1 3:85 10 Tg þ Mi MN2 ¼ pffiffiffiffiffiffiffiffi pffiffiffiffi ð 3 Σi þ 3 ΣN2 Þ2 5

cg Di, m

dθk ¼ Fk ðTs , xs , θÞ k ¼ 1, 3 3 3 , ndfq dt

ð9Þ

λg ¼ 2:66 104 Tg0:805

ð3Þ

Surface coverage equations.

k ¼ ndfq þ 1, 3 3 3 , nsrf

Nu λg Dh Sh ¼ ðcg Di, m Þ Dh

where the thermal conductivity of the gas (λg) and the mass 54,55 diffusivity (Di,m) are given by

Species mass balance in gas phase and catalyst surface (for i = 1,...,nsp).

w ∂xg, i ¼ km, i Sðxg, i xs, i Þ ¼ A ∂z

ð8Þ

0:75

ð10Þ

Here, Mi represents the molecular weight of the species in grams, and Σi represents the diffusion volume of the species (obtained from Table 11-1 in The Properties of Gases and Liquids56) . We assume asymptotic (constant) Nusselt number (Nu) for simulation purposes. (Note that Nu = Sh.) The channel cross-section is assumed to be close to a square geometry, and correspondingly a value of Nu = 3 is chosen.57 Proper prediction of warm-up behavior requires an accurate, temperature-dependent specific heat for the monolith substrate and the following functional form, which was obtained from the handbook58 (page 392) for cordierite, is being used in this work cp, sb ðTs Þ ¼ 1:071 103 þ 0:156Ts

ð7Þ 9962

3:435 107 Ts 2


Industrial & Engineering Chemistry Research Table 1. Main Reactions in TWC s. no.

reactions

heat of reaction in J/mol

Oxidation Reactions 1

CO + 0.5O2 f CO2

2.83 105

2

C3H6 + 4.5O2 f 3CO2 + 3H2O

1.93 106

3

C3H8 + 5O2 f 3CO2 + 4H2O

2.04 106

4

H2 + 0.5O2 f H2O

2.42 105

5

CO + NO f CO2 + 0.5 N2

3.73 105

6

C3H6 + 9NO f 3CO2 + 3H2O + 4.5N2

2.74 106

7

H2 + NO f H2O + 0.5N2

3.32 105

8

CO + H2O T CO2 + H2

4.12 104

9

C3H6 + 3H2O f 3CO + 6H2

+3.74 105

NO Reduction Reactions

Water-Gas and Steam Reforming Reactions

Ceria Reactions (Oxygen Storage) 10 2Ce2O3 + O2 f 4CeO2

2 105

11 Ce2O3 + NO f 2CeO2 + 0.5N2

1.90 105

12 CO + 2CeO2 f Ce2O3 + CO2

1.83 105

13 C3H6 + 12CeO2 f 6Ce2O3 + 3CO + 3H2O

4.77 105

14 C3H8 + 14CeO2 f 7Ce2O3 + 3CO + 4H2O 15 H2 + 2CeO2 f Ce2O3 + H2O

4.95 105 1.42 105

ARTICLE

metal loading. Kang et al.61 derived an analytical expression for the activity function (active site density as a function of aging or miles driven) and showed that a global reaction kinetic model coupled with this activity function is capable of describing the alteration of TWC performance with respect to the catalyst mileage, regardless of the catalyst metal content. However, it should be noted that it is possible that the reaction mechanism changes when the catalyst formulation or aging behavior is changed. Hence, in cases where we are strongly led to believe that reaction mechanism changes and a different form of global kinetics should be used, the kinetic expressions and parameters need to be re-estimated. In this work, the kinetic parameters are obtained for aged catalysts (around 100K to 120K miles) with similar washcoat/catalyst compositions, and the same set of kinetic parameters is used to predict any new vehicle data and/ or new catalyst compositions. In cases where experimental and simulation results do not match a fine-tuning of specific parameters is done. The specific reaction rate expressions for the various reactions are listed below. The rate expressions are expressed in terms of the species concentrations in the solid phase, which are related to the mole fractions as p xs, i Cs, i ¼ Rg T s Oxidation Reactions

1. 2. 3. 4.

’ TWC REACTIONS AND KINETICS The set of reactions considered in this work for modeling a TWC is listed in Table 1 along with the respective heat of reactions. For this set of reactions, the number of surface coverage parameter is 1 i.e., nsrf = 1 (with ndfq = 1) and the coverage (or the oxidation) parameter, θ, is the ratio of the oxygen storage component (e.g., Cerium) in the higher oxidation state to the total available oxygen storage capacity, as given by θ¼

CO + 0.5O2 f CO2 C3H6 + 4.5O2 f 3CO2 + 3H2O C3H8 + 5O2 f 3CO2 + 4H2O H2 + 0.5O2 f H2O ! k1 CCO CO2

R1 ¼

moles of CeO2 moles of CeO2 þ 2 moles of Ce2 O3

E1 exp Rg T s G

k2 CC3 H6 CO2 R2 ¼

The specific reaction rates or the turnover rate (Rj) for the jth reaction is expressed in mol/mol-site/s, where mol-site refers to the active sites which actually participate in the catalytic reactions.The number of active sites (aj) depends on the active metal loading level and dispersion. The specific reaction rate is multiplied by the (active) site density (aj) in mol-site/m3 to obtain the total reaction rate in mol/m3/s. The advantage of writing the rate expression as a product of the specific reaction rate and the number of active sites are the following. First, since the specific reaction rates are not dependent on the amount of catalyst, they need not (ideally) be changed when switching over from one catalyst formulation to another. Second, the aging behavior is usually incorporated through the number of active sites available, and we would not need to estimate the parameters for the specific reaction rates from scratch for a catalyst with different aging level/history. It has been demonstrated by Hauff et al.59 that only one parameter, the active site density (determined by CO chemisorption), needs to be changed in the global kinetic models as the catalyst is aged or a different catalyst loading is used. Boll et al.60 used a microkinetic model and showed that only the active site density can serve as parameter to model the varying noble

E2 exp Rg T s

!

G ! k3 CC3 H8 CO2

R3 ¼

E3 exp Rg T s G !

k4 CH2 CO2 R4 ¼

E4 exp Rg Ts G

G ¼ ð1 þ K1 CCO þ k2 CC3 H6 Þ2 ð1 þ K3 C2CO C2C3 H6 Þð1 þ K4 CNO Þ The terms in the inhibition factor G are given by Ei, G Ki ¼ ki, G exp i ¼ 1, 3 3 3 , 4 Ts 9963


Industrial & Engineering Chemistry Research NO Reduction Reactions

5. 6. 7.

CO + NO f CO2 + 0.5 N2 C3H6 + 9NO f 3CO2 + 3H2O + 4.5N2 H2 + NO f H2O + 0.5N2

ARTICLE

reaction of propene (or fast hydrocarbon) can be computed (using data from Chemical Properties Handbook62) as ΔGSRF = 3.783 105 546.1Ts 3.768 102T2s ]. Ceria Reactions (Oxygen Storage) 10. 11. 12. 13. 14. 15.

! k5 CCO CNO R5 ¼

E5 exp Rg Ts G

k6 CC3 H6 CNO R6 ¼

E6 exp Rg Ts

R10

G k7 CH2 CNO exp

R7 ¼

!

E7 Rg T s

!

R11

G R12

The inhibition factor G is given by the same expression as used for the oxidation reactions. Water-Gas and Steam Reforming Reactions

8. 9.

R13

CO + H2O T CO2 + H2 C3H6 + 3H2O f 3CO + 6H2

The reaction free energy of the Water-Gas Shift reaction is computed (using data from Chemical Properties Handbook62) as

R14

ΔGWGS ¼ 4:1034 104 þ 44:19Ts 5:553 103 Ts2 Here the subscript WGS represents Water Gas Shift (WGS) reaction. And the equilibrium constant is given by ! ΔGWGS KWGS ¼ exp Rg Ts The reaction rates are given by ! k8 E8 CH2 CCO2 R8 ¼ exp CCO CH2 O G Rg Ts KWGS ! k9 E9 R9 ¼ exp CC3 H6 CH2 O G Rg Ts Here again, the inhibition factor G is given by the same expression as mentioned in the oxidation reactions. It should be mentioned that the Steam Reforming reaction (reaction 9) of the fast hydrocarbon (SRF) is usually reversible. But it was found that, though the equilibrium favors the reverse reaction at low temperatures for the steam reforming reaction, both the forward and reverse reaction rates were lower at low temperatures. At high temperatures (>473 K), it was found that the equilibrium favors the forward reaction significantly, and hence the reverse reaction rate was low. Since the reverse reaction rate is low/small at both low as well as high temperatures, we have neglected the reverse reaction in this work. [Note: The Gibbs free energy of formation for the steam reforming

R15

2Ce2O3 + O2 f 4CeO2 Ce2O3 + NO f 2CeO2 + 0.5N2 CO + 2CeO2 f Ce2O3 + CO2 C3H6 + 12CeO2 f 6Ce2O3 + 3CO + 3H2O C3H8 + 14CeO2 f 7Ce2O3 + 3CO + 4H2O H2 + 2CeO2 f Ce2O3 + H2O ! E10 ¼ k10 exp CO2 ð1 θÞ Rg T s ! E11 ¼ k11 exp CNO ð1 θÞ Rg T s ! E12 ¼ k12 exp CCO θ Rg T s ! E13 ¼ k13 exp CC 3 H 6 θ Rg T s ! E14 ¼ k14 exp CC 3 H 8 θ Rg T s ! E15 CH 2 θ ¼ k15 exp Rg T s

The equation for the surface coverage parameter (θ) can be written as dθ ¼ FðTs , cs , θÞ dt ¼ ð4R10 þ 2R11 Þ ð2R12 þ 12R13 þ 14R14 þ 2R15 Þ By default, we take all the Cerium in the TWC to be initially in the oxidized state (more stable) so that θ = 1 at t = 0. The SI units for the various kinetic parameters are given in Table 2.

’ EXPERIMENTAL DATA FOR ESTIMATION AND VALIDATION In this work, we use engine/vehicle data for kinetic parameter estimation and validation purposes. The objective is to use a combination of evolutionary and local optimization methods to estimate the kinetic parameters (by minimizing the error between the measured and simulated data) using a few data sets and predict the remaining data sets using the estimated kinetic parameters. As mentioned earlier, the emphasis in this work is restricted to predicting the first drive cycle (125 s), particularly light-off (with more emphasis on CO and HC cold-start emissions), as this is the crucial period of cold-start emissions. Measurements from six engine/vehicle experiments performed in-house by GM-Powertrain are considered in this work. 9964



ARTICLE

Table 2. SI Units for the Various Reaction Kinetic Parameters kinetic parameter

SI units !n mol reac m3 mol site :sec mol ( 1 for i ¼ 1, 3 3 3 , 9 where n = 0 for i ¼ 10, 3 3 3 , 15

ki

Ei

J/mol

!n ( m3 1 for i ¼ 1, 2, 4 where n = 4 for i ¼ 3 mol

ki,G K

Ei,G

Table 3. Details of the Different Experimental Data data set detail, vehicle and drive cycle

abbreviation

used for

GMT 900, FTP drive cycle, set 1 GMT 900, FTP drive cycle, set 2 GMX 001, FTP drive cycle, set 1 GMX 001, FTP drive cycle, set 2 GMT 800, FTP drive cycle LY7, Euro drive cycle

GMT900-D1 GMT900-D2 GMX001-D1 GMX001-D2 GMT800 LY7

calibration calibration calibration calibration validation validation

Among the six data sets, we use four for calibration (parameter estimation) purposes and two for validation purposes. Note that among the data sets used, one of them was driven according to the European drive cycle, whereas the rest of them had the FTP (Federal Test Procedure) drive cycle. Table 5 gives the details of these data sets, and Table 4 gives the details of the catalyst (active metal) and the geometric dimensions of the converter for all the data sets. For the vehicle data considered here, the three way catalysts were close-coupled to the engine and consisted of two bricks (front and rear brick) and had slightly different catalyst formulations (in terms of the quantity of active metal loading). All the catalysts were 100K miles or more aged. It is usually observed that for an aged catalyst (because of catalyst deactivation due to thermal aging and poisoning) the catalyst dispersion near the inlet of the front brick is usually very low (almost zero), and hence we assume that the first 0.25 in. of the front brick has zero catalyst dispersion.63 This means that there is a catalyst dispersion profile along the converter length for the front brick. The rear brick is assumed to have constant dispersion value.63 Note that the catalyst dispersion value changes only aj’s and not the Rj’s. This is because the kinetics are expressed as turnover numbers (i.e., in terms of mol/mol-site-catalyst/s) and that the active site density (mol-site-catalyst per unit converter volume) is multiplied to this reaction rate. Hence, as we switch between the different data sets, Rj remains the same and only aj is different. Typical engine-out measurements have the mole fractions of the exhaust gases (CO, NO, HC, CO2, O2), the gas temperature, and the mass flow rate. However, for modeling purposes the mole fractions of H2 and H2O are also required, and these are calculated using the formulas64,65 ½H2 ¼

n ½COð½CO þ ½CO2 Þ 2 ½CO þ ½CO2 =Kw

½H2 O ¼

n ½CO2 ð½CO þ ½CO2 Þ 2 Kw ½CO þ ½CO2

The typical values of Kw and n are taken to be 0.285 and 1.85, respectively. This is a standard formula used in the literature and is obtained using Water Gas Shift Equilibrium assumption inside the engine. The tailpipe-out measurements are typically limited to CO, CO2, NO, HC, and sometimes O2. In most cases, the O2 analyzer does not produce very reliable measurements, and hence we do not use them. Also, the tailpipe-out exhaust gas temperature is not measured. In most of the data sets, the measured engine-out (EO) (converter inlet) and the measured tailpipe-out (TP) (converter outlet) were not in perfect synchronization (both in terms of time alignment and the absolute value measured) during the first few seconds. This is illustrated by the instantaneous emission measurement plots for CO and HC in Figure 1 (for GMX001-D1). This lack of synchronization seems to be a measurement issue as it cannot be explained by any physical process happening in the system. [Note: The only storage component present in the system is Ceria and its purpose is not to store and release CO/HC before light-off temperatures.] To address these measurement issues, vertical shifts have been introduced in the corresponding cumulative emission values (shown in Figure 1) so that the measured TP data matches with the EO data for the first 10 s (before any significant conversion is seen in the catalytic converter). This is to say that during the first 10 s when no significant conversion happens in the catalytic converter, the EO and TP values should match. [Note that if there were only time delays, a horizontal shift in the instantaneous data would suffice, but since the absolute magnitude of the values was also different a vertical shift in the cumulative emission values was necessary]. The vertical shift is established by visual examination and is shown in Figure 1. The numerical details of the vertical shifts for the various data sets are given in Table 5.

’ PARAMETERS TO BE ESTIMATED For the 15 reactions listed earlier in Table 1, there are a large number of kinetic parameters that could be chosen for parameter estimation purposes a total of 38 parameters (a total of 19 rate constants including the ones that appear in the inhibition term and each rate constant is determined by two parameters the pre-exponential factor and the activation energy). A blind-folded approach to determine all 38 parameters at once would be fruitless, and even a powerful and smart optimizer will not be able to arrive at the optimum and, additionally, optimizing 38 parameters at once would be computationally expensive. It is important that we use the understanding of the system to decide on the parameters that need to be estimated. The vehicle data sets used to estimate the parameters have output (or tailpipe) measurements only for NO, CO, and total hydrocarbons. But not all kinetic parameters in the model are directly (or first-order) sensitive to these measurements. Only the kinetic parameters that have some significant sensitivity to these measurements should be chosen for estimation purposes. (The ideal way to determine the sensitivities of the various parameters to the measurements is to do a systematic time-dependent sensitivity analysis.) For the current work, the parameters for estimation purposes were chosen based on the sensitivity information available in the literature and the authors’ understanding of the system. The rate constant is dependent on both the pre-exponential and activation energy, and they are not completely independent (statistically) and are correlated. Because of this correlation, simultaneous estimation of both the parameters can get difficult, unless there are controlled experiments (e.g., experiments at 9965



ARTICLE

Table 4. Catalyst and Geometric Properties of the TWC System for Various Vehicles GMT-900 front

GMX-001 rear

front

GMT-800 rear

front

rear

LY7 front

rear

length (in cm)

5.08

7.62

5.08

10.16

5.08

10.16

6.35

8.89

frontal area (in cm2)

109.7

109.7

87.75

87.75

109.7

109.7

79.57

79.57 600

cell density (cells per square inch)

600

600

600

600

600

400

600

substrate thickness (milli inches)

4.3

4.3

4.3

4.3

4.3

6.5

3.5

3.5

catalyst loading (in g/ft3)

80

15

188.8

20.8

55.2

17.5

100

20

catalyst loading ratio (Pt:Pd:Rh)

0:77.5:2.5

0:12.5:2.5

0:35.7:1

1.27:0:1

0:1:0

4:0:1

0:97:3

0:17:3

Ceria loading (g/L)

75

75

0

75

0

0

82

82

Figure 1. Lack of synchronization in measured engine out and tailpipe out emissions during the first few seconds [shown in instantaneous and cumulative plots].

different temperatures) that can delineate the effect of preexponential and activation energy. It is typically seen that when the catalyst/washcoat formulations are similar, the activation energy of the reaction remains the same and only the pre-exponential factor changes. Since the catalyst and washcoat formulations are similar across the various data sets, we keep the activation energies fixed and estimate only the pre-exponential factors, and this reduces the number of parameters to estimate to 19. The steam reforming and water gas shift reactions are fairly insensitive to the measured data (particularly during the first drive cycle) and hence we keep them constant (i.e., not estimated), and this reduces the number of parameters to 17. Among these parameters, the oxidation reaction parameters are tuned/estimated initially, followed by the NOx reduction reactions, followed by the Ceria reactions. This is further

explained in the section “Evolution of Estimated Parameters”. This sequential increase of complexity helps to get to the correct ranges of these parameters and also good initial guesses for future optimizations. It is fairly well-known that the oxidation reactions play the crucial part during the light-off period, and the NOx reduction reactions play a crucial part in reducing the NO (after the cold-start). This information is used in the initial estimating process as we increase the complexity of the optimization. We also specify lower and upper bounds to the parameters that are being estimated so that the optimizer does not go beyond these bounds. After a particular optimization run, if any of the parameter hits the upper or lower bound, this information was used to understand the system better and/or the bounds were relaxed so as to let the optimizer search a larger space. Most often 9966



ARTICLE

Table 5. Vertical Shifts Introduced in Measured Tailpipe-out Data vertical shift introduced data set

HC (kg)

CO (kg)

NO (kg)

GMT900 data set1

0.0

0.0

0.0

GMT900 data set2

7 105

0.0

0.0

GMTX001 data set1 GMTX001 data set2

1.3 104 7 105

7 104 0.0

1.11 105 3.2 106

GMT 800

0.0

0.0

0.0

when a parameter hits the bound, that parameter could be seen as one of the fairly sensitive parameters.

’ OBJECTIVE FUNCTIONS The definition of the objective function is the key to successful parameter optimization. There are various ways to define the objective function. In this work, the norm of the relative errors between cumulative predicted and cumulative measured tailpipeout emissions has been chosen as the objective function for optimization. [It should be mentioned that the exhaust gas temperature values were not chosen to be a part of the objective function as they were not measured.] Cumulative values are chosen for the objective function instead of the instantaneous values because it is found that cumulative values reflect better sensitivity to the kinetic parameters.43 By definition, cumulative values at any time instant ‘t’ contain all the information that has happened until ‘t’. So using cumulative values for the objective function can also be looked upon as giving more weights to the initial few points (particularly the light-off) in time. This causes cumulative species concentrations to be more sensitive toward kinetic parameters as catalyst light-off is mainly driven by kinetics of underlying chemical reactions. More details on the comparison in performance of the optimization process by using cumulative objective function and instantaneous objective functions can be found here.43 Choosing objective functions which are not sensitive to the parameters to be estimated can lead to ill-posed optimization/ estimation problem, and hence it is important to choose the right values and weights for the objective function. For example, trying to estimate kinetic parameters using data in the region where the reaction kinetics is very slow (or negligible close to zero conversion) or where the reactor is transport limited (i.e., kinetics is very fast close to higher or 100% conversion) would be ineffective. In typical vehicle data measurements, the confidence level for the measured data during the initial part of the FTP cycle is low and during this time, there is very little conversion (or reaction) and hence it is best to avoid these data points while calculating the objective function. A conversion efficiency based on the measured data is defined and used to automatically identify the data points that need to be included (or discarded) in calculating the objective function. Conversion efficiency is defined on the basis of the cumulative input and output and is given, for species j, at time t by ! 1 0 Z t

xg, j wMj dt B C B C 0 B C outlet ! C 100 ηj ðtÞ ¼ B1 Z t B C @ A xg, j wMj dt 0

inlet

Figure 2. Plot of the conversion efficiency and calculating the time from which the objective function needs to be evaluated.

where xg,j is the measured gas-phase mol-fraction of species j at converter inlet/outlet, w is measured instantaneous molar flow rate, and Mj is the molecular weight of species j. By definition, positive conversion efficiencies for any species mean that the measured values at the outlet are lower compared to the measured values at the inlet (i.e., the species is being consumed by reaction), and negative conversion efficiencies mean that the measured values at the outlet of the converter are higher compared to the measured values at the inlet. From the conversion efficiency plot, we define the last time instant after which the conversion efficiency remains positive as the cutoff time. Any data point before this cutoff time is discarded in the calculation of the objective function. We choose the last time instant because there are cases where the conversion efficiencies fluctuate between positive and negative values, and it is only after the last time instant the kinetics become important (because once the kinetics becomes important the conversion usually does not drop). This is illustrated by an example in Figure 2 by showing the cumulative plot for CO and the corresponding conversion efficiency plot for the data from GMT900-D1. In this figure, it is easily seen that the conversion efficiency is negative (measured tailpipe-out is higher than the measured engine-out) for the initial few seconds. The negative conversion efficiency may be due to incorrect measurements and/or some release of stored CO in the converter (before the start of the experiment, possibly from a previous experimental run). But since our simulation/model assumes that there is no stored CO to start with and that the kinetics used in the simulation consumes (and does not produce except for the 9967



ARTICLE

WGS reaction, which typically does not favor CO formation at low temperatures unless the hydrogen concentration is 23 orders of magnitude higher than CO, which is not the case with typical engine exhaust) the reactant, estimating the kinetic parameters using the data from the negative conversion points will be completely ineffective and hence it is best to ignore this data. The cutoff time is marked in the figure as t_CO. In a similar manner, we find the cutoff times for the various species (t_CO, t_HC, t_NO) for each data set. All of the above procedures are performed to ensure that only good and relevant (or valid) data points are used to evaluate the objective function.

The objective function for the ith data set is calculated using one of the formulas (depending upon whether for NO is included or not) given below vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! !2 u u 1 t_final mCO:meas mCO, sim 2 1 t_final mHC, meas mHC, sim t fi ¼ þ ni, CO t_CO ni, HC t_HC mCO, meas mHC, meas

∑

∑

ð11Þ or

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! !2 !2ffi u t_final t_final u 1 t_final mCO, meas mCO, sim 2 1 mHC, meas mHC, sim 1 mNO, meas mNO, sim fi ¼ t þ þ ni, co t_CO ni, HC t_HC ni, NO t_NO mCO, meas mHC, meas mNO, meas

∑

∑

∑

Here mj,meas and mj,sim represent measured and predicted cumulative mass of species j at converter outlet at time t. The cumulative amount of species j which has left the converter until a time instant t is defined as Z t wxg, out, j ðtÞMj dt ð13Þ mj ðtÞ ¼

ð12Þ

ð17Þ

of the objective function to a specific kinetic parameter). For example, using individual objective functions can help in identifying the parameters to which only to a specific objective function (say NO) is sensitive. Only those specific parameter(s) can then be tuned/estimated to minimize the objective function for NO. However, in a reactor system represented by a set of complex and coupled nonlinear processes, we would rarely find cases where a single kinetic parameter affects only a single objective function. It is important to scale (or normalize) the values used in the computation of the objective function. In this work, by the very definition of the objective function, it is automatically scaled (with the corresponding measured values). If the measured value did not appear in the denominator in the definition of the objective function, then each and every data set and/or species would contribute in different proportions to the objective function value. For example, NO values are of the order of 10, whereas CO values are of the order of 100 and hence CO values will mask the NO values while calculating the objective function. Though using an objective function without scaling the outputs for optimization purposes may lead to a minimum value of the objective function, it will not give a good match between experiments and measured value for some specific species and/or data sets. Hence, it is very important to have each and every value that contributes to the objective function to be of order 1 and the total objective function is also of the same order. Another important aspect of scaling is to scale the objective function with the number of data points for each species, and this is done through ni,CO, ni, HC , and n i,NO . This is important because each species will have their own cutoff time, and hence each species will contribute a different number of data points to the objective function. If the objective function is not scaled with the number of data points for each species, then the objective function will be biased toward the species having more data points compared to the rest. Hence, scaling with the number of data points is also as important as scaling the output variables, and these have been implemented in the objective function calculations.

Since it is easy to track and record the parameters and the corresponding objective function(s) during the optimization process (i.e., during the several iterations), defining individual objective functions can also help in getting sensitivity information (sensitivity

’ OPTIMIZATION ROUTINES/ALGORITHMS The software package iSIGHT (version 9.0) is used for optimization purposes. iSIGHT was preferred to MATLAB (release 2007b from Mathworks) for optimization purposes

0

where Mj and xg,out,j represent the molecular weight of species ‘j’ and mole fraction of species ‘j’ leaving the converter, respectively. In the above equations, ni,CO, ni,HC, and ni,NO represent the number of data points for CO, HC, and NO, respectively. As mentioned earlier, we consider only the first 125 s of the drive cycle (excluding, as explained before, the initial few seconds when the conversion efficiency is negative). The combined or the total objective function for all the data sets is obtained as rffiffiffiffiffiffiffiffiffi ð14Þ fi2 f ¼

∑i

The total objective function is used for minimization purposes by the optimization process. However, as we progressed further in estimating the kinetic parameters, it was found that splitting the objective function into three separate objective functions (one each for CO, HC, and NO) leads to increased efficiency of the optimizer. In such a case, all the three objective functions were minimized simultaneously. These individual objective functions are given by vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !ffi u u 1 t_final mCO, meas mCO, sim 2 fco ¼ t ni, CO t_CO mCO, meas

∑

fHC

fNO

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !ffi u u 1 t_final mHC, meas mHC, sim 2 ¼t ni, HC t_HC mHC, meas

∑

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! u u 1 t_final mNO, meas mNO, sim 2 t ¼ ni, NO t_NO mNO, meas

∑

ð15Þ

ð16Þ

9968



ARTICLE

Figure 3. Flowchart of the optimization process and the link between MATLAB, iSIGHT, and FORTRAN.

because iSIGHT offers the convenience of choosing different kinds of optimization algorithms (local and exploratory or evolutionary algorithms) and also offers the flexibility to do two or more optimization algorithms in sequence. Since iSIGHT provides the option of interfacing with MATLAB, the basic TWC model code (in FORTRAN) was modified so that the reaction kinetic parameters (or the optimization parameters) are read from an external text file, and this modified code was then linked to iSIGHT optimizer using MATLAB. The main steps (and transfer of information) that occur during each iteration are as follows: 1. iSIGHT provides the new reaction kinetic parameters to MATLAB. 2. MATLAB creates a text file of these kinetic parameters. 3. MATLAB calls the FORTRAN executable to simulate each data set. The FORTRAN executable reads the kinetic parameters from the text file. 4. The FORTRAN output from each data set is stored in MATLAB. 5. Once the simulation of all the data sets (that are included for estimation purposes) is complete, MATLAB calculates the objective function(s). 6. The calculated objective function(s) is (are) passed to iSIGHT. 7. The old set of kinetic parameters and the objective functions are recorded in an Excel sheet. 8. iSIGHT performs its internal computations and checks if the optimization termination criteria has been satisfied. If it is satisfied, the optimization process ends. If it is not satisfied, iSIGHT generates a new set of kinetic parameters, and the process from step 1 is repeated.A flowchart illustrating the above procedure and the link between iSIGHT, MATLAB, and FORTRAN is shown in Figure 3. Most of the iSIGHT optimization algorithms can broadly be classified into two major categories: (1) local methods and (2) exploratory or evolutionary methods. Local methods generally assume the parameter space is unimodal, convex, and continuous, and these methods search for a local optimum close to the initial guess provided. Exploratory methods avoid focusing only on a local region. They generally evaluate designs throughout the parameter space in search of the global optimum. The success of exploratory methods in finding an optimum in the global

parameter space is heavily dependent on the number of function evaluations and the number of parameters to be determined. The optimum obtained by exploratory methods could be further improved by using this value as an initial guess to a local method. This is because exploratory methods search the entire space, and hence the optimum obtained would not be an absolute optimum on the local scale and a local optimization method could narrow down to the absolute optimum value. In this work, the following exploratory methods are chosen: Multi-Island Genetic Algorithm, Adaptive Simulated Annealing, and Downhill Simplex Optimizer along with the following local methods: Methods of Feasible Directions, Generalized Reduced Gradient, and Sequential Quadratic Programming. iSIGHT allows to incorporate various optimization methods as a series of steps into their execution, and this is called an optimization strategy. Apart from creating an optimization strategy of our own, iSIGHT also provides an optimization method which is a blend of exploratory and local methods. This method is called the Pointer Automatic Optimizer (PAO). The PAO consists of a complementary set of optimization algorithms: linear simplex, sequential quadratic programming (SQP), downhill simplex, and genetic algorithms. PAO can efficiently solve a wide range of problems in a fully automatic manner due to a special automatic control algorithm. As the optimization proceeds, the PAO determines which algorithms are most successful and also decides on the optimal internal control parameters (step size, number of iterations, etc.). More description and working of the PAO can be obtained from the iSIGHT help manual. The optimization strategy used in this work is one of the following: (a) choose one to three of the exploratory methods followed by one to three of the local methods, (b) PAO. Though we find similar performances between the two optimization strategies, PAO helps a nonexpert user significantly. Additionally, the PAO provides the option of terminating the optimization run after specific time duration, and this feature is not available in the other optimization strategy (which terminates on reaching the termination criteria which is usually a constraint on the relative/ absolute step size for the parameters or a constraint on the relative/absolute change in the objective function). Terminating an optimization strategy after a specific time helps in obtaining the best possible results within the available time. This feature is 9969



ARTICLE

Table 6. Optimization Schedules Used for Reaction Kinetics reaction kinetics A1

optimization

Pointer Automatic Optimizer

A2

A3

same as A1

same as A1

manual tuning

algorithm optimization

A4

11 preexponential factors:

same as A1

total: 17

--

parameters k1, k2, k3, k4, k5, k6, k7, k10, k11, k12, k13,

13 preexponential factors (11 as in

k15, and k4 = 10 k1k14 = k13

A1, k4, k14) + k1G, k2G, k3G, k4G

initial values

‘base’ kinetics

objective

CO and HC (single objective function)

function include

‘base’ kinetics

A2

A3

CO, HC, and NO (single

CO, HC, and NO (three different

--

objective function)

objective functions, one each for CO, HC, and NO)

data sets used

GMT900-D1, GMT900-D2, GMX001D1, GMX001-D2

same as A1

particularly useful because it helps to analyze the results after a particular time period and rework the optimization strategy and parameters (and bounds on the parameters) to create a modified optimization run for a longer duration. It should be mentioned that the parameter estimation approach in this work is kept as generic as possible to exhaust aftertreatment devices, and a set of guidelines for parameter estimation (specifically for use in exhaust aftertreatment devices) based on the last 3 sections is presented (in the Appendix).

’ EVOLUTION OF ESTIMATED PARAMATERS This section explains the evolution of the kinetic parameters to obtain 4 major (significant) sets of reaction kinetics (including the final set of kinetics represented here as A4). The evolution of the kinetic parameters and the procedure by which these were obtained are listed in Table 6 and are explained below. Various combinations of optimization algorithms, optimization parameters, and constraints on parameters were used and analyzed before arriving on the final list of parameters, algorithms, and constraints listed in Table 6. The typical number of iterations was of the order of 5000 to 10000 iterations, and the mathematical model runs in about 10 real-time seconds (to simulate around 125 s of the drive cycle, on a PC with Intel Core duo processor with 2 GB RAM) for one data set and hence for an optimization to complete it could take anywhere between 2.5 and 5 days. The starting kinetic parameters for this work were partly taken from the literature66 and were partly obtained from laboratory experiments and is represented as the ‘base’ (or initial) set of kinetic parameters in Table 7. As mentioned earlier, the SI units for the various reaction kinetic parameters are given in Table 2. This base set was modified manually by tuning the pre-exponential factors of HC and CO oxidation reactions (reactions 1, 2, and 3) so as to get a quantitative match between model predictions and measured cumulative HC and CO concentrations at TWC outlet for GMT900-D1 data set. The modified values of reaction kinetics parameters were then used as initial values for optimization using iSIGHT optimizer.

same as A1

GMT900-D1

Optimization process using iSIGHT was started by including only HC and CO in the objective function (defined in eq 11). Table 7 shows the first improved set of reaction kinetics parameters termed as A1. This set was obtained by using PAO method treating 11 pre-exponential factors (k1, k2, k3, k5, k6, k7, k10, k11, k12, k13, and k15) as the parameters to be optimized. The 11 parameters were included in the optimization gradually so that the complexity increases gradually, and Table 7 gives the detail of the final kinetics set. To obtain this kinetics set, reactions 8 and 9 were not included in optimization process and hydrogen oxidation reaction was assumed to be 10 times faster than CO oxidation (k4 = 10k1). In addition, k14 was constrained to be equal to k13 (i.e., k14 = k13) to keep the number of optimization parameters small. This is equivalent to assuming that slow and fast HC reactions with Ceria proceed at equal rates. It should be mentioned that this constraint on k14 was later relaxed (while finding A3 and A4). Figures 4 and 5 show model predictions for the 4 data sets used in optimization. As can be observed, the model performance using A1 kinetics shows significant improvement in HC and CO predictions, but NO predictions worsened for all the 4 data sets. This is expected as NO is not included in the objective function calculations. Figures 6 and 7 present this behavior in terms of absolute percentage error between predicted and measured final cumulative emissions at the end of first drive cycle (first 125 s). In order to improve model predictions for NO, the objective function defined in eq 12 (which includes the NO) was used for optimization purposes. With the same optimization method, and the same set of parameters to be optimized, a new set of reaction kinetics A2 (refer to Table 7) is obtained. As shown in Figures 47, model performance using kinetics A2 show significant improvement in prediction for NO for all the data sets. The model performance also improved for CO and HC for some of the data sets. For the data set GMX001-D2, cumulative CO predictions for A2 and A1 were similar, but cumulative HC predictions for A2 showed an undesirable steady increase after 100 s. The significant improvement in quantitative predictions with kinetics A1 and A2 predominantly came by decreasing the rate of 9970



ARTICLE

Table 7. Sets of Reaction Kinetics Parameters values reaction no.

parameter

base

A1

A2

3.11 10

A3

2.089 10

A4

1

k1

4.21 10

E1

111450

111450

111450

111450

121450

2

k2

1.02 1015

5.88 1013

8.595 1013

1.917 1014

1.917 1015

E2

129530

129530

129530

129530

129530

3

k3

4.68 1010

4.09 1010

7.1 1010

6.404 1010

6.404 1015

E3

165160

165160

165160

165160

165160

4

k4

4.21 1013

3.111 1014

2.089 1015

1.814 1015

1.814 1015

5

E4 k5

111450 1.19 109

111450 1.301 109

111450 1.827 109

111450 6.857 109

111450 2.857 109

E5

52374

52374

52374

52374

52374

6

k6

6.2 1010

5.991 1011

2.648 1010

2.994 1010

2.994 1011

E6

90063

90063

90063

90063

80063

7

k7

6.1 109

4.155 109

4.086 109

7.88 109

7.88 1010

E7

69237

69237

69237

69237

69237

8

k8

1.8 105

1.8 105

1.8 105

1.8 105

1.8 105

9

E8 k9

56720 1.23 105

56720 1.23 105

56720 1.23 105

56720 1.23 105

56720 1.23 105

E9

81920

81920

81920

81920

81920

10

k10

9.87 101

2.37 101

2.014 l01

2.943 10°

2.943 10°

13

13

14

14

2.542 10

5.542 1013

E10

5296

5296

5296

5296

5296

11

k11

4.61 101

1.145 102

2.488 102

7.92 101

7.92 102

E11

25101

25101

25101

25101

25101

12

k12

5.05 10°

1.515 101

1.481 101

1.824 101

1.824 101

13

E12 k13

21768 2.27 103

21768 7.946 102

21768 2.441 103

21768 1.357 103

31768 1.357 101

E13

39070

39070

39070

39070

39070

14

k14

2.25 103

7.946 102

2.441 103

1.77 102

1.77 101

E14

39680

39680

39680

39680

39680

15

k15

5.05 10°

6.927 10°

1.017 101

2.845 10°

2.845 10°

E15

21768

21768

21768

21768

31768

k1,G

7.33 10°

7.33 10°

7.33 10°

4.314 10°

4.314 10°

E1,G k2,G

485 2.57 102

485 2.57 102

485 2.57 102

485 1.289 10°

485 1.289 10°

E2,G

166

166

166

166

166

k3,G

1.8 104

1.8 104

1.8 104

2.147 104

2.147 104

Parameters in Inhibition Term (G)

E3,G

10163

10163

10163

10163

10163

k4,G

3.14 106

3.14 106

3.14 106

8.699 105

8.699 105

E4,G

3685

3685

3685

3685

3685

the oxidation reactions (reactions 1 to 4) as compared to other reactions in the reaction mechanism, especially the Ceria reactions (reactions 10 to 15). However, this decrease in oxidation reaction rates resulted in increased sensitivity of model performance to Ceria reactions. It was also observed that a parametric study on the GMT900-D1 data set using A2 kinetics showed low sensitivity for CO and HC oxidation reactions. Additionally, there were unresolved issues with model prediction for GMX001-D2 as explained earlier. These strongly indicated the need for further improvement in reaction kinetics. The four pre-exponential factors in inhibition term (k1,G, k2,G, k3,G, and k4,G) along with the deactivation length for the front brick (which is usually assumed to be first 0.25 in. of the front

brick and no deactivation in the rear brick) were included in the list of optimization parameters. In addition, the number of data sets for parameter estimation was increased to 5 by including the GMT800 data set during the optimization process. The GMT800 data set was included because the catalyst does not have any Ceria in front brick, and hence, including this data set in the calculation of the objective function could help avoid the low model sensitivity toward oxidation reactions. All these modifications to the optimization schedule increased the optimization parameters to 18 (11 pre-exponential factors in reaction rates, 4 pre-exponential factors in inhibition term, and 3 deactivation lengths, one each for GMX001, GMT900, and GMT800). However, no significant improvements in model 9971



ARTICLE

Figure 4. Performance of various sets of reaction kinetics for GMT900-D1 and GMT900-D2 data sets.

performance over kinetics A2 were observed by using the modified optimization schedule, and hence the results of this run are not reported. It was found that within limits (0.2 to 0.3 in.), varying the deactivation length did not have significant effects on model predictions. Hence, the deactivation lengths were removed (in view of the above observation and more over, it is not feasible to have deactivation length as a tuning parameter as it becomes specific to the monolith used, and hence validation of the kinetics with different monolith reactors will become inconclusive) from the list of parameters to be optimized, and the data sets used in optimization were restored to the original set of 4 data sets i.e., GMT900-D1, GMT900-D2, GMX001-D1, and GMX001-D2. The constraints on k4 and k14 were also removed and were made independent optimization parameters. These

modifications reduced the number of optimization parameters to 17. Additionally, the individual objective functions for CO, NO, and HC (given by eqs 1517) were used instead of the total objective function. This meant that the PAO (Pointer Automatic Optimizer) would have three objective functions to be minimized simultaneously. This optimization schedule resulted in a new set of reaction kinetics and is listed in Table 7 as A3. The model performance using kinetics A3 improved predictions for HC and NO, but improvements in predictions for CO were not significant. However, kinetics set A3 was an improvement over A2 because model performance with kinetics A3 showed higher sensitivity toward oxidation reactions for CO and reduction reactions of NO. But the reaction kinetics required further improvement because model predictions for HC were still more sensitive toward Ceria 9972



ARTICLE

Figure 5. Performance of various sets of reaction kinetics for GMX001-D1 and GMX001-D2 data sets.

reactions. Manual tuning was performed using the GMT900-D1 data set to further tune the kinetics (CO and HC oxidation reactions, NO reduction reactions, and HC and NO reactions with Ceria) and to establish the correct sensitivities to obtain kinetics A4. It should be mentioned that the evolution of the kinetics set from Base to A3 provided sufficient understanding of the system to perform this final manual tuning.

’ VALIDATION OF ESTIMATED PARAMETERS The optimized set of reaction kinetic parameters is validated against two data sets, GMT800 and LY7, which were not used for calibration purposes. It should be pointed out that the GMT800

data set and the LY7 data set were generated for FTP drive cycle and European Drive Cycle, respectively. As is known, the European drive cycle and the FTP driving cycles are quite different in terms of the engine-out transients. GMT800 catalyst has no Ceria in the front brick, while the LY7 catalyst has Ceria in both the bricks. Essentially, the model was being tested on two vehicles with different catalyst compositions, converter properties, and also different engine power (i.e., different engine-out exhaust compositions). The estimated kinetic parameters (A3 and A4) were used to predict the experimental data without any further tuning or estimation of the kinetic parameters. Figure 8 presents model performance using kinetics A3 and the optimized set A4 for these two data sets. Model performance using kinetics 9973



Figure 6. Performance of various sets of reaction kinetics, shown in terms of percentage error between simulated and measured final cumulative emissions at the end of first driving cycle, for GMT900-D1 and GMT900-D2 data sets.

ARTICLE

Figure 7. Performance of various sets of reaction kinetics, shown in terms of percentage error between simulated and measured final cumulative emissions at the end of first driving cycle, for GMX001-D1 and GMX001-D2 data sets.

oxidation and with some importance to the NOx reduction reactions. The kinetics of the reactions involving Ceria are still not resolved, and this along with the NOx reduction reactions determine the conversion efficiencies after lightoff (during the warmed up phase). This could be a potential study for the future. b) Second, is to facilitate the use of the global kinetic model for control purposes and in 3D models. Global or microkinetic models that represent the reactions in a three-way catalyst (or exhaust aftertreatment devices in general) are quite complicated, and the number of reactions and species are significant. Though these kinetic models are very useful for evaluating various designs (using 1D models like the one shown in this paper) and to understand the system better, these can be computationally very demanding when used in 3D models. Additionally, models used for control purposes need to be very robust and fast and look-up tables or rate maps are usually preferred. Hence, there is a need to reduce these kinetic mechanisms to a smaller set and also create rate maps that can quantitatively capture the reaction behavior reasonably well. Systematic approaches for reaction mechanism reduction techniques like principle component analysis,67 sensitivity analysis,68,69 computational singular perturbation techniques,70 simulation error minimization methods,71 or repro-modeling72,73 can be used for mechanism reduction. One of the ways of repro-modeling is to create rate maps based on detailed surface reaction mechanisms, and it has been found to be effective for simpler systems like CH4 oxidation on Platinum catalysts74 and for NH3 oxidation on Platinum catalysts.75 These mechanism reduction and rate mapping approaches need to be extended to the full reaction set for TWC and other exhaust aftertreatment devices so that the global and microkinetic mechanisms developed using experimental data can be used for control purposes and in 3D models. This could be an important extension of this work.

A4 shows improvement in predictions as compared to A3 for HC and NO. Although kinetics A4 does not show any improvement over A3 in terms of prediction for CO, kinetics A4 is an improvement over A3 in terms of sensitivities of model output to the oxidation reaction.

’ SUMMARY The main contribution of this work has been in obtaining a global kinetic model (that can predict the light-off) for TWC using vehicle data. In this work, six different vehicle data sets are used to calibrate and validate the TWC global kinetic model. The emphasis in this work is restricted to predicting the first drive cycle (125 s), particularly light-off, as this is the crucial period of cold-start emissions. The model was calibrated using 4 vehicle data sets (which use the FTP drive cycle) using iSIGHT software package. The kinetic parameters of the various reactions occurring in the TWC were estimated to match the model predictions with experimental data through a combination of exploratory and local optimization methods. A systematic approach (with increasing complexity) is being used to estimate the kinetic parameters. The estimated parameters are then used to validate the model on two different vehicle data sets generated on different drive cycles (and hence different engine-out transients) with different catalyst compositions and engine power (and hence different engine out exhaust compositions). The model was able to predict the light-off reasonably well. The parameter estimation approach in this work is kept as generic as possible to exhaust aftertreatment devices, and a set of guidelines for parameter estimation (specifically for use in exhaust aftertreatment devices) is presented (in the Appendix). The work presented here can be extended in two major directions: a) First, is to focus on improving the global kinetic model. The kinetics development was focused on predicting the light-off only, and hence the focus was mainly on the 9974



ARTICLE

Figure 8. Performance of kinetics sets A3 and A4 for GMT800 and LY7 data sets.

’ APPENDIX GUIDELINES FOR PARAMETER ESTIMATION

0 Known parameters: These are parameters that have very little uncertainty and include length, area, flow rates, geometric dimensions, etc. 0 Unknown parameters: These are usually the uncertain parameters such as reaction rate parameters expressed in turnover numbers (i.e., reaction rate per mole site of the catalyst). Most often we also have to assume the number of active sites.Usually, all (or a part) of the unknown parameters are the ones that need to be estimated. Listing the parameters into these categories helps in deciding what parameters need to be estimated and how much certain (or uncertain) a particular parameter is. If any of the uncertain parameters are not included in the parameter estimation

In this section, we provide guidelines for parameter estimation (specifically for use in exhaust aftertreatment devices). The information provided below is by no means exhaustive or complete but attempts to provide insights to parameter estimation using experimental data more efficiently. While most of them are generic, some of them may specifically pertain to catalytic reactors. Parameters, Experiments, and Mathematical Model . Choosing the parameters to be estimated: The parameters of the model can be classified in two broad categories 9975


Industrial & Engineering Chemistry Research process, it is important to understand that the estimated parameters are only as good as the uncertain parameter(s). . Choose the right experiments and measurements to be done to get the parameters. Once the parameters to be estimated are decided, the right set of experiments and measurements must be chosen for parameter estimation purposes. The experiments could be one of the following: transient or steady-state data from lab/reactor experiments or vehicle data. One of the indicators which could help in choosing the correct experiments/measurements is the sensitivity of the parameter(s) to be estimated. The parameter(s) to be estimated should have a good sensitivity on the measured outputs (and this depends heavily on the operating conditions). . The model used to estimate the parameters should describe the experimental conditions as much as possible. It is important that the mathematical model used to describe the experimental conditions takes into account all the dominant physical and chemical processes. Some of the effects that are usually ignored are heat loss in the reactor and flow effects. Any additional effects that a default model would not take into account should be accounted for. While it may be (mathematically) possible to estimate the parameters of any model to any data, the test of the model lies in its prediction capabilities. If any of the dominant phenomenon(s) is (are) not accounted for, then the predictive capabilities of the model decrease drastically. Objective Function . The form of the objective function is critical in estimating the parameters. Some of the factors that need to be considered in framing the objective function are listed below: 0 Difference between experimental and simulated data can be expressed in various forms. Some of the conventional formulation of the objective functions are: 0 Sum of squares of the difference (simulatedmeasured)2 0 A logarithmic based objective function such as Log(simulated/measured) as mentioned here.76 0 Any other form that goes to “zero” when simulated and measured values are the same. 0 Objective function can be one number or vector of small size. Separate objective functions for each measured output (e.g., species concentrations or temperatures) can be defined. The database of iteration results can be used for sensitivity analysis of the parameters to the specific objective functions. It is crucial to pick the best objective function that can show sufficient sensitivity to the parameters that are being estimated. . Weights to various data points in the objective function 0 Weights should be given to data points based on the time instants (for transient experiments), run number/data set, and various observed/measured quantities. 0 The objective function can also be defined in such a way so as to provide automatic weights (to the data values) depending on the form chosen (i.e., each and every data point will be automatically weighted). 0 Sometimes, there is more confidence in the measurements of a particular experiment or a particular time period in the experiment or in the measurement of a specific quantity compared to the rest. Giving higher

ARTICLE

weight (while calculating the objective function) to the more confident measurement/data points is very helpful. . Scaling 0 It is very important to get all the Y’s (outputs) scaled to the same order or magnitude (e.g., concentration and temperatures are magnitudes apart) else the contribution of one of the outputs to the objective function will over shadow the rest of the outputs. 0 Most often, scaling does not really matter when the quantity “simulated/measured” (ratio of simulated and measured values) is used in the objective function calculation. 0 The objective function should also be scaled with the number of data points used in calculating the objective function. (For example, while calculating the objective function, if 10 data points are used for NO measurements but 20 data points are used for CO measurement, this difference may tend to overshadow the importance of NO in the objective function calculation.) . Exclusion of inappropriate data or regions which are insensitive to the parameters and exclusion of inaccurate experimental data. 0 This is crucial because trying to fit parameters in a region where the parameters are completely insensitive to the measured data could be a wasteful effort. For example, when the measured values for CO, NO, and HC at converter outlet exceeds the corresponding values at converter inlet in a vehicle test, no amount of tuning of model parameters can match the measured data as there is no physical process in the model that explains this behavior. (Refer to the section: ‘Objective Functions’ for more explanation.) Various Optimization Methods. This has been explained in detail in one of the sections in this paper. The important points are summarized here. The optimization methods can be classified into two broad categories. . Exploratory or evolutionary methods 0 Suited for searching the global space (within bounds) 0 Examples: Multi-Island Genetic Algorithm, Simulated Annealing, Downhill Simplex Optimizer . Local Methods 0 Suited for narrowing down to the optimum in a local space and are mostly gradient based methods 0 Examples: Methods of Feasible Directions, Conjugate Gradient, Sequential Quadratic Programming . Pointer Automatic Optimizer (in iSIGHT) is a combination of Local and Exploratory methods 0 Includes genetic algorithm, Nelder and Mead downhill simplex, sequential quadratic programming (NLPQL), and a linear solver Parameter Estimation Method/Process . Initial guess is crucial for local optimization methods. A good guess can usually be obtained from older attempts of optimization or from literature or from fitting a smaller set of experimental conditions (i.e., if 4 vehicle data sets need to be optimized simultaneously to get a set of parameters that will fit all, then the process could be started with optimizing one, which will be easier, and extend it to the rest). 9976


Industrial & Engineering Chemistry Research . Bounds/Constraints for the parameters: The estimated parameters should not be very far away from the physics (closer to literature numbers). Hence, bounds/constraints would be a way of incorporating the physical bounds into the mathematical optimization method. It is easily seen that smaller the space of parameters (or narrow bounds), lesser the time for optimization. . Scaling of parameters is important to avoid singularity or illconditioned matrices during optimization. All parameters that are to be estimated should be of the same order (O(1) values are usually preferred). . Estimating the right parameters from the right data: It is important to choose only those parameters that have some sensitivity to the objective function. Though it is usually seen that some specific data points would reveal information on specific kinetic parameters, sometimes there is no direct dependence of some of the parameters on the measured data. For example, the H2 reaction rate cannot be easily obtained unless H2 measurements are made. In such cases, the H2 reaction rate is obtained indirectly by comparing other measured outputs. . Increasing the complexity systematically: It is best to start with fewer parameters and slowly increase the number of parameters. The physics and chemistry of the system can help in choosing the most dominant parameters. Starting with the most significant (or sensitive) parameters and slowly adding the lesser significant parameters is found to be an efficient way to estimate the parameters. Local sensitivity analysis can be used to determine significant parameters (for local optimization methods). Runtime, Analysis, and Other Important Factors . It is suggested that all unknown (or even slightly uncertain) parameters are chosen as possible optimization parameters and the parameters to be estimated are chosen during optimization. This helps in maintaining a consolidated database of parameters versus outputs (particularly when using iSIGHT, as iSIGHT automatically builds an Excel file with the various parameters and the corresponding objective functions). . The computational time for each simulation (or each objective function evaluation) and for the entire optimization process helps in deciding the best optimization methods (Exploratory, Local, or a mix of the two). . It is important to account for any failures in the code that is being used to calculate the objective function. For some parameter values (during optimization) the code or the model could fail. One of the ways to overcome this is to assign a penalty (or higher objective function value) for the failure runs. . Among the various data points available, some of them should be used for calibration (parameter estimation) and some for validation. Trying different combinations can help in ensuring that the optimization results are not vastly different. . Analysis of optimization results: The huge database of runs generated (for various set of parameters in the parameter space) can be used to get some sensitivity information. Some parameters may hit the bounds (during local optimization methods), and this will give a hint as to whether that parameter bounds need to be relaxed. . For a new set of optimization runs, the previous optimization results can be used as an initial guess (particularly useful for local methods).

ARTICLE

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected]. Present Addresses §

Laboratory of Thermodynamics in Emerging Technologies (LTNT), Department of Mechanical Engineering (DMAVT), Swiss Federal Institute of Technology (ETH), Zurich.

’ ACKNOWLEDGMENT The authors would like to acknowledge K. Raghunathan (GM Global R&D center) for providing valuable comments on this report. The authors acknowledge Dave Belton, Dave Radke (both from GM Powertrain), and Wei Li (GM Global R&D center) for providing the vehicle data sets for the NA vehicles (GMT900, GMX001, and GMT800). The authors acknowledge Michael Keat and Jeremy Tassone (both from GM Powertrain) for sharing the vehicle data on the LY7 vehicle. ’ NOMENCLATURE aj catalyst surface area per unit volume of the monolith molsite/m3 A monolith frontal area m2 c total molar concentration mol/m3 cp,g specific heat capacity of gas J/mol/K specific heat capacity of the substrate J/kg/K cp,sb specific heat capacity of the washcoat J/kg/K cp,wc hydraulic diameter of the channel m Dh mass diffusivity of the ithth species m2/s Di,m objective function for i data set fi f combined objective function for all data scts fCO,fNO,fHC separate objective functions for CO, NO, HC, respectively volume fraction of the substrate in the monolith fsb volume fraction of the washcoat in the monolith fwc function to specify coverages s1 Fk ΔG Gibbs free energy of formation J/mol (ΔH)j heat of reaction for reaction j(

Kinetic Parameters Estimation for Three Way ... - ACS Publications

Recommend Documents