SciELO - Scientific Electronic Library Online

vol.19 issue3Nonlinear dynamic modeling of multicomponent batch distillation: a case studyConversion of n-Butane to iso-Butene on Gallium/HZSM-5 Catalysts 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.19 no.3 São Paulo July/Sept. 2002 



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


(Received: December 03, 2001 ; Accepted: June 11, 2002)



Abstract - This work presents a model to predict the behavior of velocity, gas holdup and local concentration fields in a pseudo-two-phase gas-liquid column reactor applied for thermal hydrocracking of petroleum heavy fractions. The model is based on the momentum and mass balances for the system, using an Eulerian-Eulerian approach. Using the k-e model,fluid dynamics accounts for both laminar and turbulent flows, with discrete small bubbles (hydrogen) flowing in a continuous pseudohomogeneous liquid phase (oil and catalyst particles). The petroleum is assumed to be a mixture of pseudocomponents, grouped by similar chemical structural properties, and the thermal hydrocracking is taken into account using a kinetic network based on these pseudocomponents.
Keywords: slurry bubble column, reactor, hydrocracking, heavy oil, modeling.




Three-phase and two-phase flow column reactors have many applications in industry. One of them is in the hydrocracking process, where large molecules of petroleum in the presence of excess hydrogen are broken into smaller molecules. In this process there is an increase in the amount of valuable oil subproducts, which is of great economic importance.

The modeling and simulation of such reactors are quite complex, since they need to take into account both the fluid dynamics and the reactions that occur. There are a number of works in the literature about the modeling of these reactors (Celik and Wang, 1994; Chen et al, 1995; Gasche et al, 1990; Grienberger and Hofmann, 1992; Hillmer et al, 1994; Torvik and Svendsen, 1990), but they are usually applied to a different set of reactions. The present work presents a model that considers these features, specifically applied to the hydrocracking of oils, using some reasonable approximations in order to solve the problem numerically. The oil, which is a complex mixture of many compounds, is assumed to be a mixture of a small number of pseudocomponents, so that the reaction network is based on these pseudocomponents.

A good mixture and uniform distribution of temperature is achieved with this equipment. Therefore, in this work the reactor is considered an isothermal system, where thermal hydrocracking reactions occur, with excess hydrogen (gas phase) flowing upward concurrently with the heavy oil (liquid phase) and the catalyst particles (solid phase).

The hydrocracking modeled in this work is assumed to be thermal. The catalyst is only used to promote other kind of reactions, such as the removal of heteroatoms. However, the presence of the catalyst does affect the fluid dynamics. Considering that the catalyst particles are very small and have a small terminal velocity, in this work the slurry (oil + solid catalyst) was assumed to be a pseudo-homogeneous liquid phase.



The Kinetic Network

The hydrocracking reactions considered here are based on the kinetic network proposed by Krishna and Saxena (1989), established in terms of pseudocomponents, as shown in Figure 1. A pseudocomponent contains a set of molecules that have similar chemical structures. 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; otherwise they are considered heavy. The values for the constants of reaction are given by Krishna and Saxena (1989).




In this study, the following assumptions were made:
i)  steady state process;
ii)  uniform cross sectionfor the reactor (cylinder);
iii) uniform temperature inside the reactor
iv)  pseudo-two-phase model (slurry + gas);
v)  uniform density for the liquid phase;
vi)  excess hydrogen in the gas phase and oil completely saturated with dissolved hydrogen (so that rates of reactions do not explicitly depend on the partial pressure of hydrogen).

Since the catalyst particles considered in this work are very small, the solid and liquid phases were assumed to be a pseudohomogeneous slurry phase. However, the solid phase does not enter or leave the reactor, since it remains confined within it. In this work we have not calculated the solid distribution inside the reactor. We have only considered that the solid concentration in the slurry is uniform, which means

where fs  is the average solid concentration in the slurry inside the reactor.

The density and the holdup of the slurry phase are given by

