Sensitivity Analysis and Kinetic Parameter Estimation in a Three Way

To whom correspondence should be addressed. Email: [email protected]. Tel.: +91-44-22574161. Fax: +91-44-22570509., † ... In this work, the kinetic p...
2 downloads 4 Views 307KB Size
Download PDF
Ind. Eng. Chem. Res. 2009, 48, 3779–3790

3779

Sensitivity Analysis and Kinetic Parameter Estimation in a Three Way Catalytic Converter S. Koteswara Rao,† Rayees Imam,† Karthik Ramanathan,‡ and S. Pushpavanam*,† Department of Chemical Engineering, Indian Institute of Technology, Madras, 600036, India, and India Science Laboratory, General Motors Research and DeVelopment DiVision, Bangalore 560066, India

A heterogeneous one-dimensional model is used to simulate a three way catalytic converter (TWC). The mathematical model consists of mass and energy balance equations in the gas and solid phases and accounts for 15 heterogeneous chemical reactions. Langmuir-Hinshelwood type kinetics is used to represent most of the reaction rates. In this work, the kinetic parameters of the various reactions occurring in a TWC are estimated using vehicle data based on a genetic-algorithm-based optimization. The method uses the test data to estimate the reaction kinetic parameters by minimizing the deviation between the model predictions and experimental observations. It is found that the objective function based on the cumulative mass of a species that leaves the reactor until a certain time instant is a preferred performance measure for the optimization. The sensitivity of the exit concentration of various species to various kinetic parameters was found. This is used to identify the subset of reactions which have a significant effect on the various species concentrations in the effluent stream and is used to determine the set of kinetic parameters for optimizing the reactor performance. Single parameter and multiple parameters tuning using genetic algorithm have been demonstrated successfully for a TWC application. 1. Introduction Three way catalytic converters are used to reduce the emission levels of pollutants in the exhaust gases coming out from internal combustion engines. As the emission standards of CO, hydrocarbons (HC), and NOx become more stringent, the design of TWCs is being increasingly supported by mathematical modeling. This requires good understanding of the physical (transport processes including mass and heat transfer) and chemical phenomenon (heterogeneous catalytic reactions, adsorption, desorption, and storage of chemical species). In the last 30 years, many models have appeared in the literature focusing on these physical and chemical phenomena with varying detail. Mathematical models ranging from simple one-dimensional models to more complex multidimensional models have been proposed to explain the behavior of these reactors. To understand and improve the performance of catalytic converters, it is important to develop a model incorporating the transport processes involving heat and mass transfer from the gas phase to the solid phase along with the heterogeneous reaction processes. There have been different models in the literature with varying degrees of complexity to describe the TWC (Young and Finlayson,1 Oh and Cavendish,2 Chen et al.,3 Oh. et al.,4 Dubien and Schweich5). Jirat et al.6 considered the storage of NOx and O2 on the washcoat in their model. They analyzed the effect of periodic switching between lean and rich combustion conditions. Shamim et al.7 have investigated the performance of the TWC in a monolith for different sets of chemical reactions. The model is similar to that proposed in Oh and Cavendish2 but includes the heat loss to surroundings by convection and radiation. Koci et al.8 considered diffusion in the washcoat and employed micro kinetics for the oxidation of CO, C2H2 and reduction of NOx with a mechanism consisting of 23 reaction steps. For practical applications, a one* To whom correspondence should be addressed. Email: spush@ iitm.ac.in. Tel.: +91-44-22574161. Fax: +91-44-22570509. † Indian Institute of Technology. ‡ General Motors Research and Development Division.

dimensional model capturing all the physics has been found sufficient to explain the important aspects of these reactors. These are optimal as they are computationally inexpensive and are mathematically robust. The chemical kinetics for the various heterogeneous reactions in a catalytic converter depend on several factors and have to be determined accurately for good qualitative and quantitative predictions by the model. The reaction rates are different for each washcoat formulation and cannot be determined a priori because (1) the kinetics are very sensitive to the washcoat formulation, (2) catalytic processes suffer from aging and deactivation, (3) the reaction path depends on the possible presence of species in real exhaust gas that are not present in controlled experimental studies, and (4) the kinetics is sensitive to feed composition and temperature. Most of the kinetic models employ empirical Langmuir-Hinshelwood type kinetic expressions for the chemical kinetics which involve tunable parameters. All phenomena that are not explicitly accounted for by the model are lumped into the values of these tunable parameters. The Langmuir-Hinshelwood type of kinetics used in automotive catalytic converters has its basis in the pioneering work of Voltz et al.,9 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. The first reported computer aided tuning of kinetic parameters in the literature in the context of the three way catalytic converter based on the method of conjugate gradients was by Montreuil et al.10 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/or Pd/Rh catalyst to justify the improved performance of different catalyst formulations. They used a model similar to that of Oh and Cavendish2 and also considered heat loss to surroundings. They considered 13 different reactions with kinetic expressions of the Langmuir-Hinshelwood form taking into account inhibition effects. There were approximately 95 pa-

10.1021/ie801244w CCC: $40.75  2009 American Chemical Society Published on Web 03/20/2009

3780 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009

