SciELO - Scientific Electronic Library Online

vol.17 issue1Preparation of silica with controlled pore sizes for enzyme immobilizationOn the acidity and/or basicity of USY zeolites after basic and acid treatment author indexsubject indexarticles search
Home Pagealphabetic serial listing  

Services on Demand




Related links


Brazilian Journal of Chemical Engineering

Print version ISSN 0104-6632On-line version ISSN 1678-4383

Braz. J. Chem. Eng. vol.17 n.1 São Paulo Mar. 2000 

Modelling and simulation of the hydrocracking of heavy oil fractions


E.M. Matos and R. Guirardello
Faculdade de Engenharia Química, Departamento de Processos Químicos, Universidade Estadual de Campinas
(UNICAMP), C.P. 6066, CEP 13081 - 970, Phone (+ 55) (021-19) 788-3955,
Fax (+ 55) (021-19) 788-3965, Campinas - SP, Brazil,


(Received: April 9, 1999; Accepted: September 8, 1999 )



Abstract - This work presents a model to describe the behavior of the concentration of constituents of heavy fractions from petroleum during the hydrocracking process. An approximation based on pseudocomponents or lumps is adopted due to the complexity of the mixture. The system is modeled as an isothermal tubular reactor with an axial dispersion, where the hydrogen flows upward concurrently with the oil while the solid catalyst particles stay inside the reactor in an expanded bed regime. Simulations are carried out for different values of liquid superficial velocity, reactor length and degree of mixing (Peclet number).
Keywords: reactor, simulation, hydrocracking




Hydrocracking is the phenomenon where big petroleum molecules are broken into smaller molecules with a lower boiling point in the presence of hydrogen. In this process the number of valuable oil subproducts increases, due to their formation from the residues with low commercial values naturally present in petroleum. The composition of the cracked oils is highly dependent on the extent of hydrocracking, usually at high temperatures (350 to 600 ºC) and pressures (50 to 110 atm). The increase in the weight of the marketable fraction of the petroleum compensates the high costs of the process (Mohanty et al., 1990).

The increase in the marketable fraction of petroleum is of great economic importance because this valuable scanty resource is used to produce plastics, fertilizers, rubber, disinfectants, fuels and many other products. The description of the concentration profiles in this process may be useful for future work in the field of design, optimization and control of hydroprocessing reactors.

The reactor considered in this work is an isothermal tubular reactor with excess hydrogen in (gas phase) flowing concurrently upward with heavy oil (liquid phase) while the catalyst (solid phase) remains inside the reactor in an expanded bed regime (confined bed). The hydrocracking reactions can be mainly thermal (Mosby et al., 1986) or catalytic (Krishna and Saxena, 1989), depending on the catalyst being used and the operating conditions. In thermal hydrocracking the catalyst is used for other reactions, such as hydrodemetallation and removal of other heteroatoms.



The Kinetic Network

The hydrocracking reactions considered here are based on the kinetic network proposed by Krishna and Saxena (1989), which was established in terms of pseudocomponents, as shown in Figure 1. The aromatic lump is the pseudocomponent formed by the set of all molecules that have at least one benzenic ring in their structure. The naphthene lump is the pseudocomponent formed by the set of all molecules that have at least one cyclic ring in their structure that is not a benzenic ring. The paraffin lump is the pseudocomponent formed by the set of all molecules that do not have rings in their structure. The pseudocomponents are considered light if they are formed by fractions with boiling points not greater than some cut temperature, but otherwise they are considered heavy. The values for the constants of reaction ki for several cut temperatures are given by Krishna and Saxena (1989). These authors used the data from Bennett and Bourne (1972), who reported an operating pressure of circa 137 atm and an operating temperature in the range 402.5 ° C to 411 ° C, using a type of catalyst commonly used in gasoil hydrorefining with additional acidic properties, for the processing of Kuwaiti, Iranian Light and Libyan vacuum gas oil and Kuwaiti heavy vacuum gas oil.