For laminar flows, the viscosity for the slurry phase is given by (Hillmer et al., 1994)

For turbulent flows, the k-e model was used to calculate the turbulent viscosity. The effective viscosity is given by

In this work, the value of Rp was assumed to be 1.

The dispersion coefficient can be calculated by (Hillmer et al., 1994)

Mass Balance Equations

Since we are modeling the oil fractions as pseudocomponents, all concentrations and rates of reaction will be expressed on a mass basis. Therefore, the mass rates per unit volume of liquid (not slurry) ri  satisfy the equation

where ri > 0 for formation and ri < 0  for consumption.

Gas Phase (Mainly Hydrogen)

where and NP  are the mass transfer rates of hydrogen and products between gas and slurry phases, respectively. Since in this model we are assuming a large excess of hydrogen, these two terms were neglected.

The dispersion of the gas phase (small bubbles) is caused mainly by the turbulence in the slurry phase. For laminar flows, the dispersion terms are zero. For turbulent flows, all variables and properties are time averaged.

Slurry Phase

Under the previously mentioned assumptions, the global mass balance for the slurry phase results in the following equation:

where  and NP  are the mass transfer rates of hydrogen and products between gas and slurry phases, respectively. These two terms were also neglected in this equation, because they are smaller than the others.

The mass balance for each pseudocomponent j in the slurry phase, based on the kinetic network in Figure 1, is given by

where wj,sl  is the mass fraction of j in the slurry, rj,sl  is the rate of reaction of j in the slurry and Nj is the mass transfer rate of each pseudocomponent j between gas and slurry phases. These mass transfer terms were also neglected in this equation.

It is important to note that the summation of equation (13) over all species in the slurry (dissolved hydrogen, solid catalyst, pseudocomponents) should result in equation (12), the global mass balance for the slurry phase.

It is desirable to express the concentration of j in mass of j/volume of liquid phase and not slurry, since the rates of reaction are usually expressed in the liquid phase. It can easily be shown that

and since rsl and fs (equation 1) are considered constants in this model, equation (13) can be written as

Using the kinetic network shown in Figure 1, the reaction rates for the pseudocomponents are given by

Heavy Aromatics

Light Aromatics

Heavy Naphthenes

Light Naphthenes

Heavy Paraffins

Light Paraffins

An equation similar to (13) may be written for the hydrogen in the liquid phase. Since excess hydrogen is present in the gas phase at a high pressure, flowing as small bubbles, it was assumed that mass transfer coefficients between phases were high enough so that  in the liquid phase was close to the equilibrium concentration and uniformly distributed. Therefore, the mass balance for hydrogen in the liquid may be written as

If we try to calculate the rate of consumption of hydrogen in the liquid phase, using equation (10), we get the following unrealistic result:

because we are assuming that 1 kg of component p converts into 1 kg of component q. Unless we have information about the average molecular weight of each pseudocomponent and the average stoichiometric ratio between pseudocomponents, this model does not allow calculation of the rate of consumption of hydrogen.

Momentum Balance Equations

The momentum and force balances for the system, under the previously mentioned considerations, are given by the following equations (Torvik and Svendsen, 1990), considering gravitational forces in the z direction only.

Gas Phase

Axial Component

Radial Component

Slurry Phase

Axial Component

Radial Component



Turbulent Flow

For turbulent flow, the k-e model was used to describe the turbulent viscosities in the flow field. The equations for the slurry phase are described as follows (Hillmer et al., 1994):




The turbulent viscosity for the slurry phase is calculated from

The turbulent viscosity for the gas phase is calculated from equations (6) and (7), using the turbulent viscosity for the slurry phase.

Boundary Conditions

For this work, the slurry is assumed to be leaving the reactor uniformly at the lateral position at the top of the reactor ( r = R , z = L), while the gas is assumed to be leaving the reactor at the top ( z = L , 0 < r < R). This is in fact an approximation, and the actual area for the slurry outlet is given by a small area at the lateral position. At the inlet, both slurry and gas phases enter at the bottom ( z = 0 , 0 < r < R).

