当前位置:首页 >> 能源/化工 >>

Dynamic Modeling and Simulation of Catalytic Hydrotreating


Energy & Fuels 2006, 20, 936-945

Dynamic Modeling and Simulation of Catalytic Hydrotreating Reactors
Fabian S. Mederos,? Miguel A. Rodr?guez,?,? Jorge Anche

yta,*,?,§ and Enrique Arce§ ? ?
Instituto Mexicano del Petroleo, Eje Central Lazaro Cardenas 152, Col. San Bartolo Atepehuacan, Mexico ? ? ? ? D.F. 07730, Facultad de Ingenier?a, UNAM, Ciudad UniVersitaria, Mexico D.F. 04510, and ESIQIE-IPN, ? ? UPALM, Mexico D.F. 07738 ? ReceiVed December 6, 2005. ReVised Manuscript ReceiVed March 9, 2006

This paper describes a dynamic heterogeneous one-dimensional model of trickle-bed reactors used for catalytic hydrotreating of oil fractions. The model considers the main reactions in the hydrotreating process of oil fractions: hydrodesulfurization, hydrodenitrogenation, and hydrodearomatization. The dynamic model was first validated using experimental data obtained in an isothermal pilot reactor during hydrotreating of vacuum gas oil over a commercial NiMo catalyst. Then, the model was applied to predict the dynamic behavior of a commercial hydrotreating reactor. Changes in concentration, partial pressure, and temperature profiles are obtained and discussed as a function of reactor axial position and time. The simulations obtained with the proposed dynamic model showed good agreement with experimental data.

1. Introduction Catalytic hydrotreating (HDT) is one of the most important processes in the petroleum-refining industry from technical, economic, and environmental points of view. HDT process has been used for over 60 years to reduce the polluting compound content (sulfur, nitrogen, aromatics, etc.) in fossil fuels and to fulfill the applicable legal norms of gas emissions. Nevertheless, for the next years (in mid-2006), new regulations will be applied worldwide and automotive fuels will have to meet the so-called ultralow sulfur specifications (30 ppm in gasoline and 15 ppm in diesel).1 Among other changes (i.e., new reactor design, use of new generation catalysts, etc.), these environmental restrictions will force the refineries to increase the severity in the operating conditions of the HDT plants to match these more stringent specifications.2-4 With these new reaction severity conditions, the behavior of HDT reactors will be altered; for instance, light hydrocarbons may be produced, H2 consumption will increase, the yield of liquid fuel will reduce, reactor delta-T will also increase, etc. To face this situation, refiners require either detailed experimental data or mathematical modeling tools that allow them to understand the phenomena occurring in the HDT reactor, with the main purpose of establishing the optimal operating conditions. Except for naphtha hydrodesulfurization (HDS), in which a gas-solid fixed-bed reactor is used, reactors employed in all HDT applications operate in trickle-bed flow. A trickle-bed reactor (TBR) consists of a cocurrent downward flow of gas and liquid over a fixed bed of catalyst particles (porous solids)
* To whom correspondence should be addressed. E-mail: jancheyt@imp.mx. ? Instituto Mexicano del Petroleo. ? ? UNAM. § ESIQIE-IPN. (1) Kemsley, J. C&EN 2003, October 27, 40-41. (2) Irvine, R. L.; Varraveto, D. M. Ptq Summer 1999, 37-44. (3) Lamourelle, A. P.; Nelson, D. E.; McKnight, J. Ptq Summer 2001, 51-56. (4) Palmer, R. E.; Torrisi, S. P. Ptq ReVamps & Operations 2004, February 7, 15-18.

in which different reactions occur.5 Of course, hydrodesulfurization is the most known reaction occurring in HDT processes. TBR reactors have been the subject of many investigations, and several authors have summarized them in various reviews.6,7 Modeling of TBRs for several applications has also received attention;8-18 however, in the case of TBRs for HDT of oil fractions, most of the reports describe pseudo-homogeneous and heterogeneous models on steady state,19-35 while studies on
(5) Satterfield, C. N. AIChE J. 1975, 21, 209-228. (6) Biardi, G.; Baldi, G. Catal. Today 1999, 52, 223-234. (7) Dudukovic, M. P.; Larachi, F.; Mills, P. L. Catal. ReV. 2002, 44, 123-246. (8) Hofmann, H. Int. Chem. Eng. 1977, 17, 19-28. (9) Ramachandran, P. A.; Smith, J. M. Chem. Eng. Sci. 1979, 34, 7591. (10) Warna, J.; Salmi, T. Comput. Chem. Eng. 1996, 20, 39-47. ¨ (11) Hanika, J.; Lange, R. Chem. Eng. Sci. 1996, 51, 3145-3150. (12) Valerius, G.; Zhu, X.; Hofmann, H.; Arntz, D.; Haas, T. Chem. Eng. Process. 1996, 35, 11-19. (13) Lange, R.; Gutsche, R.; Hanika, J. Chem. Eng. Sci. 1999, 54, 25692573. (14) Salmi, T.; Warna, J.; Toppinen, S.; Ronnholm, M.; Mikkola, J. P. ¨ ¨ Braz. J. Chem. Eng. 2000, 17. (15) Wilhite, B. A.; Wu, R.; Huang, X.; McCready, M. J.; Varma, A. AIChE J. 2001, 47, 2548-2556. (16) Dietz, A.; Julcour, C.; Wilhelm, A. M.; Delmas, H. Catal. Today 2003, 79, 293-305. (17) Lange, R.; Schubert, M.; Dietrich, W.; Grunewald, M. Chem. Eng. ¨ Sci. 2004, 59, 5355-5361. (18) Silveston, P. L.; Hanika, J. Can. J. Chem. Eng. 2004, 82, 11051142. (19) Mears, D. E. Chem. Eng. Sci. 1971, 26, 1361-1366. (20) Henry, H. C.; Gilbert, J. B. Ind. Eng. Chem. Process Des. DeV. 1973, 12, 328-334. (21) Dudukovic, M. P. AIChE J. 1977, 23, 940-944. (22) Froment, G. F.; Depauw, G. A.; Vanrysselberghe, V. Ind. Eng. Chem. Res. 1994, 33, 2975-2988. (23) Korsten, H.; Hoffmann, U. AIChE J. 1996, 42, 1350-1360. (24) Bhaskar, M.; Valavarasu, G.; Meenakshisundaram, A.; Balaraman, K. S. Pet. Sci. Technol. 2002, 20, 251-268. (25) Chowdhury, R.; Pedernera, E.; Reimert, R. AIChE J. 2002, 48, 126135. (26) Avraam, D. G.; Vasalos, I. A. Catal. Today 2003, 79-80, 275283. (27) Pedernera, E.; Reimert, R.; Nguyen, N. L.; van Buren, V. Catal. Today 2003, 79-80, 371-381.