rameters in the model. They used a multidimensional conjugate gradient method for tuning all the parameters including those occurring in the inhibition terms. This work was based on data obtained under controlled steady-state conditions in an experimental laboratory scale setup. The transient behavior for a real automobile exhaust system has also been analyzed by Pontikakis and Stamatelos.11 They used a two-dimensional model for the solid phase energy balance in the channel. They tuned about 12 kinetic parameters using the conjugate gradient method. They considered the oxidation of CO, H2, HC, oxygen storage, reduction of NOx, and steam reforming reactions (a total of about 10 reactions). The model is similar to that used by Oh and Cavendish2 but includes the effect of heat loss to ambient by (natural or forced) convection and radiation. 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 minima to get the lowest objective function (for physically realistic kinetic parameters). This can be achieved using an evolutionary algorithm like the genetic algorithm (GA). One of the early works in the literature to use this approach in the context (TWC) of this problem was by Glielmo and Santini.12 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 seven parameters simultaneously. Recently, Pontikakis et al.13 have given a good review of this subject. They considered the oxidation of CO and HC and reduction of NOx as well as oxygen storage reactions with Ceria. The model considered is similar to that of Pontikakis and Stamatelos.11 They used a genetic-algorithm-based method to tune around 10 parameters. The objective of this work is to develop and implement an optimization algorithm which would use transient vehicle test data from experiments to tune the reaction kinetic parameters. This is achieved by minimizing the deviation of the model prediction from experimental observation. Our emphasis is on predicting the light-off (50% CO or hydrocarbon conversion) in TWC. Most often, the data from the vehicle tests also have uncertainty, and hence we also discuss the sensitivity of the various estimated kinetic parameters on the performance of the converter. This sensitivity information gives us insight into which reactions (and kinetic parameters) have a dominant effect on the various species concentrations in the exit stream of the reactor. This also helps in choosing a minimal set of kinetic parameters for estimation purposes. In this work, genetic algorithm is used for estimating the kinetic parameters. This choice is motivated since GA uses a random set of parameters as an initial guess and does not involve computing the gradients of the objective function. This also allows us to look for physically realistic optimum solutions in a wide range of parameter space. 2. Mathematical Modeling As mentioned earlier, the behavior of a TWC has been described by a variety of models in the literature. They span the entire width of the spectrum from approximate models which are computationally fast to very detailed computationally intensive models. The TWC model can be viewed as comprising two major parts, one that focuses on the kinetics of the reaction

and the other that focuses on the momentum, heat, and mass transfer effects. 2.1. Kinetics Modeling. One of the most widely used kinetic models is of the form proposed by Voltz et al.,9 which is based on Langmuir-Hinshelwood-type kinetics. This approach can explain the experimental behavior of oxidation of CO and propylene over platinum pellets. This model includes the inhibition effect of CO, propylene, and NOx. This form of global kinetic approach then became conventional and was (and is) used for even Pd and Rh catalysts (Rh is used to provide the active sites for NOx reduction reactions in addition to serving the oxidation reactions). Some authors have modified this kinetic model by including additional production and inhibition terms. Recently researchers (Chatterjee et al.,14 Koci et al.8) have started using microkinetics in which the reaction is assumed to proceed with multiple intermediate (elementary) steps each of which follows a simple Arrhenius kind of rate expressions. But simulations using this approach are computationally intensive. In this work, we use the Langmuir-Hinshelwood form of kinetics for most of the heterogeneous reactions and a simple Arrhenius-type kinetics for the oxygen storage reactions. This choice of kinetic expressions captures important details (as mentioned in the literature) to explain and understand the real time behavior of the system. 2.2. Heat and Mass Transport Modeling. TWC models can have different levels of complexity depending on how the heat and mass transfer effects are incorporated. To determine the heat and mass transfer effects between the exhaust gas and the solid phase of the converter and to determine the exhaust gas characteristics (temperature and species concentration), the various chemical and physical phenomena have to be modeled. Following Oh and Cavendish,2 a transient one-dimensional two phase model with the corresponding heat and mass balance equations is considered to study the dynamic response of the monolith. The model neglects the radial variations of the gas phase temperature, concentration, and velocity within the individual channels. Axial diffusion of mass and heat in the gas phase is neglected. The gas phase reactions are neglected and it is assumed that the heterogeneous reactions occur only on the external surface of the catalytic wall. 2.3. Fundamental Conservation Equations. The mathematical model is very similar to that used by Oh and Cavendish2 but with a more detailed set of reactions and kinetics. The corresponding energy and mass balances are given below. Solid phase energy balance: (fsbFsbcp,sb + fwcFwccp,wc)

∂Ts ∂2Ts ) fsbλsb 2 + ∂t ∂z nrct

hS(Tg - Ts) -

∑ a (z)(∆H) R j

j

j

(1)

j)1

Fluid phase energy balance: w ∂Tg c ) hS(Ts - Tg) A p,g ∂z

(2)

Species mass balance in fluid phases and catalyst surface (for i ) 1,2, ..., nsp): w ∂xg,i ) -km,iS(xg,i - xs,i) ) A ∂z

nrct

∑ a (z)s j

ij

Rj

(3)

j)1

Here, Ts is the solid phase temperature, Tg is the fluid temperature, xg and xs are the gas and solid phase mole fractions

Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3781 Table 1. Design Parameter Values Used for Simulations parameter

front brick

Table 2. List of All Reactions Considered in TWC rear brick

