Ind. Eng. Chem. Res. 2007, 46, 1653-1656
1653
RESEARCH NOTES On the Construction of a Continuous Concentration-Reactivity Function for the Continuum Lumping Approach Jagannathan Govindhakannan† and James B. Riggs* Department of Chemical Engineering, Texas Tech UniVersity, Lubbock, Texas 79409-3121
Kinetic modeling based on continuum lumping is a powerful methodology for working with complex petroleum mixtures. The concentration of the lumped species in the reaction mixture, characterized by reactivity, is expressed as a continuous function. The continuous concentration-reactivity function can be constructed from the experimental data, which involves inverting an integral equation. A simple methodology based on a piecewise linear interpolation formulates the inversion problem as a constrained optimization problem. The proposed method is applied to two different vacuum gas oil feedstocks, and the weight fraction values calculated from the continuous concentration-reactivity function are compared with the actual weight fraction values from the true boiling point (TBP) curve. 1. Introduction The feedstocks processed in refinery conversion processes such as thermal cracking, hydrocracking, and catalytic cracking contain an infinitely large number of individual molecules that undergo complex chemical transformations. Even with modernday sophisticated analytical methods, it is impossible to identify and quantify all of the individual chemical species present in such feed mixtures. In this scenario, lumping methodologies offer significant advantages for developing kinetic models. Although the feed mixture is inherently discrete in nature, a mathematically expedient approach has been used to approximate the discrete system by a continuum one, and kinetic lumping, in this case, is performed via integration, instead of summation.1 An important aspect of continuum lumping is to invert the equation
C(t) )
∫0
∞
of this communication. Continuum lumping concepts relevant to inverting eq 1 are briefly reviewed in section 2. We shall present a simple methodology in section 3 that formulates eq 1 into a constrained optimization problem whose solution gives c(k,0). The proposed method has been applied to two different vacuum gas oil (VGO) fractions. The weight fractions calculated from c(k,0) are compared with the weight fraction values derived from the true boiling point (TBP) curve. Conclusions are presented in section 4. 2. Continuum Lumping: Theoretical Background Chou and Ho1 provided the framework for continuum lumping of nonlinear reactions. They showed the conventional procedure for lumping in a continuous mixture by
C(t) ) c(k,t)D(k) dk
(1)
to obtain a continuous function for concentration c(k,0) as a function of reactivity k for the feed mixture from a measured concentration function C(t). The continuous concentration function c(k,0) is needed in material and energy balance equations. D(k) is a species-type density function, which is independent of time, introduced to address the coordinate transformation, from species index i to reactivity k. Chou and Ho1 suggested that the most practical approach is to choose known functions for D(k) and c(k,0) with adjustable parameters and fit the C(t) data to find the best parameter values. In the special case of first-order kinetics, the continuous function c(k,0) may be obtained from C(t) using the eigenfunction method described by Ostrowsky et al.2 If the species-type density function D(k) is known, the methodology of how to obtain c(k,0) from C(t) is the subject * To whom correspondence should be addressed. E-mail: jim.riggs@ ttu.edu. † Currently with: National Center for Upgrading Technology, 1 Oil Patch Drive, Suite A202, Devon AB T9G 1A8, Canada.
∫0∞ g(k,t) dk
(2)
where g(k,t) must be interpreted as a density function, and this led to several inconsistencies. To alleviate these problems, the authors introduced a species-type density function D(k), which satisfies the following normalization condition:
1 N
∫0∞ D(k) dk ) 1
(3)
where N is the total number of species types present in the reaction mixture. Expressing g(k,t) as a product of c(k,t)D(k) in eq 2 frees the function c(k,t) from being a density function in the integral
C(t) )
∫0∞ c(k,t)D(k) dk
(4)
where c(k,t) is just the concentration of a lump with reactivity k, which can be used to formulate the material balance equation without the problems that are presented using g(k,t), as such, in eq 2. Laxminarasimhan et al.3 applied the continuum theory of lumping to develop a simulation model for hydrocracking of a
10.1021/ie0607191 CCC: $37.00 © 2007 American Chemical Society Published on Web 02/06/2007
1654
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007
VGO fraction in a plug flow reactor, and many practical applications of the continuum lumping approach were discussed by Basak et al.4 Laxminarasimhan et al.3 used a TBP profile as the characterization parameter and derived the relationship between the species-type index i and the reactivity k by
di di dθ D(k) ) ) dk dθ dk
(5)
which can be written as
xi ) ai,1c(ki-1,0) + ai,2c(ki,0) where ai,1 and ai,2 are the constant coefficients, given by
][(
[
(
(6)
where kmax and R are model parameters. The function form for D(k) satisfies the normalization criteria in eq 3. When D(k) is known, the continuous concentration function c(k,0) can be constructed from the TBP distribution using the procedure explained in section 3. 3. Construction of a Continuous Concentration Function c(k,0) The starting point for constructing a continuous concentration function for the continuum lumping model is the TBP data of the reactor feed. A TBP curve can be divided into several pseudo-components, and the weight fraction of each pseudocomponent can be assigned with a boiling range and a midboiling point. The reactivities corresponding to the normalized initial and final boiling points of a pseudo-component can be calculated using eq 6. The relationship between the weight fraction and the boiling range of the pseudo-components is itself a density function. More precisely, it is an approximation of a true density function, in the sense that a train of spikes (Dirac delta functions) are replaced by a train of rectangles while preserving the total area equal to unity. With this information, the weight fraction xi of pseudo component i, defined by reactivities ki-1 and ki, is given by
∫k
(for i ) 1, 2, ..., n)
c(k,0)D(k) dk
(7)
i-1
where n is the total number of pseudo-components derived from the TBP curve. The objective is to find the function c(k,0) (g0 for all k) so that the following normalization criteria for the material balance is satisfied:
∫0
kmax
n
c(k,0)D(k) dk )
xi ) 1 ∑ i)1
(8)
c(k,0), in the interval of ki-1 e k e ki, can be expressed using a piecewise linear interpolation form, given by
c(k,0) ) c(ki-1,0)
(
)
(
k - ki k - ki-1 + c(ki,0) ki-1 - ki ki - ki-1
)
∫kk
i
i-1
(
)
k - ki D(k) dk + ki-1 - ki k - ki-1 ki c(ki,0) k D(k) dk (10) i-1 k - k i i-1
∫
(
)
ai,2 )
[
NR (ki - ki-1)kRmax
][(
)
kiR+1 ki-1kiR R+1 R
(12a)
(
)]
(12b)
ki-1R+1 ki-1R+1 R+1 R
Equation 11 can be written for each pseudo-component present in the reactor feed, and the resulting set of n equations with n + 1 unknown variables form a linear system, as given below.
A(k)c(ki,0) ) X
[
(13)
]
The matrix A(k) and vectors c(k,0) and X are given by
a1,1 0 A(k) ) ... ... 0
a1,2 a2,1 ... ... ...
0 a2,2 ... ... ...
... 0 ... ... 0
... 0 ... 0 ... ... ... ... an,1 an,2
[]
c(k1,0) c(k2,0) l c(ki,0) ) l c(kn,0) c(kn+1,0)
[]
x1 x2 X) l l xn
(14)
The unknown variables c(ki,0) can be interpreted as the coefficients of the linear splines that form the continuous function c(k,0). Equation 13 is indeterminate and has an infinite number of solutions. Setting one of the spline coefficients to an arbitrary value might yield negative values for other spline coefficients. Thus, the calculation of spline coefficients is posed as a constrained optimization problem.
(9)
Equation 9 and the species-type distribution function D(k) are substituted into eq 7, yielding
xi ) c(ki-1,0)
)]
ki-1R + 1 kiki-1R R+1 R
and
k ) kmaxθ1/R
xi )
)
kiR+1 kiR+1 R+1 R
NR ai,1 ) (ki-1 - ki)kRmax
where θ is the normalized TBP. The authors approximated di/ dθ by the total number of species types (N) and obtained dθ/dk using a simple power-law-type equation:
ki
(11)
n
Min J(c(k,0)) )
(c(ki,0) - c(ki+1,0))2 ∑ i)1
s.t.: A(k)c(ki,0) ) X c(ki,0) g 0
(15)
where c(k,0) is a vector of dimension n + 1 and A(k) is an [n × (n + 1)] matrix. The objective function in eq 15 forces the adjacent spline coefficients to remain as close to each other as
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007 1655
Figure 1. Continuous concentration-reactivity function c(k,0) constructed from the true boiling point (TBP) curve for feedstock No. 1 with various numbers of pseudo-components.
Figure 2. Parity plot for feedstock No. 1: weight fraction from TBP versus weight fraction calculated from continuous function c(k,0) constructed with the number of pseudo-components being n ) 10.
possible, thus making smooth variations in c(k,0), with respect to reactivity k. By imposing the constraint in the form of a linear system, which arises by expressing c(k,0) as a piecewise linear interpolation function in eq 7 for each weight fraction, the condition for material balance in the continuum framework as given in eq 8 is satisfied. The proposed optimization approach has been applied to two different feedstocks: feedstock No. 15 and feedstock No. 2.6 The TBP curves generated from these feedstocks are broken into 40 pseudo-components each, and every pseudo-component is assigned a mid-boiling point and a boiling range. The values of model parameters kmax and R, estimated by Laxminarasimhan et al.3 for the two feedstocks, are used to calculate the reactivities. The resulting optimization problem is solved using sequential quadratic programming (SQP) software.7 The continuous concentration-reactivity function constructed from the TBP curve of feedstock No. 1 is shown in Figure 1 for a different number of pseudo-components. The parity plots in Figures 2-4 show how the weight fraction values calculated from the continuous function (cf. eq 7) compare with the weight fractions derived from the TBP curve. The concentrationreactivity function and the corresponding parity plots for feedstock No. 2 are shown in Figures 5-8. From the parity plots of both feedstocks, one can conclude that the fitting accuracy increases as the number of pseudo-components used to construct the continuous function c(k,0) increase. It is worth noting that if the proposed approach is used to develop a kinetic model based on continuum lumping, the parameter estimation would involve solving two nested optimization loops: one to construct c(k,0) and the other to estimate the model parameters that are relevant to reactor simulation.
Figure 3. Parity plot for feedstock No. 1: weight fraction from TBP versus weight fraction calculated from continuous function c(k,0) constructed with the number of pseudo-components being n ) 20.
Figure 4. Parity plot for feedstock No. 1: weight fraction from TBP versus weight fraction calculated from continuous function c(k,0) constructed with the number of pseudo-components being n ) 40.
Figure 5. Continuous concentration-reactivity function c(k,0) constructed from the TBP curve for feedstock No. 2 with various numbers of pseudocomponents.
The details of such an exercise applicable to the continuum lumping model for VGO hydrocracking3 are presented elsewhere.8 4. Conclusion A continuous concentration-reactivity function required for a continuum lumping approach is constructed from the true boiling point (TBP) curve, using a simple methodology based on a piecewise linear interpolation, and the resulting indeterminate linear system is solved by a constrained optimization problem. As the number of pseudo-components considered for constructing the continuous function increases, the fitting accuracy increases. The methodology has been successfully
1656
Ind. Eng. Chem. Res., Vol. 46, No. 5, 2007
Nomenclature
Figure 6. Parity plot for feedstock No. 2: weight fraction from TBP versus weight fraction calculated from continuous function c(k,0) constructed with the number of pseudo-components being n ) 14.
C(t) ) total concentration of species at time t C(0) ) total concentration of species in the feed c(k,t) ) continuous concentration-reactivity function at time t c(k,0) ) continuous concentration-reactivity function of feed c(ki,0) ) value of continuous concentration function at reactivity ki D(k) ) species-type density function g(k,t) ) density function i ) species-type index k ) reactivity kmax ) model parameter N ) total number of species types in the reaction mixture n ) number of pseudo-components xi ) weight fraction of a pseudo-component i Greek Letters θ ) normalized TBP R ) model parameter Acknowledgment This work was supported by the member companies of the Texas Tech Process Control and Optimization consortium. The authors would like to thank Professor Chang-Bock Chung (Chonnam University, South Korea) for fruitful discussions.
Figure 7. Parity plot for feedstock No. 2: weight fraction from TBP versus weight fraction calculated from continuous function c(k,0) constructed with the number of pseudo-components being n ) 20.
Figure 8. Parity plot for feedstock No. 2: weight fraction from TBP vs weight fraction calculated from continuous function c(k,0) constructed with the number of pseudo-components being n ) 40.
applied for the construction of a concentration-reactivity function for vacuum gas oil (VGO) fractions. There are two important issues that will be considered in our future work: first, whether the optimization procedure shows local minima and hence, the possibility of multiple solutions; and second, a comparison of the proposed methodology with the techniques available in the open literature.
Literature Cited (1) Chou, M. Y.; Ho, T. C. Continuum theory for lumping non-linear reactions. AIChE J. 1988, 34 (9), 1519. (2) Ostrowsky, N.; Sornette, D.; Parker, P.; Pike, E. R. Exponential sampling method for light scattering polydispersity analysis. Opt. Acta 1981, 28, 1059. (3) Laxminarasimhan, C. S.; Verma, R. P.; Ramachandran, P. A. Continuous lumping model for simulation of hydrocracking. AIChE J. 1996, 42 (9), 2645. (4) Basak, K.; Sau, M.; Manna, U.; Verma, R. Industrial hydrocracker model based on novel continuum lumping approach for optimization in petroleum refinery. Catal. Today 2004, 98, 253. (5) Bennett, R. N.; Bourne, K. H. Hydrocracking for middle distillates A study of process reactions and corresponding product yields and qualities. In Symposium on AdVances in Distillate and Residual Oil Technology, American Chemical Society, New York Meeting, August 27-September 1, 1972; pp G45-G62. (6) El-kady, F. Y. Hydrocracking of vacuum distillate fraction over bifunctional molybdenum-nickel/silica-alumina catalyst. Indian J. Technol. 1979, 17, 176. (7) Gill, P. E.; Murray, W.; Saunders, M. A.; Wright, M. H. User’s Guide for NPSOL (Version 4.0), A Fortran Package for Non-linear Programming; Systems Optimization Laboratory, Stanford, CA, 1986. (8) Govindhakannan, J.; Riggs, J. B. Continuum lumping model for hydrocracker simulation, Internal Report, Texas Tech University, Lubbock, TX, 1998.
ReceiVed for reView June 6, 2006 ReVised manuscript receiVed December 16, 2006 Accepted December 27, 2006 IE0607191