10.1021/ef050407v CCC: $33.50 ? 2006 American Chemical Society Published on Web 04/13/2006

Catalytic Hydrotreating Reactors

Energy & Fuels, Vol. 20, No. 3, 2006 937

dynamic modeling and simulation with either pseudo-homogeneous or heterogeneous models are scarce.36-39 The objective of the present contribution is to illustrate the use of a dynamic one-dimensional TBR model to predict the dynamic behavior of a commercial HDT unit. 2. Brief Literature Review 2.1. Importance of Dynamic Modeling of TBRs. The study of modeling and simulation of TBRs for HDT of oil fractions in dynamic conditions is of great interest for process and control engineers because of the following main reasons:40 (i) It becomes a tool for the design of new TBRs. (ii) Automation and optimizing control has been and continues to be introduced in the industry on an ever-increasing scale. (iii) Because the optimum operation of a reactor is often located near constraint boundaries imposed by, e.g., the strength of construction materials, catalyst deactivation, safety considerations, etc., the optimization is only feasible when suitable dynamic models of reactors are available. (iv) When predicting the dynamic behavior of the reactor during start-ups, shutdowns, or disturbances in the process, it is possible to establish the control strategies adapted to face such events. (v) To plan the necessary modifications in the operating conditions caused by changes in feed properties or product quality. (vi) In the case of HDT of heavy oils, catalyst activity declines relatively rapidly; this requires cyclic operation with dynamic changes in concentration and temperature profiles during operation. (vii) To be used as tool for training. Most of the TBR models reported in the open literature to describe the HDT of oil fractions only take into consideration the HDS reaction because of the complexity of the reaction system. Recently, a few papers have included the following reactions: hydrodearomatization (HDA), hydrodenitrogenation (HDN), hydrogenation of olefins (HGO), and hydrocracking (HDC).25-29,31,32,38 Pseudo-homogeneous models are not capable to describe the real behavior of TBRs for HDT of oil fractions because of the following reasons: (1) they do not reproduce in a suitable way the inhibition effect of hydrogen sulfide on the conversion of organic sulfur compounds; (2) mass and heat transfers between the phases present in the system are not considered; and (3) there is not distinction between the liquid-gas phase and the solid catalyst. Therefore, to obtain reliable estimations and scale up of pilot-plant data to the commercial HDT unit, it is then necessary to develop three-phase reactor models that give a more realistic performance than the pseudo-homogeneous models.14,23
(28) Bhaskar, M.; Valavarasu, G.; Sairam, B.; Balaraman, K. S.; Balu, K. Ind. Eng. Chem. Res. 2004, 43, 6654-6669. (29) Cheng, Z.; Fang, X.; Zeng, R.; Han, B.; Huang, L.; Yuan, W. Chem. Eng. Sci. 2004, 59, 5465-5472. (30) Mac?as, M. J.; Ancheyta, J. Catal. Today 2004, 98, 243-252. ? (31) Melis, S.; Erby, L.; Sassu, L.; Baratti, R. Chem. Eng. Sci. 2004, 59, 5671-5677. (32) Rodr?guez, M. A.; Ancheyta, J. Energy Fuels 2004, 18, 789-794. ? (33) Yamada, H.; Goto, S. Korean J. Chem. Eng. 2004, 21, 773-776. (34) Sertic-Bionda, K.; Gomzi, Z.; Saric, T. Chem Eng. J. 2005, 106, 105-110. (35) Stefanidis, G. D.; Bellos, G. D., Papayannakos, N. G. Fuel Process. Technol. 2005, 86, 1761-1775. (36) van Hasselt, B. W.; Lebens, P. J. M.; Calis, H. P. A.; Kapteijn, F.; Sie, S. T.; Moulijn, J. A.; van den Bleek, C. M. Chem. Eng. Sci. 1999, 54, 4791-4799. (37) Chen, J.; Ring, Z.; Dabros, T. Ind. Eng. Chem. Res. 2001, 40, 32943300. (38) Lopez, R.; Dassori, C. G. SPE Int. 2001, SPE 69499, 1-14. (39) Hastaoglu, M. A.; Jibril, B. E. Chem. Eng. Commun. 2003, 190, 151-170. (40) van Doesburg, H.; de Jong, W. A. Chem. Eng. Sci. 1976, 31, 4551.

