9 1997 Elsevier Science B. V. All rights reserved. Dynamics of Surfaces and Reaction Kinetics in Heterogeneous Catalysis G.F. Froment and K.C. Waugh, editors
571
D e v e l o p m e n t o f a c o m p u t a t i o n a l tool for the t r a n s i e n t kinetics o f c o m p l e x chemical heterogeneous reaction systems G. A. Carrillo Le Roux a, I. Bergault b, H. Delmas b and X. Joulia b aLaborat6rio de Simulag~o e Controle de Processos - DEQ - EPUSP - P.O. Box 61548, 05424970 - $5.o Paulo - Brazil - e-mail:
[email protected] bLGC - INPT / ENSIGCT - 18, Chemin de la Loge - 31078 Toulouse CEDEX 4 - France - email: Xavier.Joulia@ensigct. fr
We present a program for the modelling of the transient of complex heterogeneous chemical reaction systems. The philosophy of the program is based on the formulation of a generic model for heterogeneous reactions derived from a phenomenological description of the system, where all the reaction steps should be detailed. The challenge of dealing with the complexity of the model thus generated is addressed by the robustness of the numerical tools employed for its simulation and parameter estimation. The package is able to obtain a fit even for complex models. This fact justifies a meaningful discussion on the usual philosophy of model analysis which consists of steps of reduction and discrimination. The possible issues are illustrated with a study on the hydrogenation of acetophenone over Rhodium catalyst. 1. I N T R O D U C T I O N The usual approach for the generation of models for heterogeneous reacting systems, derived from Langrnuir-Hinshelwood Hougen-Watson (LHHW) equations, involves stages of model reduction by means of hypothesis of rate-controlling steps, quasi steady-state approximations, fast reversible reactions hypothesis and the implementation of lumping of species (1). The approach is based on model discrimination techniques, which involves fitting and comparison of a high number of candidate models, generated by the combinatorial of all pertinent hypothesis (e.g. the dehydrogenation of 1-butene to butadiene (1) which involves the discrimination of 15 candidate models, all based on elementary rate-determining steps assumptions). For complex systems this task can be very tedious and its results disappointing. However, in most cases, the necessity of this step can only be justified by the pragmatic argument that it is very difficult to fit models too complex to data. In this work, we propose an approach to the resolution of the equations derived from the modelling of heterogeneous systems and their fit to experimental data, that permits to relax the necessity of formulation of some of the restrictive hypothesis cited above. In consequence, the worth of the arguments that justified the approach based on steps of reduction and discrimination is reviewed. We present some ideas that can contribute to the development of alternative paths.
572 2. M O D E L
IMPLEMENTATION
The model is formulated in terms of mixed differential and algebraic equations. It is suitable for liquid-solid and gas-liquid-solid reactors. The main assumptions are that the liquid volume and solid mass inside the reactor are constant and that the effects of mass and heat transfer can be neglected. It was generated from balances on both the liquid and solid systems. 2.1 Material Balance in the Liquid Phase The material balance for the liquid phase is expressed as follows:
dc dt
- FIlr r ( 0 , c , | 1 7 4 2.... ,|
(1)
+ FI,er
where c is a no-vector (where nc is the total number of species) which groups the concentrations (in mole per unit volume of liquid) of all the species, in the liquid phase, 0, is a r~-vector of parameters which assembles all the kinetic and equilibrium parameters and | (with i varying from 1 to ns, the number of different sites types) is a n~+l-vector which groups the surface concentration (in mole of species per total number of active sites of type i) of all the species in the site of type i, including the surface concentration of free sites of this type, i noted | f. r is a real n~-valued function by means of which all the reaction velocities are evaluated, without distinction of whether they take place in liquid or solid phase. However, as the dimension of r is mole per unit volume of liquid phase, per unit of time, if a chemical reaction taking place on the solid phase is considered, the reaction terms are to be multiplied by the catalyst concentration in mass of catalyst per unit volume of liquid phase, such that the expression be consistent with equation (1). FI~ is a nr x n~ matrix, each of its columns corresponds to a reaction and contains the stoichiometric coefficients associated to each specie in the liquid phase. Algebraic constraints, in a number of ne, are used to represent fast reversible reactions and equilibrium terms like those arising from adsorption phenomena. They will be presented further. ~ is a n~-vector of displacement velocities. They are slack variables (2) that have the same unit as r. They can be understood as the velocity by which each specie must vary, while the system is reacting, in order to keep the algebraic equations that expresses equilibrium constraints satisfied. FI~ is a nc x n~ matrix which multiplies ~, in order to keep elemental balances satisfied. 2.2 Material Balance in the Solid Material balances are built for each type of site over the catalyst. For a given site of type i, we write down balances for each of the j = 1..... no, adsorbed species:
w C is d|at - FI~, r(O,c,O~,O 2,...,|
+ ~r
(2)
573 and for the free sites of type i we use a material balance expressed in algebraic form: nr
(3)
O j~+ O f ~ --1 j=l
where | j~ is the surface concentration of each of the adsorbed species (j = 1, ..., no, currently grouped in |
and |
f
that of the free sites, w is the catalyst concentration in mass of catalyst
per unit volume of liquid phase and C~s the concentration of active sites of type i, in mole per unit mass of catalyst. ~ r and ~ e are the stoichiometric coefficients of the species adsorbed on the site of type i corresponding to reactions and displacement velocities and are matrices of size nc x n~ and n c x nr respectively. The algebraic constraints, which can express fast reversible reactions or adsorption equilibria correspond to a set of ne equations which can be expressed by:
0=feq(0,C, Ol,O2,...,O ns)
(4)
where feq is a ne-valued real function. A complex heterogeneous reacting system can be thus described by a set of nc (1 +ns) differential equations (equations (1) and (2)) and ne+n s algebraic equations (equations (3) and (4)). It is not necessary to the user to write down the differential-algebraic equations because the package generates the corresponding system on the basis of a detailed description of the reaction mechanism and the kinetic assumptions. It is not necessary to formulate physical hypothesis in order to obtain a simple expression for the overall reaction, however, some usual physical hypothesis can be easily implemented, whose parametrization is equivalent to that of the reduced one. For instance, if it is supposed that sites can be considered in pseudo-steady state it suffices to set the parameters C~sto zero, so that the differential equations (2) are transformed into algebraic ones. If simple adsorption equilibrium is to be considered for all the species over all the sites, this would be equivalent to write down equations (4) as: 0 = c ~ - K I|
i = 1 .... ,n~,
j = l .... ,n~
(5)
Evidently, each of these approaches generates specific implementation problems in the treatment of the initialization of the EDA system and integration error control, but this was generically previewed in the package (3). 3. N U M E R I C A L
TOOLS:
FACING
COMPLEXITY
The EDA system corresponding to the model is solved by a modified version of the LSODI routine, which is based on Gear's method. The version implemented performs the solution of the EDA system concomitantly to the evaluation of the parameters' sensitivities based on the decoupled direct method (4). As a matter of fact, the simulation of the system is
574 very robust, such that it can be carried out with success for very different parameter values while the sensitivities obtained are very accurate. In the estimation procedure most of the techniques described by Farris & Law (5) (essentially the rotational discrimination) and some other features (3) are implemented. This procedure, in conjunction with the great accuracy of the sensitivities, allows a fit to be obtained for models with complex parametrization. However, in most cases a large amount of computational effort is necessary, but it does not worth to remember that this kind of resource is getting less expensive everyday. The main properties of the solution proposed by the algorithm are: it does not correspond to the set of parameters that performs the best fit of the data. If it was the best one, it would very probably be useless for prediction because of the inflation of the estimate due to a possible over-fit of the experimental data (6). it can be said to be inspired on "biased estimation" techniques (7), but it cannot be classified as an estimator. T h e solution proposed is strongly dependent on the initialisation of the parameters values and on other numerical settings, so that it cannot be considered in oneto-one relation with an experimental data set, as an estimator should be. If the model performs a good fit of the data, with parameter values physically consistent, it should be considered as a serious candidate to represent the system. It cannot be considered inadequate because of the great uncertainty of some parameters, neither because of the fact that the solution proposed is not an estimate. -
-
4. D I S C U S S I O N The objective in which the conventional approach is founded is the elucidation of the intrinsic phenomenology of the system, for which the stages of discrimination-reduction are thought to be adequate. The necessity of reduction because of the incapacity of the numerical tools and hardware and the hideous results obtained with conventional estimation procedures on models with a large number of parameters are some related arguments that help to strengthen the necessity of this approach. However, the only scientific proof proper to demonstrate that a model is the "one" that represents a system resides in showing that all the other candidate models are inadequate to represent it. Any different approach does not correspond to a proof, and should be classified as selection procedures, for which preference criteria are designed. For instance, for nested models, a great deal of criteria have been developed like information criteria (8) and likelihood ratio tests (9). In this case the paradox the selection procedure is thought to solve, consists on the fact that while complexity increases, the ability of a model to represent a system also increases. It is usual that mechanisms be assumed without serious physical considerations, an overall expression being chosen from a LHHW table. These models would have the same worth of black-box ones (10) while in the "quest for the true mechanism", the possible fit obtained can lead to the illusion that the proposed mechanism is validated, while a more serious experimental analysis would have been necessary. Our approach permits that a potential candidate model be fitted to the data, independently of its complexity. This does not means that it can fit any model to any data set, however such a feature allows that more attention could be paid to the phenomenological aspects, rather than to the mathematical manipulations necessary to deduce a single expression
575 for the overall rate. It allows also that discrimination on fundamental aspects can be carried on, as many physical hypothesis can be considered in the model without much effort. Furthermore, after a fit is obtained, the Principal Components Analysis (11) can be applied in order to suggest reduction steps a p o s t e r i o r i , which can allow the model to assume a parametrization in accordance with the usual approach. 5. A P P L I C A T I O N We studied the hydrogenation of acetophenone, where a great deal of side-reactions and production of intermediaries take place in a three-phase reactor, with liquid batch and gas continuously fed, at constant pressure and temperature. The catalyst is Rhodium (3%) over activated carbon. The solvent is Cyclohexane. Samples are taken at different instants and analysed by gas chromatography. The species for which measures are available are: acetophenone, AC, phenyl-ethanol, PE, methyl-cyclohexyl-ketone, ethyl-benzene, EB, ethylcyclohexane, cyclohexenyl-ethanol, CNE, methyl-cyclohexenyl-ketone and cyclo-hexylethanol, CE. The proposed reaction mechanism is presented on figure 1. Seventeen runs were performed at 80~ at a pressure of 25 bar. The initial concentration of commercially available species (AC, PE, MCC and EB) was varied between runs.
AC ~ 2 H 2
Y
.-.,
PE
MCNC
X'8 7 EB
EC
CNE
,v
MCC
CE
Figure 1. Reaction Mechanism for the Hydrogenation of Acetophenone An Eley-Rideal mechanism was assumed for the hydrogen: as the pressure was the same for all the runs, the data was non-informative about its adsorption. Three different approaches were used for the analysis of the data: I. In this approach it was assumed that all the species are in adsorption equilibrium and that the sites can be considered in quasi-steady state. II. Independent kinetics were supposed for the adsorption and desorption of all the species. No hypothesis was made on the concentration of sites. III.PCA was applied to eliminate some parameters from the model fitted through approach II. In consequence, it was assumed that the sites are in quasi-steady state and that AC, EC, MCC and CE are in adsorption equilibrium. It was also assumed that kinetic constants for reactions 4 and 8, and 1 and 2 are not independent and are related through a linear relation.
576 The results obtained, summarized on table 1, illustrate that the PCA is able to decrease the ill-condition of the information matrix of the parameters. The criterion presented on the table was deduced from the hypothesis that the covariance matrix of measurements is diagonal and unknown, so that they can be negative. It can be noticed that posing a model with few restrictive assumptions (II) and performing reduction a posteriori, leads to a better fit. This can be an evidence that some of the assumptions made a priori could not be proved experimentally. Table 1 Summary of the results obtained for the hydrogenation of acetophenone Approach I Approach II Number of Parameters 17 26 Final Criterion -2,129 -2.429 Condition Number of the 2 . 1 0 .5 4 . 1 0 .6 Information Matrix
Approach III 21 -2.383 2. 10 -3
6. CONCLUSIONS The great advantage of the approach here presented is that it flees the analyst from having to make simplifying assumptions since the beginning of the modelling process, which can be a dangerous approach (1). The extension of the approach to other systems, following the same philosophy, is intended.
REFERENCES 1. G.F. Froment, K. B. Bischoff, Chemical Reactor Analysis and Design, John Wiley & Sons (1990) 2. V.G. Dovi, Use of slack variables in the mathematical modelling of reaction equilibria in complex chemical kinetics, Comput. Chem. Eng., 19 (4), 489-491 (1995) 3. G.A. Carrillo Le Roux, Stratrgie d'Identification de Modrles Algrbro-Diffrrentiels. Application aux Systrmes Rractionnels Complrxes, Thrse de Dr. Ingrnieur de l'Institut National Polytechnique de Toulouse, France (1995) 4. A.M. Dunker, The decoupled direct method for calculating sensitivity coefficients in chemical kinetics, J. Chem. Phys., 81 (5), 2385-2393 (1984) 5. R.H. Farris, V.H. Law, An efficient computational technique for generalized application of maximum likelihood to improve correlation of experimental data, Comput. Chem. Engng., 3, 95-104 (1979) 6. S. Vajda, H. Rabitz, E. Walter, Y. Lecourtier, Qualitative and quantitative analysis of nonlinear chemical kinetic models, Chem. Eng. Comm., 83, 191-219 (1989) 7. D.W. Marquardt, Generalized inverses, ridge regression, biased linear estimation and nonlinear estimation, Technometrics, 12 (3), 591-612 (1970) 8. L. Ljung, System Identification, Prentice-Hall Inc., N.J. (1987) 9. D.M. Bates, D.G. Watts, Nonlinear Regression Analysis and its Application, Wiley (1980) 10.J.B. Butt, Reaction Kinetics and Reactor Design, Prentice-Hall, New Jersey (1980) 11.S. Vajda, P. Valko, T. Turanyi, Principal component analysis of kinetic models, Int. 3.. Chem. Kinet, 17, 55-81 (1985)