The boundary conditions for equations (11) ¾ (12) and (25) ¾ (28) are given by

where k = g,sl. Observe that the limits for the intervals for conditions (48) and (53) do not apply for vr,sl  at z = L and r = R.

The boundary conditions for equation (16) are given by

For laminar flows conditions (44) and (45) in the wall work fine. However, for turbulent flows the velocity gradients near the wall are very high, making it very difficult to calculate numerically the velocity profiles. So, instead of condition (44), for turbulent flows the Prandtl logarithmic law is used as a boundary condition near the wall (Launder and Spalding, 1974):

where y = R - r is the distance from the wall and tw,k is the wall stress (equation (38) for r = R), which can be related to the pressure gradient by integrating equations (25) and (27).

For the k-e model, the boundary conditions near the wall are (Xu et al., 1998)

Model Resolution

The equations were discretized by applying the finite volume method, which converted the nonlinear equations into a linearized system that was solved by LU decomposition (in each iteration). The staggered grid scheme was used to locate the variables in the control volumes. The upwind scheme for convective terms and the central difference for diffusive terms were used to avoid negative coefficients and a consequently numerical instability. The source terms were linearized in such a way as to ensure that the matrix coefficients were diagonally dominant.

Due to strong nonlinearities, the velocity (momentum) and holdup equations must be underrelaxed and a specific sequence in solving the equations must be obeyed to get convergence (Patankar, 1980).

The gas holdup was evaluated by combination of the continuity equations for the two phases, subtracting one from the other. For the pressure correction, the SIMPLEC method was used (Maliska, 1995), adding up the continuity equations for the two phases.

The model was solved in 12 x 12 grid, using C as the programming language. The criteria for convergence used in the iterative procedure was:

Except pressure, all variables are specified at the inlet of the reactor. All velocities are zero at the wall, since it is impermeable for gas and slurry. At the centerline, radial velocities and the deviates of all variables are zero. At the free surface (end of the reactor), all gradients were set to zero. The outlet velocities are calculated in such a way that the mass balance must be obeyed.

For the pressure field, all derivatives are zero at the boundary (borders), because we can not determine velocities and pressure at the same point with the grid scheme we are using. To avoid indetermination in the pressure field, one point in the grid was chosen as a reference value, since the flows are related to gradient of pressures and not their absolute values (Patankar, 1980). In this work, an arbitrary reference value of pressure equal to zero was set at the last numbered point in the grid.

About the pressure term in the momentum equations, either form (ek . P) or (ek . P) may be used in equations because both satisfy the condition that, when the corresponding momentum equations for the two phases are added, the resulting pressure gradient must be P(Celik and Wang, 1994).

The value for Cw in equations (29) and (30) was given by (Grienberger and Hofmann, 1992):

in kg/m3.s, considering a variation with radial position in order to get a better agreement between experimental and calculated values.

Some approximations were also used to solve the model. The mass transfer between gas and liquid phases was neglected, considering that the hydrogen was present in excess and assuming that the more volatile hydrocarbons remained dissolved in the oil, due to the high pressure. The oil density was assumed to be constant inside the reactor, thus neglecting some of the effects of reaction in the fluid dynamics. Also, the hydrogen density was assumed to be constant, since the pressure drop was much smaller than the absolute pressure in the reactor (around 100 atm).



Air-Water System

In order to test the performance of the model, numerical results were compared with experimental dataavailable from the literature. Figure 2 shows experimental results from an air-water system (Grienberger & Hofmann, 1992) comparing the gas holdup profile in the radial direction with numerical results for the same operating conditions. Figure 3 shows experimental results from an air-water system (Gasche et al., 1990), comparing the gas velocity profile in the radial direction with numerical results for the same operating conditions. Figure 4 shows experimental results from an air-water system (Torvik and Svendsen, 1990), comparing the liquid velocity profile in the radial direction with numerical results for the same operating conditions. In all figures, the squares represent calculated values, while circles represent experimental values. The indicated inlet velocities are the superficial velocities, given by the mass flow of each phase divided by the phase density and the area of cross section of the reactor. In order to get the actual inlet velocities, these values must be divided by the inlet holdup of each phase. In all figures, the profiles were calculated for h = 2.25 m, the middle height of the column, in order to avoid inlet and outlet effects.