length (cm) 5.08 7.62 hydraulic diameter of channel (m) 8.065 × 10-4 8.065 × 10-4 solid fraction of substrate (fsb) 0.2 0.2 solid fraction of washcoat (fwc) 0.2 0.2 catalyst loading ratio (palladium to 77.5:2.5 12.5:2.5 rhodium) 0.266 0.2503 active site density (mole-site/m3) specific heat of substrate (J/(kg K)) 1071 + 0.156TS - 3.435 × 107Ts-2 specific heat of washcoat (J/(kg K)) 1000 substrate thermal conductivity (W/(m K)) 1.5 Nusselt number/Sherwood number 3 3 density of washcoat (Fwc) kg/m 1300 density of substrate (Fsb) kg/m3 1720

of species “i”. Here, nrct 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 sij. These, as well as all other variables, are defined in the nomenclature section. The different parameter values describing the geometry and other characteristics of the monolith which are used for the simulations are given in Table 1. The boundary and initial conditions used to solve eq (1-3) are given by Tg ) Tin(t) xgi ) xin,i(t) ∂Ts )0 ∂z Ts(Z) ) Tinit(z)

@ z)0 @ z)0 @ z ) 0, L @ t)0

(4)

The equation for surface coVerage is given by dψ ) F(ψ, Ts, xs) dt

(5)

Here, ψ is the coverage parameter, and F is a function of the solid temperature, solid phase mole fractions, and the coverage parameter. The coverage (or the oxidation parameter) is the ratio of oxygen storage in the higher oxidation state to the total available oxygen storage capacity and is effectively the fraction of cerium fully oxidized. ψ)

moles of CeO2 moles of CeO2 + 2 moles of Ce2O3

It is assumed that all of the CeO2 is initially in the oxidized state (as the higher valence state is more stable) so that ψ ) 1 at t ) 0. The various reactions (15 of them) considered along with their heat of reactions are given in Table 2. Langmuir-Hinshelwood form of kinetics is used for most of the heterogeneous reactions while Arrhenius-type kinetics (with the assumption of reactions being elementary) is used for the oxygen storage reactions. The kinetic expressions/constants for the reactions considered in this work are obtained from the open literature and are listed in Table 3. The default values of the parameters used for the simulation are shown in Table 4. These values (obtained from the open literature) are indicative of the order of magnitude estimates of those parameters. These are sensitive to the catalyst and the methods used for catalyst preparation and hence have to be estimated for a particular catalyst on the basis of experimental data. 3. Genetic Algorithm The optimization algorithm used in this work is based on genetic algorithm which mimics the process of natural biological

S. no.

reactions

heat of reaction in J/mol

Oxidation Reactions 1 2 3 4

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

5 6 7

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

8 9

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

10 11 12 13 14 15

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

-2.83 × 105 -1.93 × 106 -2.04 × 106 -2.42 × 105

NO Reduction Reactions -3.73 × 105 -2.74 × 106 -3.32 × 105

Water-Gas and Steam Reforming Reactions -4.12 × 104 +3.74 × 105

Ceria Reactions (Oxygen Storage) -2 × 105 -1.90 × 105 -1.83 × 105 -4.77 × 105 -4.95 × 105 -1.42 × 105

evolution. Here a population of solutions evolves with time such that each member of the population is a solution candidate. Individuals making the population go through the processes of birth, mating, reproduction, mutation, and death. The fitness of each individual is determined and more fit individuals are chosen to get off-springs which are capable of surviving a screening process. This process is repeated over a period of generations and the surviving individuals are candidates for a solution. A good discussion of features of genetic algorithms and their applications can be found in Goldberg.15 The genetic algorithm is adopted and is applied to our system for kinetic parameter estimation. The main features/steps involved in the algorithm are as follows: 1. Lower and upper bounds on the kinetic parameters are set and this restricts our region of search for the optimal value to a restricted space. Note: Lower and upper bounds are chosen to be not more than ( 20% of the values reported in the literature. 2. The number of individuals is fixed. This number is typically around 30 to 60 and depends on the number of parameters to be tuned. The maximum number of generations for which the algorithm is processed is fixed at 200. It is found that for single parameter estimation convergence occurs in less than 15 generations while estimation of seven parameters requires about 150 generations. 3. Based on the number of individuals specified an initial population is created. This population is encoded in binary format to enable the operations of crossover and mutation. This is done using the GA toolbox of MATLAB. 4. The objective function is evaluated for the various individuals created. Fitness values are assigned from the objective function. For every individual, the front brick of the monolith is simulated, followed by simulation of the rear brick to get the outlet gas concentrations from the TWC (see Figure 1). The objective function calculated from this predicted (outlet) value and the experimental data is used to estimate fitness value. 5. A subset of these individuals is selected based on the fitness values of this generation. Individuals (also called chromosomes or strings) are subject to the processes of recombination (crossover and mutation) to produce off-springs. Crossover is applied to the selected chromosomes with a probability of 0.7 and mutation is applied randomly with a low probability of 0.03. The population for the next generation consists of the off-springs

3782 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 Table 3. Reaction Kinetics of the Various Reactions Listed in Table 2 rate expressions R1 ){k1cCOcO2 exp(-E1/(RgTs))}/G R3 ) {k3cC3H8cO2 exp(-E3/(RgTs))}/G R5 ) {k5cCOcNO exp(-E5/(RgTs))}/G R7 ) {k7cH2cNO exp(-E7/(RgTs))}/G R8 ) (k8/∆GWGS) exp(-E8/(RgTs)){cCOcH2O - cH2cCO2/KWGS} R9 ) (k9/∆GSR) exp(-E9/RgTs)){cC3H6cH2O - cH26cCO3/(KSRcH2O2cref5)} R10 ) k10 exp(-E10/(RgTs))cO2(1 - ψ) R12 ) k12 exp(-E12/(RgTs))cCOψ R14 ) k14 exp(-E14/(RgTs))cC3H8ψ where G ) (1 + K1cCO + K2cC3H6)2(1 + K3cCO2cC3H62)(1 + K4cNO) K1 ) k1,G exp(-E1,G/Ts), K2 ) k2,G exp(-E2,G/Ts), K3 ) k3,G exp(-E3,G/Ts) K4 ) k4,G exp(-E4,G/Ts), KWGS ) exp(-∆GWGS/(RgTs)), KSR ) exp(-∆GSR/(RgTs)) ∆GWGS ) -4.1034 × 104 + 44.19Ts - 5.553 × 10- 3Ts2 ∆GSR ) 3.783 × 105 - 546.1Ts - 3.768 × 10- 2Ts2 Table 4. List of Default Kinetic Parameters and Activation Energies for Reactions Listed in Table 2 Prefactor mol - reac / (mol - site)sec ( m3 / mol )n, where n ) 1 for (k1-k9) 0 for (k10-k15)