a07f01.gif (4468 bytes)


Krishna and Saxena (1989) adjusted the values of ki from data for a Kuwaiti vacuum gas oil (feed final boiling point around 621 ° C), but they showed that these values also provide results that are in good agreement for other oils such as Iranian Light and Libyan vacuum gas oil (under the same operating conditions), which indicates that the model could be used for other kinds of feedstocks. In order to test the model, these authors assumed that the constants derived for a cut temperature of 371 ° C also hold for higher temperature cuts of 454 and 538 ° C.


In this study, the reactor model is based on the following assumptions:

(a) steady state process;

(b) uniform cross section for the reactor;

(c) uniform temperature (isothermal process);

(d) uniform density for the liquid phase;

(e) flat radial velocity profile throughout the reactor (plug flow);

(f) flat radial concentration profile: concentrations only change with axial position;

(g) hydrogen in excess in the gas phase and oil completely saturated with dissolved hydrogen:
rate of reactions does not depend explicitly on partial pressure of hydrogen;
there is no need for a mass balance for hydrogen;

(h) the axial dispersion coefficient is considered to be the same for all pseudocomponents.

Mass Balance Equations

The overall mass balance for the liquid phase, using the assumption of uniform cross section and uniform density for the liquid phase, leads to a constant superficial velocity, U. The material balance equations, for all pseudocomponents are given by

a07i01.gif (642 bytes)

with the following boundary conditions (Fogler, 1992):

a07i02.gif (499 bytes)

a07i03.gif (397 bytes)

Using dimensionless variables, the material balance equations can be written as

a07i04.gif (725 bytes)

with the following boundary conditions:

a07i05.gif (710 bytes)


a07i06.gif (559 bytes)



The system of ordinary differential equations that describe the process was solved numerically using the collocation method. This method was chosen due to the smoothness of the solution of this particular model and the boundary conditions at the ends of the reactor. The polynomial trial function proposed for the concentration profiles was

a07i07.gif (807 bytes)

where Pi(x) is the dimensionless concentration Yi and x is the dimensionless axial position. Since the polynomials have nine undetermined coefficients, the function requires nine collocation points, two points satisfying the boundary conditions and seven internal points satisfying the ordinary differential equation for each individual pseudocomponent.

Although equation (4) indicates that the pseudocomponent concentrations are related, there is no need to solve (4) simultaneously for all pseudocomponents. In Figure 1, it can be seen that the concentration profiles can be calculated sequentially in the following order: heavy aromatics, light aromatics, heavy naphthenes, light naphthenes, heavy paraffins and light paraffins. The concentration profile calculated for a previous pseudocomponent is used to calculate the concentration profile for the other pseudocomponents in the sequence.

Therefore, applying the trial polynomial at the boundary points, rearranging the coefficients and then applying the resulting function in mass balace, it results in a set of residual functions given by

a07i08.gif (4938 bytes)


a07i09.gif (4489 bytes)

where Ga and B8(x) are given by for heavy aromatics:

a07i10.gif (1038 bytes)

for light aromatics:

a07i11.gif (977 bytes)

for heavy naphthenes:

a07i12.gif (1088 bytes)

for light naphthenes:

a07i13.gif (511 bytes)

a07i14.gif (983 bytes)

for heavy paraffins:

a07i15.gif (904 bytes)

for light paraffins:

a07i16.gif (203 bytes)

a07i17.gif (1298 bytes)

The inlet molar concentration, Ci0, can be calculated from the inlet weight fraction by

a07i18.gif (341 bytes)