It can be observed from the figures that in some cases the numerical results are close to the experimental data (Figures 2 and 3), while in others they are not (Figure 4). However, in all cases the physical behavior is qualitatively correct.

The fluid-dynamic variables (liquid and gas velocities, gas and liquid holdups) influence the mass balance by affecting the residence time and the mixing of the components. When the internal velocities are very high, the physical behavior of the reactor approaches the perfect mixing reactor (CSTR). However, the recirculating movement of the slurry phase leads to dead zones, so that it is important to study the fluid-dynamic profiles in order to get a better understanding of reactor behavior.

Hydrogen-Oil System

The values for the parameters used in the model are shown in Table 1. The kinetic constants were taken from Krishna and Saxena (1989) for a cut temperature of 700 °F (to distinguish heavy from light). As in Hillmer et al. (1994), the dispersion coefficients were expressed in terms of viscosity. The value for CM was taken from Torvik and Svendsen (1990), while the value for CW was calculated by equation (63). The other parameter values were arbitrarily set.



In the numerical simulations, the inlet liquid phase superficial velocity was set to 0.01 m/s, while gas phase superficial velocity at the reactor inlet was 0.08 m/s. The inlet gas holdup was 0.25, while inlet oil holdup was 0.75. The solid catalyst was assumed to be confined inside the reactor. The reactor inner diameter was 0.3 m and the total length of the reactor was 4.5 m.

The oil density was assumed to be 660 kg/m3, while the catalyst density was assumed to be 1660 kg/m3, both at the operating conditions. The slurry density depends on the value of fs (equations (1) and (2)), so that for fs equal to 0.14 the slurry density is 800 kg/m3 and for fs equal to 0.34 the slurry density is 1000 kg/m3.

The numerical results are shown in Figures 5 through 14. In all figures the slurry density was 800 kg/m3, with the exception of Figures 8 and 9, which also used slurry density 1000 kg/m3. Whenever not stated in the figures, for profiles with radial variation the fixed axial position is the middle of the column (h = 2.25 m) and for profiles with axial variation the fixed radial position is the center of the column (r = 0).















Figure 5 shows the variation in gas holdup with radial position for two different heights (axial position) inside the reactor. This figure shows that the gas is more concentrated in the center of the reactor, which is usually observed in this kind of reactor.

Figure 6 shows the variation of slurry axial velocity with radial position for two different heights (axial position) inside the reactor. This figure shows the upward movement in the center of the column and the downward movement close to the wall, which is also usually observed in this kind of reactor.

Figure 7 shows the variation in axial velocities with radial position for liquid and gas. This figure shows that gas velocity is higher than slurry velocity, due to the higher volumetric flow of hydrogen than oil.

Figure 8 shows the variation in gas holdup with radial position for two different slurry densities (800 and 1000 kg/m3). This figure shows that the gas holdup profile has only a small variation with the slurry density.

Figure 9 shows the variation in axial velocity with radial position for slurry for two different slurry densities (800 and 1000 kg/m3). This figure shows that the slurry axial velocity profile has only a small variation with the slurry density.

Both Figures 8 and 9 indicate that the slurry density have a small effect in the fluid dynamics of the column. However, one should not conclude from that that gradients in oil or slurry densities inside the reactor would not change the fluid dynamics, because in the simulated cases we have assumed an uniform slurry density through the reactor, although with different values in each simulation.

Figure 10 shows the variation in pressure with axial position (reactor height). This figure shows the expected drop in pressure with column height due mainly to the weight of the slurry.

Figure 11 shows the fluid dynamic vector field for slurry, with the upward flow in the center of the reactor and downward flow close to the wall.