In the following paragraphs, a very brief description of steadystate and dynamic models for HDT of oil fractions reported in the literature to simulate HDT units at pilot and commercial scales is presented. 2.2. Steady-State Models. Froment et al.22 developed a onedimensional heterogeneous model using rate equations of model sulfur compounds (benzothiophene, dihydrobenzothiophene, dibenzothiophene, and 4,6-dimethyldibenzothiophene) to describe diesel HDS. Korsten and Hoffmann23 have presented a three-phase reactor model for describing the HDS reaction of vacuum gas oil in a TBR operated under isothermal conditions; the reaction rate was described by a Langmuir-Hinshelwood kinetic. Bhaskar et al.24 used the model of Korsten and Hoffmann23 to simulate the HDS of an atmospheric gas oil fraction. The first attempt toward a thorough description of the HDA reaction was made by Chowdhury et al.25 Such a model accounts for three classes of aromatic compounds, namely, mono-, di-, and tri-aromatics; this model, in addition to the hydrogenation of aromatics in a diesel oil, also considers the HDS reaction rate. Avraam et al.26 developed a model for TBRs applied to HDT of light oil feedstocks with volatile compounds. Mass balances, overall two-phase flow momentum balance, and phase energy balances were described in detail, and four general reactions were taken into account: HDS, HDN, HDA, and HGO. Perdernera et al.27 investigated the influence of the compositions of the oil fractions on the conversion of sulfurous compounds. The reactor model was tested with experimental data obtained at lab scale with straight run gas oil as feed. The model used by these authors is an extension of that presented by Chowdhury et al.,25 which includes the HDS, HDA, HDN, and HDC rate equations. Bhaskar et al.28 have developed a more complete model to simulate the performance of pilot-plant and industrial TBRs applied to the HDS of diesel fractions. It employs a three-phase heterogeneous model and is based on the two-film theory. The major HDT reactions are modeled: HDS, HDN, HDA, HGO, and HDC. The model was used to simulate an industrial reactor using kinetic parameters obtained from pilot-plant experiments. Cheng et al.29 have investigated the performance of a fixed-bed reactor in concurrent and countercurrent flows to remove sulfur and aromatics in diesel fuel. The model presented by this group is a one-dimensional heterogeneous model that accounts for the HDS and HDA reactions to simulate the concentration profiles of the reactants and products in the gas, liquid, and solid phases. The effect of different catalyst particle shapes on the HDS reaction was studied by Mac?as and Ancheyta.30 They employed ? an isothermal heterogeneous reactor model, which was validated with experimental information obtained from a small HDS reactor using straight-run gas oil as feed. Melis et al.31 presented a reactor model for the interpretation of the HDA reaction during gas-oil HDT. This model only considers the HDA reaction by means of a lumped scheme for gas-oil composition and assumes that hydrogenation/dehydrogenation reactions occur according to the Langmuir-Hinshelwood mechanisms. Rodr?guez and ? Ancheyta32 have extended the model of Korsten and Hoffmann23 to include mathematical expressions for the rate of HDS, HDN, and HDA reactions. The HDS reaction was described by kinetic equations of the Langmuir-Hinshelwood type; the HDN reaction was modeled as a consecutive reaction scheme in which nonbasic compounds are hydrogenated first to basic nitrogen compounds (HDNNB), which undergo further reactions to eliminate the nitrogen atom from the molecule (HDNB); and the HDA reaction was represented by a first-order reversible

938 Energy & Fuels, Vol. 20, No. 3, 2006

Mederos et al.

Table 1. Correlations To Estimate Heat- and Mass-Transfer Coefficients, Gas Solubilities, and Properties of Oil and Gases23 parameter oil density correlation FL(P,T) ) F0 + ?FP - ?FT ?FP ) [0.167 + (16.181 × 10-0.0425F0)](P/1000) 0.01[0.299 + (263 × 10-0.0603F0)](P/1000)2 ?FT ) [0.0133 + 152.4(F0 + ?FP)-2.45](T - 520) [8.1 × 10-6 - 0.0622 × 10-0.764(F0+?FP)](T - 520)2 Hi ) VN/λiFL λH2 ) -0.559729 - 0.42947 × 10-3T + 3.07539 × 10-3(T/F20) + 1.94593 × 10-6T2 + (0.835783/F2 ) 20 λH2S ) exp(3.367 - 0.008 47T) L L L 1/2 ki aL/Di ) 7(GL/?L)0.4(?L/FLDi ) ?L ) 3.141 × 1010(T - 460)-3.444[log10(API)]a a ) 10.313[log10(T - 460)] - 36.447 DL ) 8.93 × 10-8(V0.267/V0.433)(T/?L) i L i V ) 0.285V1.048 c m 0.2896 -0.7666 VC ) 7.5214 × 10-3TMeABPd15.6 kS/DLaS ) 1.8(GL/aS?L)1/2(?L/FLDL)1/3 i i i aS ) (6/dp)(1 - ∈) jH ) (hLS/cpLuLFL)(cpL?L/kL)2/3