{

k1 ) 4.21 × 10 k4 ) 4.21 × 1014 k7 ) 6.1 × 109 k10 ) 0.987 k13 ) 2270 13

}

k2 ) 1.02 × 10 k5 ) 1.19 × 109 k8 ) 1.8 × 105 k11 ) 46.1 k14 ) 2250

15

k3 ) 4.68 × 1010 k6 ) 6.2 × 1010 k9 ) 1.23 × 105 k12 ) 5.05 k15 ) 5.05

Activation Energy (J/mol) E1 ) 111450 E4 ) 111450 E7 ) 69237 E10 ) 5296 E13 ) 39070

E2 ) 129530 E5 ) 52374 E8 ) 56720 E11 ) 25100 E14 ) 39680

E3 ) 165160 E6 ) 90063 E9 ) 81920 E12 ) 21768 E15 ) 21768

Inhibition Terms (Mixed Units) k1,G k2,G k3,G k4,G

) ) ) )

7.33 2.57 × 102 1.8 × 10-4 3.14 × 106

E1,G E2,G E3,G E4,G

) ) ) )

-485 166 -10163 3685

that are generated from the parents and the individuals that have high fitness values. As the number of generations increases the algorithm produces better fit individuals. The concentration of these individuals becomes more in a population. The process is continued till the parameter converges and the objective function reaches a minimum value. 6. During this optimization process it is possible a variable may hit the specified lower or upper bounds. This indicates that the optimal value lies outside the range specified. In this case the parameter bound is relaxed to a new value (to an extent that the new bound is not physically unrealistic).

R2 ){k2cC3H6cO2 exp(-E2/(RgTs))}/G R4 ) {k4cH2cO2 exp(-E4/(RgTs))}/G R6 ) {k6cC3H6cNO exp(-E6/(RgTs))}/G

R11 ) k11 exp(-E11/(RgTs))cNO(1 - ψ) R13 ) k13 exp(-E13/(RgTs))cC3H6ψ R15 ) k15 exp(-E15/(RgTs))cH2ψ

Some minor improvements over the standard formulation using GA have been incorporated. Instead of starting the optimization with a single random initial population, two random initial populations are used and the best halves of each population are combined get the new initial population for GA. One main advantage of selecting two random populations is that the individuals in the new population have higher fitness values which enable faster convergence of the GA. If an individual parameter or a set of parameters remains the same across successive generations (i.e., the evolution of the optimization process does not lead to change in parameters), then the optimization is terminated automatically. 4. Objective Function Formulation The experimental data used in this work has been obtained from a vehicle/engine experiment over a typical FTP (federal test procedure) cycle. The TWC had two bricks, a front and a rear brick, and a schematic of the two bricks is shown in Figure 1. The concentration of species and flow rates are measured at the inlet to the front brick and at the exit from the rear brick. The system of partial differential equations (1-5) is solved numerically using the numerical method described in Oh et al.4 Kinetic parameters estimation is done by minimizing the error between the predicted and experimental data of species concentrations. In the next two subsections, the formulation of two possible objective functions (one based on the cumulative data and the other based on instantaneous data) are defined and their performance (for parameter estimation purposes) is compared. 4.1. Instantaneous Objective Function. The time series data contains the concentration of various species at different time instants. The objective function, which has to be minimized, is calculated based on the instantaneous mole fraction of species CO, HC, and NOx measured in the outlet stream. It is calculated as the weighted sum of squares of the deviation of the instantaneous predicted concentrations from experimental values of various components. J(θ) )

Figure 1. Schematic representation of experimental setup showing two TWCs in series.

N

Ncomp

i)1

j)1

∑ ∑w

j

(

pred X meas i, j (t) - X i, j (t)

X in i, j(t)

)

2

(6)

