Mechanistic Kinetic Modeling of the Hydrocracking of Complex

Jul 31, 2007 - Vacuum Gas Oil Hydrocracking on NiMo/USY Zeolite Catalysts. .... catalytic hydroprocessing of middle-temperature coal tar to clean liqu...
6 downloads 0 Views 588KB Size
Ind. Eng. Chem. Res. 2007, 46, 5881-5897

5881

Mechanistic Kinetic Modeling of the Hydrocracking of Complex Feedstocks, such as Vacuum Gas Oils Hans Kumar and Gilbert F. Froment* Artie McFerrin Department of Chemical Engineering, Texas A&M UniVersity, College Station, Texas 77843-3122

A mechanistic kinetic model has been developed for the hydrocracking of complex feedstocks such as vacuum gas oil (VGO), based on an exhaustive computer-generated reaction network of elementary steps. The model utilizes a detailed composition of VGO, characterized by 16 different molecular classes up to C40. These classes are divided into 45 subclasses by distinguishing the isomers of a class according to the number of methyl branches, leading to 1266 groups of isomers/pure components. The frequency factors of the rate coefficients for the acid site steps are modeled using the single-event concept, and the activation energies are modeled based on the nature of the reactant and product carbenium ions. The saturation of polyaromatics is modeled according to the sequential hydrogenation of aromatic rings on the metal sites. The competitive chemisorption of aromatics and naphthenes is taken into account, based on the number of aromatic/naphthenic rings on the metal sites, and based on the corresponding gas-phase proton affinities on the acid sites. The model contains 33 feedstock and temperature-independent parameters that have been estimated from experimental data on VGO hydrocracking. The kinetic model is inserted into an adiabatic multibed trickleflow reactor model. The vapor-liquid equilibrium coefficients in this model are calculated using the PengRobinson equation of state. The model has been used to study the effect of operating conditions on the product profiles along the reactor. An analysis of the distribution of isomers of a class among its different subclasses shows the increase in VGO conversion when the content of isomers with a higher degree of branching is increased in the feedstock. This indicates the necessity of dividing a class into four subclasses. 1. Introduction Hydrocracking is used to upgrade the heavier fractions obtained from crude oil distillation, such as vacuum gas oils (VGOs), but it is also practiced on residues and other products obtained from various refining processes, such as coker gas oil, deasphalted oil, and fluidized catalytic cracking (FCC) cycle oils. Very heavy hydrocarbons, such as those extracted from tar sands and shale, can also be upgraded by hydrocracking. A typical hydrocracking feedstock consists of paraffinic, naphthenic, aromatic, and naphtheno-aromatic species, along with heteroatom impurities such as sulfur, nitrogen, and oxygen compounds. A significant amount of metals may also be present. Typically, the products obtained from fractionation of the hydrocracker effluent are defined based on their boiling-point range, such as light ends (C4-), light naphtha (C5-80 °C), heavy naphtha (80-150 °C), jet fuel/kerosene (150-290 °C), diesel fuel (290-370 °C), and the unconverted fractionator bottoms (370 °C +). Recently, environmental concerns have set the agenda for the development of hydrocracking and hydrotreatment processes, producing high-quality fuels with low aromatics, sulfur, and nitrogen content. Significant technological developments have been made, such as single-stage once through, single-stage with recycle, two-stage process, two-stage with recycle, etc., for the purposes of higher operational flexibility and maximization of the financial returns. At the same time, various catalysts have been commercialized for processing the feedstocks of vastly different compositions. The development of a fundamental kinetic model for this type of complex conversion processes has been rather limited, with most of the kinetic models still * To whom correspondence should be addressed. Tel.: 979-845-0406. Fax: 979-845-6446. E-mail address: [email protected].

based on simplified lumped schemes. The design and optimization of hydrocracking units require a detailed kinetic model that takes into account the complexity of the feedstock, while satisfying the rules of the underlying carbenium ion chemistry. A comprehensive process model with accurate hydrocracking kinetics at the elementary step level can be used to reduce expensive experimentation in pilot plants. This type of model can be used to predict optimum operating conditions and detailed product distribution for a wide range of VGO feedstocks in existing or prospective units. 2. Literature Review To study the conversion of complex feedstocks through processes such as hydrocracking and catalytic cracking, most efforts have focused on the development of lumped kinetic models in which the feedstock is divided into several lumps, based on the boiling-point range. A simplified reaction network between these lumps is set up and the rate coefficients for the global conversion of lumps are estimated from experimental data. For example, the kinetic model of Weekman and Nace1 for FCC assumes that the feedstock is converted to gasoline and a remaining fraction by only two reactions. A more-detailed lumped model for FCC comprising 10 lumps was developed by Jacob et al.2 Stangeland3 considered the feedstock as a series of 50 °F boiling range cuts, assuming that each heavier cut hydrocracks via a first-order reaction to form a series of lighter cuts. To achieve higher accuracy in the predicted product yields, more and more lumps were introduced, leading to an increasing number of parameters in the models. Nevertheless, the major fundamental limitation of the lumped kinetic models persisted: dependency of the kinetic parameters on the composition of the feedstock. As a consequence, the kinetic model must be refitted and new sets of parameters must be estimated for every different feedstock.

10.1021/ie0704290 CCC: $37.00 © 2007 American Chemical Society Published on Web 07/31/2007

5882

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007

Another approach for modeling the conversion of such feedstocks is based on continuum lumping in which the reaction mixture is considered to be a continuous mixture, with respect to the feed properties, such as boiling point and molecular weight. Similar to discrete lumping, this approach is also unable to capture the fundamental chemistry of the process, providing thrust to the development of mechanistic models. Several approaches have been developed to build mechanistic models for catalytic conversion of complex feedstocks. Liguras and Allen4,5 represented the feedstock using hundreds of pseudocomponents and developed a kinetic model for catalytic cracking, based on the nature of carbon centers in the pseudocomponents. Quann and Jaffe6 constructed the molecules by incrementing the structural units encountered in hydrocarbon molecules and expressed them using a vector notation that was used to generate the reaction network. The reaction network is based on various reaction rules, selected according to the chemistry of the process. This approach is called the structural oriented lumping (SOL). Froment and co-workers7,8 developed a mechanistic kinetic modeling approach starting from the elementary steps of carbenium ion chemistry. This approach was named “singleevent” approach. Baltanas et al.7 generated a network of elementary steps involving carbenium ions, using a computer algorithm based on the approach devised by Clymans et al.9 and Hillewaert et al.10 Vynckier et al.8 applied the single-event approach to complex feedstocks by introducing the concept of lumping coefficients. Feng et al.11 used the single-event concept to model the catalytic cracking of paraffins on a RE-Y zeolite catalyst. Martens et al.13 applied single-event kinetics for the hydrocracking of C8-C12 paraffins on Pt/USY zeolites. Kumar and Froment14 developed a single-event kinetic model combined with the Evans-Polanyi relationship for the hydrocracking of long-chain paraffins encountered in Fischer-Tropsch waxes and elucidated the effect of the relative metal to acid strength of the catalyst on the selectivities of the hydrocracking products. 3. Hydrocracking Chemistry Hydrocracking is performed on dual-function catalysts that have a cracking function and a hydrogenation-dehydrogenation function. On a zeolite-based catalyst, the feed molecules are sorbed into the pores of the zeolites from the surrounding fluid phase, which can be gas or liquid, depending on the reactor operating conditions. Different reaction pathways are available for the paraffinic, naphthenic, and aromatic species. The sorbed paraffinic and naphthenic species are chemisorbed on the metal sites and are dehydrogenated into the corresponding olefinic species that migrate from the metal sites to the nearby Brønsted acidic sites and are protonated into paraffinic and naphthenic carbenium ions. The olefinic species sorbed in the zeolite pores are assumed to be in pseudo-equilibrium with the chemisorbed olefins on the metal sites, as well as with the surface carbenium ions on the acid sites. These steps are modeled using the Langmuir adsorption isotherm. The paraffinic carbenium ions undergo hydride shift (HS), methyl shift (MS), protonated cyclopropane (PCP), and β-scission steps, as shown in Figure 1. Cyclization of paraffinic species into ring species is also possible. The naphthenic carbenium ions can undergo all the elementary steps mentioned previously in their alkyl side chain or in the ring. The PCP steps involving the ring carbons result in an expansion or contraction of the ring size. Depending on the location of the positive charge, the naphthenic carbenium ions can undergo exocyclic β-scission producing a naphthenic olefin

Figure 1. Steps involved in the hydrocracking of paraffins.

and a long alkyl chain as a paraffinic carbenium ion, or can transform through endocyclic β-scission into an unsaturated carbenium ion. The various steps in the hydrocracking of mononaphthenes are summarized in Figure 2. Some of the possibilities that can occur in the hydrocracking of tetranaphthenes are shown in Figure 3. For simplicity, the positive charge has been omitted from the species. The hydroconversion of alkyl aromatics proceeds either through the acid-site elementary steps in their side chain or through the saturation of the aromatic ring and subsequent isomerization and cracking in the resulting naphthenic ring. Dealkylation of the benzyl carbenium ion formed by the direct protonation of the aromatic ring on the acid sites is also an important step in the hydrocracking of aromatics (see Figure 4). Hydrocracking of polynuclear aromatics is quite complex and involves sequential hydrogenation of the aromatic rings, along with parallel or subsequent isomerization, cracking, and ring opening steps. 4. Reaction Network Generation Considering the large number of reaction pathways in the hydrocracking of heavy feedstocks when expressed in terms of elementary steps, a reaction network for the paraffins, naphthenes, and aromatics has been generated up to C40, using Boolean relation matrices and characterization vectors. The methodology and procedure for generating the reaction network has been described in the literature.9,15,16 In the present work, the program execution was significantly accelerated using dynamic memory allocation via linked lists for storing and searching the intermediate olefinic and ionic species. The rules for deriving the number of single events of the PCP and cracking steps, based on the symmetry of the species, have also been improved. The total number of carbenium ions involved in VGO hydrocracking up to C40 is on the order of 18 million. The total number of steps and the number of rate-determining steps are of the order of 126 and 15 million, respectively. The steep increase in the number of steps with carbon number is shown in Figure 5 for several hydrocarbon classes. 5. Molecular Characterization of Vacuum Gas Oil (VGO) The development of a mechanistic kinetic model requires knowledge of the composition of the feedstock at the molecular level. Significant developments have been made recently for a detailed characterization of heavy mixtures with the application of modern analytical techniques such as gas chromatography (GC), supercritical fluid chromatography (SFC), high-performance liquid chromatography (HPLC), field-ionization mass spectrometry (FIMS), low-voltage high-resolution mass spec-

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007 5883