and since only ratios of inlet concentrations are actually used in equations (22), (24), (26), (28) and (30), we do not need to use rl0. The concentration can also be expressed in mass basis, instead of molar basis, avoiding the need of using the molecular weight and stoichiometric coefficients for the pseudocomponents, since each one is a mixture of many different compounds. If the initial concentration of any pseudocomponent is zero, then instead of using a dimensionless concentration,Yi(x), given by (8), the weight fraction, wi, could be used in order to avoid division by zero in the expressions for B8(x). After making the necessary rearrangements in the material balance, equations (11) through (18), (19), (21), (23), (25), (27) and (29) remain the same, while equations (20), (22), (24), (26), (28) and (30) should be modified by multiplying the right-hand sides by Ci0 (the concentration in the denominator) and replacing all Ci0 and Cj0 by wi0 and wj0 where Pi(x) in the expressions of B8(x) represents wi. The internal collocation points are selected according to orthogonality conditions, which depend on interval and weight functions (Villadsen and Stewart, 1967). In this case, the orthogonal condition for a set of functions f i(x) is given by

a07i20.gif (546 bytes)

where the interval is [0,1] and the weight function is x× (1-x), since 0 and 1 are also collocation points (boundary conditions). The classes of functions that satisfy (32) are the Jacobi polynomials Gn(p,q,x) with p=3 and q=2 (Abramowitz and Stegun, 1972). Since we are using seven internal points in this work, they will be the roots of G7(3,2,x), given by (Abramowitz and Stegun, 1972)

a07i21.gif (1013 bytes)

with the roots presented in Table 1.


a07t01.gif (2492 bytes)




The model was tested numerically with a program using Pascal as a programming language. The values for the parameters used in the model are presented in Table 2. The kinetic constants were taken from Krishna and Saxena (1989) for a cut temperature of 371 ° C (to distinguish heavy from light). However, these authors did not mention the molecular weight and the stoichiometric coefficients for the lumps. Therefore, all concentrations in the equations of the model were expressed in mass basis. The other parameter values were arbitrarily set. The results are shown in Figures 2 through 11.


a07t02.gif (5497 bytes)



a07f02.gif (2971 bytes)


Figures 2 and 3 present the concentration profiles for heavy and light aromatics, respectively, for different values of reactor length L. It can be seen that the extent of reaction increases with the increase in L, which is due to the increase in residence time. Similar behavior is found for the other pseudocomponents. However, a higher net production of light aromatics can be achieved for lower values of L, since this pseudocomponent is both produced and consumed in the reaction network.


a07f03.gif (3272 bytes)


Figures 4 and 5 present the concentration profiles for heavy and light naphthenes, respectively, for different values of liquid superficial velocity U. It can be seen that the extent of reaction decreases with the increase in U, which is due to the decrease in residence time. Similar behavior is found for the other pseudocomponents. Also, a higher net production of light naphthenes can be achieved for lower values of U.


a07f04.gif (3311 bytes)



a07f05.gif (3392 bytes)


Figures 6 through 11 present the concentration profiles for all pseudocomponents, for different values of the Peclet number, Pe, with fixed L and U (which implies different values for the axial dispersion coefficient). It can be seen that the increase in the Peclet number, Pe, makes the behavior of the reactor come close to that of a plug flow reactor with no axial dispersion, as expected. The shape of the concentration profiles for each pseudocomponent depends on several factors. The concentration of heavy aromatics decreases monotonically due to the fact that this species is only consumed in the reaction network. The concentration of light aromatics and heavy naphthenes reaches its maximum between the ends of the reactor for high values of Pe, since these pseudocomponents are both produced and consumed in the reactor network. The concentration of light naphthenes increases monotonically, despite the fact that this pseudocomponent is both produced and consumed in the reaction network, because the rate of production is greater than the rate of consumption in the situations tested (L £ 2.4m). The concentration of heavy paraffins decreases monotonically, despite the fact that this pseudocomponent is both produced and consumed in the reaction network, because the rate of production is lower than the rate of consumption in the situations tested (L£ 2.4m). The concentration of light paraffins increases monotonically for any reactor length L or superficial velocity U because this pseudocomponent is only produced in the reaction network.


a07f06.gif (3127 bytes)


a07f07.gif (3193 bytes)


a07f08.gif (2973 bytes)


a07f09.gif (3185 bytes)