Henry coefficient solubility of hydrogen solubility of H2S gas-liquid mass-transfer coefficient dynamic liquid viscosity diffusivity molar volume liquid-solid mass-transfer coefficient specific surface area liquid-solid heat-transfer coefficienta

Taken from ref 41.

reaction. This model was validated with experimental information obtained during the HDT of vacuum gas oil in a pilotplant reactor operated under isothermal conditions. Yamada and Goto33 also used the model proposed by Korsten and Hoffman23 to simulate the HDS of vacuum gas oil in a TBR for both modes of cocurrent and countercurrent flows. Sertic-Bionda et al.34 reported the influence of some reaction parameters on HDS in an experimental TBR, using atmospheric gas oil and light cyclic oil from FCC as feeds. A simple reactor and kinetic models were used to simulate the HDS conversion. Recently, Stefanidis et al.35 presented a study on the improvement of the representative operating temperature from the temperature profile of an industrial adiabatic reactor, which is used to simulate the reactor performance by laboratory-scale isothermal reactors. To validate this estimated temperature, a steady-state pseudohomogeneous plug flow model with no resistance in mass and heat transfer was developed to describe the mass balance of sulfur, hydrogen sulfide, hydrogen consumption, and hydrogen as well as the heat balance in the adiabatic reactor. 2.3. Dynamic Models. In this section, the few studies reported in the open literature about dynamic modeling and simulation of HDT reactors are summarized. van Hasselt et al.36 developed a dynamic model to compare the conventional mode of operation of a cocurrent TBR with two different configurations proposed for the operation of the reactor: semi-countercurrent and countercurrent, taking as a case of the study the HDS of a vacuum gas oil. The authors reported that the application of countercurrent flow results in a significant increase of conversion. Chen and Ring37 investigated the dynamic behavior of a fixed-bed pilot-plant hydrotreater of a partially stabilized lightcoker naphtha using a pseudo-homogeneous two-dimensional reactor model. Their model considers the heat conduction in the thermowell of the reactor to predict the temperature difference between the thermowell and the catalytic bed, which, if it is too high and ignored, could cause errors in the interpretation of pilot-plant data. Lopez and Dassori38 reported a dynamic model to simulate the HDT of vacuum gas oil; however, it was applied only as a steady-state model because they omitted the accumulation term. Hastaoglu and Jibril39 generated a two-dimensional reactor model for gas-solid reactions to describe fixed-bed reactors. The model was validated through a comparison of experimental and theoretical results from naphtha HDS pilot plant.

3. Reactor Model
3.1. Dynamic Model. A dynamic heterogeneous one-dimensional reactor model was considered in this work, which is based on the three-phase steady-state reactor model reported in the literature.23,32 This model is based on the two-film theory8 and makes use of correlations reported in the literature to estimate heat- and masstransfer coefficients, gas solubilities, and properties of oil and gases under process conditions,23 which are presented in Table 1. The mass transfer of the compounds in the TBR is described with the following set of partial differential equations (PDEs). The reactor model considers that there are no reactions in the gas phase, and HDS, HDA, HDNB, and HDNNB reactions are taking place along the catalytic bed. The dynamic mass-balance equation in the catalyst bed for the gaseous compounds is pG ?pG uG ?pG i i i - CL )- kLaL i i RT ?t RT ?z Hi




where i ) H2, H2S, and NH3. The dynamic mass-balance equation in the catalyst bed for the gaseous compounds in the liquid phase is ?CL ?CL pG i i i ) -uL + kLaL - CL - kSaS(CL - CS) L i i i i i ?t ?z Hi




where i ) H2, H2S, and NH3. The model assumes that the organosulfur, organonitrogen, and aromatic compounds, as well as the liquid hydrocarbons, are nonvolatile; therefore, the dynamic mass-balance equation for the liquid compounds is ?CL ?CL i i ) -uL - kSaS(CL - CS) i i i ?t ?z (3)


where i ) S, HC, NB, NNB, and A. The components transported between liquid and solid phases are consumed or produced by the chemical reaction at the surface of the catalyst, according to the following first-order ordinary differential equations (ODEs):

- ∈)

?CS i ) kSaS(CL - CS) ( FB ηjrj(CS, ..., TS) i i i i ?t


where i ) H2, H2S, NH3, S, HC, NB, NNB, and A and j ) HDS, HDNNB, HDNB, and HDA. The “-” sign is for the reactants, and the “+” sign is for the products. The reaction rate for ammonia is rNH3 ) -rHDNB + rHDNNB.

Catalytic Hydrotreating Reactors
Table 2. Kinetic Models for HDS, HDA, HDNNB, and HDNB reaction HDS HDA HDN nonbasic nitrogen basic nitrogen kinetic model
S S rHDS ) kHDS[(CS)(CH2)0.45/(1 + KH2SCH2S)2] S G rHDA ) kf pH2CS - kr(1 - CS ) A A S rHDNNB ) kHDNNB(CNNB)1.5 S S rHDNB ) kHDNNB(CNNB)1.5 - kHDNB(CNB)1.5

Energy & Fuels, Vol. 20, No. 3, 2006 939 3.3. Model Boundary Conditions. Because the mathematical model is a system of PDEs with time and spatial coordinates as independent variables, it is necessary to define the following initial and boundary conditions: (i) Initial conditions:
G G For t ) 0, at z ) 0, pH2 ) (pH2)0