Figure 3. Reaction network of multiring naphthenes. (Shown in a simplified way in the form of molecular reactions instead of elementary steps.)

Figure 2. Elementary steps of monoring naphthenes.

troscopy (LV-HRMS), nuclear magnetic resonance (NMR) spectroscopy, etc. Boduszynski,17,18 for example, coupled HPLC and FIMS in the analysis of several petroleum crude oils. They fractionated the petroleum into 10 boiling range cuts, each of which was subjected to three HPLC extractive separations to determine the various subfractions (the saturates, monoaromatics, diaromatics, triaromatics, tetraaromatics, and pentaaromatics). Each of these subfractions was analyzed in detail by FIMS to generate the molecular weight distribution or carbon number distribution of different homologous series. This type of analysis gives the composition of a heavy petroleum fraction in a matrix form in which the rows represent different molecular classes and the columns contain the carbon number distribution for each class. Compounds in each class have a distinct hydrogen deficiency number Z (as in CnH2n+Z), where Z ) -2(R + DB -1) (here, R is the number of rings and DB is the number of double bonds in the given class). One such Z-CN representation

(where CN is the carbon number) to characterize a VGO feedstock is shown in Table 1 in which VGO is divided into 16 molecular classes up to C40. Martens and Marin19 developed a hydrocracking model based on this level of analysis, but without considering the hydrogenation of aromatics. Unfortunately, the analysis of VGO, even at such a level (the Z-CN matrix form) does not allow the development of a mechanistic kinetic model in which the model parameters could be claimed to be independent of the feedstock composition. Indeed, the isomers of a particular carbon number belonging to a class and with a given number of methyl branches reach thermodynamic equilibrium, because of fast hydride- and methyl-shift elementary steps. In the VGO characterization shown in Table 1, however, all the isomers of a given class are allocated to one group of isomers (GOI), irrespective of the number of methyl branches. Such a characterization of VGO, therefore, does not allow the calculation of the rate of conversion of less-branched isomers into more-branched isomers (and vice versa) that occurs through PCP steps. In other words, according to the scheme given in Table 1, all the isomers in a GOI with different number of methyl branches would have to be assumed at pseudo-equilibrium along the entire length of the reactor, which is far from reality. Vansina et al.20 compared the experimental product distribution with the thermodynamically calculated values in the hydrocracking of n-octane and showed that the fraction of multibranched isomers is always smaller than the thermodynamic equilibrium value. This is especially pronounced for the tribranched isomers. Similarly, Schulz and Weitkamp21 reported that, in the hydrocracking of n-dodecane on different catalysts, the ratios of the monomethyl isomers to the corresponding n-paraffins are greater than those calculated from thermodynamic equilibrium, whereas the amounts of dimethyl isomers are much smaller. This pattern has also been assumed to be applicable to the ring-containing structures, given the presence of long paraffinic side chains in these species that are undergoing a change in branching through the same reaction mechanism (i.e., PCP). Because the rate of cracking through β-scission increases with the degree of branching, allocating

5884

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007

Figure 5. Total number of elementary steps in the network generated for paraffins, mononaphthenes, dinaphthenes, and monoaromatics.

Figure 4. Elementary steps and reactions of monoring aromatics.

all the isomers per carbon number of a given class in one GOI cannot give accurate values of the global rate of conversion. Therefore, in the current model, all the classes except normal paraffins and those with four rings (naphthenic or aromatic) structures have been further divided into subclasses, based on the number of methyl branches. Normal paraffins are obviously pure components. The four ring structures are not divided into subclasses, because of their relatively lower concentrations. According to the rules set for the reaction network generation,16 paraffinic species can have a maximum of three methyl branches. Therefore, the isoparaffin class has been divided into three subclasses (i.e., monobranched, dibranched, and tribranched). The species that contain naphthenic or aromatic rings are allowed to have only one long side chain and up to three methyl branches anywhere on the ring structure or the long side chain. Consequently, for a given carbon number of any ringcontaining class, the first subclass contains the species that have just the bare ring structure (i.e., the first member of the class, e.g., benzene in the monoaromatic class), or the species that have only the long side chain without any methyl branch. This subclass is referred to as unbranched throughout the remainder

of this article. The other three subclasses contain species with one, two, and three methyl branches, denoted by monobranched, dibranched, and tribranched subclasses, respectively. According to this scheme, the 16 molecular classes lead to 45 subclasses, resulting in the characterization of VGO in terms of 1266 pure components/GOIs, as shown in Table 2. It is very difficult to characterize a heavy feedstock such as VGO to such a detailed level. Mass spectrometers/high-resolution mass spectrometers cannot distinguish among the isomers with the same chemical formula. However, techniques such as 1H and 13C NMR can provide important information about the average population of carbons and hydrogens of different nature, i.e., primary, secondary, tertiary, aromatic, R-position to aromatic ring, etc., of a fraction separated from VGO by GC/HPLC/SFC. This information can be utilized to obtain an estimate of the distribution of isomers with different degrees of branching to convert the VGO composition from 16 classes to 45 subclasses. The development of more-advanced analytical and theoretical methods would certainly be useful for obtaining this type of information. In the current model, the distribution of monobranched, dibranched, and tribranched isomers in the isoparaffin class was chosen to be 65%, 25%, and 10%, respectively. In the classes with aromatic and naphthenic rings, a distribution of 20%, 50%, 20%, and 10% was selected for the respective unbranched, monobranched, dibranched, and tribranched isomers. The distribution of the degree of branching can be significantly different in the VGOs obtained from other processing units as compared to the straight run VGOs, and must be adjusted accordingly. Section 9.4 of the present paper reports on a sensitivity analysis of the hydrocracking product yields, with respect to the distribution of isomers in the feed. 6. Kinetic Model Development As shown previously, in the hydrocracking of complex feedstocks, the number of elementary steps and reactions that occur on the acidic sites and metal sites of the catalyst is extremely large; yet every elementary step and reaction has a finite contribution toward the final product distribution. For such feedstocks, conventional kinetic modeling approaches would result in a large number of reaction rate coefficients and it would be unrealistic to try to estimate their values from the experimental data. The solution of this problem requires going beyond

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007 5885 Table 1. Z-CN Matrix Representation of a Vacuum Gas Oil (VGO), in terms of 462 Components or Groups of Isomers (GOIs)a CN

NPAR

IPAR

MNA

DNA

TNA

3 4

1 2

39

5 6 7

3 4 5

40 41 42

76 77 78

8 9 10 11 12

6 7 8 9 10

43 44 45 46 47

79 80 81 82 83

112 113 114

13 14 15

11 12 13

48 49 50

84 85 86

115 116 117

143 144

16 17 18 19 20

14 15 16 17 18

51 52 53 54 55

87 88 89 90 91

118 119 120 121 122

145 146 147 148 149

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38

56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75

92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111

number of lumps

38

37

36

QNA

MAR

DAR

TAR

QAR

NMA

NDA

NTA

DNMA

DNDA

TNMA

Liquefied Petroleum Gas (LPG)

Light Naphtha (LNAP) 193 194 Heavy Naphtha (HNAP) 195 196 197 228 198 229 199 230

309 310 311

Kerosene (KERO) 200 231 201 232 259 202 233 260

312 313 314

340 341

390 391

286 287 288

315 316 317 318 319

342 343 344 345 346

367 368 369

392 393 394 395 396

417 418 419

440 441 442

123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142

Unconverted Vacuum Gas Oil (UNCONV) 150 173 208 239 266 289 151 174 209 240 267 290 152 175 210 241 268 291 153 176 211 242 269 292 154 177 212 243 270 293 155 178 213 244 271 294 156 179 214 245 272 295 157 180 215 246 273 296 158 181 216 247 274 297 159 182 217 248 275 298 160 183 218 249 276 299 161 184 219 250 277 300 162 185 220 251 278 301 163 186 221 252 279 302 164 187 222 253 280 303 165 188 223 254 281 304 166 189 224 255 282 305 167 190 225 256 283 306 168 191 226 257 284 307 169 192 227 258 285 308

320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339

347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366

370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389

397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416

420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439

443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462

31

27

31

27

23

27

23

23

170 171 172

23

Diesel (DIES) 203 234 204 235 205 236 206 237 207 238

35

31

261 262 263 264 265

27

23

a

Abbreviations for different classes: NPAR, normal paraffins; IPAR, isoparaffins; MNA, mononaphthenes; DNA, dinaphthenes; TNA, trinaphthenes; TETNA, tetranaphthenes; MAR, monoaromatics; DAR, diaromatics; TAR, triaromatics; TETAR, tetraaromatics; NMA, naphtheno-monoaromatics; NMA, naphtheno-diaromatics; NTA, naphtheno-triaromatics; DNMA, di-naphtheno monoaromatics; DNDA, di-naphtheno-diaromatics; and TNMA: tri-naphthenomonoaromatics.