where, Xmeas is the experimentally measured species mole fraction, Xpred is the species mole fraction predicted from the model and Xin is the inlet mole fraction. Here, N represents the number of data points or time instants for which concentration measurements are available, Ncomp represents the number of components/species used in the determination of the objective

Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3783

function and can take a maximum value of 3 when all the species (CO, HC, NOx) are used. θ represents the parameter vector consisting of kinetic parameters to be estimated. Only three species CO, HC, and NOx are considered in this work as in most practical situations these are the only measured species (usually, concentrations of H2 and O2, storage capacity, solid/ gas temperatures are not measured). w1,w2,w3 corresponds to the weighting factor for CO, HC, and NOx, respectively. Depending on the species that needs to be matched with the experimental data, these weights can take a value of 0 or 1. For example, if CO and HC are the species of interest then w1 and w2 will be unity and w3 will be zero. The error is minimized (similar to a least-squares minimization) using the instantaneous objective function over the considered time interval. 4.2. Cumulative Objective Function. Another approach is to use the total (or cumulative) amount (instead of using the instantaneous concentration) of a species (CO, HC, NOx) which leaves the converter till a given time instant. The cumulative amount (at time t) of species j which has left the reactor till a time instant τ is defined as Cut, j )

∫ V(t)C τ

out,jMj

0

dt

(7)

Here V(t) is the volumetric flow rate of the exhaust gas (m3/s), Cout,j (mol/m3) and Mj is the exit concentration and the molecular weight of the jth gas species. At any given time instant, the cumulative values are calculated from model predictions as well as experimental data. This value is a macroscopic indicator of the converter performance. As before, the objective function minimizes the error between the experiments and the predictions in a typical least-squares formulation. The objective function is formulated as f(θ) )

N

Ncomp

i)1

j)1

∑ ∑w

j

(

Cumeas - Cupre Cumeas

)

2

(8) i,j

where, Cumeas i,j and Cupre i,j represent the experimental and predicted cumulative outlet of component j at time i (in kg). The model is investigated in the first drive cycle (125 s) as this is the crucial period for cold-start emissions. After this time period, in most cases, the emissions are not significant as the temperature has increased and this ensures all reactions go to completion. In this work, the time considered to compare model predictions with experimental data was 140 s and this time period captures the important transient behavior (including light off). We shall now discuss some preliminary results on performance prediction using the two objective function formulations. In Figure 2, the performance of the two objective functionss instantaneous and cumulativesis compared when the parameter k11 is being estimated. In Figure 2a the time series data of NOx concentration obtained experimentally is compared with predictions of the model for the default values of the parameters. In Figure 2b, the experimental observations are compared with the model predictions using the optimized parameter k11, obtained by minimizing the instantaneous objective function using CO, HC, and NOx concentrations. Similarly Figure 2c depicts the comparison of the experimental observations with the predictions of the model using the optimized parameter k11, obtained by minimizing the cumulative objective function. The estimated kinetic parameter in Figure 2b is 135.60 and in Figure 2c is 198.1. The parameter k11 was estimated using all the three species for the objective function calculation by considering all the weighting factors (wj, where j ) 1, 2, 3) as 1. The minimum

Figure 2. Variation of NOx concentration with time. Comparison of experiments and simulations: (a) default values of parameters; (b) parameter k11 tuned using instantaneous objective function; (c) parameter k11 tuned using cumulative objective function.

value of the instantaneous objective function obtained after the parameter converged is 7.09. Similarly for the cumulative objective function the minimum value is 7.41. (Note that these objective function values are not directly comparable). To determine which of the objective functions is effective and reliable for parameter estimation, the sensitivity of the two objective functions to changes in the parameter being estimated is determined. The sensitivity is obtained by varying the parameter around the optimal values of the selected parameter and determining the corresponding change in the objective function. To find the sensitivity of the objective functions, the kinetic parameter is given a (20% deviation from the optimal value. The objective functions (instantaneous and cumulative) are calculated and their dependency on the parameters k11 and k13 are shown in Figures 3 and 4, respectively. In these figures the y-axis represents the relative objective function, calculated as the ratio of objective function at the deviated parameter (cumulative or instantaneous) to the objective function at the optimal parameter value. The x-axis represents the ratio of the deviated parameter values to the optimal parameter value. This ratio is called the relative parameter value. The optimal or reference parameter value for each objective function is the value obtained using all three species in the objective function. For each objective function, the corresponding tuned parameters are used as the reference value. Figure 3 shows the dependence of the two relative objective functions on the parameter k11. (Here, only the kinetic parameter k11 is changed and all other parameters are maintained at the default values of Table 4). Figure 3a shows the relative objective function when all the three species are considered in the calculation of the objective function (i.e., wj ) 1 where j ) 1,

3784 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009

Figure 3. Sensitivity of the two objective function formulations near optimal k11.