Because the concentration of hydrocarbons (the main component of the feed) does not change significantly during HDT, eqs 3 and 4 with i ) HC will not be considered in the dynamic simulations. To model the commercial HDT reactors operating under nonisothermal conditions and taking into account that HDT reactions are exothermic in nature, the following energy-balance equations were used:22 (i) Liquid phase ?TL ?TL ) -uLFLcpL - hLSaS(TL - TS) LFLcpL ?t ?z (ii) Solid phase (1 - ∈)FScpS ?TS ?t ) hLSaS(TL - TS) + (5)

pG ) 0, i ) H2S and NH3 i CL ) (CL)0, i ) H2, S, NB, NNB, and A i i CL ) 0, i ) H2S and NH3 i CS ) 0, i ) H2, H2S, NH3, S, NB, NNB, and A i T ) T0 at z > 0, pG ) 0, i ) H2, H2S, and NH3 i CL ) 0, i ) H2, H2S, NH3, S, NB, NNB, and A i (11)


FBηjrj(CS, i

CS ) 0, i ) H2, H2S, NH3, S, NB, NNB, and A i ..., TS)(-?HRj) (6) T ) T0 (ii) Boundary conditions:
G G For t > 0, at z ) 0, pH2 ) (pH2)0


It is usually sufficient to consider only the liquid phase in the energy-balance equations (eqs 5 and 6), because the heat capacity of the gas phase is much lower compared with that of the liquid phase.14 3.2. Reactions Kinetic Models. One way to represent the HDS reaction is the practical and widely accepted generalized stoichiometric equation that lumps the HDS reaction of all of the sulfur compounds into a single expression υSS(liquid) + υH2H2(gas) f υHCHC(liquid) + υH2SH2S(gas) (7)

pG ) 0, i ) H2S and NH3 i CL ) (CL)0, i ) H2, S, NB, NNB, and A i i CL ) 0, i ) H2S and NH3 i CS ) 0, i ) H2, H2S, NH3, S, NB, NNB, and A i T ) T0 (13)

where υS, υH2, υHC, and υH2S are the stoichiometric coefficients of the organic sulfur compounds, hydrogen, sulfur-free hydrocarbon, and hydrogen sulfide, respectively. The HDA reaction was represented by the following first-order reversible reaction: A {k } B \



The HDN reaction was modeled as a consecutive reaction scheme in which nonbasic nitrogen (NNB) compounds are hydrogenated first to basic nitrogen (NB) compounds, which undergo further reactions to eliminate the nitrogen atom from the molecule according to the following reaction scheme: NNB 9 NB 9 HC + NH3 8 8


For commercial HDT reactors, values of partial pressures (pG) i and liquid molar concentrations (CL) of H2S and NH3 at the i entrance of the catalytic bed (z ) 0) are different to 0. 3.4. Integration Method. The PDEs describing the heat and mass transfer in the reactor were transformed into a set of firstorder ODEs by discretization in the axial direction using the backward finite difference method and leaving the independent variable time without discretize. The ODEs were then solved using a fourth-order Runge-Kutta method. The program required to calculate transport coefficients and to simultaneously solve the system of ODEs was coded in FORTRAN.



Kinetic expressions of HDS, HDA, HDNNB, and HDNB reactions are shown in Table 2.23,42,43 In the case of the HDS reaction, a Langmuir-Hinshelwood kinetic model is employed, which includes an adsorption equilibrium constant of hydrogen sulfide (KH2S) to account for the influence of the temperature described by the van’t Hoff equation23 KH2S(T) ) 41 769.8411 exp

4. Results and Discussion 4.1. Dynamic Modeling of an Isothermal HDT Pilot Plant. Kinetic parameters reported by Rodr?guez and Ancheyta32 were ? used for dynamic modeling (Table 3). These data were generated from experimental information obtained at 5.3 MPa pressure, 340-380 °C temperature, 2000 standard cubic feet per barrel H2/oil ratio, and 2.0 h-1 LHSV in a pilot plant (reactor, 143 cm in length and 2.54 cm in diameter) described elsewhere44 using a vacuum gas oil and a NiMo/Al2O3 commercial catalyst, which main properties are presented in Table 4. The catalyst was loaded into the reactor diluted by an equal volume of inert material (SiC).
(44) Marroqu?n, G.; Ancheyta, J. Appl. Catal., A 2001, 207, 407-420. ?

(2761) RT


(41) Carberry, J. Chemical and Catalytic Reaction Engineering; McGrawHill: New York, 1976. (42) Yui, S. M.; Sanford, E. C. 1985 ProceedingssRefining Department: 50th Midyear Meeting (Kansas City, MO, May 13-16, 1985); American Petroleum Institute: Washington, DC, 1985; pp 290-297. (43) Bej, S. K.; Dalai, A. K.; Adjaye, J. Energy Fuels 2001, 15, 377383.

940 Energy & Fuels, Vol. 20, No. 3, 2006
Table 3. Kinetic Parameters and Heat of Reaction Used for Dynamic Modeling32 reaction HDS HDN nonbasic nitrogen basic nitrogen HDA forward reverse

Mederos et al.

Ea (J/mol) 131 993 164 942 204 341 121 400 186 400 4.266 × 109 cm3

k0 g-1 s-1 (cm3 mol-1)0.45

?HR a (J/mol) -251 000 -64 850 -255 000

3.62 × 106 s-1 (wt %)-0.5 3.66 × 1011 s-1 (wt %)-0.5 1.041 × 105 s-1 MPa-1 8.805 × 109 s-1