Figures 12 through 17 show the variation in concentration for each one of the pseudocomponents with axial position. As expected, the only pseudocomponent that increases its concentration is the light paraffin, because it is only produced in the kinetic network. The light naphtenics and light aromatics have a peak in concentration at the bottom of the reactor. All other pseudocomponents have a decrease in concentration.



The results show that the proposed model can be used to describe the behavior of a slurry bubble column reactor, in this case applied to a process of thermal hydroconversion of heavy oils. The fluid dynamic fields were based on momentum transfer and mass balances and are used to calculate the concentrations of the chemical species involved. The kinetic model, based on a pseudocomponent approach, is general and can be applied to different oils and different operating conditions. However, the specific values used in the simulations only apply to the operating conditions used to obtain the data taken from the literature.

The model has some limitations that can be improved, such as the calculations of hydrogen consumption, mass transfer coefficients between phases and solid dispersion inside the reactor.

The kinetic network can be extended to any oils, although compositional details of the heavy oil to be processed must be known. This model can be helpful for future work in optimization and control of the petroleum hydrocracking process, which is always seeks to increase economically the marketable fraction of petroleum.



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.





Celik, I. and Wang, Y.Z., Numerical Simulation of Circulation in Gas-Liquid Column Reactors: Isothermal, Bubbly, Laminar Flow. Int. J. Multiphase Flow, vol. 20, nº 6, pp. 1053-1070, 1994.        [ Links ]

Chen, Z., Zheng, C., Feng, Y. and Hofmann, H., Modelling of Three-phase Fluidized Beds Based on Local Bubble Characteristics Measurements, Chemical Engineering Science, vol. 50, n° 2, pp. 231-236, 1995.        [ Links ]

Gasche, H.-E., Edinger, C., Kömpel, H. and Hofmann, H., A Fluid-dynamically Based Model of Bubble Column Reactors, Chemical Engineering & Technology, 13, pp. 341-349, 1990.        [ Links ]

Grienberger, J. and Hofmann, H., Investigations and Modeling of Bubble Columns, Chemical Engineering Science, 47 (9-11), pp. 2215-2220, 1992.        [ Links ]

Hillmer, G., Weismantel, L. and Hofmann, H., Investigations and Modelling of Slurry Bubble Columns. Chemical Engineering Science, vol. 49, n° 6, pp. 837-843, 1994.        [ Links ]

Krishna, R. and Saxena, A. K., Use of an Axial-dispersion Model for Kinetic Description of Hydrocracking. Chemical Engineering Science, vol. 44, n° 3, pp. 703-712, 1989.        [ Links ]

Launder, B.E. and Spalding D.B., The Numerical Computation of Turbulent Flows. Comput. Methods Appl. Mech. Eng., 3, 269-289, 1974.        [ Links ]

Maliska, C.R., Transferência de Calor e Mecânica dos Fluidos Computacional. Livros Tecnicos e Científicos Editora, 1995.        [ Links ]

Mosby, J.F., Buttke, R.D., Cox, J. and Nikolaides, C., Process Characterization of Expanded-bed Reactors in Series, Chemical Engineering Science, vol. 41, n° 4, pp. 989-995, 1986.        [ Links ]

Patankar, S.V., Numerical Heat Transfer and Fluid Flow, Hemisphere Publishing Corporation, pp. 67, 1980.        [ Links ]

Torvik, R. and Svendsen, H.F., Modelling of Slurry Reactors ¾ A Fundamental Approach, Chemical Engineering Science, vol. 45, nº 8, pp. 2325-2332, 1990.        [ Links ]

Xu, Z.G., Gotham, D.H.T., Collins, M.W., Coney, J.E.R., Sheppard, C.G.W. and Merdjani, S., Validation of Turbulence Models in a Simulated Air-conditioning Unit. International Journal for Numerical Methods in Fluids, vol. 26, pp. 199-215, 1998.        [ Links ]

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