Modelling and simulation of the hydrocracking of heavy oil fractions

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,

E-mail: [e-mail]

(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

INTRODUCTION

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.

MATHEMATICAL MODEL

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 k_i 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.

Krishna and Saxena (1989) adjusted the values of k_i 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.

Considerations

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

with the following boundary conditions (Fogler, 1992):

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

with the following boundary conditions:

where

NUMERICAL RESOLUTION

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

where P_i(x) is the dimensionless concentration Y_i 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

where

where G_a and B₈(x) are given by for heavy aromatics:

for light aromatics:

for heavy naphthenes:

for light naphthenes:

for heavy paraffins:

for light paraffins:

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

and since only ratios of inlet concentrations are actually used in equations (22), (24), (26), (28) and (30), we do not need to use r_l⁰. 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,Y_i(x), given by (8), the weight fraction, w_i, could be used in order to avoid division by zero in the expressions for B_8(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 C_i⁰ (the concentration in the denominator) and replacing all C_i⁰ and C_j⁰ by w_i⁰ and w_j⁰ where P_i(x) in the expressions of B₈(x) represents w_i. 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

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 G_n(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 G₇(3,2,x), given by (Abramowitz and Stegun, 1972)

with the roots presented in Table 1.

NUMERICAL RESULTS

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.

Figures 2

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.

Figures 6 through 11 present the concentration profiles for all pseudocomponents, for different values of the Peclet number, P_e, 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, P_e, 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 P_e, 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.

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 x₁ = 0.1, x₂ = 0.2, x₃ = 0.3, x₄ = 0.5, x₅ = 0.7, x₆ = 0.8 and x₇ = 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.

DISCUSSION AND CONCLUSIONS

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.

ACKNOWLEDGMENTS

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.

NOMENCLATURE

Ah heavy aromatics

Al light aromatics

C_i concentration of lump i, mol/cm³

C_i⁰ inlet concentration of lump i, mol/cm³

D axial dispersion coefficient, m²/h

k_i constant of reaction for path i, 1/h

L reactor length, m

M_i 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

w_i weight fraction of lump i

w_i⁰ inlet weight fraction of lump i

x dimensionless axial position

z axial position, m

a_ij parameter obtained from network

e₁ liquid hold-up

r_i⁰ inlet liquid density, g/cm³

Y_i dimensionless concentration for lump i paraffins

REFERENCES

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).

Publication Dates

• Publication in this collection
  24 Apr 2000
• Date of issue
  Mar 2000

History

• Accepted
  08 Sept 1999
• Received
  09 Apr 1999