Taken from ref 45. Table 4. Main Properties of HDT Feed and Catalyst property feed API gravity molecular weight mean average boiling point (°C) sulfur (wt %) total nitrogen (wppm) basic nitrogen (wppm) total aromatics (wt %) catalyst equivalent diameter (mm) specific surface area (m2/g) pore volume (cm3/g) mean pore diameter (?) molybdenum content (wt %) nickel content (wt %) bulk density (g/cm3) 2.54 175 0.56 127 10.7 2.9 0.8163 22 441.9 476 2.009 1284 518 41.9 value

Dynamic simulations using the reactor model described in section 3 were carried out to observe the behavior of different product properties and operating conditions with time. Because in this case the reactor is operated isothermally, only the dynamic mass-balance equations were considered. Figures 1 and 2 show how sulfur, nitrogen (basic and nonbasic), and total aromatic contents in the product at the exit of the reactor change as a function of time. The four concentration profiles shown in these figures are quite similar. A small amount of hydrotreated product is detected at the exit of the reactor at about 250 s (0.07 h), which corresponds to the mean residence time given by the interstitial velocity of the liquid phase (u′L ) uL/ L); after that, concentrations start increasing and finally the steady-state is reached at 2300 s (0.64 h). Experimental values are also shown in Figures 1 and 2 (symbol “O”). The concentration of aromatics starts increasing earlier at about 150 s. This is because, different from the HDS and HDN kinetic models, HDA takes into account the hydrogen

Figure 2. Concentrations of different compounds (s, simulated; O, experimental) at the outlet of the catalytic bed (z ) 31.58 cm) as function of time (Pilot reactor, 380 °C; 5.3 MPa; uL ) 1.75 × 10-2 cm/s; uG ) 0.28 cm/s).

Figure 1. Concentration of organic sulfur compound (s, simulated; O, experimental) at the outlet of the catalytic bed (z ) 31.58 cm) as function of time (Pilot reactor, 380 °C; 5.3 MPa; uL ) 1.75 × 10-2 cm/s; uG ) 0.28 cm/s).

partial pressure in the kinetic model and the gas velocity is about 16 times higher than liquid velocity (uG/uL ) 16). The solution of the steady-state reactor model was also obtained solving the corresponding equations,32 and identical concentration profiles were found in both cases. When compar-

Catalytic Hydrotreating Reactors

Energy & Fuels, Vol. 20, No. 3, 2006 941

Figure 3. H2 partial pressures and concentration profiles down through the catalytic bed at different times. ([) 60 s, (???) 500 s, (×) 1000 s, and (s) 2300 s.

Figure 4. H2S partial pressures and concentration profiles down through the catalytic bed at different times. ([) 60 s, (???) 500 s, (×) 1000 s, and (s) 2300 s.

ing the predicted concentrations with experimental values, sulfur, nonbasic nitrogen, and aromatics are estimated very well (absolute errors of 0.6, 0.4, and 1.3%, respectively), while basic nitrogen showed the highest deviation (absolute error of 4.8%). It is possible that, in the kinetic model for basic nitrogen (rHDNB in Table 2), in our case NB disappearance may not follow a reaction order of 1.5 as found by Bej et al.43 However, because of insufficient experimental data, this cannot be confirmed. Partial pressure and concentration profiles of hydrogen and hydrogen sulfide along the reactor are shown in Figures 3 and 4, at times very close to the beginning of the operation (60 s), intermediate times (500 and 1000 s), and when the steady state is reached (2300 s). In Figure 3, it is observed that because of the high gas velocity inside the reactor the curves of the hydrogen partial pressure at 500, 1000, and 2300 s are overlapping. It is also seen from Figure 4 that a small amount of H2S production is present only at the first section of the reactor at 60 s and, at 2300 s, the observed H2S profile reported by previous papers along the reactor is found.22-24,26,32,33,38 The overall shapes of molar concentration profiles of H2 and H2S are determined by the balance between the reaction rate and mass transfer. The H2 concentration decreases while the H2S concentration increases at the beginning of the reactor because of the high reaction rate of the catalyst bed in this zone. This behavior can be clearly seen at different times. Figure 5 presents the dynamic simulated liquid molar concentration profiles of sulfur, nonbasic nitrogen, basic nitrogen, and aromatics along the reactor at several times ranging from 60 to 2300 s. In the graph of sulfur, when the model predicts the steady state at 2300 s (s), it is observed that the concentrations in both the liquid and solid phases are not equal along the reactor.

The results at the steady state obtained from dynamic simulations using the model developed in this work ([, ???, ×, and s) agree reasonably with experimental concentrations (O) of sulfur, nonbasic nitrogen, basic nitrogen, and aromatics at the exit of the reactor. Thus, the model is demonstrated to accurately predict experimental data obtained at the isothermal reactor mode of operation. 4.2. Dynamic Simulation of a Commercial HDT Reactor. After the dynamic reactor model was validated to reproduce the experimental data obtained in the isothermal reactor at the pilot-plant scale, it was applied to predict the expected dynamic behavior of a commercial HDT reactor. The following characteristics of the reactor were taken from the literature:32 internal diameter, 3.048 m; total length, 9.1 m; and catalytic bed length, 8.534 m. Because of the nonisothermal operation of the commercial reactor, the energy balance given by eqs 5 and 6 were solved simultaneously with mass-balance equations. Other reaction conditions are: pressure of 54 kg/cm2, LHSV of 2.66 h-1, and H2/oil ratio of 2000 ft3/barrel. Additionally, a highpurity hydrogen stream with the following molar composition was used for dynamic simulation: 81.63% hydrogen, 3.06% hydrogen sulfide, and 15.31% lights. As an example, Figure 6 illustrates the dynamic simulation of the HDS reaction using eqs 5 and 6, in which results of the dynamic simulation of the isothermal reactor are also shown for comparison. We can see that the steady-state in the commercial unit is reached most fast than that found in the pilot reactor. This is not surprising because operating conditions, particularly LHSV and uG/uL ratio, are greater. There are some differences in the sulfur concentration profiles at certain period of time before the steady-state operation is reached in pilot reactor and commercial reactor, which can be attributed to the

942 Energy & Fuels, Vol. 20, No. 3, 2006

Mederos et al.

Figure 6. Sulfur concentration and temperature (s and ???, simulated; s, commercial reactor; ???, pilot reactor; O, experimental) at the oulet of the catalytic bed as function of time. (Commercial reactor, 380 °C; 5.3 MPa; z ) 853.44 cm; uL ) 0.63 cm/s; uG ) 10.27 cm/s. Pilot reactor, 380 °C; 5.3 MPa; z ) 31.58 cm; uL ) 1.75 × 10-2 cm/s; uG ) 0.28 cm/s).