2, 3), and Figure 3 panels b-d depict the objective function when only CO (w1 ) 1, w2, w3 ) 0), HC (w2 ) 1, w1, w3 ) 0), NOx (w3 ) 1, w1, w2 ) 0) are considered individually in the objective function. Each panel shows the variation in cumulative and instantaneous objective function varies as the kinetic parameter is varied. It can be seen that the cumulative objective function of total (when all species are considered) and for CO and NO is very sensitive to parameter k11, that is, small changes in parameter result in significant changes in objective function. The relative objective function for HC does not change much (as shown in Figure 3c) and hence it can be concluded that varying the parameter k11 does not significantly affect the HC concentrations in the exit stream. Figure 4 shows a similar set of plots as in Figure 3 for the kinetic parameter k13. Here, the total relative objective function is calculated using only CO and HC is given in Figure 4a, that is, w1 and w2 are taken to be 1 and w3 is taken to be 0. Figure 4 panels b and c depict the objective function when only CO (w1 ) 1, w2, w3 ) 0), HC (w2 ) 1, w1, w3 ) 0), are considered individually in the objective function. These results show that the objective function f(θ) is more sensitive than the objective function J(θ) and it can be concluded that the cumulative objective function is more sensitive to the kinetic parameters compared to the instantaneous objective function. Also, the definition of the cumulative objective function intrinsically gives more weight to the initial few seconds (when the reaction starts to occur). This is because the cumulative values at any time t contain the information over all the past times and any data point in the future will have all the information from the past (starting at t ) 0). When this cumulative value is used at various time instants to calculate the objective function, the information from the past is being accounted for indirectly at all the time instants (at which the objective function is calculated). Also, the instantaneous measured values (Xmeas) can be very low (close to and sometimes equal to zero) after the initial light-off (once the converter is in the transport limited regime). This can lead to numerical errors in calculating the objective function whereas the cumulative objective function does not face this problem as the cumulative

Figure 4. Sensitivity of the two objective function formulations near optimal k13.

values do not become close-to-zero. Hence, in this work, unless mentioned otherwise the cumulative objective function is used to calculate parametric sensitivity and optimize kinetic parameters. 5. Sensitivity of Parameters The chosen heterogeneous kinetic model contains many kinetic parameters. However, not all parameters have a significant influence (or similar influence) on the behavior of the system. In this section, the sensitivity of different kinetic parameters on the exit concentrations of various species is determined. This enables us to determine which set of parameters from the entire set has to be tuned/estimated (when the model predictions of various species concentrations do not agree with the experimental results). The sensitivity of the parameters to the system behavior is obtained by giving a small change in the parameter and determining the corresponding change in the cumulative objective function. This is called local sensitivity analysis as the analysis is restricted to the case where the changes/variations in the parameter are small. The information obtained from parametric sensitivity analysis is very useful for optimization and parameter estimation as it helps us determine which parameters to choose in a multiparameter setting. To determine and understand the sensitivity of the parameters, we show the dependence of the relative objective function (defined as the ratio of objective function at a perturbed parameter value to the objective function at a default parameter value) on the relative parameter value (defined as the ratio of

Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3785

Figure 5. Sensitivity of various cumulative objective functions to variations in k11.

Figure 7. Sensitivity of various cumulative objective functions to variations in k13. Table 5. Sensitivity Levels of the Parametersa sensitivity parameter

CO

HC

NOx

total

k1 k2 k3 k4 k5 k6 k7 k8 k9 k10 k11 k12 k13 k14 k15

high poor poor poor medium poor poor poor poor poor high high poor poor poor

poor poor poor poor poor poor poor poor poor medium poor poor high poor poor

poor poor poor poor medium poor poor poor poor medium high poor poor poor poor

medium poor poor poor medium poor poor poor poor medium high medium medium poor poor

a

Figure 6. Sensitivity of various cumulative objective functions to variations in k12.

changed/perturbed parameter value to the default parameter value). In Figure 5, the solid line with circles in the graph shows the relative objective function for CO (i.e., ratio of Fco for a perturbed parameter value to Fco at the default value) as a function of the relative parameter k11 (i.e., the ratio of changed/ perturbed parameter value to the default value). The small dashed line represents the relative objective function for HC, the dashed line with square represents the relative objective function of NOx, and the total relative objective function (including NOx, CO, and HC) variations are represented by dark dashed line. It can be observed that there is nearly 30% change in the relative objective function (for CO) and around 35% change in the objective function of NOx for a 20% change/ deviation in the parameter. There is no significant change in the relative objective function based on HC. From these results we can conclude that the variations in parameter k11 result in significant changes in CO and NOx concentrations, that is, both CO and NOx concentrations are sensitive to k11 and that this parameter does not show any effect on the HCs. Figure 6 shows the variation of contribution of each species to the relative objective function for changes in parameter k12. It is easily observed that the parameter k12 has a significant effect on CO concentration compared to HCs and total objective function (w1, w2 ) 1 and w3 ) 0). There is only a slight variation in the total (CO and HC) relative objective function. This can

Note: Definition of high, medium and poor sensitivity provided in the text.

be explained as being due to the higher numerical contribution of HC in the objective function, that is, HC contributes more to the objective function than CO. Also a change in k12 has only a small effect on the NOx concentration and this is not shown. From Figure 7 it can be observed that the parameter k13 has a significant effect on HC concentration and no effect on CO. The sensitivity of the various species concentrations to other kinetic parameters have been determined and are depicted qualitatively in Table 5. The table illustrates how each parameter affects the concentration of a species in the exit gas stream. In the table, “high”, “medium”, and “poor” represent the degree of sensitivity to the parameter. This is defined on the basis of the percentage variation of the objective function of a particular species observed when a percentage change in a parameter is given around an optimal value. The sensitivity is characterized as high (H) if the percentage variation of the relative objective function is greater than 20% for a 20% change in the parameter, and is characterized as medium (M) if the change in the relative objective function is higher than 10% and it is characterized as poor (P) if the change in the relative objective function is less than 10%. It should be mentioned that the observations made in this section and the earlier section are specific and may be limited to the catalyst formulation used in the work since the kinetics