a07f10.gif (3093 bytes)


a07f11.gif (3242 bytes)


Equations (11) through (30) were solved using the orthogonal collocation points presented in Table 1. The resulting linear systems were all solved using the Gaussian elimination method. In order to test the influence of the collocation points, these equations were also solved using points arbitrarily set at x1 = 0.1, x2 = 0.2, x3 = 0.3, x4 = 0.5, x5 = 0.7, x6 = 0.8 and x7 = 0.9. There was no noticeable change in the results presented in Figures 2 through 11, which indicates that polynomials of degree 8 (equation (10)) are adequate to represent the concentration profiles for the situations tested, independently of the set of collocation points used.



The results indicate that the approach is easy to apply and it can be useful for future work in modelling and simulation of thermal hydrocracking of heavy fractions from petroleum. The approach is general, since it can be applied to different oils and under different operating conditions. However, it would be interesting to have data fitted for the specific oil being processed. Also, the values for the kinetic constants used in the simulations in this work only apply to the operating temperature and pressure and the catalyst used by Bennett and Bourne (1972). Although the reaction network can be extended to other oils, the values of the constants at the temperature and pressure of operation, the kind of catalyst used, and if the hydrocracking is catalytic or thermal all should be known. The usefulness of the approach also depends on the reliability of the methods of detailed compositional analysis of the heavy oil to be processed, since knowledge of the inlet weight fractions and molecular weights of the lumps must be accurate.



The authors would like to thank FAPESP, Fundação de Amparo à Pesquisa do Estado de São Paulo, for the funding that supported this research.



Ah heavy aromatics

Al light aromatics

Ci concentration of lump i, mol/cm3

Ci0 inlet concentration of lump i, mol/cm3

D axial dispersion coefficient, m2/h

ki constant of reaction for path i, 1/h

L reactor length, m

Mi molecular weight of lump i

Nh  heavy naphthenes

Nl   light naphthenes

Pe  Peclet number

Ph  heavy paraffins

Pl   light parafins

U  liquid superficial velocity, m/h

wi weight fraction of lump i

wi0 inlet weight fraction of lump  i

x  dimensionless axial position

z  axial position, m

aij parameter obtained from network

e1 liquid hold-up

ri0  inlet liquid density, g/cm3

Yi dimensionless concentration for lump i paraffins



Abramowitz, M. and Stegun, I. A., eds., Handbook of Mathematical Functions, M. Dover Publications, Inc., New York (1972).

Bennett, R.N. and Bourne, K.H., Hydrocracking for Middle Distillate – A Study of Process Reactions and Corresponding Product Yields and Qualities, Symposium on Advances in Distillate and Residual Oil Technology, New York, pp. G45-G62 (1972).

Fogler, S. H., Elements of Chemical reactor Engineering, Second Edition, Prentice Hall, 1992.

Matos, E.M., Modelagem e Simulação da Cinética Química do Tratamento das Frações Pesadas do Petróleo, Campinas - SP, Master’s thesis, Chemical Engineering, Universidade Estadual de Campinas (1997).

Mohanty, S., Saraf, D.N., and Kunzru, D., Hydrocracking: A Review. Fuel, vol. 69, December (1990).

Mosby, J.F., Buttke, R.D., Cox, J.A., and Nikolaides, C., Process Characterization of Expanded-Reactors in Series, Chemical Engineering Science, vol. 41, no. 4, pp. 989-995 (1986).

Krishna, R. and Saxena, A.K., Use of an Axial-dispersion Model for Kinetic Description of Hydrocracking. Chemical Engineering Science, vol. 44, no. 3, pp. 703-712 (1989).

Villadsen, J.V. and Stewart, W.E., Solution of Boundary-value Problems by Orthogonal Collocation, Chemical Engineering Science, vol. 22, pp. 1483-1501 (1967).

Creative Commons License All the contents of this journal, except where otherwise noted, is licensed under a Creative Commons Attribution License