constant temperature in the former and the variation of temperature in the latter. Figure 7 shows the predicted sulfur, nonbasic nitrogen, basic nitrogen, and aromatics contents in the product at different times for an inlet reactor temperature of 380 °C. The dynamic simulation was carried out at the same reaction conditions than those employed for simulation of the pilot reactor except LHSV (2.66 h-1). The values of concentrations reported at the exit of the isothermal pilot reactor are included in this figure (symbol “O”). Pilot-plant experimental concentrations values are in general higher than those predicted for the commercial reactor because of the higher average catalytic bed temperature observed in the nonisothermal commercial operation. The energy balance given by eqs 5 and 6 considers the heat generation on the solid surface and transfer of this heat to the liquid phase. Sometimes to simplify the model, one can lump these two processes in only one equation for a pseudo-phase. When doing this, all of the complexity introduced during the mass balances seems to be inconsistent with the consideration of a pseudo-homogeneous solid-liquid phase. To corroborate this, we have also simulated the commercial HDT reactor with the following pseudo-homogeneous heat balance:41

?T ?t

) -uLFLcpL

?T ?z


∑(FBηjrj(CS, ..., T)(-?HRj)) i j


Figure 5. Sulfur, nonbasic nitrogen, basic nitrogen, and aromatics concentration profiles down through the catalytic bed at different times. ([) 60 s, (???) 500 s, (×) 1000 s, and (s) 2300 s.

The dynamic simulation of the reactor using this latter equation gives very similar results as those presented in Figure

Catalytic Hydrotreating Reactors

Energy & Fuels, Vol. 20, No. 3, 2006 943

Figure 8. Temperature profiles down through the catalyst bed of the commercial reactor with time (a) using eq 14 and (b) using eqs 5 and 6. ([) 60 s, (???) 500 s, (×) 1000 s, and (s) 2000 s.

Figure 7. Concentration profiles of sulfur, nonbasic nitrogen, basic nitrogen, and aromatics down through the catalyst bed of the commercial reactor with time. (O) experimental pilot-plant data, ([, ???, ×, and s) predicted values for commercial reactor. ([) 60 s, (???) 500 s, (×) 1000 s, and (s) 2000 s.

6. Concentration profiles are practically the same, and the only imperceptible difference is in the temperature profile, in which the minimum temperature at the exit of the reactor is reached first by the pseudo-homogeneous energy balance, 500 versus 550 s. The difference in predictions using pseudo-homogeneous (eq 14) and heterogeneous (eqs 5 and 6) heat balances is more noticeable in parts a and b of Figure 8, respectively, in which the dynamic predicted axial profile of the reaction temperature in the commercial reactor is presented. The major discrepancy between these plots is the intermediate temperature profiles before reaching the steady state. In the profile at 500 s (???), we can see that an increase in the reaction temperature at the initial part of the reactor (40 and 35% of the catalyst bed for parts a and b of Figure 8, respectively), because of the removal of contaminants occurring in this zone, generates the behavior known as “wrong way”, which is reported in the literature to cause a decrement in the downstream temperature. This wrongway behavior occurs because a rapid heating of the feed increases the conversion in the entrance section of the reactor and therefore decreases the concentration in the cold outlet section of the reactor, which leads to a transient temperature diminution.46 Similar results are found at 1000 s, but in this case, the decrease in temperature is observed at about 80 and 75% of the catalyst bed for parts a and b of Figure 8, respectively. For times longer than 1000 s, the wrong-way behavior was not encountered.
(45) Tarhan, M. O. Catalytic Reactor Design; McGraw-Hill: New York, 1983. (46) Pinjala, V.; Chen, Y. C.; Luss, D. AIChE J. 1988, 34, 1663-1672.

944 Energy & Fuels, Vol. 20, No. 3, 2006

Mederos et al. B ) Saturated hydrocarbon cpj ) Specific heat capacity of j phase (J g-1 K-1) Cji ) Molar concentration of compound i in the j phase (mol cm-3) dp ) Particle diameter (cm) d15 ) Liquid density at 15 °C (g cm-3) d15.6 ) Specific gravity at 15.6 °C Dji ) Molecular diffusivity of compound i in the j phase (cm2 s-1) Ea ) Activation energy (J mol-1) GL ) Liquid superficial mass velocity (g cm-2 s-1) hLS ) Heat-transfer coefficient for liquid film surrounding the catalyst particle (J s-1 cm-2 K-1) Hi ) Henry’s law constant for compound i (MPa cm3 mol-1) HDA ) Hydrodearomatization reaction HDN ) Hydrodenitrogenation reaction HDS ) Hydrodesulfurization reaction jH ) j factor for heat transfer kf ) Forward HDA rate constant (s-1 MPa-1) kj ) Apparent j ()HDS, HDNB, and HDNNB) reaction rate constant, see Table 3 kL ) Thermal conductivity of liquid phase (J s-1 cm-1 K-1) kr ) Reverse HDA rate constant (s-1) k0 ) Frequency factor, see Table 3 kji ) Mass-transfer coefficient of compound i at the interface j (cm s-1) KH2S ) Adsorption equilibrium constant for H2S (cm3 mol-1) Lb ) Length of catalyst bed (cm) pji ) Partial pressure of compound i in the j phase (MPa) P ) Reactor total pressure (psia) rj ) Reaction rate j (mol cm-3 s-1, for j ) HDS mol g-3 s-1) R ) Universal gas constant (J mol-1 K-1) t ) time (s) T ) Temperature (K) TMeABP ) Mean average boiling point (°R) uj ) Superficial velocity of j phase (cm s-1) z ) Axial coordinate (cm) Greek Symbols ?HRj ) Heat of reaction j (J mol-1) ?FT ) Temperature correction of liquid density (lb ft-3) ?FP ) Pressure dependence of liquid density (lb ft-3) ∈ ) Bed void fraction j ) Hold up of j phase p ) Particle porosity ηj ) Catalyst effectiveness factor for reaction j λi ) Solubility coefficient of the compound i (Nl kg-1 MPa-1) FB ) Catalyst bulk density (g cm-3) Fj ) Density at process conditions of j phase (lb ft-3) F0 ) Liquid density at standard conditions (15.6 °C; 101.3 kPa) (lb ft-3) F20 ) Liquid density at 20 °C (g cm-3) ?L ) Absolute viscosity of the liquid (mPa s) νc ) Critical specific volume of the gaseous compounds (cm3 mol-1) νi ) Molar volume of solute i at its normal boiling temperature (cm3 mol-1) νL ) Molar volume of solvent liquid at its normal boiling temperature (cm3 mol-1) νN ) Molar gas volume at standard conditions (Nl mol-1) νCm ) Critical specific volume (ft3 lbm-1) ) Fractional volume of the catalyst bed diluted by inert particles Subscripts 0 ) Reactor inlet condition A ) Aromatics G ) Gas phase HC ) Desulfurized or denitrogeneted hydrocarbon H2 ) Hydrogen