3786 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009

Figure 8. Flow sheet depicting the genetic algorithm and the link between the simulation and the optimization units of the numerical algorithm.

is sensitive to the washcoat formulations and the gas phase concentrations. These results are also highly dependent on the initial or default set of kinetic parameters that are used. 6. Parameter Tuning From the sensitivity analysis, we determine/recognize (from Table 5) which parameters are more sensitive to the objective function and these parameters are selected for tuning (for optimization purposes). The algorithm for tuning of the parameters is coupled with the simulation of the system. The coupling between the optimization algorithm and the simulation of PDEs is depicted clearly in Figure 8. The flow sheet describes two main routines. One contains the main program for the genetic algorithm and the other subroutine is for the objective function calculation (obtained by solving the nonlinear governing equations of the model). While tuning or estimating kinetic parameters only that subset of parameters (selected for estimation) is allowed to change, and all other parameters are maintained at the default parameter values given in Table 4. We first present results of a single parameter tuning, then perform a three parameter tuning and finally perform a seven parameter tuning and study the effect on performance prediction. Tuning k11. Figure 9 panels a-c show the cumulative amount of CO, HC, and NOx leaving the reactor as a function of time using the default values of kinetic parameters given as a thin dashed line. It is seen that the cumulative amount of a species leaving a reactor reaches a constant value as time increases. This is because after the light off has occurred the concentration of the various species in the exit stream becomes almost zero (i.e., 100% conversion) and no further emission occurs. The values obtained from experiments are shown as a thick dashed line. The default values of the parameters do not predict the reactor performance accurately for any of the species. The solid line in Figure 9a,b,c shows the predicted cumulative output of CO, HC, and NOx, respectively, after tuning the parameter k11 using concentration of all three species in the

objective function. The optimal value of the estimated parameter is 198.08. A significant improvement in the predictions of CO, NOx at this optimal value is observed, but there is no significant change observed in the prediction of HC (see Figure 9b). The initial objective function value taking all three compounds (CO, HC, and NOx) into account is 86.0. The final optimal objective function value after the parameter estimated is 10.58, and this is due to a significant improvement in the predictions of CO and NOx. From the plots shown in Figure 9 it can be concluded that the parameter k11 affects both CO and NO concentrations significantly and does not affect HC concentration. These are also consistent with the sensitivity results presented in Figure 5. Three-Parameter Tuning (k11, k12, k13). The parameters selected for tuning are the pre-exponential factors of the Ceria oxidation and reduction reactions, that is, k11, k12, k13 corresponding to the reaction of Ceria with NO, CO, and propene (HC), respectively. These parameters are known to have a significant effect on the concentration of CO, HCs and NOx when the exhaust gas is rich (i.e., under less oxygen or reducing conditions). The comparison plots for the various compounds are given in Figure 10. The estimated values of the tuned parameters provide a much better match with the experimental data than the default values. There is significant improvement in the predictions of all three species when compared to that obtained from a single parameter tuning. This is to be expected as the sensitivity analysis shows that k11 affects both CO and NOx, k12 affects only CO, and k13 affects HC concentration. In tuning k11, k12, and k13 the optimal value of the objective function obtained is a low value of 2.45. Sensitivity analysis has helped in determining the right choice of a minimum set of parameters enabling us to avoid extensive computation. Seven-Parameter Tuning. As the number of tuning parameters increases it is possible to obtain a better (and move toward

Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3787

Figure 9. Comparison of cumulative mass of various species leaving the converter (parameter tuned: k11).

Figure 10. Comparison of cumulative mass of various species leaving the converter (parameter tuned: k11, k12, k13).

a global) optimal value of the objective function. Seven parameter tuning was carried out to illustrate this and this resulted in a better match of the predictions with the experimental data for the corresponding optimal values of the parameters. The seven parameters are chosen such that the parameters show high sensitivity to the objective function. The parameters that are selected for tuning with the correspond-

ing optimal values are given in Table 6. The plots of predictions of various species concentrations before tuning with those obtained after tuning along with the experimental data are shown in Figure 11. It is seen that we obtain good predictions for the emission of all the compounds CO, HC, and NO. The optimal value of the objective function obtained while tuning these parameters is 0.36. In Table 6, as expected, the optimal objective

3788 Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 Table 6. Optimal Values of the Parameters number of parameters tuned

tuning parameters

single parameter tuning

k11

198.08

10.6

three parameter tuning

k11 k12 k13

250.0a 7.537 0.354 × 103

2.45

seven parameter tuning

k1 k5 k10 k11 k12 k13 k14

3.183 × 1013a 0.421 × 109 3.34 × 10-2 236.21a 1.947 0.856 × 103 0.172 × 103

0.36

a

after tuning

objective function

Values that hit the bounds or are nearer to the bounds.