the modeling of the rate equations by applying the modeling approach to the rate parameters also, starting from transitionstate theory. 6.1. Modeling of the Frequency Factors. Single-event kinetics has been introduced to model the frequency factors of the elementary steps that occur on acidic sites of the catalyst. In the single-event kinetics, the effect of molecular structure on the frequency factor of an elementary step is described with the help of transition-state theory and statistical thermodynamics, leading to the following equation for the rate coefficient of an elementary step:7

k)

( )( ) ( ) ( ) σRgl kBT ∆Sˆ °† -∆H°† exp exp R RT σ† h

(1)

gl

In the aforementioned equation, the effect of changes in symmetry in going from reactant to activated complex on the frequency factor of an elementary step has been factored out as

the ratio of the global symmetry numbers of the corresponding species. This ratio is defined as the number of single events (ne), so that the frequency factor for the elementary step can be written as

A ) n eA ˜

(2)

where A ˜ is the single-event frequency factor, expressed as

A ˜)

( ) ( ) kBT ∆Sˆ °† exp h R

(3)

The single-event frequency factor A ˜ is independent of the structures of the reactant and activated complex for all elementary steps that belong to the same type of step. As a result, the application of this concept requires only one independent parameter (namely, A ˜ ) to model the frequency factors of all elementary steps of that type. Consequently, for the six types of elementary steps that occur on the acidic sites (i.e., PCP,

5886

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007

Table 2. Z-CN Matrix Representation of a VGO, in terms of 1266 Components/GOIsa carbon number, CN

Subclass Subclass Subclass Subclass Subclass Subclass Subclass Subclass Subclass Subclass Subclass 1 2 3 4 5 6 7 8 ‚‚‚ ‚‚‚ ‚‚‚ 43 44 45 NPAR MBP DBP TBP MNA0 MNA1 MNA2 MNA3 ‚‚‚ ‚‚‚ ‚‚‚ DNMA3 DNDA TNMA

3 4 5 6 7 8 9 10 11 12 l 18 l l 38 39 40

1 2 3 4 5 6 7 8 9 10 l l l l 36 37 38

39 40 41 42 43 44 45 46 47 l l l l 73 74 75

76 77 78 79 80 81 82 l l l l 108 109 110

CN No. of lumps cumulative

NPAR 38 38

MBP 37 75

DBP 35 110

a

111 112 113 114 115 116 l l l l 142 143 144

145 146 147 148 149 150 151 152 l l l l 178 179 180

181 182 183 184 185 186 l l l l 212 213 214

215 216 217 218 219 l l l l 245 246 247

248 249 250 251 l l l l 277 278 279

TBP 34 144

MNA0 36 180

MNA1 34 214

MNA2 33 247

MNA3 32 279

‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ l l l l ‚‚‚ ‚‚‚ ‚‚‚

‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ l l l l ‚‚‚ ‚‚‚ ‚‚‚

‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ ‚‚‚ l l l l ‚‚‚ ‚‚‚ ‚‚‚

l l l l 1218 1219 1220

‚‚‚ ‚‚‚ ‚‚‚ DNMA3 ‚‚‚ ‚‚‚ ‚‚‚ 23 ‚‚‚ ‚‚‚ ‚‚‚ 1220

1221 l l 1241 1242 1243

1244 l l 1264 1265 1266

DNDA 23 1243

TNMA 23 1266

Numbers 0 to 3 in the second row represent the number of methyl branches in that subclass.

acyclic β-scission, endocyclic β-scission, exocyclic β-scission, dealkylation, and cyclization), only six single-event frequency factors are used as model parameters. The ne values are estimated for all elementary steps during the reaction network generation and stored with the respective elementary steps. 6.2. Modeling of the Activation Energies. The activation energies for the acid site elementary steps are modeled based on the nature of the reactant and the product carbenium ions. It is considered that the difference in the energy levels of the reactant and product carbenium ion are dependent only on their nature, i.e., secondary or tertiary (primary ions that are very unstable do not need to be considered). Therefore, four distinct activation energiessnamely, E(s;s), E(s;t), E(t;s) and E(t;t)sare required for the four subtypes of a given type of elementary step. A more-rigorous treatment for the modeling of the activation energies is presented by Kumar and Froment14 for the hydrocracking of pure paraffinic feedstocks. It is based on the application of the Evans-Polanyi relationship, which explicitly accounts for the differences in the energy levels of the reactants and the products starting from their heats of formation. A reduction in the total number of parameters required for modeling the activation energies is performed by observing the similarity in structural transformation in different types of β-scission elementary steps: acyclic, exocyclic, and endocyclic steps. Figure 2 shows that, in all the acyclic β-scission steps (i.e., when β-scission occurs in a paraffin or in the side chain of the ring species), the double bond would be formed in the paraffinic fragment or in the side chain of the naphthenic fragment. In all the exocyclic β-scission steps, however, the double bond would be formed in the naphthenic ring. The presence of the double bond in the naphthenic ring causes an additional ring strain. This leads to an average heat of reaction for exocyclic steps of any given subtype (for example, (s;s)), slightly higher than the average heat of reaction of the corresponding subtype of acyclic β-scission steps. Considering the application of the Evans-Polanyi relationship, the activation energies for the exocyclic steps would be higher than those of acyclic steps by an amount of R|∆Hexo-β - ∆Hacyc-β|, provided the intrinsic activation energy (E0) for these two types of steps is the same. This difference is not dependent on the nature of

the reactant and product carbenium ions. Therefore, the activation energies for the four subtypes of exocyclic steps were obtained from the activation energies of the corresponding acyclic steps via the introduction of only one additional parameter, ∆Eexo-β, i.e., Eexo-β(m;n) - Eacyc-β(m;n) (where m and n represent secondary or tertiary). A similar approach has been followed for the calculation of the activation energies of endocyclic β-scission steps. Only one parameter, ∆Eendo-β, was introduced to account for the difference in the activation energies of a given subtype of endocyclic steps from the corresponding subtype of acyclic steps, i.e., Eendo-β(m;n) - Eacyc-β(m;n). Cyclization steps are the reverse of endocyclic β-scission. In the absence of the equilibrium coefficients of the endocyclic steps at the surface of the catalyst, four different activation energies (Ecyc(s;s), Ecyc(s;t), Ecyc(t;s), and Ecyc(t;t)) are used for the four subtypes of cyclization steps. The carbenium ions involved in the dealkylation of aromatics are the ring-protonated aromatic species (the “benzenium” ions). Because of the strong delocalization of the positive charge on the aromatic ring, the energy levels of benzenium ions cannot be distinguished on the basis of their secondary or tertiary nature. This is why, depending on the nature of the paraffinic carbenium ions formed in the dealkylation steps, only two activation energies (Edealk(s) and Edealk(t)) are used for all the dealkylation steps. 6.3. Development of Rate Equations. Figure 6 depicts physical and chemical phenomena that occurs inside the catalyst particle during the hydrocracking of VGO. Paraffins, naphthenes, and aromatic (PNA) species are sorbed from the liquid phase into the zeolite pores. Depending on the difference in the reactive intermediate, the elementary steps and reactions that occur on the acid/metal sites of the catalyst are classified as follows (shown by different blocks in Figure 6): (a) PCP, acyclic β-scission, exocyclic β-scission, endocyclic β-scission: proceed through olefinic intermediates, which are protonated into carbenium ions for these elementary steps to occur. (b) Cyclization: proceeds through diolefinic intermediates, which are protonated into olefinic carbenium ions before the cyclization step.

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007 5887

Figure 6. Physical and chemical phenomena inside the catalyst particle.

(c) Dealkylation: for these elementary steps to occur, the aromatic species sorbed inside the zeolite pores are directly protonated on the ring. (d) Saturation of aromatics occurs on the metal sites of the catalyst itself. 6.3.1. PCP, Acyclic/Exo/Endo β-Scission Steps. The concentration of a paraffinic, naphthenic, aromatic, or naphthenoaromatic species Si sorbed from the liquid phase into the zeolite pores is obtained from the Langmuir sorption isotherm:

[Si] )

HSiCSliqi DL

(4)

species Si and the corresponding olefinic species SOij are the same, the concentration of the olefinic species inside the zeolite pores is given by

[SOij] )

∑i

KL,SiCSliqi

(5)

For these elementary steps, the species Si is chemisorbed on the metal sites of the catalyst as Smi and dehydrogenated into the corresponding olefinic species, SOmij. Naphthenic species can dehydrogenate either in the side chain or in the ring, whereas aromatic species can dehydrogenate in the side chain only.

Si a SOmij + Hliq 2

(6)

It is assumed that the catalyst has sufficient metal/acid activity to bring these hydrogenation/dehydrogenation steps to quasiequilibrium. The resulting olefinic species are desorbed from the metal sites. Assuming that the Henry coefficients for the

ij

DLCHliq2

(7)

Protonation/deprotonation steps have been shown to be potentially much faster than the other steps on the acidic sites and are assumed to be in quasi-equilibrium. The concentration of the surface carbenium ions, [SR+ik], is then obtained from

where

DL ) 1 +

liq HSiCSliqi KDH,(S iaSO )

[SR+ik]

)

Kpr(SO aSR+ )[SOij]Ct ij

1+

∑ij

ik

(8) Kpr(SO aSR+ )[SOij] ij

ik

