Chemical Engineering Science 62 (2007) 4510 – 4528 www.elsevier.com/locate/ces A new generic approach for the modeling of fluid catalytic cracking (FCC) riser reactor Raj Kumar Gupta a , Vineet Kumar a,∗ , V.K. Srivastava b a Department of Chemical Engineering, Thapar University, Patiala 147 004, India b Department of Chemical Engineering, Indian Institute of Technology, Delhi 110 016, India Received 17 April 2006; received in revised form 4 November 2006; accepted 14 May 2007 Available online 21 May 2007 Abstract A new kinetic model for the fluid catalytic cracking (FCC) riser is developed. An elementary reaction scheme, for the FCC, based on cracking of a large number of lumps in the form of narrow boiling pseudocomponents is proposed. The kinetic parameters are estimated using a semi-empirical approach based on normal probability distribution. The correlation proposed for the kinetic parameters’ estimation contains four parameters that depend on the feed characteristics, catalyst activity, and coke forming tendency of the feed. This approach eliminates the need of determining a large number of rate constants required for conventional lumped models. The model seems to be more versatile than existing models and opens up a new dimension for making generic models suitable for the analysis and control studies of FCC units. The model also incorporates catalyst deactivation and two-phase flow in the riser reactor. Predictions of the model compare well with the yield pattern of industrial scale plant data reported in literature. 䉷 2007 Elsevier Ltd. All rights reserved. Keywords: Fluid catalytic cracking; Mathematical modeling; Simulation; Cracking kinetics; Pseudocomponents; Computational fluid dynamics 1. Introduction Fluid catalytic cracking (FCC) is a process in which the heavy hydrocarbon molecules are converted into lighter molecules. The hydrocarbon feed enters a transport bed tubular reactor (riser) through feed atomizing nozzles and comes in contact with the hot catalyst coming from the regenerator. The feed gets vaporized and cracks down to the lighter molecules as it travels upwards along with the catalyst. As a result of cracking, the velocity of the vapors increases along the riser height. Coke, the byproduct of cracking reactions, gets deposited on the catalyst surface thus causing the catalyst to loose its activity. The cracked hydrocarbon vapors are separated from the deactivated catalyst in a separator; the vapors adsorbed onto the surface of the catalyst are also stripped off using steam in the catalyst stripper. The cracked hydrocarbon vapors are sent to the main distillation column for further ∗ Corresponding author. Tel.: +91 175 2393063. E-mail address: vikumar@tiet.ac.in (V. Kumar). 0009-2509/$ - see front matter 䉷 2007 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2007.05.009 separation into various fractions, and the deactivated catalyst flows into the regenerator. In the regenerator, the coke deposited on the catalyst surface is burned off to regenerate the catalyst. The catalyst also becomes hot during the regeneration process. This hot-regenerated catalyst is recycled back to the riser reactor. Thus the catalyst acts as a heat carrier also and provides the heat required for endothermic cracking reactions in the riser reactor as well as the heat required for the vaporization of feed. Detailed modeling of the riser reactor is a challenging task for theoretical investigators not only due to complex hydrodynamics and the fact that there are thousands of unknown hydrocarbons in the FCC feed but also because of the involvement of different types of reactions taking place simultaneously. It is believed that catalytic cracking begins with the formation of + carbenium ion (R1–CH+ 2 –R2 or R–CH3 ) by the interaction of olefin molecules with the acidic site on the catalyst followed by the beta scission of the carbenium ion. In the beta scission reaction, the  bond of carbenium ion breaks to form an olefin and a new carbenium ion. The carbenium ion formed by beta scission can undergo further cracking reaction. The olefin can also be cracked further after being converted to a R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 carbenium ion through hydrogen addition. Thus, large hydrocarbon molecules can be cracked repeatedly producing successively smaller hydrocarbons. However, as the hydrocarbon chain becomes smaller, cracking rates become slower. At the same time, a short, less reactive ion can transfer its charge to a larger, more reactive molecule by hydrogen transfer that makes the cracking reaction of larger molecules still faster (Wilson, 1997). Such detailed chemistry of catalytic cracking coupled with a large number of unknown compounds present in the feedstock is very difficult to be used in the mathematical modeling of an industrial scale FCC riser reactor because of the analytical and computational limitations. The traditional and global approach of modeling of cracking kinetics is based on lumping of compounds. Mathematical models dealing with riser kinetics can be categorized into two main types. In one category, the lumps are made on the basis of boiling range of feed stocks and corresponding products in the reaction system. This kind of model has an increasing trend in the number of lumps of the cracked components. The other approach is that in which the lumps are made on the basis of molecular structure, characteristics of hydrocarbon group composition in reaction system. This category of models emphasizes on more detailed description of the feedstock (Wang et al., 2005). These both categories of models do not include chemical data such as type of reaction and reaction stoichiometry. The number of kinetic constants in these models increases very rapidly with the number of lumps. All these models assume that FCC feed and products are made of a certain number of lumps, and kinetic parameters for these lumps are estimated empirically considering the conversion of one lump to the other. In both of these categories, however, reaction kinetics being considered is that of ‘conversion’ of one lump to another and not the ‘cracking’ of an individual lump. As the reactivity of the feed depends on the composition of the feed, and catalyst composition and reactivity is important in determining conversion and product selectivity, the values of kinetic constants obtained by the above discussed models depend on the particular pair of feedstock and catalyst for which they are obtained. Hence, these kinetic constants cannot be used for different feed and catalyst pairs. Other modeling schemes include, models based upon reactions in continuous mixtures (Aris, 1989), structure oriented lumping (Quam and Jaffe, 1992), and ‘single-events’ cracking (Feng et al., 1993). Nevertheless, the application of these models to catalytic cracking of industrial feedstocks (vacuum gas oil), is not realized because of the analytical complexities and computational limitations. Liguras and Allen (1989a) proposed a lumped kinetic model so as to utilize the pure components cracking data for the catalytic cracking of oil mixtures. The authors in their subsequent work (Liguras and Allen, 1989b) divided the petroleum feedstock into a number of pseudocomponents. These pseudocomponents were characterized by grouping the feedstock components into compound classes and selecting a set of representative compounds in a compound class and then assigning the concentrations to the representative compounds. This pseudocomponents characterization method 4511 requires extensive use of analytical procedures and suffers from the same drawbacks that are cited above. A more detailed compilation of various lumping schemes, and their use in advanced modeling and control of FCC unit are presented in our recent work (Gupta et al., 2005). Even today, in the advanced models of the FCC riser (models considering various aspects of reactor modeling) three lump or four lump kinetic schemes are being used by the investigators to avoid the mathematical complexities and load of computation that will be there if more number of lumps are considered in the kinetic scheme as the number of cracking constants increases rapidly with the number of lumps. With this in view, in the present work, a new approach of kinetic scheme, considering a large number of lumps (the lightest being methane and the heaviest containing all compounds of the feedstock having boiling point close to its end-point), for the FCC riser is introduced. These lumps are characterized using constant Watson characterization factor, which is an indicator of the feedstock composition. Characterization was made in such a way that all physico-chemical properties of the lumps are known and they can be treated as hypothetical pure components. The proposed model considers that each lump on cracking gives two other lumps in one single reaction step. The proposed model falls under the category of models in which lumps are formed on the basis of boiling point, but in this approach, each individual lump is considered as a pure component with known physicochemical properties. A separate coke lump is also considered and it is assumed that when one mole of a lump cracks down it gives one mole each of two other lumps and the balance material gives the coke. The proposed model also incorporates two-phase flow and catalyst deactivation. Since a new cracking reaction mechanism is introduced, a new semi-empirical approach based on normal probability distribution is also proposed to estimate the cracking reactions’ rate constants. The proposed correlation contains four parameters that depend on the feed characteristics, catalyst activity, and coke forming tendency of the feed. This eliminates the need of estimating the rate constants by regression analysis. Stangeland (1974) had proposed a kinetic model, for the prediction of hydrocracker yields, in which the feedstock was divided into a series of 50 ◦ F boiling range cuts without characterizing these cuts. In the present work we have characterized the pseudocomponents as discussed in Appendix A. 2. Riser model The riser is modeled as a vertical tube comprising of a number of equal sized compartments (or volume elements) of circular cross section. Volume elements are designated by symbol j (j = 1, 2, . . . , NC ) and numbering of the volume elements is done from bottom (inlet) to top (outlet). Each volume element is assumed to contain two phases (i) solid phase (catalyst and coke) and (ii) gas phase (vapors of feed and product hydrocarbon, and steam). In one volume element, each phase is assumed to be well mixed so that heat and mass transfer resistances can be ignored. Model equations are written for both the phases in each j th volume element for all pseudocomponents 4512 R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 PCi (i = 1, 2, . . . , N ). Cracking reactions’ rate constants in each volume element are evaluated at the local temperature of the j th volume element and are subsequently used for calculating the change in molar concentration of each component and the heats of reactions. This approach of finite volume method (FVM) is widely used in computational fluid dynamics (CFD) and has been elaborated elsewhere for membrane separation process by Kumar and Upadhyay (2000). It is assumed that known amount of reactants enter the j th volume element, cracking reactions take place in the element for a period equal to the residence time of the hydrocarbons entering the element. The concentration of the reacting mixture at the outlet of j th volume element serve in defining the feed for the (j + 1)th volume element. To start computation from volume element number one, i.e., the inlet of the riser reactor, the known process parameters at the inlet of the riser are used. To simplify this FVM based approach for solving material and energy balance equations in FCC riser reactor, following commonly used assumptions were made. 2.1. Assumptions • At the riser inlet, hydrocarbon feed comes in contact with the hot catalyst coming from the regenerator and instantly vaporizes (taking away latent heat and sensible heat from the hot catalyst). The vapor thus formed moves upwards in thermal equilibrium with the catalyst (Corella and Frances, 1991; Martin et al., 1992; Fligner et al., 1994; Ali et al., 1997; Derouin et al., 1997). • There is no loss of heat from the riser and the temperature of the reaction mixture (hydrocarbon vapors and catalyst) falls only because of the endothermicity of the cracking reactions (Corella and Frances, 1991; Ali et al., 1997; Theologos et al., 1999; Gupta and Subba Rao, 2001). • Ideal gas law is assumed to hold while calculating gas phase velocity variation on account of molar expansion due to cracking and gas phase temperature (Gupta and Subba Rao, 2001). • Catalyst particles are assumed to move as clusters to account for the observed high slip velocities (Gupta and Subba Rao, 2001). • Heat and mass transfer resistances are assumed as negligible (Corella and Frances, 1991; Martin et al., 1992; Ali et al., 1997; Derouin et al., 1997; Theologos et al., 1999). • Both phases are assumed in plug flow condition hence back mixing in both phases is neglected. 2.2. Kinetic modeling Most of the kinetic schemes currently being used for the modeling of FCC riser reactor are based on the specified number of lumps (four lump and 10 lump schemes being the most common). In the proposed model, any suitably large number of lumps can be considered. However, in the present work number of lumps, N, is taken as 50. Each of these lumps is characterized on the basis of average boiling point and specific gravity, and treated as hypothetical pure component with all physical, thermodynamic, and critical properties known. A brief description for generating these pseudocomponents is given in Appendix A. For the generation of these lumps, the feedstock is divided into 12 lumps and average boiling point of each lump is determined by area averaging of the true boiling point (TBP) curve of the feedstock. Thus one individual lump is a collection of all compounds having boiling point close to each other. The temperature range between boiling point of the lightest lump of the feedstock and the boiling point of propane is divided into 31 equal intervals which were considered to be boiling points of 31 lumps of liquid products to be formed after cracking. Molecular weight of each lump is also calculated using empirical correlation proposed by Edmister and Lee (1984) assuming constant Watson characterization factor for product and feed lumps. Constituents of dry and wet gas are also considered as seven lumps consisting of methane, ethane, propane, butene, butane, pentene, and pentane. Properties of these seven lumps were taken equal to that of pure components due to the fact that empirical correlations fail to predict properties of very light hydrocarbons and also due to the fact that gaseous products are rather known mixtures of these hydrocarbons. Thus we have lumps of seven pure components and 43 hypothetical components. Henceforth, all these components are collectively called pseudocomponents and abbreviated as PCs. Apart form these 50 components, a separate coke lump is also considered which is formed as a byproduct of the cracking reactions. The pseudocomponent based approach for design and simulation of crude distillation unit is highly successful. Sufficient empirical correlations are available in literature to predict almost all physico-chemical, thermodynamic, and critical properties of pseudocomponents. Vapor–liquid equilibrium of these pseudocomponents are predicted fairly accurately without resolving these components in aromatics, olefins, naphthenes, and paraffins. For applying this approach to FCC modeling, it is assumed that one mole of a pseudocomponent on cracking gives 1 mole each of two other smaller pseudocomponents and some amount of coke may also form. A schematic diagram of the reaction mechanism is given in Fig. 1. There are total N blocks in each column of the diagram and each block in a column represents one pseudocomponent (blocks in a row represents same pseudocomponent). Pseudocomponents are numbered in increasing order of normal boiling point, which also ensures increasing order of molecular weight. According to the proposed scheme, there are several possible ways through which one particular pseudocomponent (say ith pseudocomponent, PCi ) can crack down to give a pair of pseudocomponents PCm and PCn along with some amount of coke as cracking byproduct according to the following pseudoreaction mechanism: ki,m,n PCi → PCm + PCn + i,m,n , (1) where i, m, and n are pseudocomponents’ numbers, i,m,n is the amount of coke formed (kg) when one kmol of ith pseudocomponent cracks to produce one kmol each of mth and nth pseudocomponents. 4513 R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 Reactant Feed Product 1 + Product 2 Mol. Wt. MW1 MW2 MW3 1 2 3 1 2 3 1 2 3 MWm m m m MWn n n n MWi i { MWN-11 MWN PCi ki , m ,n PC m + PC n + coke N-11 N-11 N-11 N N N + coke Amount of coke (gram) =MWi - (MWm+ MWn) Fig. 1. Schematic diagram of reaction mechanism. The proposed pseudoreaction mechanism is different from the conventional reaction mechanism. In Eq. (1) PCi is the ith pseudocomponent consisting of several hundreds of compounds having molecular weight close to MWi . After cracking, each individual molecule of PCi gives at least two product molecules of molecular weight less than the cracking molecule. There is a wide range of possibilities in which one particular molecule of PCi can crack. In one case the products may be of widely different molecular weights (i.e., one molecule is of small molecular weight and other of high molecular weight close to the original molecule). In other case both the product molecules can be of almost equal molecular weight. Also, amongst these reactions, some may be coke forming and some may not be coke forming reactions. Eq. (1) represents the collection of only those reactions in which each compound of molecular weight MWi gives two molecules of molecular weights MWm and MWn along with the associated coke formed in this process. Here it should be noted that for incorporating all possible cracking reactions of ith pseudocomponent, value of m varies from 1 to i as no product can be heavier than the reactant. Similarly, value of n ranges from 1 to m only as the products PCm and PCn are interchangeable. Therefore, in the proposed mechanism, there are (i + 1) × i/2 possible ways in which an ith pseudocomponent can crack. For the law of mass conservation, it is required that mass of reactant should be equal to the mass of products in one reaction step. Therefore, when MWi kg of PCi cracks to give MWm kg of PCm and MWn kg of PCn , rest of the material gives coke (i,m,n ). Thus, the value of i,m,n can be calculated by taking difference in molar masses of reactant and product hydrocarbons as follows: i,m,n = MWi − (MWm + MWn ). (2) In Eq. (1) there are (i + 1) × i/2 possible ways in which cracking reaction for one pseudocomponent PCi can be written. Only two of such possible reactions are shown by solid line in Fig. 1. The dashed line represents the same reaction, as the values of m and n in Eq. (1) are interchangeable. Thus the total number of reactions that could be written for the complete set of reactions becomes N × (N − 1) × (N − 2)/3 (for N = 50, total number of reactions become 39 200). However, all of these reactions are not feasible. The feasibility of a reaction is found using the stoichiometry of the cracking reaction (Eq. (1)). Only those reactions are considered feasible for which the value of i,m,n calculated by Eq. (2) is either zero or positive. This natural constraint greatly reduces the number of feasible cracking reactions but still the number is too large (more than 10 000 in the present case with N = 50). To handle such a large number of reactions in riser reactor, special considerations for the estimation of reaction rate constants and technique for solving material and energy balance equations are required. In the present work, a new semiempirical scheme for the estimation of rate constants is developed. This scheme makes the kinetic model more versatile. In the beginning, six parameters were introduced to adjust more than 10 000 reaction rate constants (ki,m,n ) needed to explain complete reaction mechanism in a typical FCC riser reactor. Later, it was observed that only four of these parameters are significant. However, for the pedagogical point of view, this correlation is discussed in detail in Appendix B in the same sequence as it was developed. The final form of this correlation is given by the following equation:  2 ki,m,n = [k0 MWi ]e−1 (MWm −MWn )  −i,m,n /2  − e−MWi −E0 MWi /RT e ×e . 1 − e−MWi (3) In the above equation, parameters to be estimated by using experimental data are k0 , , E0 , , 1 , and 2 . Later it was observed that two parameters  (correlating frequency factor in terms of molecular weight of pseudocomponent) and 1 (an indicator of cracking tendency of pseudocomponents from 4514 R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 middle or from sides) are insignificant, thus the correlation reduces to the following four parameter equation:  −i,m,n /2  − e−MWi e  ki,m,n = k0 e−E0 MWi /RT . (4) 1 − e−MWi Therefore, the parameters required to be estimated are k0 , E0 , , and 2 . All these parameters are constant for a particular pair of feed and catalyst. Parameter 2 correlates the coke forming tendency of the feed and parameter  correlates activation energy of individual pseudocomponent in terms of its molecular weight. The theoretical kinetic scheme developed in the present work assumes that all the cracking reactions are of first order. This is in agreement with the reported cracking rate of gas lump and gasoline lump. Many researchers have reported the cracking of heavier fractions of feed such as gas oil lump to be of second order due to deeper cracking action. In the present case feed is broken into large number of narrow boiling lumps (pseudocomponent), therefore, first order cracking can be assumed safely even for heavier lumps. Thus for first order reactions, rate of disappearance of ith pseudocomponent due to cracking in j th volume element through one reaction as indicated in Eq. (1) is given by ri,m,n = j · ki,m,n · Ci,j · (Mcat · t Cat j ), (5) equations. Using higher order numerical techniques, such as Runge–Kutta method, is also not feasible for solving these ODEs due to excessive computational load. Therefore, an FVM based approach was used for computing material balances. As discussed earlier, in this FVM approach it is assumed that the riser reactor is made up of a large number of volume elements of circular cross-section and of very small height, placed one over the other. Volume elements are numbered from bottom to top, i.e., from inlet to outlet (Fig. 2). Both gas phase and solid phase are moving through each volume element in upward direction, with different velocity. Since the height of a typical j th volume element zj is very small, therefore the residence-time for the gas phase, tj , and the residence-time of the solid phase, t Catj , is also very small. Applying concepts of finite volume approach, it is assumed that all reactions take place for tj time at a constant rate determined by the prevailing temperature, pressure, and concentration at the inlet of j th volume element. After elapse of time tj the temperature and concentration of the stream (which is now outgoing stream from j th element) are determined by the following material balance equations. Material balance over j th volume element for the gas phase: Rate of mass in from the (j − 1)th element − rate of mass out from j th element where the concentration of ith pseudocomponent at the inlet of j th volume element is given by Ci,j = Pi,j . (ug,j · g,j )Ar (6) B +1 , B + exp(A · Ccj ) or N  i=1 Pi,j in Eq. (6) is kmol of ith component entering per second into the j th volume element. Ar , ug,j , and g,j are area of cross-section of the riser, local gas velocity, and gas phase volume fraction in the j th volume element. t Cat j is the residence time of catalyst in the j th volume element, hence Mcat · t Cat j in Eq. (5) is the mass of catalyst present in j th volume element. Non-selective deactivation of catalyst, because of coke deposition, is assumed. The activity coefficient (j ) depends on the coke concentration on the catalyst. Pitault et al. (1995) proposed the following correlation for the estimation of activity coefficient: j = = rate of mass converted to coke in j th element (7) The values for deactivation constants A and B reported by the authors are 4.29 and 10.4, respectively; the same are used in this work. Ccj is the concentration of coke on catalyst surface (wt%). 2.3. Material balance Calculation of material balances with several thousands of parallel reactions is practically impossible by analytical solution of ordinary differential equations (ODE) representing rate =  Pi,j −1 MWi − N  m i   i=1 m=1 n=1 N  (8) Pi,j MWi i=1 ri,m,n i,m,n  . (9) for all feasible reactions Similarly, material balance over j th volume element for solid phase: Rate of mass out from j th element − rate of mass in from the (j − 1)th element = rate of mass accumulated in j th element = rate of mass of coke formed or (Mcat + Mcokej ) − (Mcat + Mcokej −1 ) N i m     = ri,m,n i,m,n i=1 m=1 n=1 (10) , (11) for all feasible reactions where Mcokej is the cumulative mass of coke formed upto j th volume element. For the first volume element (j = 1), Mcokej −1 is the mass of coke on regenerated catalyst. In Eqs. (9) and (11) summation for i is done from 1 to N. However, summation for m is done from 1 to i only because no product species can have molecular weight greater than the reactant species, whereas summation over n is done from 1 to m only as the products PCm and PCn are interchangeable. 4515 R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 Fig. 2. A typical volume element in the riser reactor. 2.4. Heat balance It is assumed that the hydrocarbon feed, after coming into contact with hot catalyst from the regenerator, vaporizes completely in the first volume element of the riser. Afterwards the gas phase is in thermal equilibrium with the catalyst all along the reactor. At the entrance (the first volume element of the riser reactor) the equilibrium temperature, Tin , is calculated by taking into account the specific heat of catalyst and the latent heat of vaporization of the feed. The pseudocomponents are chemically stable hypothetical components, and according to the proposed reaction mechanism one pseudocomponent produces two other chemically stable pseudocomponents. This indicates that the proposed reaction mechanism includes all possible reactions taking place during cracking (such as hydrogen-transfer, condensation, isomerization, etc.) that leads to the formation of chemically stable products in the form of pseudocomponents. Therefore, the overall heat of reaction (Hr ), of Eq. (1) can be estimated by finding the difference between the heat of combustion of all products and the heat of combustion of the reactant. Heat of combustion of pseudocomponents, H comb, is estimated by the following equations in terms of the API gravity of the hydrocarbon: H combi = 133.976 APIi + 41170 for APIi 25, (12) H combi = − 0.4017 API2i + 57.859 APIi + 1953.14 ln(APIi ) + 37037.09 for 25APIi 50, (13) H combi = 55.824 APIi + 43775.32 for APIi 50. (14) Eqs. (12)–(14) were obtained by curve fitting of graphical data proposed by Maxwell (1968). Curve fitting was done in such a way that there is no discontinuity in the heats of combustion values predicted by these three equations. Thus for the cracking of ith pseudocomponent, giving mth and nth pseudocomponents, heat of reaction becomes H r i,m,n = i,m,n · H coke + (MWm H combm + MWn H combn − MWi H combi ). (15) Thus the energy balance equation for the j th volume element can be written as  Mcat · Cpcat + Mcokej −1 · Cp coke + Mst · Cp st + N  i=1 =  Pi,j −1 MWi · Cp i (Tj −1 − Tj )   i m n ri,m,n · H r i,m,n  (16) . for all feasible reactions Rearranging the above energy balance equation, we get Tj = Tj −1 −    ( i m n ri,m,n · H r i,m,n )for (Mcat · Cp cat + Mcokej −1 · Cp coke + Mst all feasible reactions  · Cp st + N i=1 Pi,j −1 MWi · Cp i ) . (17) 4516 R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 The temperature (Tj ) thus calculated was used for the estimation of kinetic parameters in the next volume element and other temperature dependent properties of the pseudocomponents. Re = 2.5. Hydrodynamics Eq. (19) can be solved as a difference equation and the value of solid phase velocity in the next volume element can be calculated with the initial condition of solid phase velocity at the riser entrance calculated by the following equation: In the proposed riser model two phases (cluster phase and gas phase) are considered. Cluster phase includes the loosely held particles of catalyst and the coke. The cluster phase and gas phase hold up vary along the riser height. Solid particles spend more time in the riser than hydrocarbon vapor due to slip between the two phases. The slip velocities observed in the riser are higher than the terminal settling velocity of a single particle. The reason for the higher slip velocities is attributed to particles moving in clusters (Subba Rao, 1986; Fligner et al., 1994; Horio and Kuroki, 1994). The clusters are agglomerates of loosely held particles (Fligner et al., 1994). Cluster voidage ( c ) is assumed to be 0.5, in line with the two-phase theory of fluidization. As proposed by Tsuo and Gidaspow (1990), solid phase momentum balance along the riser height may be written as d( c c uc uc ) = Cf (ug − uc ) + 2fs c dz 2 c uc − c c g. (18) Assuming the change in the mass of solids along the riser height as negligible (this assumption is valid as the change in the mass of solid phase from riser entrance to riser outlet due to coke deposition is typically less than 1%, hence in a volume element of the riser the change in the mass of solids will be negligibly small) solid phase continuity equation can be written as 2fs c c u2c duc = Cf (ug − uc ) + − c c uc dz D c c g. (19) The frictional force per unit volume at the gas particle interphase due to differing phase velocities can be calculated by the following expression (Markatos and Shinghal, 1982): F = 0.5CD AP g |ug − uc |(ug − uc ) = Cf (ug − uc ), (20) where AP is total projected area of particles per unit volume, CD is interphase friction coefficient between the two phases. Projected area per unit volume can be calculated based on equivalent spherical diameter as AP = 1.5 c /dc . (21) In the above equations dc is cluster diameter, 6.0 × 10−3 m (Fligner et al., 1994), and c is cluster density. The cluster density can be approximated by the following expression: c= p (1 − c) + g c p (1 − c ). (22) The empirical correlations for CD used by Arastoopour and Gidaspow (1979) are uc0 = g g |ug − uc |dc g . Mcat , c Ar c0 (25) (26) where uc0 and c0 are the values of solid phase velocity and solid phase volume fraction at the entrance of first volume element. The value of cluster volume fraction for the next volume element is calculated by the equation cj = Mcat . c Ar ucj (27) Gas phase volume fraction is obtained using the relation c + g = 1. (28) Having obtained the values of the solid phase velocity and cluster volume fraction for the second volume element, the total pressure drop for the first volume element can be calculated. The pressure drop is assumed to be composed of four main components (Pugsley and Berruti, 1996):      dP dP dP dP dP = + + + , dz total dz s dz acc dz f s dz fg (29) where (dP /dz)s is the pressure drop due to the hydrostatic head of the solids, (dP /dz)acc is the pressure drop due to solids acceleration, and (dP /dz)f s and (dP /dz)fg are the pressure drops due to solids friction (defined as the frictional force per unit volume between solids and wall) and gas friction (defined as frictional force per unit volume between the gas and the solids), respectively. These components can be calculated using the following relations:  dP = c gc , (30) dz s  dP c u2c , (31) = c dz acc 2 z  2fs c c u2c dP , (32) = dz f s D  dP dz fg = 2 g ug fg D . (33) Blasius friction factor given by the following empirical equation is used as the gas friction factor fg = 0.316Re−1/4 . (34) 24 CD = (1 + 0.15Re0.687 ) for Re < 1000, Re (23) Konno’s correlation is used to calculate the solids friction factor CD = 0.44 (24) fs = 0.0025u−1 c . for Re1000, (35) R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 The pressure in the next volume element is obtained by the following relation:  dP Pj = Pj −1 − z. (36) dz total The gas phase density is calculated from the ideal gas law as    Pj N i=1 yi,j MWi (37) gj = RTj and the gas phase velocity is calculated by    Mst + N i=1 Pi,j MWi . ugj = Ar g,j gj (38) Residence time of catalyst in volume element j is t Catj = (z)j . uc,j (39) 2.6. Model solution Model simulation is done on a P-IV computer, which took less than 5 min for giving the simulation results. Although we can take variable height of volume elements, however, in the present work the height of each volume element of the riser was kept 10 mm. Further decrease in the height of the volume element had no appreciable effect on the results. 3. Results and discussions The material balance equations were combined with reaction kinetics and the hydrodynamic model equations to obtain the moles of each pseudocomponent coming out of any volume element j (=1 to Nc) of the riser reactor. Thus, the model could predict the yield pattern along the riser height. Model validation is done as three case studies by using industrial data reported in the literature. Results of the proposed model can be obtained by adjusting six parameters as discussed earlier. Numerical values of these parameters were obtained separately for each case study. Results of the simulator are being discussed in the following case studies. Due to the unavailability of the feed TBP data for each case, boiling point characteristic of feed (simulated distillation, SD) reported by Pekediz et al. (1997) is used (Table A1) for all the case studies. Various other parameters common to all cases are given in Table 1 and plant data used in these cases are presented in Table 2. 3.1. Case study 1 Industrial FCC plant data reported by Ali et al. (1997), presented in Table 2, was used in this case study. In this case five industrial data at the riser outlet—gasoline yield, gas yield, unconverted hydrocarbon, coke yield, and riser outlet temperature—are available to obtain the rate constant parameters. 4517 Model results were obtained by adjusting four rate constant parameters and two tuning parameters of Eq. (3) in such a way that the deviation in predicted and actual data, at the riser outlet, is minimum. To achieve this, these parameters (k0 , , E0 , , 1 , and 2 ) were determined by line-search technique followed by golden-section method. To estimate the value of these parameters, an iterative method was adopted in which one parameter was varied at a time (keeping other five constant). Starting with an initial guess value of all six parameters, a linesearch was made (by changing one parameter with a constant increment/decrement) to find a condition where the sum of absolute value of deviation in the predicted and experimental values of product yield and temperature at riser outlet are minimum. Then the value of this parameter was fine tuned between two consecutive values considered during line-search by the method of golden-section. Since the system of equations under consideration is highly nonlinear, on-line graphical observation of reported and predicted data was made almost after each iteration. The simulator results, yield pattern and temperature profile, are compared with the industrial data at the riser outlet in Figs. 3 and 4, respectively. The values of six parameters are k0 = 0.01,  = 0.01, E0 = 1540,  = 0.43, 1 = 0.0, 2 = 17.0. Fig. 4 indicates that as the hydrocarbons and catalyst mixture travel upwards, the temperature inside the FCC riser reactor decreases because of the endothermic cracking reactions. The catalyst temperature at the inlet of the riser (960 K) falls sharply to 880 K because sensible heat of catalyst coming from the regenerator is utilized in providing heat for raising the sensible heat of the feed, for vaporizing the feed, and for further heating of the vaporized feed. Afterwards, within first 10 m height inside the riser reactor the temperature drops from 880 to 790 K, as most of the cracking takes place within first 10 m of the riser height. The temperature at the outlet of the riser is 774 K (Fig. 4). The decrease in reaction mixture’s temperature and catalyst activity along the riser height cause a decline in the reaction rate, hence the temperature gradient falls appreciably with the increasing riser height. Sensitivity analysis: Sensitivity analysis of these parameters is presented in Figs. 5–7. Sensitivity analysis is the process of varying the parameters over a wide range about the mean value and recording the relative change in predicted gas, gasoline, and coke yields (Figs. 5–7). The sensitivity of one parameter relative to other is also demonstrated in these figures. Such analysis is useful in cases when the experimental data are few in number and statistical analysis cannot be applied to predict confidence interval. A detailed discussion about the sensitivity analysis is given by Saltelli (2000). Figs. 5–7 show that in the proposed model, gas yield and coke yield, are strong functions of E0 , , and 2 , and gasoline yield is very sensitive to E0 , and  and moderately sensitive to k0 , and 2 , whereas,  is almost insensitive parameter. It is evident from these figures that  has virtually no effect on any product yield. This indicates that the frequency factor (k0,i ) is independent of molecular weight. Also, the value of parameter 1 is zero. Therefore only four parameters may be used to match the model results with the industrial data, with 1 = 0, and  = 0. Hence for all the 4518 R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 Table 1 Parameters used for the simulation of riser reactor Parameter Value Source Heat of combustion of coke Molecular weight of coke Volume fraction of clusters at inlet Specific heat of catalyst Specific heat of steam Mass flow rate of steam Feed temperature at the riser inlet Latent heat of feed vaporization Catalyst particle density Catalyst particle diameter Specific gravity of feed Cluster diameter −32950 kJ/kg 12 kg/kmol 0.5 1.15 kJ/kg K 2.15 kJ/kg K 1.33 kg/s 494 K 96 kJ/kg 1200 kg/m3 75 m 0.9292 g/cm3 6 mm Austin (1984) Arbel et al. (1995) Gupta and Subba Rao (2001) Ali et al. (1997) Blasetti and de Lasa (1997) Blasetti and de Lasa (1997) Ali et al. (1997) Gupta and Subba Rao (2001) Gupta and Subba Rao (2001) Gupta and Subba Rao (2001) Pekediz et al. (1997) Fligner et al. (1994) Table 2 Plant data used for simulation of riser reactor Riser height Riser diameter Riser pressure Catalyst temperature Feed rate Feed temperature C/O ratio Ali et al. (1997) (Case study—1) Derouin et al. (1997) (Case study—2) Theologos and Markatosa (1993) (Case study—3) 33 m 0.8 m 2.9 atm 960 K 20 kg/s 496 K 7.2 (32 m) 1.0 m 3.15 atm (960 K) 85 kg/s 650 K 5.53 50 m 1.24 m (2.9 atm) 1025 K 17.5 kg/s 568 K 8.0 Data given in the parentheses are those used in the present work in place of data either not reported or reported in ranges. a Literature data. 1000 120 Plant data Gas yield Gasoline yield Coke yield Unconverted 100 Model prediction Gas yield Gasoline yield Coke yield Unconverted Plant data (Source: Ali et al., 1997) 950 Riser temperature (K) 80 Yields (wt%) Model predictions 60 40 900 850 20 800 0 750 0 5 10 15 20 Riser height (m) 25 30 35 0 5 10 15 20 Riser height (m) 25 30 35 Fig. 3. Case study 1, comparison with the data reported by Ali et al. (1997). Fig. 4. Axial temperature profile along the riser height. subsequent simulations only four parameters (k0 , E0 , , and 2 ) are used. Case study 1 with four parameters: The results for case study 1 using the four parameters are presented in Figs. 8–11. The values of the four parameters are k0 =0.01, E0 =1540, =0.43, and 2 = 17.0. The product yield profiles given in Figs. 3 and 8 match very closely. Also, the temperature profiles obtained in Figs. 4 and 9 are similar. 4519 R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 Fig. 7. Sensitivity analysis for coke yields. Fig. 5. Sensitivity analysis for gas yields. 120 Plant data 100 Yields (wt%) 80 Model prediction Gas yield Gasoline yield Coke yield Gas yield Gasoline yield Coke yield Unconverted Unconverted 60 40 20 0 0 5 10 15 20 25 30 35 Riser height (m) Fig. 6. Sensitivity analysis for gasoline yields. The activity of the catalyst also decreases rapidly as byproduct (coke) of the cracking reactions gets deposited on the catalyst surface (Fig. 10). Fig. 11 shows an initial decline in the gas velocity because of the sharp increase in the gas void fraction due to increase in the moles of the gas as a result of cracking. After this initial decline, the gas velocity starts increasing as the cracking reactions along the riser height continues to increase the moles of the gas causing a continuous decline in the gas phase density. The initial sharp increase in the catalyst velocity is due to the sharp fall in the solid volume fraction and drag exerted by the gas. After this initial sharp increase the Fig. 8. Case study 1, comparison with the data reported by Ali et al. (1997) with four parameters. catalyst velocity keeps on increasing gradually all along the riser height. Also, high values of slip factor are predicted in riser entry zone which gradually decreases along the riser height and finally reaches at 1.8 (Fig. 11). Catalyst volume fraction falls from about 0.5 to 0.1 in first few meters of riser height (Fig. 12). This sharp decline can be attributed to the fact that most of the cracking takes place within first few meters of the riser height. Furthermore, the catalyst volume fraction along the riser height is plotted for two initial values (0.5 and 0.3). 4520 R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 0.55 14 Plant data 12 0.50 Model prediction 2 950 10 0.45 Velocity (m/s) Riser temperature (K) Gas phase molar flux (kmol/m s) 1000 900 8 0.40 6 0.35 850 4 Gas phase velocity Catalyst velocity 2 800 0.30 Gas phase molar flux 0 0.25 0 5 750 10 15 20 25 30 35 Riser height (m) 0 5 10 15 20 25 30 35 Fig. 11. Predicted catalyst and gas velocity profiles and gas phase molar flux along the riser height. Riser height (m) Fig. 9. Axial temperature profile along the riser height with four parameters. 0.50 1.2 Initial cluster volume fraction = 0.5 0.45 Initial cluster volume fraction = 0.3 1.0 Catalyst volume fraction 0.40 Catalyst activity 0.8 0.6 0.4 0.35 0.30 0.25 0.20 0.15 0.10 0.2 0.05 0.00 0.0 0 0 5 10 15 20 25 30 35 5 10 15 20 25 30 35 5 Riser height (m) Riser height (m) Fig. 12. Predicted catalyst volume fraction along the riser height. Fig. 10. Predicted catalyst activity along the riser height. The plots for both the values are almost similar because this value gets adjusted very quickly at the riser entrance itself (in the first 2 m of the riser height itself the catalyst volume fraction value reaches 0.12 for both the cases) and hence the yield profiles remain unaffected. Since plant data for the product yields were available at the riser outlet only, few more comparisons were made in subsequent case studies. Even without changing the values of the parameters the results of the simulation were encouraging. However, for better comparison, the parameters were adjusted further in both the following case studies. It is observed that the parameters k0 , and 2 are the two parameters those were needed to be adjusted for different cases. 3.2. Case study 2 In this case FCC plant data (Table 2) reported by Derouin et al. (1997) was used to compare the simulator predictions. Authors have reported the product data for gasoline yield and conversion at different positions along the riser height (Fig. 13). R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 ous transport equations. Riser operating conditions given in Table 2 (Theologos and Markatos, 1993) were used for the simulation. A comparison of the gasoline yield from the model presented in this work, with the gasoline yield from 3-D, twophase flow, heat transfer and reaction model of Theologos and Markatos (1993) is made. Only two parameters k0 , and 2 were changed appreciably. The final values of the parameters were k0 = 0.0012, E0 = 1540,  = 0.4, and 2 = 13.0. Comparison of results predicted by present approach and those reported by Theologos and Markatos (1993) is given in Fig. 14. Although the model of Theologos and Markatos (1993) is mathematically complex, its yield prediction profiles are similar to the models that use constant values for the gas and solid velocities and four to six lumps. The present model overpredicts the product yields in the first few meters of riser height. This may be attributed to the fact that in this region the heat and mass transfer resistances are not negligible. 80 Conversion or yield (wt%) 70 60 50 40 Model gasoline 30 Plant gasoline Model conversion 20 Plant conversion 10 0 0 5 10 15 20 25 30 4521 35 4. Conclusion Riser height (m) Fig. 13. Case study 2, comparison with the plant data reported by Derouin et al. (1997). 100 Theologos and Markatos model prediction Present model prediction Conversion or yields (wt%) 80 Conversion 60 Gasoline 40 Gas+coke 20 0 0 10 20 30 40 50 60 Riser height (m) Fig. 14. Case study 3, comparison with the simulator data reported by Theologos and Markatos (1993). The results of simulation with k0 = 0.045, E0 = 1540,  = 0.43, and 2 = 17.0, are presented in Fig. 13. Model predictions for the gasoline yield along the riser height matches satisfactorily with the plant data. 3.3. Case study 3 The objective of this case study is to compare the results from the present work with other models based on rigor- A new technique for modeling the FCC riser has been developed. The model incorporated a more realistic kinetic scheme for the cracking reactions, and a new correlation to evaluate Arrhenius type reaction rate constants. The rate constant parameters can easily be obtained for each combination of feed and catalyst. Although there is significant variation in the yield pattern of different case studies, activation energy parameter (E0 ) remained same for all cases. However, to account for different characteristics of the feedstock and catalyst, only frequency factor parameter (k0 ) and feed coking tendency parameter 2 were required to be adjusted to compare the yield patterns of different case studies with the model results. The proposed model is capable of predicting overall conversion, products yields, temperature, and catalyst activity along the riser height. The model results are in close agreement with the industrial data reported in the literature and the data predicted by other simulators. The predictions of the FCC riser reactor model are dependent on the values of cracking reactions’ rate constants, which can easily be obtained with the help of proposed kinetic model for different characteristics of the feedstock, type of catalyst, activity of catalyst, and operating parameters. Therefore, it seems to be more appropriate to use these rate constant parameters obtained for a pair of feedstock and catalyst in place of using the kinetic constants from the literature which are obtained for a different combination of feedstock and catalyst by regression analysis. Further, this detailed kinetic model can be easily used for the other advanced studies (such as control and optimization) of FCC modeling. Notation Ar Cc Ci,j cross-sectional area of riser, m2 coke concentration on catalyst surface, wt% concentration of ith pseudocomponent in j th volume element, kmol/m3 4522 Cpcat Cpi Cpmix,j Cpst Ei H coke H comb Hr k0,i ki,m,n Mcat Mcokej Mst MWi N NC p Pi,j ri,m,n R tj t Cat j Tin Tj u uc ug R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 specific heat of catalyst, kJ/kg specific heat of ith pseudocomponent, kJ/kg specific heat of the gas and solid mixture in j th volume element, kJ/kg specific heat of steam, kJ/kg activation energy for cracking of ith pseudocomponent giving kmax,i , kJ/kmol heat of combustion of coke, kJ/kg heat of combustion of pseudocomponents, kJ/kg heat of reaction, kJ/kmol frequency factor for cracking of ith pseudocomponent giving kmax,i m3 /(kg Cat s) rate constant for the cracking of ith pseudocomponent to produce mth and nth pseudocomponents, m3 /(kg Cat s) mass flow rate of catalyst, kg/s mass flow rate of coke at the outlet of j th volume element, kg/s mass flow rate of steam, kg/s molecular weight of ith component, kg/kmol total number of components, including pure components and pseudocomponents total number of hypothetical volume elements in the riser pressure, atm molar flow rate of ith component,PCi , through j th volume element, kmol/s rate of disappearance of ith pseudocomponent giving mth, and nth pseudocomponents, kmol/s gas constant, atm m3 /(kmol K) residence time of gas phase in j th volume element, s residence time of catalyst in j th volume element, s temperature of reaction mixture in the riser inlet, K temperature of reaction mixture leaving j th volume element, K superficial gas velocity, m/s cluster velocity, m/s actual gas velocity, m/s Subscripts i, m, n j ith, mth, and nth component j th volume element in the riser starting from the bottom Greek letters i,m,n g,j   mass of coke formed when 1 kmol of pseudocomponent PCi cracks to give 1 kmol each of PCm and PCn , kg coke/kmol PCi volume fraction of gas in j th volume element exponent of molecular weight for frequency factor exponent of molecular weight for activation energy cat coke g,j 1 , 2 j density of catalyst, kg/m3 density of coke, kg/m3 density of gas phase in j thvolume element tunable parameters catalyst activity coefficient Appendix A. Petroleum fractions are mixtures of innumerable components which are difficult to be identified individually. However, Watson characterization factor can be treated as an indicator of the composition of various groups of compounds (such as paraffin, olefin, naphthene, aromatic, etc.) present in the petroleum fraction. Watson and Nelson (1933) made a remarkable ob1/3 servation that the factor KW (=Tb /sg), known as Watson characterization factor, is closer to 12 for paraffins and olefins, approximately 10 for aromatics, and between 11 and 12 for naphthenes when the normal boiling point of the component, Tb , is in Rankin and sg is the specific gravity at 60◦ /60 ◦ F. The characterization factor of the mixture of hydrocarbons is given by KW = MeABP1/3 /sg, where MeABP is the mean average boiling point of the mixture (API Data Book, Chapter 2, Characterization of Hydrocarbons, 1976). Using this, Miquel and Castells (1993) proposed a method along with a computer program (Miquel and Castells, 1994) that can represent an oil fraction by an equivalent mixture of small number of hypothetical components or pseudocomponents. To use this approach, atmospheric TBP distillation curve and the entire fraction density is required. This method assumes that if the difference in final boiling point (FBP) and initial boiling point (IBP) of a petroleum oil is not too high (i.e., < 300 K) then the Watson characterization factor of any narrow-boiling fraction (boiling range between 15 and 25 K) of this oil remains equal to that of original petroleum oil. In the present case, due to unavailability of TBP curve for the FCC feed, SD curve reported by Pekediz et al. (1997) was used. The SD curve was first converted to ASTM-D86 curve and then to TBP curve by the correlation proposed by Daubert (1994). The two-step conversion of SD data to TBP data is given in Table A1. To generate pseudocomponents, the TBP curve of the feed was divided into 12 parts, out of which four were of 5 vol% each and eight of 10 vol% each (shown as vertical bars in Fig. A1 ). These vertical bars represent 12 pseudocomponents of the feed. The boiling point of each individual pseudocomponent was determined by area-averaging of the TBP curve (clearly visible in Fig. A1). Considering constant Watson characterization factor, specific gravity of each pseudocomponent was determined by the equation 1/3 T sg = 1.21644 b KW (where Tb is in K). (A.1) Molecular weights of these pseudocomponents were then calculated by the following equation proposed by Edmister and Lee (1984), which requires knowledge of boiling point and the R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 4523 Table A1 Distillation data of hydrocarbon feed (Pekediz et al., 1997) Vol. % distilled (wt%) SD (K) ASTM-D86a (K) TBPa (K) IBP 10 30 50 70 90 FBP 532 587 621 650 683 730 800 585.5 615.5 632.8 652.3 680.5 721.4 756.6 558.9 604.3 636.4 665.3 700.6 743.9 808.1 Density of feed (at 15 ◦ C) = 929.20 kg/m3 . a Estimated using correlation proposed by Daubert (1994). The lightest seven components were taken as pure components which are the major constituent of gases. Initially, volume fraction of all these seven pure components and 31 pseudocomponents, which are not present in feed, were taken zero. Thus total 50 components (seven pure components and 43 pseudocomponents) were considered in the present approach for the simulation of FCC riser reactor. After determining normal boiling point, specific gravity, and molecular weight of all pseudocomponents, heat capacities were determined using the correlations of Kesler and Lee (1976). Pseudocomponents thus generated are listed in Table A2. Predicted concentrations of pseudocomponents in the product stream are given in Fig. A2. Also, various product streams, viz., gas, gasoline, LCO, and residue are marked on the basis of boiling points. Appendix B. Fig. A1. Pseudocomponents generated from feed TBP. specific gravity of a hydrocarbon fraction: MW = 204.38 · e(0.00218·Tb ) · e(−3.07·sg) · Tb0.118 · sg1.88 . (A.2) Having known values of the volume fraction, boiling point, specific gravity, and molecular weight, each individual bars of Fig. A1 can be treated as a pure component (of course, hypothetical pure component or pseudocomponent). To make use of the Eqs. (A.1) and (A.2), an iterative method has to be adopted as the value of Watson characterization factor is not known beforehand. Miquel and Castells (1993, 1994) have explained this iterative approach in detail. After breaking the FCCU feed into 12 pseudocomponents, and determining the exact value of Watson characterization factor, properties of other 31 pseudocomponents were also determined by using Eqs. (A.1) and (A.2). Boiling points of these 31 pseudocomponents were taken at equal intervals between the boiling point of n–pentane and the boiling point of first pseudocomponent of FCCU feed (570 K in the present case). The reaction rate constant, k, is normally determined by Arrhenius equation in terms of frequency factor, k0 , and energy of activation E. In most of the cases of FCC kinetic modeling, these parameters are determined empirically, using experimental data. In the present case, however, all data available in literature relevant to calculate the rate constants are for lumped reaction mechanism which, in fact are reaction rate data for conversion of one lump to other and not for cracking of one lump giving two other lumps. Due to lack of experimental data, a purely hypothetical correlation for predicting Arrhenius type rate constant is being used. The proposed scheme can be perfected in future after performing more and more experimental work (may be in different laboratories). Hypothetical scheme: It is well observed fact that almost any physical, thermodynamic, or transport properties of hydrocarbons of a particular group have similar behavior, and properties of these hydrocarbons can be well correlated empirically (Daubert, 1998). In the present case, we are dealing with pseudocomponents, which do not fall under a particular group of hydrocarbon, but are mixtures of large number of hydrocarbons of almost equal boiling point but widely different properties. In fact, pseudocomponents are neither paraffin, nor olefin, and not aromatic either. However, the average characteristic 4524 R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 Table A2 Properties of pseudocomponents Component ID Component name Boiling point (K) Molecular weight H comb (kJ/kmol)a PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16 PC17 PC18 PC19 PC20 PC21 PC22 PC23 PC24 PC25 PC26 PC27 PC28 PC29 PC30 PC31 PC32 PC33 PC34 PC35 PC36 PC37 PC38 PC39 PC40 PC41 PC42 PC43 PC44 PC45 PC46 PC47 PC48 PC49 PC50 Methane Ethane Propane Butene Butane Pentene Pentane Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent Pseudocomponent 111.65 184.50 231.09 266.90 272.65 303.11 309.21 317.37 325.54 333.70 341.86 350.03 358.19 366.35 374.52 382.68 390.84 399.01 407.17 415.33 423.50 431.66 439.83 447.99 456.15 464.32 472.48 480.64 488.81 496.97 505.13 513.30 521.46 529.62 537.79 545.95 554.11 562.28 570.44 593.12 612.48 628.50 643.76 658.25 674.30 691.91 711.55 733.25 760.13 792.20 16.043 30.070 44.097 56.108 58.124 70.135 72.150 88.563 91.443 94.402 97.443 100.569 103.781 107.083 110.478 113.969 117.558 121.249 125.045 128.948 132.964 137.094 141.343 145.713 150.210 154.835 159.595 164.491 169.530 174.714 180.049 185.538 191.187 197.000 202.982 209.138 215.474 221.995 228.706 248.399 266.496 282.446 298.497 314.565 333.352 355.231 381.318 412.311 454.185 509.673 62764.79 58622.94 51983.86 50464.83 49960.31 49073.02 48952.84 47318.37 47226.62 47137.89 47052.01 46968.81 46888.17 46809.94 46734.00 46660.24 46588.54 46517.73 46445.89 46373.59 46300.89 46227.80 46154.34 46080.57 46006.47 45932.07 45857.39 45782.42 45707.18 45631.67 45555.89 45479.85 45403.53 45326.93 45250.05 45172.88 45095.39 45017.59 44939.46 44720.49 44531.20 44353.16 44187.51 44034.96 43871.15 43697.31 43510.14 43311.36 43075.64 42808.56 PC1 to PC7 are pure components constituting gas; PC8 to PC29 constitute gasoline fraction; PC30 to PC41 constitute light cycle oil fraction; PC42 to PC50 constitute residual fraction whereas feed contains PC39 to PC50 . a Heat of combustion values are calculated by using Eqs. (12), (13), and (14). of these mixtures of hydrocarbons with almost equal boiling point is characterized by Watson characterization factor KW (as discussed in Appendix A). Therefore it can safely be assumed that the over all cracking behavior of all pseudocomponents should follow similar trend, as all the pseudocomponents are generated with exactly the same value of KW . The Watson char- acterization factor is an indicator of an average characteristic in terms of paraffinicity, as well as aromaticity of the hydrocarbon mixtures (pseudocomponents). Therefore, these pseudocomponents are treated just as pure components with specific characteristics, represented by KW , which helps in determining all their physico–chemical properties. R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 4525 800 Residue Boiling Temperature (K) 700 Feed 600 LC O 500 Fig. B1. Schematic representation of cracking of a pseudocomponent: (a) cracking from middle of the molecule (without coke formation), (b) cracking from side of the molecule (without coke formation), (c) cracking from middle of the molecule (along with coke formation), (d) cracking from side of the molecule (along with coke formation). Gasoline 400 300 Gas 200 100 0 20 40 60 80 100 Liquid Volume % (in Product) Fig. A2. Mass fraction of pseudocomponents in the product. In the absence of experimental data, the rate constant of a particular reaction can be considered as the probability of reaction to take place. A higher probability of the reaction will correspond to higher cracking rate and hence a higher rate constant. According to the proposed reaction scheme (Eq. (1)) a pseudocomponent PCi , cracks to give two other pseudocomponents PCm , and PCn and some amount of coke () is also formed. For the cracking of one specific PCi (i.e., for one fixed value of i between 1 to N) there are a large number of parallel ‘feasible’ reactions taking place through which two components PCm and PCn are formed such that m and n can lie between 1 to i only (Fig. 1). However, the rate constants of all these feasible reactions (with different values of i, m, and n) must be different from one another. There are only three possible ways through which variation in the magnitude of cracking rate constant of ith pseudocomponent with changing molecular weight of PCm and PCn can occur, (i) the rate constant of a cracking reaction is maximum when molecular weights of PCm and PCn are almost equal (Fig. B1(a)), (ii) the rate constant is highest when molecular weights of PCm and PCn are widely apart (Fig. B1(b)), and (iii) the rate constant is almost constant for all values of m and n. These three possibilities can be expressed in terms of probability distribution as, (i) probability of cracking of a pseudocomponent from its middle is highest, (ii) cracking from its sides is more probable than from the middle, and (iii) cracking from anywhere is equally probable. It was later observed that cracking reactions follow the last case in which rate constant is almost equal for all possible combinations of cracking of a particular pseudocomponent. However, for the pedagogical point of view, following discussion explains how the present model was developed. First of all we consider the first case, i.e., when cracking from middle of the molecule has highest probability. Pitault et al. (1994) has also supported this assumption. Here it should be noted that ‘middle of the molecule’ means the molecular weight of PCm and PCn are equal, since there is no differentiation between ring-chain and straight chain molecules of pseudocomponents. Thus in this first case, cracking of pseudocomponents from the middle has highest probability, therefore, numerical value of the rate constant should also be highest when molecular weights of the two product pseudocomponents are equal. It implies that for the cracking of ith pseudocomponent (PCi ) giving two other pseudocomponents PCm and PCn, the rate constant (ki,m,n ) is maximum (kmax,i ) when molecular weight of mth and nth components are equal. In the present scheme, however, this is possible only when m=n since there are no two pseudocomponents having equal molecular weight. For cases when m = n (Fig. B1(b)), a function f (x) is defined to predict ki,m,n in such a way that the function value approaches kmax,i when x (=MWm − MWn ) tends to zero, and f (x) is less than kmax,i for all |x| > 0. This function could be any even function with a maxima at x = 0. However, in the present case, it is assumed that the probability of cracking of a pseudocomponent may follow a normal distribution. Therefore the function f (x), the normal distribution function in standard form, becomes 1 2 2 f (x) = √ e−x /2 . 2 (B.1) At x = 0, this function has maximum value which corresponds to the rate constant of the most probable cracking reaction, kmax,i . Therefore 1 f (0) = √ = kmax,i . 2 Hence 2 2 = 1 (kmax,i )2 (B.2) . (B.3) Substituting Eqs. (B.2) and (B.3) in (B.1) we get the normal distribution function f (x) = kmax,i e−x 2 (kmax,i )2 . (B.4) Here it should be noted that x is the difference in the molecular weights of the two product components and it is possible that the value of x is same for different combination of PCm and PCn . Another parameter in Eq. (1) is the amount of coke formed during cracking reaction. The point to ponder is, for a fixed value of x whether the cracking reaction rate shall increase, decrease, or remain unchanged when coke formation increases? 4526 R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 Table B1 Kinetic data for the cracking reactions reported by Arbel et al. (1995) Cracking reaction HFO to LFO HFO to gasoline HFO to coke LFO to gasoline LFO to coke Gasoline to coke Activation energy (kJ/mol) Frequency factor (h−1 ) Rate at 538 ◦ C (h−1 ) Molecular weight of cracking lump 60.7086 23.0274 73.269 23.0274 73.269 41.868 1.422 × 107 1760.4 3380.4 712.8 2707.2 356.4 172.8 380 380 380 255 255 120 1.026 × 105 3.704 × 107 8.215 × 104 1.852 × 107 8.555 × 104 HFO: heavy fuel oil; LFO: Light fuel oil. Table B2 Kinetic data for the cracking reactions reported by Lee et al. (1989) Cracking reaction Activation energy (kJ/mol) Frequency factor (h−1 ) Rate at 538 ◦ C (h−1 ) Molecular weight of cracking lump Gas oil to gasoline Gas oil to gas Gas oil to coke Gasoline to gas Gasoline to coke 68.2495 89.2164 64.5750 52.7184 115.4580 7.978 × 105 4.549 × 106 3.765 × 104 3.255 × 103 7.957 × 101 39.364 9.749 6.012 2.470 1.364 380 380 380 120 120 The formation of coke can be depicted as in Fig. B1(c) and (d), and the mass of coke formed (in kg), when one kmol of reactant (PCi ) cracks, can be given by i,m,n = MWi − (MWm + MWn ). (B.5) It appears from Table B.1 and Table B.2 that the rate constant for the conversion of a lump to coke is significantly lower than the rate constants for the conversion of same lump to other hydrocarbon lumps. It indicates that the formation of coke is not favored during the cracking reaction. From this fact an inference can be drawn that the probability (or the rate constant) of a cracking reaction decreases when coke formation (i,m,n ) is increased. Since we are considering cracking of pseudocomponents, which itself are mixtures of large number of actual compounds, it seems that the fractional decrease in rate constant with increasing coke formation should be a continuous function say g(), where  is the amount of coke formation (i,m,n ). Thus the general correlation for the rate constant can be expressed as the product of these two functions as ki,m,n = f (x)g(). (B.6) e−i,m,n − e−MWi , 1 − e−MWi such that, g(i,m,n ) = 1 at i,m,n = 0 (i.e., no coke formation) and g(i,m,n ) = 0 at i,m,n = MWi (i.e., all hydrocarbon mass is converted to coke). 2 2 ki,m,n = kmax,i e− (MWm −MWn ) (kmax,i ) e−(MWi −(MWm +MWn )) − e−MWi × . 1 − e−MWi (B.8) A schematic three dimensional surface of ki,m,n as functions of coke formation (i,m,n ) and x (the difference between MWm and MWn ) for the cracking of a typical pseudocomponent is shown in Fig. B2. It can be observed that at different values of coke formation, the grid-line parallel to x-axis follows the normal distribution function of different curvature, and the function has a maxima at x = 0. Up to this point we have considered only hypothetical correlations, without any experimental verification. Eq. (B.8) is quite rigid to accommodate experimental data as the surface shown in Fig. B2 can change only by changing kmax,i . Therefore, at this stage, we introduced two tunable parameters 1 and 2 which can be used to adjust the span (variance) of the normal distribution function and the curvature of the decay function, respectively, such that 2 In the present work, we hypothesize that the function g() is an exponential decay function of the form g(i,m,n ) = Combining Eqs. (B.4), (B.6) and (B.7), we get (B.7) ki,m,n = [kmax,i · e−1 [(MWm −MWn )·kmax,i ] ] × e−(MWi −(MWm +MWn ))/2 − e−MWi 1.0 − e−MWi . (B.9) The parameter 2 is an indicator of coking tendency and its value will depend on the nature of the feed. Parameter 1 is an indicator of cracking tendency of pseudocomponents from middle or from sides. By increasing 2 decay of exponential function g() is reduced. The thick solid line on the grid surface (Fig. B2) is elevated vertically upward to thick dashed line when 2 is increased. Hence by increasing 2 reaction rate constant can be increased even with higher coke formation. On the other hand, R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 4527 the following relationship with molecular weight of cracking lump:  k0,i = k0 MWi , (B.11) Ei = E0 MWi . (B.12) Therefore, the parameter kmax,i used in Eq. (B.9) can be expressed in terms of the frequency factor (k0,i ) and the energy of activation (Ei ) with four rate constant parameters k0 , , E0 , and . After the sensitivity analysis, as discussed in the main body of the text, it was observed that the parameter  is insignificant; therefore it can be eliminated from the list of six parameters. References Fig. B2. Schematic representation of the probability distribution function f (x) (on the vertical plane) and the reaction rate constant ki,m,n as a function of coke formation and the difference in molecular weights of the products. by decreasing value of 1 the function f (x) (hence the entire surface of Fig. B2) can be made more flat. If the value of 1 is reduced to zero, the intersection of the surface with ki,m,n –x plane becomes a straight line. By further decreasing the value of 1 the surface becomes concave upward. In Fig. B2, the intersection of such surface with ki,m,n –x plane is shown by a thin dashed line. Thus if 1 is positive, the cracking behavior of pseudocomponents is as explained for case (i), whereas, negative value of 1 leads to case(ii), i.e., the rate constant is highest when molecular weights of PCm and PCn are widely apart. For the third case, when rate constant is independent of molecular weights MWm and MWn , 1 is zero. Thus by introducing these two tuning parameters, 1 and 2 , Eq. (B.9) can be used to predict rate constant for all possible schemes. With change in temperature, it is expected that there should not be any change in the probability distribution for the cracking of a particular pseudocomponent i.e., shape of the normal distribution curve remains unchanged at all temperatures. At the same time, the reaction rate constant ki,m,n has to follow the Arrhenius equation. This is possible only if kmax,i follows the equation kmax,i = k0,i e−Ei /RT , (B.10) where the energy of activation (Ei ) and frequency factor (k0,i ) are the characteristics of ith pseudocomponent (PCi ), therefore numerical values of these parameters should be different for different pseudocomponents. In other words, Ei and k0,i should be a function of some characteristic property of ith pseudocomponent. One such property could be the molecular weight (MWi ). To predict such a relation, data of Tables A1 and A2 can be used. Although, these Tables represent data for conversion reaction (not a cracking reaction), even then, a wild guess can be made from data given in these tables. Overall, it appears that parameters Ei and k0,i have an increasing trend with increasing molecular weight of reactant lump. These trends may have Ali, H., Rohani, S., Corriou, J.P., 1997. Modeling and control of a risertype fluid catalytic cracking (FCC) unit. Transactions of the Institution of Chemical Engineers 75, 410–412. Arastoopour, H., Gidaspow, D., 1979. Vertical pneumatic conveying using four hydrodynamic models. Industrial and Engineering Chemistry Research 18, 123–130. Arbel, A., Huang, Z., Rinard, I.H., Shinnar, R., Sapre, A.V., 1995. Dynamics and control of fluidized catalytic crackers. Modeling of the current generation FCC’s. Industrial and Engineering Chemistry Research 34, 1228–1243. Aris, R., 1989. Reactions in continuous mixtures. A.I.Ch.E. Journal 35, 539–548. Austin, G.T., 1984. Shreve’s Chemical Process Industries. fifth ed. McGrawHill, Singapore. Blasetti, A., de Lasa, H., 1997. FCC riser unit operated in the heat-transfer mode: kinetic modeling. Industrial and Engineering Chemistry Research 36, 3223–3229. Corella, J., Frances, E., 1991. Analysis of the riser reactor of a fluid catalytic cracking unit: model based on kinetics of cracking and deactivation from laboratory tests. In: Occelli, M.L. (Ed.), Fluid Catalytic Cracking—II: Concepts in Catalyst Design, ACS Symposium Series, vol. 452. American Chemical Society, Washington, pp. 165–182. Daubert, T.E., 1994. Petroleum fraction distillation interconversions. Hydrocarbon Processing 73, 75–78. Daubert, T.E., 1998. Evaluated equation forms for correlating thermodynamic and transport properties with temperature. Industrial and Engineering Chemistry Research 37, 3260–3267. Derouin, C., Nevicato, D., Forissier, M., Wild, G., Bernard, J.R., 1997. Hydrodynamics of riser units and their impact on FCC operation. Industrial and Engineering Chemistry Research 36, 4504–4515. Edmister, W.C., Lee, B.I., 1984. Applied Hydrocarbon Thermodynamics, vol. 1, second ed. Gulf Publishing, Houston. Feng, W., Vynckier, E., Froment, G.F., 1993. Single event kinetics of catalytic cracking. Industrial and Engineering Chemistry Research 32, 2997–3005. Fligner, M., Schipper, P.H., Sapre, A.V., Krambeck, F.J., 1994. Two phase cluster model in riser reactors: impact of radial density distribution on yields. Chemical Engineering Science 49, 5813–5818. Gupta, A., Subba Rao, D., 2001. Model for the performance of a fluid catalytic cracking (FCC) riser reactor: Effect of feed atomization. Chemical Engineering Science 56, 4489–4503. Gupta, R., Kumar, V., Srivastava, V.K., 2005. Modeling and simulation of fluid catalytic cracking unit. Reviews in Chemical Engineering 21 (2), 95–131. Horio, M., Kuroki, H., 1994. Three dimensional flow visualization of dilutely dispersed solids in bubbling and circulating fluidized beds. Chemical Engineering Science 49, 2413–2421. Kesler, M.G., Lee, B.I., 1976. Improve prediction of enthalpy of fractions. Hydrocarbon Processing 55, 153–158. 4528 R.K. Gupta et al. / Chemical Engineering Science 62 (2007) 4510 – 4528 Kumar, V., Upadhyay, S.N., 2000. Computer simulation of membrane processes: ultrafiltration and dialysis units. Computers and Chemical Engineering 23, 1713–1724. Lee, L., Chen, Y., Huang, T., Pan, W., 1989. Four lump kinetic model for fluid catalytic cracking process. Canadian Journal of Chemical Engineering 67, 615–619. Liguras, D.K., Allen, D.T., 1989a. Structural models for catalytic cracking. 1. Model compound reactions. Industrial and Engineering Chemistry Research 28, 665–673. Liguras, D.K., Allen, D.T., 1989b. Structural models for catalytic cracking. 2. Reactions of simulated oil mixtures. Industrial and Engineering Chemistry Research 28, 674–683. Markatos, N.C., Shinghal, A.K., 1982. Numerical analysis of one-dimensional, two phase flow in a vertical cylindrical passage. Advance Engineering Software 4, 99. Martin, M.P., Derouin, P.T., Forissier, M., Wild, G., Bernard, J.R., 1992. Catalytic cracking in riser reactors. Core-annulus and elbow effects. Chemical Engineering Science 47, 2319–2324. Maxwell, J.B., 1968. Data Book on Hydrocarbons. Krieger, New York. Miquel, J., Castells, F., 1993. Easy characterization of petroleum fractions. Hydrocarbon Processing 72, 101–105. Miquel, J., Castells, F., 1994. Easy characterization of petroleum fractions. Hydrocarbon Processing 73, 99–103. Pekediz, A., Kraemer, D., Blasetti, A., de Lasa, H., 1997. Heats of catalytic cracking. Determination in a riser simulator reactor. Industrial and Engineering Chemistry Research 36, 4516–4522. Pitault, I., Nevicato, D., Frossier, M., Bernard, J.R., 1994. Kinetic model based on a molecular description for catalytic cracking of vacuum gas oil. Chemical Engineering Science 49, 4249–4262. Pitault, I., Forissier, M., Bernard, J.R., 1995. Determination de constantes cinetiques du craquage catalytique par la modelisation du test de microactivite (MAT). The Canadian Journal of Chemical Engineering 73, 498–503. Pugsley, T.S., Berruti, F., 1996. A predictive hydrodynamic model for circulating fluidized bed risers. Powder Technology 89, 57–69. Quam, R.J., Jaffe, S.B., 1992. Structure-oriented lumping: describing the chemistry of complex hydrocarbon mixtures. Industrial and Engineering Chemistry Research 31, 2483–2497. Saltelli, A., 2000. In: Chan, K., Scott, E.M. (Eds.), Sensitivity Analysis. Wiley, New York. Stangeland, B.E., 1974. A kinetic model for the prediction of hydrocracker yields. Industrial and Engineering Chemistry Process Design & Development 13 (1), 71–76. Subba Rao, D., 1986. Clusters and lean phase behavior. Powder Technology 46, 101–107. Theologos, K.N., Markatos, N.C., 1993. Advanced modeling of fluid catalytic cracking riser-type reactors. A.I.Ch.E. Journal 39, 1007–1017. Theologos, K.N., Lygeros, A.I., Markatos, N.C., 1999. Feedstock atomization effects on FCC riser reactors selectivity. Chemical Engineering Science 54, 5617–5625. Tsuo, Y.P., Gidaspow, D., 1990. Computation of flow patterns in circulating fluidized beds. A.I.Ch.E. Journal 36, 885–896. Wang, L., Yang, B., Wang, Z., 2005. Lumps and kinetics for the secondary reactions in catalytically cracked gasoline. Chemical Engineering Journal 109, 1–9. Watson, K.M., Nelson, E.F., 1933. Symposium on physical properties of hydrocarbon mixtures—improved methods for approximating critical and thermal properties of petroleum fractions. Industrial and Engineering Chemistry 25, 880–887. Wilson, J.W., 1997. Fluid Catalytic Cracking Technology and Operation. PennWell Publishing Co, Tulsa, Oklahoma.