Figure 9. Variation with time and axial position of the temperature in the catalyst bed of the commercial reactor.

Figure 9 summarizes the simulated dynamic temperature along the catalyst bed as a function of time for the commercial reactor (using eqs 5 and 6 for the energy balance). It is clearly observed that the temperature in the first part of the reactor remains more or less without changes with time, whereas the temperature near the outlet of the reactor changes significantly with time. At short times, the already mentioned wrong-way behavior is observed, and at long times, the reactor temperature stabilizes. Conclusions A TBR model for simulating the dynamic behavior of pilot and commercial HDT reactors was developed in this work. The model is heterogeneous one-dimensional and includes the most important HDT reactions (HDS, HDN, and HDA). The dynamic model is capable of simulating isothermal and nonisothermal modes of operation. Experimental data of HDT of vacuum gas oil were employed for model validation, and the model simulation results agree reasonably well with experimental ones. The solution of the model considering only the steady-state equations gave the same results as the dynamic model when this one gets the steady-state solution (by running the dynamic model for a long time). The steady state was found to be almost the same for both reactor scales, isothermal pilot reactor and commercial reactor. The steady state in the pilot reactor scale was obtained at 2300 s, whereas for the commercial reactor, it was obtained at 2000 s. Concentrations at the exit of the pilot reactor were measured when this steady state was obtained. For short times (<1000 s), the known wrong-way behavior was observed when simulating temperature axial profiles for the commercial reactor. This phenomenon caused a decrease in the reactor temperature in the final part of the reactor. At times close or higher than the steady state of simulations, results for this behavior were not observed.
Acknowledgment. The authors thank Instituto Mexicano del Petroleo for its financial support. M. A. R. also thanks CONACyT ? for financial support.

a ) Dimensionless number of Glaso’s correlation aj ) Specific surface area at the interface j (cm-1) A ) Aromatic compound API ) API gravity

Catalytic Hydrotreating Reactors H2S ) Hydrogen sulfide j ) Reaction (HDS, HDNB, HDNNB, or HDA) L ) Liquid phase or gas - liquid interface NB ) Basic nitrogen NNB ) Nonbasic nitrogen NH3 ) Ammonia S ) Organic sulfur compound, solid phase or liquid-solid interface Superscripts

Energy & Fuels, Vol. 20, No. 3, 2006 945

G ) Gas phase L ) Liquid phase or gas-liquid interface S ) Solid phase or liquid-solid interface