In the hydrocracking of pure paraffins, the concentrations of the olefinic intermediates are several orders of magnitude smaller than those of the paraffinic species;12 therefore, the concentration of surface carbenium ions is negligible, with respect to the total acid site concentration. On the other hand, in VGO hydrocracking, the presence of significant amounts of aromatics and their preferential protonation on the acid sites leads to significant site coverage. The procedure for the estimation of fractional site coverage is described in section 6.5. Based on the concentration of the carbenium ion, [SR+ik], the rate of an elementary step SR+ik f SR+uV is given by

5888

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007

HSiCSliqi

liq rω(m;n) ) ne,ikuVk˜ω(m;n)KDH,(S K + )Ct iaSO ) pr(SOijaSR ij

DA ) 1 +

DLDACHliq2

ik

∑ij Kpr(S

[SOij]

(9) (10)

+ OijaSRik)

where ω represents the type of the elementary step. 6.3.2. Cyclization Steps. To model the rate of cyclization steps, a reaction mechanism is proposed that proceeds through the sequential formation of monoolefin and diolefin species, followed by the protonation on the acidic sites producing olefinic carbenium ions. Cyclization of these carbenium ions results in the formation of an extra naphthenic ring, as shown in Figure 7. Sorbed concentrations of the monoolefin SOij and diolefin SDOij can be obtained from

[SOij] )

[SDOij] )

liq H Cliq KDH(S iaSO ) Si Si ij

liq Kliq H Cliq KDH(S iaSO ) DH(SO /SDO ) Si Si ij

Table 3. Classification of Aromatics, Based on the Number of Moles of Hydrogen Required for Hydrogenationa

(11)

DLCHliq2 ij

Figure 7. Mechanism for the cyclization steps.

ij

(12)

DLCHliq2 2

Using the protonation equilibrium coefficient of SDOij, the rate + is given by of cyclization of the carbenium ion SRik

rcyc(m;n) ) liq Kliq K ne,ikuVk˜cyc(m;n)KDH(S + )C t iaSOij) DH(SiaSDOij) pr(SDO aSR ij

ik

HSiCSliqi DLDACHliq2 2 (13)

6.3.3. Dealkylation of Aromatics. The surface concentration of the benzenium ions is given in terms of the liquid-phase concentration of the aromatic specie, SAi, by

a The symbol “H” within a ring structure indicates the ring to be hydrogenated.

liquid-phase concentrations so that

[SA+ik] ) Kpr(SA aSA+ )Ct i

ik

HSA CSliqA i

i

(14)

DLDA

rhyd )

rdealk(m) ) ne,ikuVk˜dealk(m)Kpr(SA aSA+ )Ct i

ik

(

comp K HSA CSliqA (CHliq2)n khyd,nH 2 C,SA

The rate of dealkylation steps is then given by

i

HSA CSliqAi i

DLDA

(15)

where m is the nature of the product paraffinic carbenium ion. 6.3.4. Saturation of Aromatics. The present model considers aromatic species with up to four rings. Several experimental studies and review articles have been published on the hydrogenation of aromatic model compounds such as benzene, toluene, alkyl benzenes, tetraline, naphthalene, biphenyl, phenanthrene, pyrene, chrysene, flourene and fluoranthene, etc.22-42 Korre et al.32,34 established a relationship between molecular structures and the hydrogenation reactivities in heavy oil hydroprocessing via the elucidation of the rate controlling reaction pathways of one-, two-, three-, and four fused ring aromatics. A modified version of the rate expression used by Korre has been employed in the current model. The hydrogen partial pressure was replaced by the liquid-phase molar concentration of hydrogen and the liquid-phase concentration of hydrocarbons by their sorbed concentrations. The sorbed concentrations are then expressed in terms for the observable

[

DL +

∑i KC,S

i

i

CSliqN

i

liq Keqm

)

HSA CSliqA + KC,SN HSN CSliqN + KC,SP HSP CSliqP

Ai

i

i

i

i

i

i

i

i

]

(16)

comp is the composite rate In the aforementioned equation, khyd,nH 2 coefficient for the hydrogenation of aromatics, given by the product of the rate coefficient of the surface controlling reaction and the total concentration of the active metal sites (i.e., khyd,nH2cm). The subscript nH2 relates to the consumption of n moles of hydrogen per mole of aromatics. Based on the observations of Korre et al.32,34 and the structures of the aromatic species considered in the reaction network generation in this work, only two hydrogenation rate coefficients have been used for all the possible aromatic hydrogenation reactions: one for the reactions in which three moles of hydrogen are consumed and another for the reactions in which two moles of hydrogen are consumed. Table 3 shows all the possible hydrogenation modes for onering to four-ring aromatic species that are considered in the

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007 5889

reaction network of aromatics. The “H” in the ring represents the ring to be hydrogenated. All the molecular classes of subgroup 1 have been assumed to have the same hydrogenation rate coefficient, irrespective of the length and number of alkyl substituents and the number of saturated rings. The same is true for the classes of subgroup 2. The difference in the reactivities of the molecules of each subgroup is taken into account by the differences in their sorption equilibrium coefficients in the zeolite pores and their chemisorption equilibrium coefficients on the metal sites. The values of the liquid-phase hydrogenation equilibrium liq ) are obtained from the gas-phase hydrogenacoefficient (Keqm gas ), the molar concentration of tion equilibrium coefficient (Keqm the liquid phase, and the vapor-liquid equilibrium partition coefficients (KVLE) of the hydrogen, aromatics, and naphthenic species by means of eq 17: liq Keqm

)

gas Keqm

( )( ) Cliq total

3

KSVLE N i

KHVLE 2

(17)

KSVLE A i

gas The Keqm terms were calculated using Benson’s group contribution method,43 and the VLE partition coefficients using the Peng-Robinson equation of state.44 The chemisorption equilibrium coefficients on the metal sites are modeled using the correlation developed by Korre and Klein.30 6.4. Development of Global Rate Expressions. It was shown, in section 5, that 1266 pure components/GOIs were required for VGOs up to C40. All the isomers in a particular GOI belong to one molecular class and have the same number of methyl branches. Using the computer-generated reaction network, the contributions of all the individual elementary steps are added up to yield the global rate of conversion of a GOI g to another GOI h, as

(g;h) Rω(m;n)

)



( )

liq ne,ikuVk˜ω(m;n)KDH,(S K + )Ct iaSOij) pr(SO aSR ij

ik

HSiCSliqi

DLDACHliq2

(18)

The aforementioned equation represents the conversion through PCP, acyclic, exocyclic, or endocyclic β-scission steps. In eq 18, the concentration of species Si is obtained from its equilibrium distribution in GOI g:

CSliqg CSliqi ) ySeqm i

(19)

The same value of Henry’s coefficient is used for all the isomers of a GOI, so that the denominator term for sorption is given by

[

DL ) 1 +

Nlumps



g)1

KL,gCSliqg

]

(20)

To estimate the value of DA, the same value of the equilibrium coefficient is considered for the protonation of all olefinic species that originate from the isomers of a GOI: Nsat_lumps+Naro_lumps

Nsat_lumps

DA ) 1 +



g)1

Kpr,Og[SOg] +



g)Nsat_lumps+1

Kpr,Ag[SAg] (21)

where [SOg] is the sorbed concentration of the olefins formed

from the isomers of GOI g in the case of paraffinic and naphthenic GOIs and [SAg] is the sorbed concentration of aromatic GOI g. In previous work on the single-event kinetic modeling of hydrocracking, the protonation equilibrium coefficient was estimated by involving a reference olefin SOr, so that the protonation step SOij a SR+ik was written as

SOij a SOr a SR+ik

(22)

The protonation equilibrium coefficient is given by

˜ isom(SO aSO )K ˜ pr(SO am) K ˜ pr(SO aSR+ ) ) K ij

ij

ik

r

(23)

r

It was shown, in the current work, that the involvement of the reference olefins is not necessary to obtain the protonation equilibrium coefficients.45 Therefore, irrespective of the olefinic species SOij, only two single-event protonation equilibrium coefficients, K ˜ pr(SOijas) and K ˜ pr(SOijat), are required in the model. The final rate equation is obtained by combining eqs 18 and 19: (g;h) Rω(m;n)

)

(g;h) comp LCω(m;n) k˜ω(m;n)

( )( ) HgCSliqg

Cliq total

DLDACHliq2 KHVLE 2

(24)

comp where k˜ω(m;n) is the composite single-event rate coefficient, which is given by

comp k˜ω(m;n) ) k˜ω(m;n)K ˜ pr(SO am)Ct

(25)

ij

(g;h) and LCω(m;n) by

(g;h) LCω(m;n)

)

∑ ne,ikuV

( ) σSglOij

+ σSglRik

gas KDH,(S yeqm iaSO ) Si ij

(26)

The LCs are functions of temperature only and are not dependent on the feedstock composition. Their values, calculated for temperature intervals of 5 °C and stored in the computer memory, are utilized during the integration of the ordinary differential equations (ODEs) in the adiabatic reactor simulations. 6.5. Competitive Chemisorption on the Acidic Sites. Accounting for the strong competitive chemisorption of aromatics on the acidic sites requires estimation of the denominator term DA that appears in the rate expressions. The evaluation of DA, in turn, requires the dehydrogenation equilibrium coefficients of the paraffins and naphthenes and the protonation equilibrium coefficients for olefins and aromatics. Representative values of the dehydrogenation equilibrium coefficients were estimated using Benson’s group contribution method. It was assumed that the aromatics are primarily protonated in the ring. The gas-phase heats of protonation for olefins and aromatics are derived from the corresponding proton affinities (PAs) that have been compiled by Hunter and Lias.46 A single value of PA was considered for all the GOIs of a particular class, assuming that its value is invariant with respect to the carbon number (CN). Because the PAs of all the classes considered in the model are not available, certain simplifying assumptions were made,45 based on the trends from the available data. The heats of protonation at the surface of the catalyst are obtained using the heats of stabilization of carbenium ions from the gas phase to the catalyst site. Because of the relatively low

