ARTICLE pubs.acs.org/EF
Modeling the Hydrocracking Kinetics of Atmospheric Residue in Hydrotreating Processes by the Continuous Lumping Approach Haitham M. S. Lababidi† and Faisal S. AlHumaidan‡ †
Chemical Engineering Department, College of Engineering and Petroleum, Kuwait University, Post Office Box 5969, Safat 13060, Kuwait ‡ Petroleum Research and Studies Center, Kuwait Institute for Scientific Research, Post Office Box 24885, Safat 13109, Kuwait ABSTRACT: The primary objective of this work is to study the hydrocracking associated with the hydrotreatment of atmospheric residue (AR) feedstock and develop a kinetic model describing the undergoing cracking reactions. Experimental data were obtained for three types of conventional hydrotreating catalysts [hydrodesulfurization (HDS), hydrodemetalization (HDM), and hydrodenitrogenation (HDN)] at three space velocities and three operating temperatures. The developed kinetic cracking model is based on the continuous lumping approach, which assumes a continuous concentration and reactivity of the reacting mixtures. Species concentrations were represented as an integro-differential equation with only five modeling parameters. The developed continuous lumping models predicted the concentration profile of the complete true boiling point (TBP) range with reasonably high accuracy. Analysis of experimental results indicated that upgrading of residual fractions is achieved through both catalytic and thermal conversion, where the extent of each depends upon the catalyst type and operating conditions. The yield of the distillate fraction was found to be the highest for the HDN catalyst followed by HDS and HDM catalysts.
1. INTRODUCTION The heavy consumption of conventional crude oil and the depletion of its reserves have increased the demand of using less desirable heavy crude oil and residues. A number of hydroprocessing technologies are available for upgrading residual materials. One of the most successful and famous residue upgrading technologies today is the atmospheric residue desulfurization (ARDS) process.1,2 The ARDS process was developed back in 1965, and its main objective at that time was to produce desulfurized fuel oil. In recent years, the process has significantly improved, and its objectives today include the removal of feedstock impurities (e.g., S, N, Ni, and V), saturation of unsaturated hydrocarbon compounds (e.g., olefins and aromatics), and hydrocracking the large hydrocarbon molecules to produce the desired fuels and lubes.35 The ARDS process allows for the residue to be partially converted into a light fraction through a sequence of stages in which the characteristics of the product distillates are adjusted and improved. Researchers have extensively investigated the ARDS process, and considerable effort has been directed toward the development of kinetic models that describe the various hydrotreating reactions and catalyst deactivation.611 However, the hydrocracking reactions for atmospheric residue (AR) feedstock have received fewer attention, and therefore, very few models are available for the hydrocracking of heavy feedstock. In a previous publication, AlHumaidan and co-workers12 studied the effect of the catalyst type and reaction severity on the hydrocracking reactions of atmospheric AR feedstock and developed kinetic models that describe the basic characteristics of the hydrocracking reactions using the discrete lumping approach. In this work, the main objective is to use the continuous lumping approach to develop a kinetic model for the mild hydrocracking reactions accompanying the catalytic hydrotreatment of AR feedstock. r 2011 American Chemical Society
2. BACKGROUND Modeling of complex hydrocracking kinetics is necessary for petroleum refining industries, because it allows the refiners to predict the product yields at different operating conditions. Such prediction has a significant impact on process optimization, unit design, and catalyst selection. There are two distinct classes of kinetic models: the lumped empirical models and the detailed molecular models. The detailed molecular models have a high predictive power, as compared to the lumped empirical models, because they account for all possible reactions and mechanisms. However, their application to heavy feedstock is difficult because of the complexity of the mixture and the number of reactions involved.13 The lumped empirical models, on the other hand, are commonly used in modeling the hydrocracking of heavy petroleum fractions. Lumped empirical models have mainly two approaches: the discrete lumping approach and continuous lumping approach.14,15 The discrete lumping is a simplified approach, in which the complex hydrocracking chemistry and kinetics are viewed as a set of model compounds or pseudo-components. In other words, the chemically similar species of the complex mixture are combined or lumped together and treated as pseudo-components. The selection of pseudo-components can be based on product slate, true boiling point (TBP), carbon number, or molecular weight. Discrete lumping models are not characterized by very accurate predictive powers. Nevertheless, their predictive performance is quite sufficient for many applications. The success of discrete lumping models lies in their ease of application and incorporation into reactor models, considering the limited Received: January 27, 2011 Revised: April 18, 2011 Published: April 22, 2011 1939
dx.doi.org/10.1021/ef200153p | Energy Fuels 2011, 25, 1939–1949
Energy & Fuels
ARTICLE
Figure 1. Process flow diagram of the pilot plant.
number of reactions and rate parameters involved. Their disadvantages, on the other hand, are mainly the rapid increase in the number of kinetic parameters as the number of lumps increases and the time-consuming adjustment of rate parameters for each feedstock.12 Various discrete lumping hydrocracking kinetic models were developed in the past decades. Some of these models were developed for the hydrocracking of gas oil and vacuum gas oil1618 while others were developed for the hydrocracking of heavier residual fractions.12,1922 In contrast to discrete lumping, continuous lumping considers the reactive mixtures to form a continuous mixture with respect to its species type, boiling point, molecular weight, etc. Continuous lumping models exhibit better predictive capabilities of product yield than discrete lumping and, therefore, are considered to be a step toward the recognition of the chemical and physical diversity of heavy hydrocarbon feedstock.23 The application of the continuous lumping approach in crude oil cracking has significantly developed over the past few decades. Aris24 provided a continuum description of the first-order reaction in thermal cracking. Weekman25 established the principle of parallel reactions in a continuous mixture to the cracking of crude oil. Chou and Ho23 proposed a continuous lumping approach for the complex mixture whose components undergo nonlinear reactions and demonstrated their lumping procedure on two types of reactions: irreversible nth-order and reversible first-order reactions. Further research in this area was carried out, and significant progress has been achieved.14,2628 Laxminarasimhan and co-workers14 applied continuous lumping to describe the hydrocracking of vacuum gas oil and proposed a model formulation to determine the concentration distribution of the reaction mixture at any given residence time. They tested the developed model on two sets of pilot-plant data published in the literature. Basak et al.26 simulated the hydrocracking kinetics and heat effects in a commercial hydrocracker using the contimuous lumping approach. Their hydrocracker model was capable of predicting the yields, reactor-bed temperature profile, hydrogen consumption, pressure drop, and catalyst deactivation rate. Moreover,
Elizalde et al.28 reported a continuous lumping kinetic model for the hydrocracking of heavy crude oil at moderate conditions. Although continuous lumping models exhibited better predictive capabilities of product yields than discrete lumping, they still lack the incorporation of the fundamental reaction chemistry. Some researchers tried to include as much molecular information as possible to the continuous lumping approach. Quann and Jaffe29 introduced the fundamental reaction chemistry into a lumped model by considering each hydrocarbon as an assembly of a limited number of structural building units. This promising approach that orients the lumping strategy toward the molecular structure was named “structure-oriented lumping”.30,31
3. EXPERIMENTAL SECTION 3.1. Apparatus. The experimental runs were conducted in a microreactor pilot plant. The unit has the ability to operate at elevated temperature and pressure and the flexibility of operating at different liquid hourly space velocities. The microreactor is a flange-type reactor, 63.5 cm in length and 1.2 cm internal diameter. The preheated feed and hydrogen are introduced to the reactor in a down co-current flow mode. The reactor is equipped with a built-in furnace that consists of five independently controlled heating zones. Each heating zone is equipped with a K-type thermocouple to measure and control the reactor skin temperatures. The reactor is also equipped with a thermowell with five K-type thermocouple bundles to monitor the bed temperature. The process flow diagram of the pilot plant is illustrated in Figure 1, while detailed designed specifications of the unit are given by AlHumaidan.20 3.2. Experimental Run. Three sets of experimental runs were conducted, using the same feedstock and three different types of hydrotreating catalysts. The feedstock was Kuwait Export Crude AR (KEC-AR), whose properties and distillation ASTM D-1160 are given in Tables 1 and 2, respectively. Each experimental set included separate runs at three different temperatures (T = 370, 390, and 410 C) and three liquid hourly space velocities (LHSV = 0.5, 1, and 2 h1), which are equivalent to changing the feed flow rates. The operating procedures and sampling schedule for a single catalyst type is illustrated in Figure 2. For 1940
dx.doi.org/10.1021/ef200153p |Energy Fuels 2011, 25, 1939–1949
Energy & Fuels
ARTICLE
each operating condition, two samples of the hydrocracked product were obtained and analyzed. A stabilization period of 36 h was allowed after each change in operating conditions. The total duration of a complete run is 18 days. After catalyst loading, the catalyst presulfiding procedure was carried out, followed by operating the reactor at fixed conditions for some time to ensure stable (middle-of-run) operation. All experimental runs were conducted at a fixed system pressure of 120 bar and a H2/oil ratio of 680. Process conversion was determined by high-temperature simulated distillation (AC SimDist HT-750) that reports the boiling point
Table 1. Properties of KEC-AR property
unit
value
American Petroleum Institute (API) gravity
degree
12.84
sulfur
wt %
4.4
nitrogen
ppm wt
3100 (0.31 wt %)
vanadium
ppm wt
63
nickel
ppm wt
21
iron
ppm wt
5
sodium
ppm wt
2
carbon hydrogen
wt % wt %
83.7 12.2
asphaltenes
wt %
4.5
Conradson carbon
wt %
11.81
kinematic viscosity at 50 C
cSt
765.1
kinematic viscosity at 100 C
cSt
53.77
density at 15 C
g/cm3
0.9774
density at 65 C
g/cm3
0.9442
pour point flash point
C C
21 187
Table 2. Distillation ASTM D1160 for KEC-AR volume % distilled
temperature (C)
initial boiling point (IBP)
367
5 10
391 394
20
438
30
474
40
512
50
560
65
583
final boiling point (FBP)
589
total volume recovery (%)
57.1
distribution up to 750 C. In addition, all samples were analyzed for sulfur, nickel, vanadium, and nitrogen contents. An Oxford 3000 sulfur analyzer was used to determine the sulfur content, while metal contents (Ni and V) were analyzed by inductively coupled plasma (ICP) equipment manufactured by Varian. Moreover, an Antek 7000 elemental analyzer was used for the determination of the total nitrogen content. The SimDist analysis is performed with a capillary chromatographic (GC) column and is based on the assumption that individual components of a sample escape from the column in the order of their boiling point. The method is well-known for its accuracy, repeatability, and rapidity (a complete simulated distillation analysis is performed in 40 min). The correlation between retention time and distillation temperature is normally established through a calibration mixture that covers the complete boiling range of the analyzed sample. The sample should be diluted with carbon disulfide (CS2) at a concentration of 1% by weight prior to its injection. The column temperature is gradually raised, and the area under the chromatogram is recorded through the run. The boiling temperatures are assigned to the time axis using a calibration curve. From these data, the boiling point distribution is obtained. Maximum error between the TBP profiles of the two samples taken at the same operating conditions was only 1.6%, while the average error for all samples was 0.75%. On the other hand, the average error in the sulfur concentration between the two samples was 2.77%. 3.3. Catalysts. Three commercial hydrotreating catalysts were used in this study: desulfurization catalyst (HDS), demetalation catalyst (HDM), and denitrogenation catalyst (HDN). They are the same catalyst types used in the ARDS units of a local refinery. Before the experiments were started, the catalysts were presulfided to create the essential surface required for optimum activity. The catalyst presulfiding practice reduces the coke formation tendency by passivating the highly active acidic oxide sites.32 To enhance the oil distribution and reduce the channeling effect, the catalyst bed was diluted with inert particles (carborandum). More details about catalyst characterization, catalyst loading, and activation are given by AlHumaidan.20 For residual oil hydroprocessing, pilot-plant test runs that extended for 1015 months showed that the catalysts undergo rapid deactivation by coke deposition during the first few hours after presulfiding. After this early deactivation, the main mechanism for deactivation is metal deposition, which is a slow process that extends for a few months.33
4. HYDROCRACKING CONTINUOUS KINETIC MODEL The continuous lumping model developed in this work is based on the formulation reported by Laxminarasimhan and coworkers.14 The main concept in this formulation is that the rate constant of hydrocracking is assumed as a monotonic function of the TPB, and the mass balance equations are reformulated by considering the reaction rate constant as a continuous variable.
Figure 2. Operating procedure and sampling schedule for the experimental runs. 1941
dx.doi.org/10.1021/ef200153p |Energy Fuels 2011, 25, 1939–1949
Energy & Fuels
ARTICLE
4.1. From Discrete to Continuous Lumping. The main difference between the discrete and continuous lumping approaches is the way in which kinetic lumping is performed. In discrete lumping, the kinetic lumping is performed by summation, while in continuous lumping, it is performed through integration. In discrete lumping, a nth-order hydrocracking model for Nl lumps, undergoing Nl series and parallel reactions, may be written as
ri ¼
dci ðtÞ ¼ Ki c i n þ dt
i
∑
ðγji kji cj n Þ,
j¼1 j 6¼ i
ci jt ¼ 0 ¼ cif
ð1Þ
ði ¼ 1, :::, Nl Þ
where ci(t) is the concentration of the i-type species and cif is its feed concentration. kij is the rate constant or reactivity of the itype species cracked to j-type species. γji values are the stoichiometric coefficients that account for the increase in weight that occurs when hydrogen is consumed by the lumped hydrocarbon components.34 Ki is the sum of rate constants for the parallel reactions, represented as Ki ¼
Nl
∑
j¼i þ 1
kij
ði ¼ 1, :::, Nl 1Þ
ð2Þ
AlHumaidan and co-workers12 used the discrete lumping approach in modeling hydrocracking of AR feedstock. Discrete lumping models were developed for 4, 5, and 10 lumps and assuming first and second order reaction schemes. Such models provided a perfect fit to the experimental data and predicted the product yields of the proposed lumps accurately. Simplicity is the main advantage of the discrete lumping approach. However, increasing the number of lumps increases the number of kinetic parameters numerously. For a fixed reaction order, the number of parameters (rate constants) that need to be estimated is 1 /2Nl (Nl 1). Hence, for 4, 5, and 10 lumps, the number of parameters is 6, 10, and 45, respectively. To avoid defining a large number of k values, it would be more practical to define a continuous function relating the reactivity of the i-type species to one of its properties, such as molecular weight or TBP. Equivalently, a continuous concentration function, c(k,t), may be defined to express the concentration of species with reactivity k at time t. Thus, the total reactant concentration can be represented by the following weighted integral: Z ¥ cðk, tÞDðkÞdk ð3Þ CðtÞ ¼ 0
Equation 3 approximates the sum over i-type species, by an integral over reactivity k. Moving from discrete lumps to continuous reactivity may be viewed as a transformation from an “i coordinate” to a “k coordinate”. This transformation is achieved by introducing the species-type distribution function D(k), defined for i-type species transformed into species of reactivity ki as Δi di Dðki Þ ¼ ¼ Δki dk
ð4Þ
number of species types with rate constants between k and k þ dk. For large number of species, eq 4 can be expressed as Z 1 ¥ Dðki Þdk ¼ 1 ð5Þ N 0 4.2. Continuous Model Formulation. During the course of the hydrocracking of a particular feed, the TBP distribution curve of the reaction mixture changes continuously with the changing residence time of the reaction mixture in the reactor. As the residence time of the reaction increases, more of the heavier components in the reaction mixture are converted into lighter components. As a result, the distribution curve of the mixture changes with the increased length of residence time. The TBP curve is converted into a distribution function describing the weight fraction of any component as a function of normalized TBP, which may be defined as
θ¼
TBP TBPmin TBPmax TBPmin
ð6Þ
TBPmin and TBPmax correspond to the lowest and highest boiling points of the reaction mixtures, respectively. It is required now to define the transformation function from the normalized TBP, θ, and the species reactivity, k. This transformation can be described as a monotonic power law type14 k ¼ θ1=R kmax
ð7Þ
where kmax is the reactivity (rate constant) of the species with the highest TBP, namely, θ = 1, and R is a model parameter. The species-type distribution function, D(k), can be rewritten as di di dθ ¼ ð8Þ dk dθ dk Assuming that the species indices, i values, are equally spaced, then the term di/dθ can be approximated to N, where N is the total number of species (pseudo-components) in the mixture (N f ¥). In addition, the dθ/dk term can be obtained by differentiating eq 7. Hence, D(k) becomes DðkÞ ¼
DðkÞ ¼ NR
kR 1 kRmax
ð9Þ
which satisfies the normalization function defined by eq 5. Thus, coordinate transformation resulted in defining a continuous species-type distribution function, D(k), which is a function of species reactivity k and the parameter R. The rate of cracking reaction, defining the rate of disappearance of species with reactivity k, can be formulated as14 Z kmax dcðk, tÞ ¼ k cðk, tÞ þ pðk, KÞ K cðK, tÞ DðKÞ dK dt k ð10Þ where p(k,K) is a yield distribution function that determines the amount of formation of species of reactivity k from the cracking of the species with reactivity K. This distribution function can be defined as
While D(ki) is defined only at the discrete points of ki, it can be treated as a continuous function of k as the number of lumps, N, approaches infinity. In view of that, D(k)dk will represent the
pðk, KÞ ¼ 1942
1 pffiffiffiffiffiffiðA B þ CÞ S0 2π
ð11Þ
dx.doi.org/10.1021/ef200153p |Energy Fuels 2011, 25, 1939–1949
Energy & Fuels where
ARTICLE
8 !2 9 a0 < ðk=KÞ 0:5 = A ¼ exp : ; a1 n o B ¼ exp ð0:5=a1 Þ2
and
k C ¼ δ 1 K
ð12Þ
ð13Þ
ð14Þ
The parameters a0 and a1 defined in eqs 12 and 13 are model parameters that determine the location of the maximum reactivity in the interval k∈(0,K). δ is another model parameter, which accounts for the fact that p(k,K) should be a finite small quantity when k = 0. The term B, defined by eq 13, comes from the condition p(k,K) = 0 when k = K. The parameter S0 defined in eq 11 is estimated as follows: Z K 1 ðA B þ CÞDðkÞdk ð15Þ S0 ¼ pffiffiffiffiffiffi 2π 0 4.3. Calculation Procedure. The continuous lumping model formulated above consists of five parameters: kmax, R, a0, a1, and δ. These parameters can be estimated by fitting the experimental data to the model defined above using nonlinear regression. This involves solving N integro-differential equations derived from eq 10. To obtain the desirable accuracy, a large number of species (pseudo-components) should be considered (N f ¥). Accordingly, the TBP experimental data (in the range from TBPmin to TBPmax) was divided into 50 equally spaced divisions. Nonlinear regression and numerical solution of the integrodifferential equations were performed using MATLAB programming language. The calculation procedure is summarized by the flowchart shown in Figure 3. The MATLAB regression function “lsqcurvefit” minimizes the error between the experimental and calculated values of C(k,t) by searching for the best values of the five model parameters. For each iteration, the system of 50 integro-differential equations are solved using a program developed as a m-file. Returned regression results consist of the best fitted model parameters (kmax, R, a0, a1, and δ) and the predicted concentrations of the lumps. The best fit is quantified by determining the coefficient of determination, R2, evaluated as follows:
R2 ¼
SST SSE SST
ð16Þ
where SSE ¼
N
N
exp exp 2 2 ðCi Ccalc ∑ i Þ and SST ¼ ∑ ðCi CÞ i¼1 i¼1
ð17Þ
calc are the experimental and calculated concentrations Cexp i and Ci of the lump i, respectively, while C is the average of the experimental concentrations.
5. EXPERIMENTAL RESULTS The AR feedstock and the cracked products are lumped as pseudo-components, termed as naphtha (N), middle distillate (D), vacuum gas oil (G), and vacuum residue (V), with boiling
Figure 3. Calculation procedure for nonlinear regression integrated with the solution of the system of integro-differential equations.
ranges C1160 C, 160360 C, 360525 C, and þ525 C, respectively. These fractions are usually obtained from ARD units in refineries. The yield of the naphtha fraction is normally very low, while the middle distillate fraction yields diesel and jet fuel, which are normally the most significant hydrocracked products of residue hydrotreating. The gas oil fraction is an important cut because it is used as a feedstock for fluid catalytic cracking and hydrocracking, which are mainly operated for the production of transportation fuels. The demand for the vacuum residue fraction is very low because of its deteriorating effects on the catalyst of downstream processes, such as fluid catalytic cracking. The AR feed consists of no naphtha, 7.37% distillate, 38.66% gas oil, and 53.97% vacuum residue. Hydrocracking conversion results were determined using the SimDist analysis data. Table 3 lists the yields of the four lumps for the three types of catalyst and different LHSVs and temperatures. Effects of temperature, residence time, and catalyst type are illustrated in the SimDist TBP plots shown in Figures 4, 5, and 6, respectively. 1943
dx.doi.org/10.1021/ef200153p |Energy Fuels 2011, 25, 1939–1949
Energy & Fuels
ARTICLE
Table 3. Conversion Analysis Results concentration (wt %) LHSV (h1)
temperature (K)
370
lump
390
410
0 9.81
11.15
11.2
G
44.18
40.16
41.77
R
46.01
48.69
47.03
1.0
390
410
370
0.5
390
410
0
0
HDN
D
0.4
0
0
D
11.35
13.43
14.21
G
44.27
41.91
44.83
R N
44.38 0.72
44.26 0.92
40.96 0.9
D
14.13
16.73
18.19
G
46.03
43.85
46.89
R
39.12
38.5
34.02
N 370
HDS
N
N 2.0
HDM
0
0
Figure 4. Effect of the reaction temperature on conversion for the HDS catalyst at LHSV = 1 h1.
0.69
D
10.15
11.79
12.01
G
44.13
41.23
42.46
R N
45.72 0
46.98 0.71
44.84 1.1
D
12.81
16.58
15.48
G
45.26
46.98
44.46
R
41.93
35.73
38.96
N
1.28
1.54
1.29
D
17.55
20.41
21.94
G
46.5
45.02
47.05
R N
34.67 0
33.03 0
29.72 0.77
D
11.6
14.08
13.82
G
45.08
42.28
43.89
R
43.32
43.64
41.52
N
0.58
0.93
D
14.64
17.5
18.71
G
45.4
43.95
46.15
R N
39.38 1.47
37.62 2.15
34.23 2.56
D
20.48
25.9
27.7
G
47.31
46.15
46.38
R
30.74
25.8
23.36
Figure 5. Effect of the residence time on conversion for the HDS catalyst at T = 410 C.
0.91
Experimental results showed (see Table 3) that conversion of the AR fraction to lighter fractions varied from 9.8% for HDS catalyst at 370 C and 2 h1, to the maximum conversion of 56.7% for HDN catalyst at 410 C and 0.5 h1. This achieved conversion is considerably high, taking into account that hydrocracking is a secondary objective of the hydrotreating processes. The results also indicated clearly that the cracked fractions are converted mainly to the middle distillate cut. For the lowest severity, the concentration of the middle distillate in the cracked product increased by around 50%, while at the highest severities, the concentration of the middle distillate almost tripled (275%). On the other hand, the increase in the concentration of the VGO fractions varied between 4 and 20%. Moreover, the naphtha cut appeared only at medium to high severities with a maximum concentration of 2.56% (Table 3). 5.1. Effect of the Temperature. The effect of the temperature on the conversion of AR and generation of lighter fractions are
Figure 6. Effect of the catalyst type on conversion for T = 410 C and LHSV = 1 h1.
illustrated by the shift of the TBP curve to the left of the feed TBP curve, as shown in Figure 4, for HDS catalyst type. Similar trends were also obtained for the HDN and HDM catalyst types. The TBP curves clearly show how the thermal effect enhances the yields of the distillate and VGO. The results also showed that the effect of the temperature on conversion is not linear; that is, more impact is shown on conversions above 390 C. For example, the production of the distillate fraction has a faster rate 1944
dx.doi.org/10.1021/ef200153p |Energy Fuels 2011, 25, 1939–1949
Energy & Fuels
ARTICLE
in the range of higher temperatures (390410 C), as compared to its generation in the lower temperature zone (370390 C). This observation may possibly be explained as that in the lower temperature zone; the conversion is mainly catalytic (associated with the removal of heteroatoms), whereas in the higher range, it is both catalytic and thermal. Despite the fact that operating the hydrotreater at high temperatures would increase the yield of the middle distillate, it is not favorable because this will surely accelerate the deactivation rate of the catalyst. Refiners always try to balance between conversion and the rate of catalyst deactivation to achieve the optimum revenue. 5.2. Effect of the LHSV. Decreasing the LHSV increases conversion of the residual fraction to lighter fractions, because more residence time is allowed in the reactor. The results showed that the effect of LHSV is limited at a low reaction temperature (370 C) and quite apparent at higher temperatures (390 and 410 C). The highest yield of the distillate fraction has been achieved at the highest reaction severity (Table 3). However, the generation of the gas oil fraction is poorly influenced by LHSV. This observation can be either related to the fact that gas oil is Table 4. Model Parameters of the Continuous Lumping Hydrocracking Kinetic Model catalyst
model
type
parameter
HDM
HDS
HDN
370 C
390 C
consecutively converted to distillate or to the higher conversion selectivity from the residue fraction to distillate. In fact, the balance between support acidity and hydrogenation activity in the catalyst influences product selectivity, and a high concentration of molybdenum active sites allowed the catalyst to be more selective toward distillate production. The effect of LHSV at elevated temperatures was noticeable for all catalyst types, but it was more evident for HDS and HDN. These results indicate that the increase in the conversion rate with the contact time at elevated temperatures is ruled by not only the thermal reaction but also the catalytic reaction to a large extent. 5.3. Effect of the Catalyst Type. Hydrotreating catalysts are the backbone of hydroprocessing reactions. The characteristics of each catalyst in terms of shape, size, chemical composition, and pore structure mainly depend upon the type of reactions that it should perform. During the removal of sulfur, nickel, vanadium, and nitrogen atoms through hydrotreating operations, some conversion in the boiling range of the feedstock takes place. This conversion occurs because the molecules that contain the heteroatoms are broken down to lighter fractions, thus shifting the molecular weight and the boiling range of the feedstock. This partial conversion mainly depends upon the severity of the process and the acidity of the catalyst, which promote the
410 C
R
0.17574
0.39099
0.67758
a0
8.57152
12.1660
7.93510
a1
1.27207
1.56998
1.89825
kmax
6.27824
2.56788
1.34991
δ
7.110 104
6.987 104
1.757 103
R
1.15189
0.79752
0.89111
a0 a1
4.24511 2.44525
7.15175 2.09843
5.63812 1.92192
kmax
0.21231
0.89075
0.96206
δ
4.477 103
1.793 103
3.374 103
R
0.86510
0.84412
0.73726
a0
9.09330
6.45107
7.42454
a1
2.29719
1.88166
1.72660
kmax
0.58493
0.72199
1.44046
δ
1.100 103
2.081 103
2.442 103
Figure 7. Comparing experimental and predicted concentration profiles for the HDS catalyst at 410 C.
Table 5. Modeling and Testing Performance of the Continuous Lumping Modelsa R2 modeling catalyst type
HDM
temperature (C)
0.5 h
2.0 h
370
0.9985
390
0.9982
410 370 HDS
HDN
390
2
0.9995 0.9997
testing
modeling
testing
1.0 h
overall
0.5 h
2.0 h
0.9977
0.9915
0.9958
1.47 102
1.88 102
3.60 102
2.49 102
0.9986
0.9915
0.9961
1.65 102
1.45 102
3.60 102
2.44 102
2
1.27 10
2
4.44 10
2
2.85 102
1.28 10
2
2.46 10
2
1.66 102
1.77 10
2
2
2.46 102
0.9990 0.9988 0.9980
0.9874 0.9957 0.9903
0.9948 0.9980 0.9960
1.74 10
3
8.07 10
3
6.74 10
1.0 h
overall
3.82 10
410 370
0.9991 0.9996
0.9980 0.9989
0.9912 0.9962
0.9962 0.9982
1.14 10 8.01 103
1.73 10 1.30 102
3.63 10 2.39 102
2.41 102 1.63 102
390
0.9983
0.9985
0.9882
0.9950
1.57 102
1.54 102
4.29 102
2.78 102
2
2
2
4.05 102
410 a
0.9980
RMSE
0.9964
0.9988
0.9728
0.9894
2.32 10
2
2
2
1.36 10
6.48 10
0.5
R is defined by eq 16. RMSE = (SSE/N) , where N is the number of data points and SSE is defined by eq 17. 1945
dx.doi.org/10.1021/ef200153p |Energy Fuels 2011, 25, 1939–1949
Energy & Fuels
ARTICLE
cracking reaction. The difficulty of extracting heteroatoms varies significantly. The most difficult to be extracted is nitrogen. Therefore, the HDN catalyst has higher acidity and more active metals than HDS and HDM. The acidity allows the HDN catalyst to perform more catalytic hydrocracking on the molecules, while the active metals perform the extraction of nitrogen from the cracked molecules. The extraction of sulfur is relatively easier than nitrogen but more difficult compared to metals. Therefore, it is expected to have more cracking in HDS and HDN reactions than Table 6. Prediction Performance of the Continuous Lumping Modelsa kerosene
HGO
catalyst type temperature (C) 0.5 h 1.0 h 2.0 h 0.5 h 1.0 h 2.0 h
HDM
370
5.33
390
0.53 35.03 44.2
410 370 HDS
HDN a
390
1.82 38.4
2.38 15.27 18.45
410
10.9
29.4
370
19.8
390 410
24.82 19.04 11.75 30.6
0.21
8.65 5.93
4.76 11.18 7.62
3.32 4.04 11.22 1.85
1.35 11.34 0.98
5.38 13.1
Percent relative errors are reported.
3.63
3.96 0.86 11.2
5.02 0.91 8.11 1.22
0.24
4.87 1.36
6.76 2.91
7.75 1.39
7.81 4.83 9.46 0.48 5.9 5.44 13.18 0.36
the HDM reaction. The HDM catalyst has the lowest acidity and lowest active metals. Therefore, it is expected to have less catalytic hydrocracking in the HDM catalyst. This conclusion is supported by the experimental results, which showed that the HDN and HDS conversions are comparable, while the gap is wider for HDM, as illustrated in Figure 6. The cracking behavior of the three catalyst types may be explained by considering the diffusion mechanism of large hydrocarbon molecules into catalyst channels and pores.35 At low reaction severities, cracking of molecules is mainly controlled by the catalytic effect because the thermal effect is normally very limited. In addition, for small-pore-size catalysts, such as HDS and HDN, the effect of catalytic hydrocracking is also limited because the diffusion of the large asphaltene molecules into the catalyst channels is highly resistant. In these circumstances, the HDM catalyst is favored because it has large pores that allow for easier diffusion of asphaltene molecules. Consequently, asphaltene molecules would diffuse inside the pores and crack. It is important to note here that the pore distribution of the HDM catalyst varies in size to control the rate of catalytic cracking and to allow some thermal cracking to occur. This is a very important aspect in the design of the HDM catalyst because an excessive rate of catalytic cracking can rapidly drop the catalyst activity. The extent of both thermal and catalytic hydrocracking starts to increase as the temperature increases. At higher temperatures, the impact of the limited diffusion rate starts to diminish as the thermal cracking breaks down the asphaltene to smaller molecules, which can diffuse into small pore sizes. This clearly explains the
Figure 8. Species reactivity as a function of normalized TBP (θ), for different catalyst types and operating temperatures. 1946
dx.doi.org/10.1021/ef200153p |Energy Fuels 2011, 25, 1939–1949
Energy & Fuels
ARTICLE
Figure 9. Normalized species reactivity as a function of normalized TBP (θ), for different catalyst types and operating temperatures.
conversion data listed in Table 3. At low temperatures, the HDM catalyst shows superior behavior in cracking heavy molecules because of its large pore size that exerts less resistance against diffusing molecules. As reported by Calemma et al.,36 the balance between the acidity of the support and the hydrogenation activity of the metal is of primary importance in determining the selectivity of cracking products. Panariti et al.37 also indicated that the distillate yield is significantly affected by the concentration of acidic sites. Table 3 illustrates that the yield of the distillate fraction is the highest for the HDN catalyst, which has the highest acidity, followed by HDS and HDM catalysts.
6. CONTINUOUS LUMPING RESULTS Experimental data were used to determine the model parameters for the continuous lumping hydrocracking kinetic model described above. The five model parameters are R, a0, a1, kmax, and δ. Nonlinear optimization was used in estimating the values of these parameters, for the three types of catalyst (HDM, HDS, and HDN) and at three operating temperatures (370, 390, and 410 C). Yield experimental data for hydrocracking AR feedstock is available at three residence times (1/LHSVs); 0.5, 1.0, and 2.0 h. The strategy followed in developing the continuous lumping models is to estimate the model parameters using TBP distribution curves for cracking occurring at 0.5 and 2.0 h. The third data set, 1.0 h, has been kept for testing the developed models.
Moreover, the prediction performance of the continuous models will be validated further using the testing lumps [kerosene and heavy gas oil (HGO)]. Values of the five model parameters are listed in Table 4 for different catalyst types and operating conditions. Data fitting errors are listed in Table 5. Both the coefficient of determination (R2) and the root-mean-squared errors (RMSEs) demonstrate that the proposed continuous lumping model fitted the experimental data perfectly. The average overall R2 is 0.9955 (average RMSE is 2.53 102), which indicates excellent fits and very small errors in correlating the experimental data. Table 5 also compares the modeling errors with the testing errors. The average testing R2 is 0.9894, while the average modeling R2 is 0.9985. Such results illustrate the effectiveness of the continuous model in simulating the hydrocracking reactions. In addition, the outstanding performance in predicting the concentration profiles at 1.0 h residence time proves that the model has succeeded in capturing the effect of the residence time on the extent of the cracking reactions. The main advantage of the continuous lumping model is that it provides a continuous description of the species concentrations with respect to the normalized TBP, θ, and as a function of the reaction time. The effectiveness of the developed continuous models in simulating the concentration profiles is demonstrated in Figure 7, for the HDS catalyst at T = 410 C. The plots illustrate clearly that the proposed continuous lumping model succeeded in predicting the species concentration profiles as a function of the normalized TBP, θ. The model also showed 1947
dx.doi.org/10.1021/ef200153p |Energy Fuels 2011, 25, 1939–1949
Energy & Fuels excellent simulation capabilities for the effect of the residence time on the cracking yields. The prediction performance of the developed continuous lumping models will be further validated against the kerosene (205260 C) and HGO (315425 C) testing lumps, representing the light and heavy cuts, respectively. Percent relative errors resulting from comparing the predicted and experimental weight fractions of the testing lumps are listed in Table 6. In general, the reported errors indicate that the developed continuous lumping models predicted the concentrations of the testing lumps with reasonable accuracy. The models demonstrated much better performance in predicting the concentrations of the heavy testing lump. Average errors for kerosene and HGO are 14 and 4.75%, respectively. The maximum error for the HGO lump is 13.18%, while the maximum error for the kerosene lump is as high as 44.2%. The main reason for the relatively higher error values in predicting the concentrations of the light lump is that the weight fractions for this lump are quite small. In fact, the kerosene lumps starts from zero concentration (at 205 C) for most of the data sets. The HDM model predicted the concentration of the kerosene lump at 390 C and 2.0 h as 0.024 wt %, while the experimental concentration is 0.0166 wt %. Such prediction accuracy may be regarded as reasonably acceptable, while the same resulted in 44.2% error. The significance of the parameters used in developing the continuous lumping models will be studied using R and kmax, which describe the functional relation between the species reactivity (rate constants) and the normalized temperature, θ, as represented by eq 7. The reactivity profiles were found to be dependent upon the catalyst type and operating temperature. The HDM catalyst showed different behavior, as illustrated Figure 8. In fact, the HDM catalyst showed the highest kmax value at 370 C (Table 4). This means that, for this type of catalyst, cracking of the heavy species is much faster (easier) than the light species. As shown by Figure 8a, this difference in reactivity becomes less with an increasing temperature. The discrepancies of the HDM reactivity profiles compared to those of HDS and HDN are more evident in the plots shown in Figure 9. These plots represent the normalized reactivity of the species (k/kmax) versus the normalized TBP (θ). Two main observations may be derived from these plots: (i) For the HDM catalyst, the normalized reactivity increases with an increasing reaction temperature. However, for the HDS and HDN catalysts, the normalized reactivity decreases with an increasing temperature. (ii) The HDM catalyst shows more sensitivity to the operating temperature. The plots illustrate clearly that the highest temperature effect is shown for the HDM catalyst (Figure 9a), followed by the HDS catalyst (Figure 9b), while the HDN catalyst exhibits the least effect. These observations support the fact that thermal cracking dominates in the case of the HDM catalyst, while for the HDS and HDN catalysts, catalytic cracking is the main driving force.
7. CONCLUSION A kinetic model has been developed to describe the hydrocracking reactions associated with the hydrotreatment of AR. The developed kinetic model is based on the continuous lumping approach, which provides a continuous description of the species concentrations with respect to the normalized TBP. The continuous lumping model is very effective in simulating the hydrocracking reactions of heavy petroleum fractions and the yield
ARTICLE
concentration profiles. The effectiveness of the model was demonstrated for three hydrotreating catalysts at different operating severities. The prediction performance of the developed continuous lumping models was validated against experimental testing lumps that represent the light and heavy cuts. Percent relative errors resulting from comparing the predicted and experimental weight fractions of the testing lumps indicate reasonable accuracy. Reactivity profiles support the fact that thermal cracking dominates in the case of the HDM catalyst, while catalytic cracking is the main driving force in the case of HDS and HDN catalysts.
’ AUTHOR INFORMATION Corresponding Author
*Telephone: þ965-2498-5778. Fax: þ965-2484-6147. E-mail:
[email protected].
’ NOMENCLATURE A, B, and C = model parameters a0 = model parameter a1 = model parameter ci(t) = concentration of the i-type species C(t) = total reactant concentration cif = feed concentration c(k,t) = continuous concentration function D(k) = species-type distribution function kij = rate constant or reactivity of the i-type species cracked to j-type species (h1) Ki = sum of rate constants for the parallel reactions (h1) kmax = reactivity (rate constant) of the species with the highest TBP, i.e., for θ = 1 N = total number of species (pseudo-components) in the mixture (N f ¥) Nl = number of discrete lumps p(k,K) = yield distribution function S0 = model parameter TBP = true boiling point temperature (K) TBPmin = lowest boiling point temperature of the reaction mixtures (K) TBPmax = highest boiling point temperature of the reaction mixtures (K) R = model parameter δ = model parameter γij = stoichiometric coefficients θ = normalized TBP temperature ’ REFERENCES (1) Abbas, A. K.; Chaudhuri, S. R.; Bhattacharya, S. Performance optimization of atmospheric residue desulfurization units. Proceedings of the 8th Annual SaudiJapanese Symposium; Dhahran, Saudi Arabia, Nov 2930, 1998; pp 175189. (2) Speight, J. G. The Desulfurization of Heavy Oils and Residues; Marcel Dekker: New York, 2000. (3) Takahashi, T.; Higashi, H.; Kai, T. Development of a new hydrodemetallization catalyst for deep desulfurization of atmospheric residue and the effect of reaction temperature on catalyst deactivation. Catal. Today 2005, 104 (1), 76–85. (4) Al-Adwani, H. A. H.; Lababidi, H. M. S.; Alatiqi, I. M.; Al-Dafferi, F. S. Optimization study of residuum hydrotreating processes. Can. J. Chem. Eng. 2005, 83 (2), 281–290. 1948
dx.doi.org/10.1021/ef200153p |Energy Fuels 2011, 25, 1939–1949
Energy & Fuels (5) Marafi, A.; Maruyama, F.; Stanislaus, A.; Kam, E. Multicatalyst system testing methodology for upgrading residual oils. Ind. Eng. Chem. Res. 2008, 47 (3), 724–741. (6) Kodama, S.; Nitta, H.; Takatsuka, T.; Yokoyama, T. Simulation of residue hydrodesulfurization reaction based on catalyst deactivation model. J. Jpn. Pet. Inst. 1980, 23 (5), 310–320. (7) Lababidi, H. M. S.; Shaban, H. I.; Al-Radwan, S.; Alper, E. Simulation of an atmospheric residue desulfurization unit by quasisteady state modeling. Chem. Eng. Technol. 1998, 21 (2), 193–200. (8) Oyekunle, L. O.; Kalejaiye, B. O. Kinetic modeling of hydrodesulphurization of residual oils. I. Power law model. Pet. Sci. Technol. 2003, 21 (910), 1475–1488. (9) Kam, E.; Al-Shamali, M.; Juraidan, M.; Qabazard, H. A hydroprocessing multicatalyst deactivation and reactor performance modelpilot-plant life test applications. Energy Fuels 2005, 19 (3), 753–764. (10) Alvarez, A.; Ancheyta, J. Modeling residue hydroprocessing in a multi-fixed-bed reactor system. Appl. Catal., A 2008, 351 (2), 148–158. (11) Ma, C.; Weng, H. Study on residual oil HDS process with mechanism model and ANN model. China Pet. Process. Petrochem. Technol. 2009, 1, 39–44. (12) AlHumaidan, F. S.; Al-Adwani, H. A. H.; Lababidi, H. M. S. Hydrocracking of atmospheric residue feedstock in hydrotreating processes. Kuwait J. Sci. Eng. 2010, 37 (1B), 129–159. (13) Baltanas, M. A.; Froment, G. F. Computer generation of reaction networks and calculation of product distribution in the hydroisomerisation and hydrocracking of paraffins on Pt-containing bifunctional catalysts. Comput. Chem. Eng. 1985, 9 (1), 71–77. (14) Laxminarasimhan, C. S.; Verma, R. P.; Ramachandran, P. A. Continuous lumping model for simulation of hydrocracking. AIChE J. 1996, 42 (9), 2645–2653. (15) Ancheyta, J.; Sanchez, S.; Rodríguez, M. A. Kinetic modeling of hydrocracking of heavy oil fractions: A review. Catal. Today 2005, 109 (14), 76–92. (16) Krishna, R.; Saxena, A. Use of axial dispersion model for kinetic description of hydrocracking. Chem. Eng. Sci. 1989, 44 (703), 177–206. (17) Balasubramanian, P.; Pushpavanam, S. Model discrimination in hydrocracking of vacuum gas oil using discrete lumped kinetics. Fuel 2008, 87 (89), 1660–1672. (18) Sadighi, S.; Ahmad, A.; Rashidzadeh, M. 4-Lump kinetic model for vacuum gas oil hydrocracker involving hydrogen consumption. Korean J. Chem. Eng. 2010, 27 (4), 1099–1108. (19) Aoyagi, K.; McCaffrey, W.; Gray, M. Kinetics of hydrocracking and hydrotreating of coker and oil sands gas oils. Pet. Sci. Technol. 2003, 21 (5), 997–1015. (20) AlHumaidan, F. S. Modeling hydrocracking kinetics of atmospheric residue by discrete and continuous lumping. M.Sc. Thesis, College of Graduate Studies, Kuwait University, Kuwait City, Kuwait, 2004. (21) Sanchez, S.; Ancheyta, J. Effect of pressure on the kinetics of moderate hydrocracking of Maya crude oil. Energy Fuels 2007, 21 (2), 653–661. (22) Krishna, R.; Balaramanian, P. Analytical solution for discrete lumped kinetic equations in hydrocracking of heavier petroleum fractions. Ind. Eng. Chem. Res. 2009, 48 (14), 6608–6617. (23) Chou, M. Y.; Ho, T. C. Continuum theory for lumping of catalytic cracking process. AIChE J. 1992, 34 (9), 1519–1527. (24) Aris, R. Reaction in continuous mixtures. AIChE J. 1989, 35 (4), 539–548. (25) Weekman, V. W. Lumps, models and kinetic in practice. Chem. Eng. Prog. Monogr. Ser. 1979, 75 (11), 3. (26) Basak, K.; Sau, M.; Manna, U.; Verma, R. Industrial hydrocracker model based on novel continuum lumping approach for optimization in petroleum refinery. Catal. Today 2004, 98 (12), 253–264. (27) Lababidi, H. M. S.; AlHumaidan, F. S.; Al-Adwani, H. A. H. Modeling hydrocracking kinetics of atmospheric residue feedstock. ACS National Meeting Book of Abstracts, 2005.
ARTICLE
(28) Elizalde, I.; Rodríguez, M. A.; Ancheyta, J. Application of continuous kinetic lumping modeling to moderate hydrocracking of heavy oil. Appl. Catal., A 2009, 365 (2), 237–242. (29) Quann, R. J.; Jaffe, S. B. Building useful models of complex reaction systems in petroleum refining. Chem. Eng. Sci. 1996, 51 (10), 1615–1631. (30) Martens, G. G.; Marin, G. B. Kinetics for hydrocracking based on structural classes: Model development and application. AIChE J. 2001, 47 (7), 1607–1621. (31) Froment, G. F. Single effect kinetic modeling of complex catalytic processes. Catal. Rev.Sci. Eng. 2005, 47, 83–124. (32) Absi-Halabi, M.; Stanislaus, A.; Qamara, A. Effect of presulfiding on the activity and the deactivation of hydrotreating catalysts in processing Kuwait vacuum residue. Stud. Surf. Sci. Catal. 1995, 243–251. (33) Marafi, A.; Stansilaus, A.; Furimsky, E. Kinetics and modeling of petroleum residues hydroprocessing. Catal. Rev.Sci. Eng. 2010, 52, 204–324. (34) Squires, R. G.; Reklaitis, G. V.; Yeh, N. C.; Mosby, J. F.; Karimi, I. A.; Anderson, P. K. Purdue-industry computer simulation modules. Chem. Eng. Educ. 1991, 98–101. (35) Nomura, M.; Akagi, K.; Murata, S.; Matsui, H. Hydrocracking of polycyclic aromatic compounds using zeolite catalysts. Catal. Today 1996, 29 (14), 235–240. (36) Calemma, V.; Peratello, S.; Perego, C. Hydroisomerization and hydrocracking of long chain n-alkanes on Pt/amorphous SiO2Al2O3 catalyst. Appl. Catal., A 2000, 190 (12), 207–218. (37) Panariti, N.; Del Bianco, A.; Del Piero, G.; Marchionna, M. Petroleum residue upgrading with dispersed catalysts: Part 1. Catalyst activity and selectivity. Appl. Catal., A 2000, 204 (2), 203–213.
1949
dx.doi.org/10.1021/ef200153p |Energy Fuels 2011, 25, 1939–1949