function value decreases with an increase in the number of parameters being tuned. The predictive capability of the model with the estimated or optimized kinetic parameters is shown in Figure 12. The data used for this purpose was from the same engine/vehicle (and hence same catalysts) but had different inlet/operating conditions. Figure 12 panels a,b,c show the cumulative emissions of CO, HC, and NOx, respectively. Each shows the measured cumulative emissions along with model prediction using the default parameters (listed in Table 4) and model prediction using the estimated/new parameters (listed in Table 6) using seven parameter tuning. This result confirms that the new/estimated set of parameters give a very good prediction of the CO, HC, and NOx concentrations even when the inlet conditions are varied. 7. Discussion and Conclusions A systematic methodology has been proposed and studied to tune the kinetic parameters for a three way catalytic converter, which is used as an after treatment device for automotive exhaust. This methodology has been applied to data obtained from an engine/vehicle application. An optimization technique based on genetic algorithm is used to estimate the kinetic parameters by minimizing the error between experimental and predicted values of concentrations of different species. Two objective functions formulations were studied and it is found that the cumulative objective function formulation is more sensitive compared to the instantaneous objective function formulation to changes in kinetic parameters, and this formulation is used for parameter tuning purposes. Single parameter as well as multiple parameter tuning has been successfully demonstrated for experimental data using genetic algorithm. Tuning seven parameters gave a global optimal value of the objective function. It is also shown that the estimated parameters have good predictive capability. It is observed that as the number of tuning parameters increases, the number of generations required for convergence increases and there is a greater probability of reaching a global optimum. The sensitivity of the exit species concentrations to various tuned parameters has been analyzed. It has been observed that CO and NOx are sensitive to the parameter k11. Similarly k12 influences CO more than HC and NOx, and the parameter k13 affects only HC. Note that these observations are specific and may be limited to this catalyst formulation, initial set of kinetic parameters, and experimental data since the kinetics are sensitive to the washcoat formulations and the gas phase concentrations. But the algorithm and the sensitivity analysis methodologies

Figure 11. Comparison of cumulative mass of various species leaving the converter for seven parameter tuning.

discussed in this work are fairly generic and can be applied to any system. So the quantitative results from this work should be viewed upon as a representative case to illustrate the capabilities of the numerical algorithm and the sensitivity analysis methodologies. One of the main contributions of this work is in identifying which of the parameters significantly affect the exit concentra-

Ind. Eng. Chem. Res., Vol. 48, No. 8, 2009 3789

parameters are included as parameters to be estimated, the optimization problem is ill-conditioned. By avoiding (or not including) the insensitive parameters for parameter estimation purposes, the minimization (or optimization) problem will be well-posed. Acknowledgment The authors acknowledge helpful discussions with Raghunathan of the India Science Laboratory, General Motors Research and Development (GM R&D) Center. We would also like to acknowledge Ed Bissett (of the GM R&D center) for the GM fortran code for emission control simulations, which has been used in this work. We also greatly acknowledge the GM R&D Center for their financial support. The authors would also like to thank the reviewers for their suggestions, which helped clarify certain points and also helped improve the focus and readability of the text. Nomenclature

Figure 12. Predictive capability of the model using the estimated parameters: comparison of cumulative mass of various species leaving the converter.

tion of the various species. This is useful in a multiparameter context. For example, if there is mismatch observed in the model predictions of the concentrations of a species (with the corresponding experimental data) the sensitivity information can be used to determine which of the kinetic parameters needs to be estimated. These would be the parameters which have a high sensitivity to the species concentration. When insensitive

aj(z) ) active site density for jth reaction, mol-site(cat)/m3 A ) monolith frontal area, m2 cp,g ) (molar) specific heat capacity of gas, J/mol/K cp,sb ) specific heat capacity of substrate, J/kg/K cp,wc ) specific heat capacity of washcoat, J/kg/K ci ) concentration of the species i (appears in rate equations), mol/ m3 Cui ) cumulative weight of the species, kg Ei ) activation energy (appears in rate equations), J/mol fsb ) volume fraction of substrate, unitless fwc ) volume fraction of washcoat, unitless f(θ) ) cumulative objective function, unitless G ) inhibition factor, unitless h ) heat transfer coefficient between gas and solid, J/m2/s/K J(θ) ) instantaneous objective function, unitless km,i ) mass transfer coefficient for ith species, mol/m2/s Kw ) equilibrium constant for water-gas shift reaction, unitless L ) length of the converter, m Mj ) molecular weight of jth gas species, kg/mol N ) number of data points (depends on length of time data), unitless Ncomp ) number of component species (CO, HC, NOx), unitless Rj ) specific reaction rate for jth reaction or turn-over rate, mol/ mol-site(cat)/s Rg ) universal gas constant, J/mol/K sij ) stoichiometric coefficient of ith species in jth reaction, unitless S ) surface area of solid substrate per converter volume, m2/m3 t ) time, s Tg ) gas temperature, K Tin ) inlet gas temperature, K Tinit ) initial solid temperature, K Ts ) solid temperature, K V(t) ) volume flow rate (at inlet), m3/s w ) molar flow rate of inlet gas, mol/s wi,j ) weight factor, unitless xg,i ) mole fractions of ith species in gas phase, unitless xin,i ) inlet mole fractions of ith species in gas phase, unitless xs,i ) mole fractions of ith species in solid phase, unitless Z ) axial coordinate, m Greek Symbols ∆Hj ) heat of reaction for reaction j (