5890

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007

importance of the heats of stabilization of paraffinic and naphthenic carbenium ions in the presence of large amount of aromatics, a constant value was taken from the literature.47 In the absence of information on the heats of stabilization of different aromatic-ring-containing ions, a single value, ∆q(S+ A ), was accepted. This value was considered as a parameter in the model. The entropy of protonation was estimated using statistical thermodynamics, based on the loss of various degrees of freedom (DOFs) in the protonation step. As the protonation occurs, the olefinic and aromatic species from the sorbed phase are linked to the acidic sites and, thus, a difference in entropy between the sorbed state and the protonated state is needed. As a molecule loses some of its rotational or translational DOFs, it gains some extra vibrational DOFs. It has been assumed that the difference in the sum of the rotational and vibrational DOFs of the molecule from the sorbed state and the protonated state is negligible. Further assuming that the sorption step results in the loss of one translation DOF and that the protonated species has no translational motion, two translational DOFs are taken as the entropy of the protonation step. 6.6. Sorption in the Zeolite Pores. Denayer et al.48 performed experiments on the sorption of alkanes, aromatics, and some organic molecules in the liquid and in the dense vapor phase on FAU zeolites. They reported the values of the partition coefficients for normal paraffins from C5 to C16, some light isoparaffins, benzene, toluene, and mesitylene in Na-USY zeolite at room temperature, with n-octane or methanol as the mobile phases. With the nonpolar n-octane as the mobile liquid phase, the partition coefficients of paraffins were determined to be very weak functions of the CN of the adsorbate. With the polar methanol as mobile liquid phase, the partition coefficients of paraffins increased with CN. Accounting for the presence of some polar compounds in the hydrocracking of VGO, it is assumed that the partition coefficients only slightly increase with CN. Three different base values of partition coefficients were selected for paraffins, naphthenes, and aromatics (i.e., Kparaf,base, KNap,base, and KAro,base, respectively). An increase in the partition coefficient by 8% per CN was considered to account for the effect of chain length. The increase in the value of the partition coefficients in the ring compounds from its base value is related to the number of saturated acyclic carbons in the molecule. For example, a diaromatic C15 molecule contains five saturated acyclic carbons, so that its partition coefficient is calculated from the base value from

KDAR,C15 ) KAro,base × (1.08)5

(27)

Initial estimates of the base values of the partition coefficients were taken from the literature.48 The percentage increase per CN and the base values were later adjusted during the parameter estimation, using the VGO hydrocracking experimental data. The Henry coefficients were calculated from the partition coefficients by dividing them by the catalyst bulk density. 7. Reactor Model The global rate expressions are plugged into a model for three-phase multibed tubular reactors operated adiabatically with interstage hydrogen quenching. The particles are assumed to be completely wetted and the gas and liquid phases are assumed to be in plug flow, under the trickle-flow regime. The interphase mass-transfer flux is described using the two-film theory.49 The continuity equations for the gas- and liquid-phase components/ GOIs and hydrogen are written as follows:

(

gas CSgas 1 dFSg g ) -kO,SgaV C,VLE - CSliqg Ω dz KSg

)

(for g ) 1, 2, ..., Nlumps) (28)

(

)

liq CSgas 1 dFSg g ) kO,SgaV C,VLE - CSliqg + RSForm g,net Ω dz KSg

(for g ) 1, 2, ..., Nlumps) (29)

In the aforementioned equations, the overall mass-transfer coefficient is given by

1 1 1 ) + kO,Sg k KC,VLE kL G Sg

(30)

The value of the liquid-side mass-transfer coefficient (kLav) was calculated from a correlation given by Sato50 and the gas-side mass-transfer coefficient (kGav) was calculated from Reiss’ correlation.51 The gas-liquid interfacial area (av) was calculated is the from the correlation derived by Charpentier.52 KSC,VLE g concentration equilibrium partition coefficient relating the equilibrium concentrations of GOI Sg in the gas and liquid ) CSgas /CSliqg . This is expressed in terms of the phases, KSC,VLE g g , as KSC,VLE ) (ySg/ true equilibrium partition coefficient, KSVLE g g gas liq liq VLE gas VLE xSg)(Ctotal/Ctotal) ) KSg (Ctotal/Ctotal). The KSg values are calculated using the Peng-Robinson equation of state. The liquid-phase molar concentrations in eqs 28 and 29 are calculated using the density of the liquid phase. The latter is obtained from the local density of the individual GOIs, as estimated by the Hankinson and Thomson (HBT) method and the liquid-phase mole fractions of the GOIs. This approach efficiently accounts for the volume expansion of the hydrocarbon mixture along the reactor. A significant amount of heat is released during the hydrogenation of aromatics/naphtheno-aromatics and the hydrocracking reactions. The liquid-phase temperature profile along the reactor beds is obtained from

(∑

Nlumps

dTliq dz

)

g)1

liq Form (-∆Hf,g )RSg,net

-

Nlumps

dFSgas g

g)1

dz



Nlumps



g)1

)

λS g

(31)

gas liq (FSgas CP,S + FLliqg CP,S ) g g g

In the numerator of eq 31, the first term represents the total heat of reaction, calculated using the heats of formation of individual GOIs/pure components obtained from the Benson’s group contribution method. The second term accounts for the heat consumed in the vaporization of hydrocarbons. The heattransfer resistance between the gas phase and the liquid phase was assumed to be negligible. The temperature of the solid phase is obtained from Nlumps

T solid ) Tliq +

∑ g)1

liq Form (-∆Hf,g )RSg,net

aLShLS

(32)

in which the solid-liquid heat transfer coefficient (hLS) has been estimated using the correlation of Whitaker.53

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007 5891

Figure 8. Comparison of the experimental and predicted product composition: (a) carbon number (CN) distribution (C3 to C40) of the reactor effluents; (b) CN distribution (C3 to C40) of the paraffins; (c) CN distribution (C5 to C40) of the mononaphthenes; (d) CN distribution (C10 to C40) of the dinaphthenes; (e) weight percentage of normal paraffins, isoparaffins, one-ring to four-ring naphthenes and one-ring to four-ring aromatics; and (f) weight percent of commercial fractions (i.e., LPG, light naphtha, heavy naphtha, kerosene, diesel, and the unconverted VGO).

The set of eqs 28, 29, and 31 consists of 2535 () 1267 × 2 + 1) stiff ODEs. It is solved via the backward differentiation method, using an approximate analytical Jacobian to accelerate the integration. The initial conditions for the gas and liquid molar flow rates are obtained by a VLE flash calculation at the reactor inlet temperature and pressure. A maximum allowable temperature rise (∆Tmax) is specified for each reactor bed. The integration is performed until the reactor temperature reaches a temperature of Tin + ∆Tmax, after which a cold stream of hydrogen is injected, to quench the reaction mixture. 8. Parameter Estimation The approach described here has reduced the number of independent rate parameters for the gigantic network encountered here to only 33 temperature and feedstock composition independent parameters. These were estimated from detailed experimental data on the three-phase hydrocracking of VGO in a bench-scale reactor and commercial/literature information. The product composition was available in terms of the weight percentage of normal and isoparaffins, one-ring to four-ring naphthenes, and one-ring to four-ring aromatics classes. An estimate of the CN distribution was available for normal paraffin, isoparaffin, mononaphthene, dinaphthene, and trinaphthene classes. The parameters were estimated using LevenbergMarquardt’s optimization algorithm by minimizing the sum of square of the residuals between the experimental and model predicted responses. The optimized parameters related to steps on the acid sites have been determined to satisfy the rules of carbenium ion chemistry. The activation energies corresponding to the saturation of aromatics are in the range of 25-50 kJ/ mol, as reported in the literature.30,54 The relative magnitude of

the parameters related to the sorption in the zeolite pores and the chemisorption on the metal sites are in accordance with their physicochemical nature. Considering the complexity of the reaction network and the relatively small number of parameters, the model provides a very good agreement with the experimental data, as shown by the parity plots of Figure 8. 9. Reactor Simulations The commercial reactor geometry was taken from Pacheco and Dassori,55 and the composition of the feedstock VGO was taken from Moustafa and Froment.56 The actual feedstock composition used as the input for the reactor simulations is derived from this composition, as discussed in section 5. The conversion of VGO is based on the cracking of C20+ feed into C20 and lighter products and is calculated by

X)

(

)

feed prod WC20+ - WC20+ feed WC20+

× 100

(33)

Simulations have been conducted at two different reactor inlet temperatures (350 and 370 °C), both at a total pressure of 150 bar and a H2-to-hydrocarbon molar ratio of 15.0 (∼5000 scfb). (Note: scfb ) standard cubic feet per barrel.) The temperature increase per bed was limited to 15 °C. Consequently, the number and length of the beds were dependent on the total amount of heat released by the hydrocracking reactions and the heat release profile corresponding to the selected operating conditions. The model can also be used to simulate a hydrocracking unit with a given number of beds and pre-specified bed lengths to obtain the temperature increase, the conversion, and the product profiles

5892

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007

Figure 9. Effect of reactor inlet temperature on the (a) temperature profile of the solid and liquid phases. (b) Weight percent conversion of VGO. (c) Profiles of the gas and liquid-phase hydrogen flux. (d) Profiles of the weight-percent vaporization of hydrocarbons. (Figure legend: solid lines, 350 °C; dotted lines, 370 °C.)

