Abstract: The aim of this study is to obtain a steady state model that can simulate an industrial Residue Fluid Catalytic Cracking (RFCC) unit. This unit is one of the technologies for producing more gasoline from residue. The yield of gasoline in RFCC strongly depends on certain process variables. In this work, an RFCC model is developed by combining the material and energy balance equations with a 7-lump kinetic model and a modified two-dimensional hydrodynamic model. Simulation has been performed based on the data from an operating unit at Khartoum Refinery Company (KRC) and the results are reported. Optimum values of process variables for a required cracking efficiency, such as space velocity, catalyst to oil ratio and catalyst circulation rate, are also reported.
INTRODUCTION
Catalytic cracking is the most important and widely used refinery process for converting heavy oils into more valuable gasoline and lighter products. The catalytic cracking processes produce carbon (coke) which remains on the catalyst particle and rapidly lowers its activity. To maintain the catalyst activity at a useful level, it is necessary to regenerate the catalyst by burning off this coke with air. The cracking reaction is endothermic and the regeneration reaction is exothermic. Recent units are designed to use the regeneration heat to supply that needed for the reaction and to heat the feed up to reaction temperature (Gary and Handwerk, 1984). This is accomplished by fluidization, where the fluidized catalyst carries the regeneration heat to the reactor.
Several studies have been conducted for FCC modeling, simulation, optimization and control, but very few of them consider the RFCC (Xu et al., 2006; Gauthier et al., 2000; Fernandes et al., 2007). Many of the models in these studies describe the riser reactor with one dimensional mass, energy and chemical species balances. These models are normally simple to formulate and to solve. More detailed models, are those that consider the riser reactor in two or three dimensions. These models are based on phenomenological concepts and use a simultaneous solution of the conservation equations of mass, momentum, energy and species for both the gas and particulate phases (Gao et al., 1999; Mathiesen et al., 1999). Therefore, Souzaa et al. (2003), combined a 2-D fluid flow field with a 6- lumps kinetic model and used two energy equations (catalyst and gas oil) to simulate the gas oil cracking process inside the riser reactor. Das et al. (2003), performed the three-dimensional simulation of an industrial-scale fluid catalytic cracking riser reactor using a novel density-based solution algorithm. The particle-level fluctuations are modeled in the framework of the kinetic theory of granular flow. The reactor model includes separate continuity equations for the components in the bulk gas and inside the solid phase.
Berry et al. (2004), modified the two-dimensional hydrodynamic model to make it predictive by incorporating the slip factor for the calculation of the cross-sectionally averaged voidage. The model has been coupled with the four-lump kinetic model of Gianetto et al. (1994), to predict the effect of operating conditions on profiles of conversion, yield, temperature and pressure in the riser. With regard to reaction and kinetics, Xu et al. (2006) proposed a 7-lump kinetic model to describe residual oil catalytic cracking, in which products especially coke were lumped separately for accurate prediction. Because in recent studies, kinetics was developed accounting for coke formation leading to catalyst deactivation. The reactor block is modeled as a combination of an ideal Plug Flow Reactor (PFR) and a Continuously Stirred Tank Reactor (CSTR).
In the present work, the 7-lump kinetic model is utilized to model the cracking reactions in the riser to determine the products yield with the required conversion and other process variables. These are the base information required for the dynamic model. The result is compared with the plant data from Khartoum Refinery Company (KRC) in Sudan and with Xu et al. (2006) results. Also the hydrodynamic model of Malcus and Pugsley (2002), which is modified by Berry et al. (2004), is solved in this work for the riser operating conditions from KRC. These models have to be combined with the steady state material and energy balances of RFCC unit with residual feedstock from Sudanese crude distilling unit.
Process description: To process the heavy, viscous residual feedstock, that produce relatively high amounts of coke, the design incorporates an advanced feed injection system, a unique regeneration strategy and a catalyst transfer system that produces extremely stable catalyst circulation. Recent technology advances have been made in the areas of riser termination, reactant vapor quench and Mix Temperature Control (MTC) (Meyers and Hunt, 2003). The goal of RFCC unit is to maximize gasoline production. Figure 1 shows the schematic diagram of the reaction - regeneration section (area of our study in this unit), for which the mathematical model is being solved.
The pre-heated feed is injected into the riser bottom, atomized into small droplets and contacted with the high temperature regenerated catalyst from the second regenerator. It vaporizes immediately, being cracked into products. Cycle oil is injected into four nozzles above the fresh feed nozzles in order to optimize product distribution by raising the catalyst to oil ratio. In the middle of the riser, nozzles for slurry recycle and naphtha are provided.
At the riser out let, the product vapor is separated from the catalyst by the Ramshorn disengager and four single stage cyclones. The stripped catalyst passes through the spent catalyst slide valve and enters into the catalyst distributor of the first regenerator, where partial combustion of the coke take place at the dense bed. The semi-regenerated catalyst from the first regenerator passes through the standpipe and slide valve and then lifted by compressed air to the second regenerator. The regenerated catalyst from the second regenerator flows into the hopper and then into the regenerated catalyst standpipe to the bottom of the riser, thus forming a continuous catalyst circulation loop. Catalyst cooler can be installed by the first regenerator to keep balance of the reactor regenerator system and to improve the operating flexibility.
MATHEMATICAL MODELS
Reaction kinetics: For kinetic modeling, the complexity of charge stocks make it extremely infeasible to characterize and describe at a molecular level, so attempts have been made to lump large numbers of chemical compounds into a kinetic species to describe the complex reactions in the process (Xu et al., 2006; Wei and Kuo, 1969). The simple models that describe the cracking kinetics with 3, 4 or 5 lumps have the advantage that just a few kinetic constants must be estimated for each feedstock, but depending on the simplicity of the model, the key FCC products cannot be predicted separately. More sophisticated models, normally with more than 10 lumps, have basically two advantages; a single group of estimated kinetic constants can be used for various feedstock and all the most important FCC products can be predicted separately.
| Fig. 1: | Process diagram of the residue fluid catalytic cracking unit RFCC, reaction section, at Khartoum Refinery Company (KRC) |
| Fig. 2: | The 7-lump reactions scheme for RFCC (2) |
The disadvantages of these models are that a large number of kinetic constants must be estimated and as each lump represents a differential equation in the mathematical model, the complexity of the numerical solution may increase exponentially (Souzaa et al., 2003).
Recently kinetics was developed accounting for coke formation leading to catalyst deactivation. Therefore, in the present work the 7-lump kinetic model, proposed by Xu et al. (2006), has been chosen. In this model coke, dry gas and LPG were lumped separately. The probability of the cracking reaction from gasoline to coke and gases to coke can be neglected since the kinetic constants for these reactions are many orders of magnitude smaller than that for the others, which simplifies the kinetic parameters determination. So, the reaction scheme for RFCC is shown in Fig. 2.
According to the principle of invariant response, lumps are formed by their boiling ranges having similar behavior. Also these lumps are not separated by chemical type with the parameters depending on feedstock and catalyst properties, which is the limitation of this model (Xu et al., 2006). The j lumps are Lump 1 for residue R, lump 2 for heavy fuel oil HFO, lump 3 for light fuel oil LFO, lump 4 for Gasoline G, lump 5 for liquefied petroleum gas LPG, lump 6 for dry gas and lump 7 for coke C.
An irreversible pseudo first order reaction was accepted for all reactions in this model. The reaction rate of a pseudospecies j, is assumed to be a function of the molar concentration of species j(ρaj), the mass density of catalyst relative to the gas volume (ρc/ε) and reaction constant (k'j ), as in Eq. 1.
| (1) |
The reaction constant k'j decayed with time as:
| (2) |
The heavy, inert aromatics adsorption function f(CArh) is described by:
| (3) |
Catalyst deactivation was accounted for by a hyperbolic function Ø(tC) as:
| (4) |
Nitrogen poisoning function f(N) is defined as:
| (5) |
With regard to high catalyst to oil ratio, the nitrogen poisoning deactivation can be neglected because of its insignificance. Taking the ideal gas assumption:
| (6) |
The average molecular weight of all lumps,
| (7) |
Therefore, the reaction rate in matrix notation can be expressed by:
| (8) |
where, the reaction rate vector is:
The composition vector consisting of j lumped species is a = [a1 a2 a3 a4 a5 a6 a7] and the rate constant matrix defined as in the appendix.
The reactor block is considered as a combination of two different reactors, according to its hydrodynamic regime, the riser is a near ideal plug-flow displacement of the catalyst and reaction mixture, while the main vessel (separator) is an ideal mixing CSTR. Therefore, the kinetic model is incorporated into an isothermal plug flow reactor model by the continuity equation, with interparticle diffusion and axial dispersion neglected Xu et al., 2006, as in Eq. 9:
| (9) |
For the steady state one dimensional equation:
| (10) |
Combining this with equation 8 to get the lumps concentration profiles with the axial dimensionless distance, Z = z/L, along the riser:
| (11) |
where, Swh = Gvε/ρcL is the true weight hourly space velocity.
The concentration profile obtained above assumes the continuity equation of each lump is well represented by the steady state equation only. However, the dynamic continuity equation must be used to get the actual concentration profile for each lump. Thus a dynamic, two dimensional continuity equation is proposed to re-optimize the kinetic parameters:
| (12) |
By this equation the concentration profile of all the lumps can be optimized to get the correct values of concentration at any point in the riser axially and radially.
Riser hydrodynamic model: Several different riser hydrodynamic models have been presented in the fluidization literature (Malcus and Pugsley, 2002; Godfroy et al., 1999). Many of these models are empirical or semi-empirical in nature and are based on experimental observations from cold-flow laboratory and pilot-scale risers. One of these is the recently developed model of Malcus and Pugsley (2002), which Berry et al. (2004) have used in their simulations. They have shown that the Boltzmann function of Eq. 13 leads to good predictions of radial profiles of the solids volume concentration at operating conditions typical of FCCU risers.
| (13) |
Berry et al. (2004) made the model of Malcus and Pugsley, predictive by using the slip factor correlation of Patience et al. (1992), in Eq. 14:
| (14) |
The average particle velocity in the riser, Vp, is defined by:
| (15) |
Combining Eq. 14 and 15 gives:
| (16) |
To complete their hydrodynamic model, Malcus and Pugsley, performed a multiple nonlinear regression on a large set of data collected from the open literature to obtain the following expression for predicting the mass flux at the centerline of the riser.
| (17) |
Combining this with the solids concentration at the riser centerline, obtained from Eq. 13, to get:
| (18) |
A second Boltzmann function is then used to describe the radial solids profile of the solids velocity:
| (19) |
Xo stands for the center of the Boltzmann function described as:
| (20) |
Where:
| (21) |
| (22) |
Berry et al. (2004) calculated the radial solids velocity profile using an iterative procedure. They set a value for the wall solids velocity (at r = R), then calculated the profile Vp,r from the above equations. They integrated this profile to get the average solids velocity and compared this value with the average solids velocity obtained from mass balance. If the mass balance was not satisfied, the value of k2 has to be updated and the calculation has to be performed again. They repeated this iterative procedure until the integrated value of the solids velocity was in agreement with the average solids velocity according to the mass balance.
However, this procedure did not take into account the axial velocity, Vp,z in calculating the average velocity. Whereas, in 2-dimensional modeling Vp,z must be included in estimating the average velocity. In order to accomplish this, momentum balance equations for both phases are proposed to get the other profiles Vp,z, vr ,vz.
Process variables, material and energy balances: In calculating the flow rate of streams in the reactor and two regenerators, the following definations and equations (Gary and Handwerk, 1984; Nelson, 1958 and Fei et al., 1997) were used:
| • | The catalytic cracking efficiency (η) depend on how much gasoline is yielded from the unit as in Eq. 23: |
| (23) |
| • | The catalyst activity = conversion/(100 - conversion) |
| • | The severity factor for highly wax-stock with characterization factor equal to 12.4 was correlated from Nelson (1958), with coefficient of determination 0.99666 as: |
| (24) |
Then some of the operating variables like space velocity and residence time can be calculated from:
| • | Space velocity = catalyst to oil ratio/severity factor |
| • | Catalyst residence time=1/[space velocity* catalyst/oil ratio] |
| • | Also delta coke is a very popular index and when increased can cause significant rise in regenerator temperature, ultimately reducing the amount of residual feed which can be processed. Therefore in order to maximize the amount of residual content in RFCC feed, a low-delta coke catalyst must be employed (must be less than 1.2), which is defined as: |
| (25) |
where, CSC is Carbon on spent catalyst and CRC is Carbon on regenerated catalyst. And as:
| (26) |
All these definitions and equations were solved simultaneously to calculate the flow rate of the streams that was needed in the steady state material and energy balance calculations where:
Sum of flow rates of input streams = Sum of flow rates of output steams
Reactor duty = Enthalpy of the input streams-Enthalpy of the output streams
RESULTS AND DISCUSSION
To solve the above mathematical models, a program in MATLAB was developed. The chosen 7-lump kinetic model was solved giving the concentration profiles of the j components (lumps), as shown in Fig. 3.
From this figure, it is observed that most of the oil cracking occurs in the feed injection zone, this means that oil fraction cracking is an instantaneous reaction.
| Fig. 3: | Axial lumps concentration profiles along riser |
| Fig. 4: | Radial solids velocity profile and mass flux profile in the riser |
The products yield increase rapidly in the first one-third of the riser height and then vary gently. The heavy fuel oil yield reaches a maximum at about 20% of the riser height and then declines due to overcracking reactions. The products yield at the end of the riser length Z = 1, is consistent with the design value for this yield. A very good agreement is shown with Xu et al. (2006), because if the two heavy products are combined together and treated as the unreacted products lump, the same trend is obtained.
Also in solving the hydrodynamic model, the radial solid velocity profile is predicted as Berry et al. (2004) have done by initially setting the wall solid velocity equal to the terminal settling velocity of a single particle. Then, from the simulator the actual wall solid velocity is fixed. So, the hydrodynamic model was solved giving the following radial solids velocity profile and solids mass flux profile, as shown in Fig. 4. It is assumed no axial changes appears in both velocity and flux, with the average mass flux = 598 kg m-2 sec-1 and the superficial gas velocity = 6.6 m sec-1.
There is a good agreement between these results and that of Berry et al. (2004), who used in their model prediction and validation, commercial data from three publications in the open literature for two different commercial risers.
The percent difference between in and out streams balances the reactor, first and second regenerators and gives the optimum values for the process variables. Changing input variables and saving the output gives the relation between process variables and some of the input parameters. It can be noticed from Fig. 5 and 6, if the task is lower residence time or higher space velocity in the reactor then the unit can easily be designed with higher cracking efficiency, but this will be at the cost of severity and % conversion, where both of them decrease with the increasing efficiency. The effect of efficiency on catalyst/oil ratio and delta coke is not big. Therefore, the best cracking efficiency ranges between 75-85%.
| Fig. 5: | Changes with cracking efficiency |
| Fig. 6: | Changes with cracking efficiency |
| Fig. 7: | Changes with gasoline yield |
Yield calculations depend on crude and residue specifications so changing gasoline yield will affect the energy balance slightly. Figures 7-9 shows how process variables change with the change in gasoline yield. Increasing gasoline yield will increase % conversion which is obviously true. But an increase in conversion does not necessarily mean an increase in gasoline yield because coke and gases yield increase at the expense of gasoline yield. This agrees with Gary and Handwerk, (1984) and Fei et al. (1997). Also increasing gasoline yield will increase catalyst/oil ratio and severity. This increase is sharp and not preferable above 0.6 gasoline yield. Besides, at this region the energy is not balanced any more. Delta coke and residence time increase with the increase in gasoline yield till an optimum point then decrease.
| Fig. 8: | Changes with gasoline yield |
| Fig. 9: | Changes with gasoline yield |
| Table 1: | Process variables |
The decrease becomes sharp after 0.55 gasoline yield and lower delta coke gives higher Carbon on Regenerated Catalyst (CRC) which leads to lower catalyst activity. Also, space velocity is too small at higher yield. This means that gasoline yield above 0.55 is dangerous.
Obtaining fixed and accurate values for the above process variables is very important in the design of the reactor and two stage regenerators, so with an objective of achieving 80% cracking efficiency the values obtained is in Table 1.
CONCLUSION
It has been observed that optimum cracking efficiency is 75-85% and gasoline yield up to 0.55. Designing with gasoline yield above this limit is not preferred. It is recommended that, as a future work, these steady state models have to be expanded to the overall dynamic model for the full riser simulation and optimization.
ACKNOWLEDGMENTS
Authors are grateful to Universiti Technologi PETRONAS (UTP) for providing all support to undertake this research work. Also authors would like to acknowledge the staff at Khartoum Refinery Company (KRC) in Sudan, for their resources.
NOMENCLATURE
| a | : | Composition vector consisting of j lumped species |
| aj | : | Concentration of lump j (moles j/g gas) |
| Carh | : | wt.% of aromatics in residual oil |
| D | : | Riser diameter (m) |
| Fr | : | Froude number = uo/(gD)0.5 |
| Frt | : | Froude number based on the terminal settling velocity = Ut/(gD)0.5 |
| Gv | : | The superficial mass flow rate of gas mixture (g/(cm2.h)) |
| g | : | Gravitational constant 9.81 (m sec-2) |
| Gs | : | Riser solids mass flux (kg m-2 sec-1) |
| Gs,cl | : | Solids mass flux at the centerline of the riser (kg m-2 sec-1) |
| L | : | Riser length (m) |
| k1 | : | Parameter in the hydrodynamic model defined by Eq. 21 |
| k2 | : | Parameter in the hydrodynamic model defined by Eq. 22 |
| kj | : | Invariant reaction rate constant |
| k'j | : | Reaction rate constant decayed with time |
| K | : | Rate constant matrix |
| Kh | : | The heavy aromatic ring adsorption coefficient, (wt. % CArh)-1, Kh = 0.128 (Xu et al., 2006) |
| Kh | : | Basic nitrogen poisoning adsorption coefficient (g basic N/g catalyst)-1 |
| R | : | Reaction rate in matrix notation |
| R | : | Riser radius (m) |
| r | : | Radial position in the riser (m) |
| rj | : | Kinetic expression of species j |
| t | : | Time (sec) |
| tc | : | Catalyst residencetime |
| uo | : | Riser superficial gas velocity (m sec-1) |
| Vp | : | Average catalyst velocity (m sec-1) |
| Vp,cl | : | Solid velocity at the centerline of the riser (m sec-1) |
| Vp,r | : | Radial solid velocity profile (m sec-1) |
| z | : | Axial position in the riser (m) |
| Z | : | Axial dimensionless distance along the riser |
Greek symbols
| εavg | : | Average voidage |
| εr | : | Voidage at radial position r |
| β | : | Catalyst deactivation constant = 162.15 (2) |
| γ | : | Catalyst deactivation constant = 0.76 (2) |
| n | : | Slip factor |
| ρ | : | Density of the catalyst (kg/m3) |
| ρs | : | Gas mixture density |
| η | : | Catalytic cracking efficiency |
APPENDIX: RATE CONSTANT MATRIX