in each bed. To integrate the set of ODEs, the model requires ∼100-120 seconds on a PC with a 2 GHz Intel Centrino Duo processor and 1 GB RAM. 9.1. Effect of Inlet Temperature on the Overall Reactor Performance. Figure 9a shows the temperature profiles of the solid and liquid phases along the reactor. For Tin ) 350 °C, the total adiabatic temperature increase in the reactor is 40.8 °C, and three beds are required to limit the ∆T per bed within the imposed range. The first bed is shorter than the second one, because of the higher rate of heat release, which is due to higher rates of the exothermic aromatics hydrogenation and cracking reactions. For Tin ) 370 °C, on the other hand, the total adiabatic temperature increase in the reactor is 58.8 °C and four beds are required. The average ∆T between the solid and the liquid phases is ∼0.4 °C at Tin ) 370 °C and ∼0.2 °C at Tin ) 350 °C. The profile of VGO conversion along the beds is shown in Figure 9b, indicating 88.6% total conversion for Tin ) 370 °C, compared to only 56.7% conversion for Tin ) 350 °C. A conversion increase of ∼26% has been reported in the literature57 with a temperature increase of 22 °C. Figure 9c shows the flux of hydrogen in the gas and liquid phases along the reactor. The gas-phase hydrogen flux (blue lines) decreases along each bed, because of the transfer of hydrogen from the gas to the liquid phase. A step increase in the gas-phase hydrogen flux occurs at the bed exits, because of the injection of quench hydrogen. The hydrogen flux in the liquid phase is dependent on the solubility of hydrogen in the liquid phase, the rate of hydrogen transfer from the gas to the liquid phase, and the total amount of hydrocarbons present in the liquid phase. Figure 9d shows that the weight-percent vaporization of

hydrocarbons increases very rapidly at higher temperature, because of higher amounts of lighter hydrocarbons formed by faster cracking. 9.2. Yields of Various Commercial Fractions. Figure 10 shows that, at Tin ) 350 °C, the diesel fraction (DIES) increases in the first bed and remains almost steady in the last two beds, ultimately resulting in a diesel yield higher than those of the other commercial fractions. At Tin ) 370 °C, the yield of the diesel fraction increases quite fast in the first bed, maintains almost a steady value in the second bed, and then sharply decreases in the last two beds. The cracking of diesel range hydrocarbons dominates in the third bed, causing a sharp increase in the heavy naphtha (HNAP) yield. At 370 °C, the yield of heavy naphtha is highest among all the fractions, indicating a shift in the hydrocracking from middle distillate to the naphtha mode. The yield of kerosene (KERO) increases faster in the first bed at 370 °C than at 350 °C but reaches almost the same value at the exit of the reactor. The light naphtha (LNAP) and LPG yields increase monotonically along the reactor and are favored by higher temperatures. 9.3. Yields of Individual Classes in Various Commercial Fractions. Figures 11a-f show the yields of several classes in various commercial fractions along the reactor at 370 °C. For LPG (C3-C4), light naphtha (LNAP, C5-C7), and heavy naphtha (HNAP, C8-C12) fractions (Figures 11a, 11b, and 11c, respectively), the yields of almost all the hydrocarbon classes increase with reactor length. In LPG, the yield of isobutane is always higher than that of n-butane, which is consistent with the industrial observations. In LNAP, isoparaffins also exhibit the highest yield, followed by n-paraffins and mononaphthenes.

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007 5893

Figure 10. Effect of reactor inlet temperature on the evolution of various commercial fractions. (Figure legend: solid lines, 350 °C; dotted lines, 370 °C.)

Figure 11. Evolution along the reactor of different components/classes of (a) LPG, (b) light naphtha (LNAP), (c) heavy naphtha (HNAP), (d) kerosene (KERO), (e) diesel (DIES), and (f) unconverted VGO (UNCONV). (Inlet temperature ) 370 °C.)

In HNAP, mononaphthenes constitute more than 50% of the heavy naphtha fraction. The rate of increase of mononaphthenes is much higher, in comparison to any other hydrocarbon class. This can be attributed to the large amount of mononaphthenes in the feedstock. When subjected to hydrocracking, they primarily produce mononaphthenes of the HNAP range. At 370 °C, the cyclization of paraffins and the saturation of monoaro-

matics are also important for the production of the mononaphthenic species of this range. The profile of monoaromatics in the last three beds of the reactor indicates that the rate at which monoaromatics of the HNAP range are formed from the cracking of heavier monoaromatics is similar to the rate at which they are hydrogenated to mononaphthenes, resulting in an almost flat profile of monoaromatics in these beds. As in LPG and

5894

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007

Table 4. Distribution of Isoparaffins in the Feed for Set I Simulations Distribution (mol %) monoditribranched branched branched base case case 1 case 2

65 100 0

25 0 100

10 0 0

VGO conversion (wt %) 56.7 53.2 62.4

Table 5. Distribution of Isomers in the Feed within the Ring-Containing Classes for Set II Simulations unbranched base case case 3 case 4 case 5

20 100 0 0

Distribution (mol %) monodibranched branched 50 0 100 0

20 0 0 100

tribranched

VGO conversion (wt %)

10 0 0 0

56.7 23.7 54.5 74.2

LNAP, the ratio of iso-paraffin to normal paraffin exceeds 1, which is a general feature of hydrocracking. The evolution of paraffins, naphthenes, aromatics, and naphtheno-aromatics in the kerosene (KERO, C13-C15), diesel (DIES, C16-C20), and unconverted VGO (UNCONV, C21-C33) fractions are shown in Figures 11d, 11e, and 11f, respectively. The yields of paraffins and naphthenes rapidly increase initially for both KERO and DIES, reach a maximum, and then decrease in the last part of the reactor. The decrease in the yields of paraffins and naphthenes in the DIES fraction starts before the KERO fraction, indicating an increase of the apparent hydrocracking rate coefficient with CN. For both fractions, the yield of aromatics decreases along the reactor, because of ring saturation. Given the larger amount of aromatics in the DIES fraction, its rate of decrease is significantly higher, compared to that of the KERO fraction. The amount of paraffins and naphthenes in the UNCONV monotonically decreases along the reactor. 9.4. Dependence of the Product Yields on the Composition of the Various Hydrocarbon Classes. The VGO composition taken from Moustafa and Froment56 has been converted to a more-detailed form, in which each hydrocarbon class is divided into subclasses based on the number of methyl branches. Two sets of simulations were performed using the distribution of isomers shown in Tables 4 and 5. In Set I, only the distribution of isoparaffins has been altered, while the distribution of all other ring-containing classes was kept at their base values. In Set II, the distribution of isoparaffins was kept at the base value and those of ring-containing classes was varied. For different cases, the VGO conversions obtained from the reactor simulations are shown in Tables 4 and 5 and are plotted in Figure 12. All the simulations reported in this section are performed at 350 °C, 150 bar, and a H2-to-hydrocarbon molar ratio of 15. The total amount of isoparaffins in the VGO feed is only 12.2 wt %, so that the effect of the variation in the distribution of isoparaffins on the VGO conversion while the distribution of all other classes is kept at the base values is relatively small (see Table 4). Changing the distribution of ring-containing classes causes a marked difference in the VGO conversion, which was attributed to the high content of ring species in the feed (see Table 5). Figure 13 shows the evolution of the molepercent composition of normal, monobranched, dibranched, and tribranched isomers in the total paraffin fraction along the length of the reactor for Set I simulations. When tribranched paraffins are present in the feed (e.g., the base case), their amount decreases rapidly without any corresponding increase in the amount of dibranched isomers, which is indicative of fast

Figure 12. Effect of the distribution of isomers in the feed on the VGO conversion (for different simulation cases of Set I and Set II).

Figure 13. Molar distribution of the normal, monobranched, dibranched, and tribranched isomers in the paraffin fraction versus reactor length for Set I simulations. (Figure legend: blue, NPA; green, MBP; red, DBP; and cyan, TBP. Solid lines, base case; dotted lines, case 1; and dashed lines, case 2.)

cracking. Dibranched paraffins are also rapidly converted by cracking when present in large concentrations in the feed (e.g., case 2). In both cases, increasing the percentage of isomers with a higher degree of branching increases the total conversion (see Figure 12), because the number of favorable pathways available for cracking steps increases as the number of methyl branches increases. In both case 1 and case 2, the final composition of isomers approaches the base case values, within a range of (5%. In all cases, the amount of isomers with different degrees of branching tend toward asymptotic values, but these are not the thermodynamic equilibrium values. If all the paraffinic isomers (NPA, MBP, DBP, and TBP; see Table 1) of a given CN were able to reach thermodynamic equilibrium among them, the total amount of tribranched isomer fraction would have been the highest, with decreasing amounts of dibranched, monobranched, and normal paraffins, in that order. This can be related to the increase in the number of isomers with higher number of methyl branches and a lower Gibbs free energy of formation of isomers with more methyl branches. Although the curves shown in Figure 13 represent the distribution of isomers for the entire CN range

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007 5895

(C4 to C33) and not just for one particular CN, it can be concluded that the final distribution is very far from the thermodynamic equilibrium distribution. The simulations of Set II for the ring-containing classes lead to the same conclusions. The deviation of the hydrocracking product yields from thermodynamic equilibrium is well-documented.20,21 The aforementioned observations illustrate the need to divide the isomers of a class of a given CN into four subclasses, based on the number of methyl branches. Without this subdivision, the predicted VGO conversion would be incorrect. 10. Conclusions The hydrocracking of vacuum gas oil (VGO) is a three-phase catalytic process (gas-liquid-solid) in which a broad variety of hydrocarbon classes are transformed through a huge network of reactions of very different types into an extremely wide product spectrum of commercial relevance. Its simulation requires a comprehensive kinetic model, reflecting the complexity of the pattern of gradual conversion of feedstock and intermediates into products. In contrast with most of the previous modeling attempts, in which drastic and unrealistic simplifications were introduced, the present paper retains the full details of the feed composition and reaction network. The kinetics of the hydrogenation steps of the unsaturated ring components on the metal sites of the catalyst were systematized, based on their molecular structure. The kinetics of the steps on the acidic sites were written in an even more fundamental way, i.e., in terms of elementary steps of carbenium ion chemistry. Although the reaction network thus generated for each reacting species becomes gigantic, it is precisely this fundamental approach, combined with the modeling of the rate parameters, that allows the number of independent parameters in the kinetic model to be reduced to only 33. This is a number that permits reliable parameter estimation from a judiciously chosen set of experimental data and commercial information. Because of their fundamental nature, these parameters are invariant, with respect to the feedstock, thus saving costly pilot experimentation. Combined with a reactor model comprising 1267 conservation equations per fluid phase, it became possible to simulate the hydrocracking process and predict its product distribution, not just in terms of commercial fractions, but on a much more detailed level, even up to single components. The model will serve the design, operation, and optimization of hydrocracking units, as well as the feedstock selection and the development of targeted hydrocracking catalysts. Acknowledgment The authors are grateful to Dr. Rayford G. Anthony for his interest and the financial support for this project provided from the C. D. Holland Professorship. Nomenclature aV ) gas-liquid interfacial area [mi2/mr3] aLS ) liquid-solid interfacial area [mi2/mr3] A ˜ ω ) single-event frequency factor for ω-type elementary steps [1/h] cm ) total concentration of the metal sites [kmol/kgcat] gas ) gas phase heat capacity of GOI Sg [kJ/(mol K)] CP,S g liq CP,S ) liquid phase heat capacity of GOI Sg [kJ/(mol K)] g ) total concentration of the gas phase [kmol/m3] Cgas total Cliq ) total concentration of the liquid phase [kmol/m3] total Ct ) total concentration of the acid sites [kmol/kgcat]

CHliq2 ) concentration of hydrogen in the liquid phase [kmol/m3] CSliqi ) concentration of species Si in the liquid phase [kmol/m3] CSliqg ) concentration of GOI Sg in the liquid phase [kmol/m3] CSgas ) concentration of GOI Sg in the gas phase [kmol/m3] g Eω(m;n) ) activation energy for the steps of the ω(m;n) subtype [kJ/mol] ) gas-phase flow rate of GOI Sg [kmol/h] FSgas g FSliqg ) liquid-phase flow rate of GOI Sg [kmol/h] GOI ) group of isomers h ) Planck’s constant [kJ h] hLS ) solid-liquid heat transfer coefficient [kJ/(m2 K h)] ∆H°† ) enthalpy of activation [kJ/mol] ∆Hexo-β ) heat of the exocyclic β-scission elementary steps [kJ/mol] ∆Hacyc-β ) heat of the acyclic β-scission elementary steps [kJ/ mol] ∆Hliq f,g ) heat of formation of GOI Sg in the liquid phase [kJ/ mol] ∆Hpr(SOijam) ) enthalpy of protonation of olefinic species SOij HSi ) Henry coefficient for the sorption of species Si [m3/kgcat] kB ) Bolztman’s constant [kJ/(molecule °C)] k˜ω(m;n) ) single-event rate coefficient for the elementary step of type ω(m;n) [1/s] kO,Sg ) overall mass-transfer coefficient for GOI Sg [mi/h] kG ) gas-phase mass-transfer coefficient [mi/h] kL ) liquid-phase mass-transfer coefficient [mi/h] comp khyd,nH ) composite rate coefficient for the hydrogenation of 2 monoaromatics [(kmol/(kgcat/h))(m3/kmol)n] KHVLE ) vapor liquid equilibrium partition coefficient for 2 hydrogen ) concentration based vapor liquid equilibrium coefKSC,VLE g ficient for GOI Sg KC,SA ) chemisorption equilibrium coefficient for aromatic i species Ai at the metal site [m3/kmol] liq KDH(SiaSO ) ) equilibrium coefficient for the given overall ij dehydrogenation reaction in the liquid phase [kmol/m3] gas KDH(S ) equilibrium coefficient for the given overall iaSOij) dehydrogenation reaction in the gas phase [kmol/m3] KL,Si ) sorption equilibrium coefficient of species Si, [m3/kmol] KL,Sg ) sorption equilibrium coefficient of GOI Sg, [m3/kmol] Kpr ) protonation equilibrium coefficient ne ) number of single events Nlumps ) total number of GOIs/pure components Nsat_lumps ) number of saturated GOIs/pure components Naro_lumps ) number of aromatic GOIs/pure components R ) universal gas constant [kJ/(mol K)] (g;h) ) global rate of conversion of GOI g to GOI h through Rω(m;n) ω(m;n) type of elementary steps [kmol/(kgcat h)] ) net rate of formation of GOI Sg on acid sites [kmol/ RSForm g,net (mr3 h)] [Si] ) sorbed concentration of species Si in the zeolite pores [kmol/kgcat] [SR+ik] ) concentration of surface carbenium ions SR+ik [kmol/ kgcat] ∆Sˆ pr ) intrinsic entropy of protonation [kJ/(mol K)] ∆Sˆ °† ) intrinsic entropy of activation [kJ/(mol K)] Tliq ) liquid-phase temperature [K] Tsolid ) solid-phase temperature [K] feed ) weight of C20+ hydrocarbons in the feed [kg] WC20+ prod WC20+ ) weight of C20+ hydrocarbons in the products [kg] xi ) mole fraction of component i in the liquid phase xSg ) mole fraction of GOI Sg in the liquid phase

5896

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007

X ) weight percent conversion of VGO ) equilibrium mole fraction of species i in the correySeqm i sponding GOI ySg ) mole fraction of GOI Sg in the gas phase z ) axial reactor coordinate [mr] Greek Symbols R ) transfer function Ω ) reactor cross-sectional area [mr2] λSg ) representative latent heat of vaporization of GOI Sg [kJ/ mol] σigl ) global symmetry number of species i Literature Cited (1) Weekman, V. W.; Nace, D. M. Kinetics of Catalytic Cracking Selectivity in Fixed, Moving, and Fluid Bed Reactors. AIChE J. 1970, 16 (3), 397. (2) Jacob, S. M.; Gross, B.; Voltz, S. E.; Weekman, V. W. Lumping and Reaction Scheme for Catalytic Cracking. AIChE J. 1976, 22 (4), 701. (3) Stangeland, B. E. Kinetic-Model for Prediction of Hydrocracker Yields. Ind. Eng. Chem. Process Des. DeV. 1974, 13 (1), 71. (4) Liguras, D. K.; Allen, D. T. Structural Models for Catalytic Cracking. 1. Model-Compound Reactions. Ind. Eng. Chem. Res. 1989, 28 (6), 665. (5) Liguras, D. K.; Allen, D. T. Structural Models for Catalytic Cracking. 2. Reactions of Simulated Oil Mixtures. Ind. Eng. Chem. Res. 1989, 28 (6), 674. (6) Quann, R. J.; Jaffe, S. B. Structure-Oriented Lumping-Describing the Chemistry of Complex Hydrocarbon Mixtures. Ind. Eng. Chem. Res. 1992, 31 (11), 2483. (7) Baltanas, M. A.; Vanraemdonck, K. K.; Froment, G. F.; Mohedas, S. R. Fundamental Kinetic Modeling of Hydroisomerization and Hydrocracking on Noble-Metal-Loaded Faujasites. 1. Rate Parameters for Hydroisomerization. Ind. Eng. Chem. Res. 1989, 28 (7), 899. (8) Vynckier, E.; Froment, G. F. Modeling of the Kinetics of Complex Processes Based upon Elementary Steps; Elsevier Science Publishers BV: Amsterdam, The Netherlands, 1991; Vol. 131. (9) Clymans, P. J.; Froment, G. F. Computer-Generation of Reaction Paths and Rate-Equations in the Thermal-Cracking of Normal and Branched Paraffins. Comput. Chem. Eng. 1984, 8 (2), 137. (10) Hillewaert, L. P.; Dierickx, J. L.; Froment, G. F. ComputerGeneration of Reaction Schemes and Rate-Equations for Thermal-Cracking. AIChE J. 1988, 34 (1), 17. (11) Feng, W.; Vynckier, E.; Froment, G. F. Single-Event Kinetics of Catalytic Cracking. Ind. Eng. Chem. Res. 1993, 32 (12), 2997. (12) Svoboda, G. D.; Vynckier, E.; Debrabandere, B.; Froment, G. F. Single-Event Rate Parameters for Paraffin Hydrocracking Oil a Pt/Us-Y Zeolite. Ind. Eng. Chem. Res. 1995, 34 (11), 3793. (13) Martens, G. G.; Marin, G. B.; Martens, J. A.; Jacobs, P. A.; Baroni, G. V. A fundamental kinetic model for hydrocracking of C-8 to C-12 alkanes on Pt/US-Y zeolites. J. Catal. 2000, 195 (2), 253. (14) Kumar, H.; Froment, G. F. A Generalized Mechanistic Kinetic Model for the Hydroisomerization and Hydrocracking of Long-Chain Paraffins. Ind. Eng. Chem. Res. 2007, 46 (12), 4075. (15) Baltanas, M. A.; Froment, G. F. Computer-Generation of Reaction Networks and Calculation of Product Distributions in the Hydroisomerization and Hydrocracking of Paraffins on Pt-Containing Bifunctional Catalysts. Comput. Chem. Eng. 1985, 9 (1), 71. (16) Govindhakannan, J. Modeling of a Hydrogenated Vacuum Gas Oil Hydrocracker. Dissertation, Texas Tech University, Lubbock, TX, 2003. (17) Boduszynski, M. M. Characterization of Heavy Crude Components. Prepr.sAm. Chem. Soc., DiV. Pet. Chem. 1985, 30, 626. (18) Boduszynski, M. M. Composition of Heavy Petroleums. 2. Molecular Characterization. Energy Fuels 1988, 2 (5), 597. (19) Martens, G. G.; Marin, G. B. Kinetics for hydrocracking based on structural classes: Model development and application. AIChE J. 2001, 47 (7), 1607. (20) Vansina, H.; Baltanas, M. A.; Froment, G. F. Hydroisomerization and Hydrocracking. 4. Product Distribution from n-Octane and 2,2,4Trimethylpentane. Ind. Eng. Chem. Product Res. DeV. 1983, 22 (4), 526. (21) Schulz, H. F.; Weitkamp, J. H. Zeolite Catalysts-Hydrocracking and Hydroisomerization of n-Dodecane. Ind. Eng. Chem. Product Res. DeV. 1972, 11 (1), 46. (22) Girgis, M. J.; Gates, B. C. Reactivities, Reaction Networks, and Kinetics in High-Pressure Catalytic Hydroprocessing. Ind. Eng. Chem. Res. 1991, 30 (9), 2021.

(23) Lin, S. D.; Vannice, M. A. Hydrogenation of Aromatic-Hydrocarbons over Supported Pt Catalysts. 1. Benzene Hydrogenation. J. Catal. 1993, 143 (2), 539. (24) Lin, S. D.; Vannice, M. A. Hydrogenation of Aromatic-Hydrocarbons over Supported Pt Catalysts. 2. Toluene Hydrogenation. J. Catal. 1993, 143 (2), 554. (25) Lin, S. D.; Vannice, M. A. Hydrogenation of Aromatic-Hydrocarbons over Supported Pt Catalysts. 3. Reaction Models for Metal-Surfaces and Acidic Sites on Oxide Supports. J. Catal. 1993, 143 (2), 563. (26) Lin, S. D.; Vannice, M. A.; Herrmann, J. M.; Wang, D.; Apestequia, C.; Duprez, D.; Figueras, F.; Conner, W. C.; Kiperman, S. L.; Hall, W. K.; Blackmond, D. G.; Grunert, W.; Butt, J. B.; Schulz, H. Toluene Hydrogenation over Supported Platinum Catalysts. Stud. Surf. Sci. Catal. 1993, 75, 861. (27) Girgis, M. J.; Gates, B. C. Catalytic Hydroprocessing of Simulated Heavy Coal Liquids. 2. Reaction Networks of Aromatic-Hydrocarbons and Sulfur and Oxygen Heterocyclic-Compounds. Ind. Eng. Chem. Res. 1994, 33 (10), 2301. (28) Girgis, M. J.; Gates, B. C. Catalytic Hydroprocessing of Simulated Heavy Coal Liquids. 1. Reactivities of Aromatic-Hydrocarbons and Sulfur and Oxygen Heterocyclic-Compounds. Ind. Eng. Chem. Res. 1994, 33 (5), 1098. (29) Klein, M. T.; Korre, S. C.; Read, C. J.; Russell, C. L. Hydrocracking of Model Polynuclear Aromatics: Reaction Pathways, Kinetics, and Structure/Reactivity Correlations. Prepr.sAm. Chem. Soc., DiV. Pet. Chem. 1993, 38 (1), 26. (30) Korre, S. C.; Klein, M. T. Development of temperature-independent quantitative structure/reactivity relationships for metal- and acid-catalyzed reactions. Catal. Today 1996, 31 (1-2), 79. (31) Korre, S. C.; Klein, M. T. Effect of Temperature on Quantitative Structure/Reactivity Relationships for Hydrocracking Polynuclear Aromatics. Prepr.sAm. Chem. Soc., DiV. Pet. Chem. 1995, 40 (4), 210. (32) Korre, S. C.; Klein, M. T.; Quann, R. J. Polynuclear Aromatic Hydrocarbons Hydrogenation. 1. Experimental Reaction Pathways and Kinetics. Ind. Eng. Chem. Res. 1995, 34 (1), 101. (33) Korre, S. C.; Klein, M. T.; Quann, R. J. Hydrocracking of polynuclear aromatic hydrocarbons. Development of rate laws through inhibition studies. Ind. Eng. Chem. Res. 1997, 36 (6), 2041. (34) Korre, S. C.; Neurock, M.; Klein, M. T.; Quann, R. J. Hydrogenation of Polynuclear Aromatic-Hydrocarbons. 2. Quantitative Structure/ Reactivity Correlations. Chem. Eng. Sci. 1994, 49 (24A), 4191. (35) Rahier, H.; Denayer, J. F.; Van Mele, B. Low-temperature synthesized aluminosilicate glassessPart IV. Modulated DSC study on the effect of particle size of metakaolinite on the production of inorganic polymer glasses. J. Mater. Sci. 2003, 38 (14), 3131. (36) Rautanen, P. A.; Aittamaa, J. R.; Krause, A. O. I. Liquid phase hydrogenation of tetralin on Ni/Al2O3. Chem. Eng. Sci. 2001, 56 (4), 1247. (37) Rautanen, P. A.; Aittamaa, J. R.; Krause, K. O. I. Solvent effect in liquid-phase hydrogenation of toluene. Ind. Eng. Chem. Res. 2000, 39 (11), 4032. (38) Rautanen, P. A.; Lylykangas, M. S.; Aittamaa, J. R.; Krause, A. O. I. Liquid-phase hydrogenation of naphthalene and tetralin on Ni/Al2O3: Kinetic modeling. Ind. Eng. Chem. Res. 2002, 41 (24), 5966. (39) Lylykangas, M. S.; Rautanen, P. A.; Krause, A. O. I. Liquid-phase hydrogenation kinetics of multicomponent axomatic mixtures on Ni/Al2O3. Ind. Eng. Chem. Res. 2002, 41 (23), 5632. (40) Lylykangas, M. S.; Rautanen, P. A.; Krause, A. O. I. Hydrogenation and deactivation kinetics in the liquid-phase hydrogenation of isooctenes on Pt/Al2O3. Ind. Eng. Chem. Res. 2004, 43 (7), 1641. (41) Singh, U. K.; Vannice, M. A. Kinetics of liquid-phase hydrogenation reactions over supported metal catalystssa review. Appl. Catal., A 2001, 213 (1), 1. (42) Singh, U. K.; Vannice, M. A. Erratum to “Kinetics of liquid-phase hydrogenation reactions over supported metal catalystsa review”. Appl. Catal., A 2003, 255 (2), 361. (43) Cohen, N.; Benson, S. W. Estimation of Heats of Formation of Organic-Compounds by Additivity Methods. Chem. ReV. 1993, 93 (7), 2419. (44) Peng, D.; Robinson, D. B. New 2-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15 (1), 59. (45) Kumar, H. Mechanistic Kinetic Modeling of the Hydrocracking of Complex Feedstocks. Dissertation, Texas A&M University, College Station, TX, 2006. (46) Hunter, E. P. L.; Lias, S. G. Evaluated gas phase basicities and proton affinities of molecules: An update. J. Phys. Chem. Ref. Data 1998, 27 (3), 413. (47) Dumesic, J. A.; Rudd, D. F.; Aparicio, L. M.; Rekoske, J. E.; Trevino, A. A. The Microkinetics of Hetrogeneous Catalysis; American Chemical Society: Washington, DC, 1993.

Ind. Eng. Chem. Res., Vol. 46, No. 18, 2007 5897 (48) Denayer, J. F.; Bouyermaouen, A.; Baron, G. V. Adsorption of alkanes and other organic molecules in liquid phase and in the dense vapor phase: Influence of polarity, zeolite topology, and external fluid density and pressure. Ind. Eng. Chem. Res. 1998, 37 (9), 3691. (49) Froment, G. F.; Depauw, G. A.; Vanrysselberghe, V. Kinetic Modeling and Reactor Simulation in Hydrodesulfurization of Oil Fractions. Ind. Eng. Chem. Res. 1994, 33 (12), 2975. (50) Sato, Y.; Hirose, H.; Takahashi, F.; Toda, M. Proc., Pac. Chem. Eng. Congr., 1st 1972, 187. (51) Reiss, L. P. Cocurrent Gas-Liquid Contacting in Packed Columns. Ind. Eng. Chem. Process Des. DeV. 1967, 6 (4), 486. (52) Charpentier, J. C. Recent progress in two phase gas-liquid mass transfer in packed beds. Chem. Eng. J. 1976, 11, 161. (53) Whitaker, S. Forced convection heat transfer correlations for flow in pipes, past flat plates, single cylinders, single spheres, and for flow in packed beds and tube bundles. AIChE J. 1972, 18 (2), 361.

(54) Stanislaus, A.; Cooper, B. H. Aromatic Hydrogenation Catalysiss A Review. Catal. ReV.sSci. Eng. 1994, 36 (1), 75. (55) Pacheco, M. A.; Dassori, C. G. Hydrocracking: An improved kinetic model and reactor modeling. Chem. Eng. Commun. 2002, 189 (12), 1684. (56) Moustafa, T. M.; Froment, G. F. Kinetic modeling of coke formation and deactivation in the catalytic cracking of vacuum gas oil. Ind. Eng. Chem. Res. 2003, 42, (1), 14. (57) Weismantel, G. E. Petroleum Processing Handbook; Marcel Dekker: New York, 1992.

ReceiVed for reView March 25, 2007 ReVised manuscript receiVed May 14, 2007 Accepted May 22, 2007 IE0704290