Molecular-Level Modeling and Simulation of ... - ACS Publications

Dec 7, 2015 - Natural Resources Canada, CanmetENERGY-Devon, One Oil Patch Drive, Devon, Alberta T9G 1A8, Canada. ABSTRACT: Advancing the ...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/EF

Molecular-Level Modeling and Simulation of Vacuum Gas Oil Hydrocracking Anton Alvarez-Majmutov, Jinwen Chen,* and Rafal Gieleciak Natural Resources Canada, CanmetENERGY-Devon, One Oil Patch Drive, Devon, Alberta T9G 1A8, Canada ABSTRACT: Advancing the capabilities of refinery process models requires fundamental knowledge of hydrocarbon composition and processing behavior at the molecular level. In common practice, however, the main obstacle to reaching such a level of understanding is the difficulty to characterize the molecular composition of petroleum and its derived products with current analytical methods. A different approach is through the use of hydrocarbon composition modeling techniques to derive the molecular make up of petroleum fractions, thus enabling the development of molecular reaction models. The purpose of this study is to illustrate the application of this concept to model and simulate the vacuum gas oil hydrocracking process at the molecular level. At first the analytical characterization of the feed sample is transformed into a computational mixture of hydrocarbon molecules that is consistent with the chemistry of the actual oil sample. This molecular representation is then used as input to model the chemical transformations occurring in the hydrocracking reactor. The reaction network is organized in terms of reaction families, and reactivity parameters are modeled with quantitative structure/reactivity correlations. The developed model is tuned using experimental data obtained from a bench-scale hydrocracking reactor. Simulations showed that the model reproduces the product distribution by boiling range and hydrocarbon type, relevant product properties (e.g., API gravity), and process parameters such as hydrogen consumption and hydrocarbon vaporization, over a wide range of operating conditions.

1. INTRODUCTION Hydrocracking is a well-established and flexible refinery process for converting heavy fractions of petroleum into high quality distillates without any waste product. Commercial units can process a wide range of feedstocks including straight-run vacuum gas oil, coker heavy gas oil, and deasphalted oil. It is expected that hydrocracking technologies will be important options to produce transportation fuels as a result of the constantly increasing availability of low-value crude oils with abundant content of heavy fractions. Understanding hydrocracking chemistry is crucial for process development and optimization. Ideally, process models should be able to predict detailed product distribution, yields of commercial fractions, and physical and chemical properties of interest (e.g., density, boiling point distribution, viscosity, etc.).1 From a practical perspective, however, it is difficult to meet all these requirements, primarily due to the complexity of petroleum composition and its analysis. Hydrocracker feedstocks (∼343−550 °C) are large-scale systems consisting of a massive number of reacting molecules (104−105), and at the same time, detecting and quantifying all these species is not possible with conventional analytical methods. Consequently, hydrocracking reactions are typically modeled in terms of the macroscopic properties of the reacting mixture. In such an approach, the hydrocarbon mixture is considered as a limited set of chemical lumps (100−101), often defined by boiling point range or some other measurable properties.2,3 Since each lump contains a broad group of molecular classes, there is inevitably a considerable loss of information on molecular behavior. Models of such empirical nature are thus very limited in capabilities and insensitive to variations in feedstock composition.4 Molecular-level models, in contrast to conventional models, offer a realistic description of the process through the Published 2015 by the American Chemical Society

incorporation of fundamental knowledge of molecular structures and reaction pathways and mechanisms. Adopting this approach to simulating hydrocarbon conversion processes has been enabled by implementing petroleum composition models to overcome the limitations of analytical characterization methods.5,6 The function of composition models is to produce a mixture of representative hydrocarbon molecules that is consistent with the experimental properties of the real sample. This representation serves not only as an efficient means of organizing feedstock composition but also as a suitable computing architecture to model chemical transformations. While a few studies provide a fundamental treatment of gas oil hydrocracking reactions,7,8 detailed feedstock composition simulation, integrated with molecular reaction models, has not been extensively explored. Recently, Pereira de Oliveira et al.9,10 demonstrated the application of the composition-reaction modeling approach to other complex refinery processes such as light cycle oil hydrotreating and residue hydroconversion. A distinctive feature of their work is the use of the kinetic Monte Carlo algorithm11 to handle large reaction networks, as an efficient alternative to traditional mathematical methods for simulation. In our previous and ongoing studies,12,13 composition models for middle and heavy distillate fractions were developed and validated with advanced characterization data, to provide a basis for molecular-level reaction modeling of refining processes. The purpose of this study was to develop a Received: September 14, 2015 Revised: November 19, 2015 Published: December 7, 2015 138

DOI: 10.1021/acs.energyfuels.5b02084 Energy Fuels 2016, 30, 138−148

Energy & Fuels

Article

Figure 1. VGO composition simulation algorithm.

tive hydrocarbon molecules that is chemically consistent with the actual oil sample. During the second stage, the molecular transformations occurring in the hydrocracking reactor are simulated using the set of molecules generated in the first stage as a starting point. The details of each step are discussed below. 3.1. Feedstock composition modeling. Description and representation of the VGO feed composition are based on the methodology reported in our previous publications.12,13 The VGO feed is considered to be composed of long chain paraffins, cycloparaffins, mono-, di-, and polynuclear aromatics, and organosulfur and organonitrogen compounds. The molecular spectrum of the VGO is broken down into basic structural blocks, such as alkyl chains, benzene rings, cyclohexane rings, and heterocyclic rings (e.g., thiophene, pyrrole, and pyridine). It is assumed that the frequency of occurrence of these blocks follows well-known probability distribution functions (PDFs), such as gamma distribution, exponential distribution, and histograms. Molecules are constructed by arranging these blocks according to coherent chemical rules. This process is executed by Monte Carlo sampling of the statistical distributions assigned to each structural feature. Figure 1 illustrates the VGO composition simulation algorithm. The “virtual” feed mixture consists of a set of 1 × 104 molecules built by repetitive molecular assemblies. Physical and chemical properties of each molecule are determined either directly from its structure (molecular weight, chemical formula, carbon types, etc.) or by group contributions (density and boiling point). Bulk properties of the hydrocarbon mixture are estimated from the individual properties of the generated molecules by applying mixing rules. The calculated properties are then compared with the analytical characterization through a least-squares objective function. The mixture is optimized by adjusting the parameters of the PDFs until minimum difference between the experimental data and the simulation is achieved. This optimization process is carried out with the simulated annealing method. The optimized hydrocarbon mixture at this stage of the simulation process consists of an equimolar set of representative molecules. The last step in this sequence is tuning the molecular concentrations by applying the principle of maximum entropy. The simulated VGO composition is then used as basis to model hydrocracking reactions.

composition-reaction model for the simulation of vacuum gas oil hydrocracking at molecular level. Bench-scale tests were conducted to generate experimental data for model tuning. The effect of vapor−liquid equilibrium on hydrocracking reactions was also addressed in this work.

2. EXPERIMENTAL SECTION Hydrocracking experiments were conducted in a bench-scale fixed-bed reactor (length 100 cm; internal diameter 9 mm) operating in up-flow mode. Reactor temperature was controlled with a four-zone heating furnace and multiple thermocouples placed along the reactor tube. Gas and liquid products coming out of the reactor were separated in a high-pressure separator (HPS). The liquid product was passed through an atmospheric stabilizer column to further remove any dissolved gases. The gases from the HPS and stabilizer column were combined and analyzed by an online gas chromatograph (HP5880 GC). A commercial zeolite-based hydrocracking catalyst was used for the tests. The middle section of the reactor was loaded with 10 mL of catalyst extrudates diluted with 10 mL of 0.2 mm glass beads for minimum axial dispersion. The remaining reactor space above and below the catalyst bed was packed with different sizes of glass beads to ensure proper gas and liquid distribution. The fresh catalyst was in situ activated with hydrogen sulfide (H2S) according to the procedure provided by the catalyst vendor, and equilibrated by processing the feed at 395 °C for about 200 h. The feedstock used in this study was vacuum gas oil (VGO) from Canadian oil sands bitumen. The VGO feed was hydrotreated prior to hydrocracking to meet the sulfur and nitrogen specifications for the zeolite catalyst. Hydrocracking tests were carried out under steadystate catalyst activity at the following operating conditions: liquid hourly space velocity (LHSV) 0.75−3.0 h−1, temperature 380−400 °C, pressure 11 MPa, and H2/oil ratio 800 NL/L. After completing the experiments catalyst activity was monitored with a check-back test at base conditions (temperature 395 °C and LHSV 1.5 h−1). Standard analytical methods were used to characterize the feedstock and products.

3. MODELING APPROACH In this study, the hydrocracking process is modeled in two sequential stages: 1) feedstock composition modeling and 2) reaction modeling. In the composition modeling stage, the analytical characterization of the feedstock is computationally transformed into a virtual mixture of thousands of representa139

DOI: 10.1021/acs.energyfuels.5b02084 Energy Fuels 2016, 30, 138−148

Energy & Fuels

Article

Table 1. Hydrocracking Reaction Families

molecular detail.15 Moreover, hydrocracking reactions have a recursive nature and therefore can be organized into a limited set of reaction families:16 aromatic ring hydrogenation, paraffin isomerization, paraffin cracking, ring isomerization, ring opening, and ring dealkylation. Table 1 presents representative examples of each reaction family and their specific chemical rules. Other reactions such as hydrodesulfurization and hydrodenitrogenation can also take place during hydrocracking but are omitted for the purposes of this work. The characteristics of each reaction family are described below. Hydrogenation of aromatic rings is considered as the first step in the aromatics hydrocracking pathway. These reactions are known to proceed sequentially ring by ring and are reversible under typical hydroprocessing conditions. Analysis of model compound reactivity patterns suggests that aromatic hydrogenations fall into three main categories:17 benzenic

3.2. Reaction modeling. Hydrocracking reactions are known to proceed through dual function catalysis. The acid function of a hydrocracking catalyst promotes cracking and isomerization while the metal function is responsible for hydrogenation and dehydrogenation. The acid centers are located on the catalyst support material, which is typically an amorphous silica−alumina or crystalline zeolites (X and Y zeolites). The metal phase of the catalyst is constituted of noble metals such as Pd and Pt or metal sulfides of NiMo and NiW.14 Hydrocracking reaction modeling can be approached at different levels. Mechanistic-level modeling offers the most detailed description of hydrocracking;8 however, these models are complex to implement due to the massive number of elementary steps and reaction intermediates required to build the reaction network. Alternatively, hydrocracking can be represented in terms of observable reactions without losing 140

DOI: 10.1021/acs.energyfuels.5b02084 Energy Fuels 2016, 30, 138−148

Energy & Fuels

Article

catalyst pellet are minor mainly due to the wide pore structure of the zeolite catalyst.25,26 Implementing eqs 1 and 2 for complex hydrocarbon feedstocks implies that thousands of individual rate, adsorption, and equilibrium parameters need to be derived from experimental data. This task would be completely impractical if not impossible. An alternative approach to significantly reduce these parameters to a manageable number is through quantitative structure/reactivity correlations (QSRC). QSRCs are essentially semiempirical correlations between reaction parameters and molecular properties.27,28 Application of this concept to the hydrocracking reactions of polynuclear aromatics was first introduced in a series of papers by Korre et al.17,18,29 and Korre and Klein.24 These authors developed a number of QSRCs for rate, adsorption, and equilibrium parameters on metal and acid centers from the study of model hydrocarbon compounds. These correlations enable the simulation of hydrocracking of complex feedstocks at molecular level by organizing the large reaction network into a limited number of reaction families and establishing reactivity patterns based on molecular structure. The set of QSRCs describing reactivity parameters of hydrocracking transformations is presented below. Adsorption on metal centers:

(saturation of a single aromatic ring), naphthalenic (saturation of a terminal aromatic ring), and phenanthrenic (saturation of a middle aromatic ring). Benzenic hydrogenations are the slowest, naphthalenic hydrogenations are intermediate, whereas phenanthrenic hydrogenations are the fastest. Hydrogenation on terminal rings is of particular importance as it enables further acid-catalyzed reactions. Acid center transformations of aromatic cores are believed to occur in a sequential manner.18,19 In the first step, a 6membered saturated terminal ring is contracted into a 5membered ring with a methyl group substituent. This is followed by opening of the 5-membered ring into the corresponding alkyl-aromatic. Dealkylation is the final step and consists of the complete removal of an alkyl branch from the parent ring structure. A requirement for dealkylation is that the leaving alkyl branch must have at least 3 carbons.16 Cracking on a side alkyl chain of an aromatic core and ring closure are also possible to a lesser extent,20 but are not considered in this study. Hydrocracking of naphthenes is considered to follow the same sequence.15,21 Reactions of paraffins have been thoroughly documented in the literature.22,23 The main reaction pathways comprise a sequence of skeletal isomerizations and cracking. Initially, the paraffin molecule undergoes successive isomerizations to produce a variety of mono- and multibranched structures. The resulting multibranched structures are cracked through a βscission reaction to produce smaller hydrocarbon fragments. Specific cracking mechanisms are available for each structural configuration (Types A, B1, and B2),22 as shown in Table 1. Type A involves cracking of tribranched isoparaffins into 2 isoparaffin molecules, whereas Types B1 and B2 apply to the cracking of dibranched isoparaffins into one n-paraffin and one isoparaffin. The smallest possible hydrocarbon fragments produced by A- and B-type cracking are C4 and C3, respectively. Type A has been determined to be significantly faster than Types B1 and B2 and therefore is considered as the preferential cracking route. 3.2.1. Rate equations. In this study, hydrocarbon transformations occurring on metal and acid sites are described by the following Langmuir−Hinshelwood−Hougen−Watson (LHHW) rate equations:24

rij =

riq =

⎛ n ksMijKadMi⎜C LiPHij2 − ⎝

(3)

RT

ln KeqMij = aeqM nij +

⎛ 1 ⎞ ⎜ − beqM ⎟ΔHrij ⎝ RT ⎠

(4)

Rate of aromatic hydrogenations on metal centers: ⎛ ΔHrij ⎞ ⎟⎟ ln ksMij = asM + bsM nij + csM ⎜⎜ ⎝ RT ⎠

(5)

Adsorption on acid centers: (badANAR i + cadANSCi)

ln KadA i = aadA +

⎞ ⎟ KeqMij ⎠

(6)

RT

Chemical equilibrium of ring isomerization and ring opening on acid centers: ⎛ 1 ⎞⎟ ln KeqA iq = aeqA + ⎜beqA − ΔHriq ⎝ RT ⎠

(1)

⎞ ⎟ KeqAiq ⎠ C Lq

1 + ∑i KadA iC Li

(badM NAR i + cadMNSCi)

Chemical equilibrium of aromatic hydrogenations on metal centers:

C Lj

1 + ∑i KadMiC Li ⎛ ksA iqKadA i⎜C Li − ⎝

ln KadMi = aadM +

(7)

Rate of ring isomerization: ln ksA iq = asRI +

(2)

bsRI c ⎞ ⎛ + ⎜1 − sRI ⎟ΔHC +iq ⎝ RT RT ⎠

(8)

Rate of ring opening:

where rij is the rate of a metal center reaction of compound i into product j and riq is the rate of an acid center reaction of compound i into product q; ksMij and ksAiq represent the surface rate parameters for metal and acid center transformations, respectively; CLi, CLj, and CLq are component concentrations in the liquid phase; PH2 is the hydrogen partial pressure and nij is the hydrogen stoichiometry of reaction ij; KeqMij and KeqAiq are the equilibrium constants for metal and acid center transformations, respectively; KadMi and KadAi are component adsorption coefficients on metal and acid sites, respectively. An implicit assumption is that diffusional limitations inside the

ln ksA iq = asRO +

bsRO c ⎞ ⎛ + ⎜1 − sRO ⎟ΔHC +iq ⎝ RT RT ⎠

(9)

Rate of ring dealkylation: ln ksA iq = asRD +

bsRD ΔHC +iq RT

(10)

In eqs 3 through 10, a, b, and c are QSRC parameters, NARi is the number of aromatic rings in the structure of compound i, NSCi is the number of saturated carbons in the structure of 141

DOI: 10.1021/acs.energyfuels.5b02084 Energy Fuels 2016, 30, 138−148

Energy & Fuels

Article

compound i, R is the universal gas constant, T is reaction temperature, ΔHrij and ΔHriq are the heats of reaction ij and iq, respectively, and ΔHC+iq is the heat of formation of the carbenium ion intermediate of reaction iq. The structural descriptors NARi and NSCi are obtained directly from inspecting the reactant’s structure. The thermodynamic descriptors ΔHr and ΔHC+ are estimated through semiempirical quantum chemistry calculations using AMPAC v. 9.2 software (PM3 method). Concerning the transformations of paraffins, previous studies30,31 concluded that hydrocracking and hydroisomerization reactivity increases with paraffin chain length due to the growing number of possible elementary reactions and intermediates. The effect of the number of carbon atoms of a paraffin molecule on its reactivity can be approximated through structural proportionality factors, as a measure of the contribution of all the potential transformations:32,33 k A iq = (NPCi − aiq)kRef

iq

boiling range and hydrocarbon type at various operating conditions. 3.2.2. Vapor−liquid equilibrium modeling. It is well documented that petroleum feedstocks vaporize under hydroprocessing conditions.39−41 The extent of vaporization is intensified in hydrocracking processes due to the conversion of heavy molecules into lighter and more volatile materials. This phenomenon is known to disturb reactor conditions by changing phase composition, properties, and velocities. Such conditions directly affect reaction rates since components in the vapor phase cannot access the active sites of the catalyst. Therefore, vapor−liquid equilibrium (VLE) is an important factor to consider in the modeling of hydrocracking processes. Assuming that gas and liquid phases are at equilibrium at any position of the reactor, the distribution of component i in the two phases can be expressed as follows:

CGi = C LiKVLEi (11)

(13)

where CGi is the gas phase concentration of component i, KVLEi is the VLE partition coefficient, and CG and CL are the total concentrations of the gas and liquid phases, respectively. The dependence of phase concentrations and VLE partition coefficients on reactor temperature, pressure, and H2/oil ratio is modeled with the Peng−Robinson equation of state (EoS). The key element of thermodynamic modeling of hydrogenhydrocarbon systems is obtaining the binary interaction coefficients of the Peng−Robinson EoS. To this end, in previous work39−41 flash experiments were conducted to generate VLE data for a variety of petroleum distillates, covering a wide range of temperatures and pressures. Interaction coefficients between hydrogen and hydrocarbons were then obtained by regression of the experimental data. Figure 2 illustrates the predicted and experimental vaporization of a middle and heavy distillate in the typical range of

where kAiq is the combined rate parameter (including the surface rate and adsorption) of reaction iq, NPCi is the number of carbon atoms of paraffinic compound i, aiq represents the number of carbon atoms that do not participate in reaction iq, and kRef iq is the reference rate parameter of reaction iq. Parameter aiq is fixed for each reaction regardless of paraffin chain length: aiq = 3 for isomerization of n-paraffins to monobranched paraffins, aiq = 4 for isomerization of monobranched to dibranched paraffins, aiq = 5 for isomerization of dibranched to tribranched paraffins, aiq = 7 for cracking Type A, and aiq = 6 for cracking Type B1 and B2. It follows that NPCi − aiq is the structural proportionality factor that reflects the number of reactive paraffinic carbons. The temperature dependence of paraffin rate parameters is given by the Arrhenius equation. Equilibrium ratios for paraffin isomerization reactions are calculated using the Benson’s group contribution method.34 Adsorption of paraffins is expressed as a function of a single fitting parameter aadP as follows:

ln KadA i = aadP

CG CL

(12)

It is noted that eq 12 has the form of the adsorption QSRCs (eqs 3 and 6) but without the molecular structure contribution term. The choice for this simplified equation is based on experimental evidence35,36 showing that adsorption coefficients of paraffins in liquid phase systems are practically invariant with respect to their chain length. The model comprises 24 QSRC coefficients for aromatic transformations on acid and metal centers, 6 paraffin rate parameters, and 1 paraffin adsorption coefficient. Values for the adsorption and chemical equilibrium QSRC parameters are taken directly from Korre et al.18 and Korre and Klein,24 respectively. The number of paraffin rate parameters is reduced by relating the rates of cracking Types A, B1, and B2 according to the rate ratios proposed by Möller et al.33 It is assumed that paraffin isomerization and cracking reactions have the same activation energy23,37 and its value is 31 kcal/mol.38 This strategy reduces the number of fitting parameters to 11. Initial values for rate QSRC coefficients are obtained from Klein et al.15 Estimates for paraffin rate and adsorption parameters are established using experimental data on hydrocracking of pure C16 and C28 n-paraffins.38 The values of these parameters are fine-tuned using pilot plant data on product distribution by

Figure 2. Vaporization profiles of LCO (P = 7 MPa and H2/oil = 940 NL/L) and VGO (P = 9 MPa and H2/oil = 700 NL/L).

hydroprocessing temperatures.39,41 The middle distillate is a light cycle oil (LCO) (density, 0.9360 g/cm3; boiling range, 128−431 °C; aromatics content, 83.4 wt %), whereas the heavy distillate is a raw vacuum gas oil (density, 0.8709 g/cm3; boiling range, 180−590 °C; aromatics content, 26.1 wt %). It is observed that the VLE behavior of the two feedstocks is wellpredicted with this approach, considering their differences in chemical composition and origin. 142

DOI: 10.1021/acs.energyfuels.5b02084 Energy Fuels 2016, 30, 138−148

Energy & Fuels

Article

Figure 3. Hydrocracking simulation algorithm.

3.2.3. Reactor simulation. The algorithm for the simulation of the hydrocracking process is shown in Figure 3. The set of molecules from the feed composition modeling step, operation conditions, and reaction rules represent the basis of the model input. Adsorption, equilibrium, and rate parameters for the hydrocracking reaction families are estimated with the QSRCs. Simulations start with an initial flash calculation to obtain the liquid phase concentrations at the reactor inlet. The flash calculations are performed with an in-house program designed to predict VLE in hydroprocessing reactors.42,43 Hydrocracking transformations in the space-time domain are then executed with the kinetic Monte Carlo (KMC) algorithm, as described by Pereira et al.9,10 Physical and chemical properties of the product are updated after each transformation. Computations are repeated starting from the flash calculation until reaching the final time.

Table 2. Experimental and Simulated Properties of the VGO Feed

4. RESULTS AND DISCUSSION 4.1. Simulation of the VGO feed. The composition of the VGO feedstock was simulated using a set of 1 × 10 4 computational molecules. Experimental and simulated properties of the feed mixture are given in Table 2. Boiling point distribution is presented in Figure 4. Results show that the set of computational molecules is consistent with the VGO in terms of liquid density, elemental composition, hydrocarbon types, 13C groups, and boiling point distribution. Table 2 shows that the feed sample is largely composed of saturates (64.6 wt %), consisting mainly of cycloparaffins (61.3 wt %) and a small amount of paraffins (3.2 wt %). The aromatic fraction (35.4 wt %) contains primarily mono- and diaromatics (18.6 and 9.1 wt %, respectively) plus small quantities of polyaromatics. Sulfur (0.09 wt %) and nitrogen (129 wppm) contents are low due to

Property

Experimental

Simulated

Density at 15.6 °C, g/cm3 API gravity C, wt % H, wt % Atomic H/C ratio S, wt % N, wppm SARA analysis Saturates, wt % Aromatics + Polars, wt % Hydrocarbon types Paraffins, wt % Cycloparaffins, wt % Monoaromatics, wt % Diaromatics, wt % Triaromatics, wt % Tetraaromatics, wt % Pentaaromatics, wt % 13 C Spectrum Aliphatic CH3, wt % Aliphatic CH2, wt % Aliphatic CH, wt % Aromatic CH, wt % Substituted aromatic C, wt % Bridge aromatic C, wt %

0.9115 23.7 87.48 12.60 1.72 0.09 129

0.9136 23.4 87.29 12.60 1.72 0.10 98

64.6 35.4

66.6 33.4

3.2 61.3 18.6 9.1 4.5 1.9 1.3

3.9 62.7 17.1 8.3 4.6 2.4 1.0

20.2 45.4 18.6 6.8 7.8 1.2

17.8 45.7 21.7 5.4 7.1 2.3

prior hydrotreating. Figure 4 also shows that the VGO feed contains a portion (25.8 wt %) of lighter material boiling below 343 °C as a result of mild hydrocracking reactions that occurred 143

DOI: 10.1021/acs.energyfuels.5b02084 Energy Fuels 2016, 30, 138−148

Energy & Fuels

Article

Figure 7. 343 °C+ conversion as a function of temperature with LHSV = 1.5 h−1, P = 11 MPa and H2/oil =800 NL/L.

Figure 4. Boiling point distribution for the VGO feed.

during the hydrotreating step. Figure 5 presents the overall carbon number distribution resulting from the simulation. It is

Figure 8. Product distribution as a function of LHSV with T = 390 °C, P = 11 MPa, and H2/oil = 800 NL/L. Figure 5. Simulated carbon number distribution for the VGO feed.

found that the generated molecular structures exhibit a gamma distribution within the C10−C40 carbon number range. 4.2. Simulation of the hydrocracking reactor. Reactor simulations were performed in the range of operating conditions used during the experiments. Comparison of the simulation results with the experimental data are reported in Figures 6 through 12. For the most part, the simulations compare well with the experimental data and the trends in process behavior are captured adequately.

Figure 9. Product distribution as a function of temperature with LHSV = 1.5 h−1, P = 11 MPa, and H2/oil = 800 NL/L.

Figure 6 shows how the conversion of material boiling above 343 °C evolves as a function of LHSV. It is observed that conversion increases by decreasing LHSV (i.e., longer contact time with the catalyst bed). In the typical industrial range of LHSV (0.75−1.5 h−1),14 conversion for this feedstock is found to be about 40−70 wt %. Conversion is also improved significantly by increasing reactor temperature, as evidenced by Figure 7. It is noted that conversion varies substantially between 22 and 65 wt % in a relatively narrow range of temperatures (380−400 °C).

Figure 6. 343 °C+ conversion as a function of LHSV with T = 390 °C, P = 11 MPa, and H2/oil = 800 NL/L. 144

DOI: 10.1021/acs.energyfuels.5b02084 Energy Fuels 2016, 30, 138−148

Energy & Fuels

Article

hydrocarbons boiling in the range of LGO are intermediates that are further cracked into lighter materials. Naphtha yield is found to be strongly favored by higher process severity. For instance, naphtha yield increases from 6 wt % at LHSV of 3.0 h−1 to 40 wt % at LHSV 0.75 h−1. Gas yield is below 2 wt % in all cases, showing that there is no overcracking of liquid products. Figures 10 and 11 illustrate the effect of LHSV and temperature, respectively, on hydrocarbon type distribution (paraffins, naphthenes, and aromatics). In essence, the hydrocracked product becomes more paraffinic and less aromatic at lower LHSVs and higher temperatures. This is an expected trend, since hydrogenation and cracking reactions reduce aromatics content while increasing the amount of paraffinic hydrocarbons. In this case, operating conditions have a major effect on both paraffins and aromatics contents. For example, at temperatures of 395−400 °C, aromatics content is brought down to about 20−22 wt %, whereas paraffins content is increased to 22−29 wt %. Also noted is that the profiles of paraffins are parallel to that of naphtha (Figures 8 and 9). This is a consistent trend, since cracked paraffinic molecules tend to concentrate in the naphtha boiling range. Naphthenes (cycloparaffins) exhibit little change (52−61 wt %) with respect to LHSV and temperature. This suggests that the rate of naphthenic hydrocarbon formation by full saturation of aromatic cores is comparable to the rate at which they are cracked into smaller naphthenic cores and paraffins. The molecular approach developed in this work also enables the prediction of important process parameters such as hydrogen consumption and API gravity of the liquid product, as shown respectively in Figures 12 and 13. As expected,

Figure 10. Hydrocarbon type distribution as a function of LHSV with T = 390 °C, P = 11 MPa, and H2/oil = 800 NL/L.

Figure 11. Hydrocarbon type distribution as a function of temperature with LHSV = 1.5 h−1, P = 11 MPa, and H2/oil = 800 NL/L.

Figure 13. API gravity of liquid product as a function of 343 °C+ conversion.

Figure 12. Hydrogen consumption as a function of 343 °C+ conversion.

hydrogen consumption and API gravity increase as more material boiling above 343 °C is cracked into distillates. Simulated values for hydrogen consumption are determined from the stoichiometry of the hydrogenation and cracking reactions, whereas API gravity is estimated from the individual densities of the molecules in the liquid product. Simulated profiles of various process variables along the catalyst bed are exhibited in Figures 14−16. The evolution of saturated hydrocarbons (n-paraffins, isoparaffins, and mono-, di-, tri-, and tetracycloparaffins) is given in Figure 14. Isoparaffin concentration increases steadily due to dealkylation reactions and isomerization of n-paraffins, reaching almost 25 wt % at the outlet of the reactor. n-Paraffins exhibit a relatively

The product yield distribution at different LHSVs and temperatures is presented in Figures 8 and 9, respectively. Definition of hydrocracking products is based on boiling point ranges as follows: gas (C4−), naphtha (IBP−204 °C), light gas oil (LGO) (204−343 °C), and VGO (343 °C+). Experimental values are obtained from distillation data of liquid products and gas analyses, while the simulated values result from grouping the molecules in the product by boiling range. It is seen that, in general, cracking of VGO into products, such as naphtha, LGO, and gas, increases at lower LHSVs and higher temperatures. LGO yield does not change much (30−40 wt %) over the entire range of LHSV and temperature, indicating that 145

DOI: 10.1021/acs.energyfuels.5b02084 Energy Fuels 2016, 30, 138−148

Energy & Fuels

Article

cycloparaffins. Lastly, monocycloparaffins concentration progressively increases along the reactor, showing no evident signs of further cracking into paraffins. The profiles of molecular weight and API gravity of the liquid product along the catalyst bed are presented in Figure 15. From the simulation it is clear that hydrocracking reactions considerably reduce the average molecular weight of the hydrocarbon mixture from 300 g/mol at the reactor inlet to about 190 g/mol at the outlet. At the same time, API gravity is improved from 23.4°API to 35.9°API. Figure 16 shows hydrocarbon vaporization profiles along the reactor at 380 and 395 °C. It is noted that vapor−liquid equilibrium is a strong function of process temperature. At the inlet of the reactor, about 15 and 20 wt % of the feed are released into the vapor phase at temperatures of 380 and 395 °C, respectively. Vaporization increases along the reactor as the hydrocarbon stream becomes more volatile due to hydrocracking, with such an effect being more pronounced at higher temperature. In this case, vaporization at the outlet of the reactor can be as high as 50 wt % at 395 °C and about 27 wt % at 380 °C.

Figure 14. Saturated hydrocarbon profiles at T = 395 °C, LHSV = 1.5 h−1, P = 11 MPa, and H2/oil = 800 NL/L.

5. CONCLUSIONS This study was conducted to develop an integrated composition−reaction model to simulate the vacuum gas oil hydrocracking process at the molecular level. In a first step, the composition model transforms the analytical characterization of the feed into a computational mixture of representative hydrocarbon molecules. Subsequently, the reaction model uses this molecular representation to simulate acid- and metal-catalyzed transformations occurring in the hydrocracking reactor. Reaction sets are systematized through the application of quantitative structure/reactivity correlations. A series of bench-scale hydrocracking experiments with a vacuum gas oil feedstock were conducted to collect data on processing behavior. Simulations revealed that the model is capable of generating detailed product composition and closely reproducing major process trends over a wide range of operating conditions. Moreover, the model can track the evolution of important parameters, such as hydrogen consumption, detailed hydrocarbon type distributions, vapor−liquid equilibrium, average hydrocarbon molecular weight, and physical properties of interest (e.g., API gravity of the liquid product.) Therefore, the model constitutes a robust computational tool for process development, operation, and optimization of existing units to improve their efficiency and minimize environmental impact.

Figure 15. Molecular weight and API gravity of liquid product profiles at T = 395 °C, LHSV = 1.5 h−1, P = 11 MPa, and H2/oil = 800 NL/L.

■

Figure 16. Vaporization profiles at different temperatures with LHSV = 1.5 h−1, P = 11 MPa, and H2/oil = 800 NL/L.

AUTHOR INFORMATION

Corresponding Author

*Tel.: 780-987-8763; fax: 780-987-5349; email: jinwen.chen@ canada.ca.

flat profile, meaning that their generation is counterbalanced by isomerization into isoparaffins. Cycloparaffin distribution evolves as a result of sequential cracking and to a lesser extent, hydrogenation of aromatics. It is observed that tetracycloparaffins are rapidly consumed along the reactor. This causes tricycloparaffins concentration to increase in the first 10% of the reactor (relative reactor length of 0.1), but subsequently there is a gradual decrease in tricycloparaffins concentration, indicating a shift toward the formation of dicycloparaffins. A similar but less pronounced trend is observed for dicycloparaffins, where the rate of formation from tricycloparaffins appears to compete with their conversion rate into mono-

Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Partial funding for this study was provided by Natural Resources Canada and government of Canada’s interdepartmental Program of Energy Research and Development (PERD). Comments and suggestions from Dr. Kirk Michaelian on revising the manuscript are greatly appreciated. The authors are grateful to the pilot plant and analytical lab staff at CanmetENERGY-Devon. 146

DOI: 10.1021/acs.energyfuels.5b02084 Energy Fuels 2016, 30, 138−148

Energy & Fuels

■

Article

NOMENCLATURE aadA, badA, cadA = Acid center adsorption parameters defined in eq 6 aadM, badM, cadM = Metal center adsorption parameters defined in eq 3 aadP = Paraffins adsorption parameter defined in eq 12 aeqA, beqA = Acid center equilibrium parameters defined in eq 7 aeqM, beqM = Metal center equilibrium parameters defined in eq 4 aiq = Paraffins rate parameter defined in eq 11 asM, bsM, csM = Hydrogenation surface rate parameters defined in eq 5 asRD, bsRD = Ring dealkylation surface rate parameters defined in eq 10 asRI, bsRI, csRI = Ring isomerization surface rate parameters defined in eq 8 asRO, bsRO, csRO = Ring opening surface rate parameters defined in eq 9 CG = Total concentration in the gas phase, mol/L CGi = Component concentration in the gas phase, mol/L CL = Total concentration in the liquid phase, mol/L CLi, CLj, CLq = Component concentrations in the liquid phase, mol/L KadAi = Component acid center adsorption coefficient, L/mol KadMi = Component metal center adsorption coefficient, L/ mol kAiq = Combined rate coefficient of paraffins acid center reactions, L/gcat·s KeqAiq = Equilibrium constant of acid center reactions, molj/ moli KeqMij = Equilibrium constant of metal center reactions, molj/ moli·atmnij kRef iq = Reference rate coefficient of paraffins acid center reactions, L/gcat·s ksAiq = Surface rate coefficient of acid center reactions, mol/ gcat·s ksMij = Surface rate coefficient of metal center reactions, mol/ gcat·s·atmnij KVLEi = Component vapor−liquid equilibrium partition coefficient nij = Hydrogen stoichiometry of hydrogenation reactions NARi = Number of aromatic rings NPCi = Number of paraffinic carbon atoms NSCi = Number of saturated carbons P = Reactor pressure, MPa PH2 = Hydrogen partial pressure, atm R = Universal gas constant, kcal/mol·K rij = Rate of metal center reactions, mol/gcat·s riq = Rate of acid center reactions, mol/gcat·s T = Reactor temperature, K

■

LHHW = Langmuir−Hinshelwood−Hougen−Watson LHSV = Liquid hourly space velocity LGO = Light gas oil LCO = Light cycle oil PDF = Probability distribution function QSRC = Quantitative structure/reactivity correlation VGO = Vacuum gas oil VLE = Vapor−liquid equilibrium

REFERENCES

(1) Quann, R. J. Modeling the chemistry of complex petroleum mixtures. Environ. Health Persp. 1998, 106, 1441−1448. (2) Ho, T. C. Kinetic modeling of large scale systems. Catal. Rev.: Sci. Eng. 2008, 50, 287−378. (3) Stangeland, B. E. Kinetic-model for prediction of hydrocracker yields. Ind. Eng. Chem. Process Des. Dev. 1974, 13, 71−76. (4) Froment, G. F. On fundamental kinetic equations for chemical reactions and processes. Curr. Opin. Chem. Eng. 2014, 5, 1−6. (5) Neurock, M.; Libanati, C.; Nigam, A.; Klein, M. T. Monte Carlo simulation of complex reaction systems: Molecular structure and reactivity in modelling heavy oils. Chem. Eng. Sci. 1990, 45, 2083− 2088. (6) Neurock, M.; Nigam, A.; Trauth, D.; Klein, M. T. Molecular representation of complex hydrocarbon feedstocks through efficient characterization and stochastic algorithms. Chem. Eng. Sci. 1994, 49, 4153−4177. (7) Martens, G. G.; Marin, G. B. Kinetics for hydrocracking based on structural classes: Model development and application. AIChE J. 2001, 47, 1607−1622. (8) Kumar, H.; Froment, G. F. Mechanistic modeling of the hydrocracking of complex feedstocks, such as vacuum gas oils. Ind. Eng. Chem. Res. 2007, 46, 5881−5897. (9) Pereira de Oliveira, L.; Verstraete, J. J.; Kolb, M. A Monte Carlo modeling methodology for the simulation of hydrotreating processes. Chem. Eng. J. 2012, 207−208, 94−102. (10) de Oliveira, L. P.; Verstraete, J. J.; Kolb, M. Simulating vacuum residue hydroconversion by means of Monte-Carlo techniques. Catal. Today 2014, 220−222, 208−220. (11) Gillespie, D. T. A general method for numerically simulating the stochastic time evolution of coupled chemical reactions. J. Comput. Phys. 1976, 22, 403−434. (12) Alvarez-Majmutov, A.; Chen, J.; Gieleciak, R.; Hager, R.; Heshka, N.; Salmon, S. Deriving the molecular composition of middle distillates by integrating statistical modeling with advanced hydrocarbon characterization. Energy Fuels 2014, 28, 7385−7393. (13) Alvarez-Majmutov, A.; Gieleciak, R.; Chen, J. Deriving the molecular composition of vacuum distillates by integrating statistical modeling with detailed hydrocarbon characterization. Energy Fuels 2015,10.1021/acs.energyfuels.5b02082 (14) Robinson, P. R.; Dolbear, G. E. Hydrotreating and hydrocracking: fundamentals. In Practical advances in petroleum processing Vol. 1; Hsu, C. S., Robinson, P. R., Eds.; Springer: New York, 2006; Chapter 7. (15) Klein, M. T.; Hou, G.; Bertolacini, R. J.; Broadbelt, L. J.; Kumar, A. Molecular modeling in heavy hydrocarbon conversions; CRC Press/ Taylor & Francis: Boca Raton, FL, 2006. (16) Quann, R. J.; Jaffe, S. B. Structure−oriented lumping: Describing the chemistry of complex hydrocarbon mixtures. Ind. Eng. Chem. Res. 1992, 31, 2483−2497. (17) 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, 101−117. (18) 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, 2041−2050. (19) Du, H.; Fairbridge, C.; Yang, H.; Ring, Z. The chemistry of selective ring-opening catalysts. Appl. Catal., A 2005, 294, 1−21.

Greek Symbols

ΔHC+iq = Heat of formation of carbenium ion intermediates of acid center reactions, kcal/mol ΔHrij = Heat of formation of metal center reactions, kcal/ mol ΔHriq = Heat of formation of acid center reactions, kcal/mol Abbreviations

EoS = Equation of state HPS = High pressure separator IBP = Initial boiling point KMC = Kinetic Monte Carlo 147

DOI: 10.1021/acs.energyfuels.5b02084 Energy Fuels 2016, 30, 138−148

Energy & Fuels

Article

(20) Russell, C. L.; Klein, M. T.; Quann, R. J.; Trewella, J. Catalytic hydrocracking reaction pathways, kinetics and mechanisms of nalkylbenzenes. Energy Fuels 1994, 8, 1394−1400. (21) McVicker, G. B.; Daage, M.; Touvelle, M. S.; Hudson, C. W.; Klein, D. P.; Baird, W. C.; Cook, B. R.; Chen, J. G.; Hantzer, S.; Vaughan, D. E. W.; Ellis, E. S.; Feeley, O. C. Selective ring opening of naphthenic molecules. J. Catal. 2002, 210, 137−148. (22) Martens, J. A.; Jacobs, P. A.; Weitkamp, J. Attempts to rationalize the distribution of hydrocracked products. I Qualitative description of the primary hydrocracking models of long chain paraffins in open zeolites. Appl. Catal. 1986, 20, 239−281. (23) Froment, G. F. Kinetics of the hydroisomerization and hydrocracking of paraffins on a platinum containing bifunctional Yzeolite. Catal. Today 1987, 1, 455−473. (24) Korre, S. C.; Klein, M. T. Development of temperatureindependent quantitative structure/reactivity relationships for metaland acid-catalyzed reactions. Catal. Today 1996, 31, 79−91. (25) 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, 4075−4090. (26) Landau, R. N.; Korre, S. C.; Neurock, N.; Klein, M. T.; Quann, R. J. Hydrocracking of heavy oils: Development of structure/reactivity correlations for kinetics. Am. Chem. Soc. Div. Fuel Chem. Prepr. 1992, 37, 1871−1877. (27) Mochida, I.; Yoneda, Y. Linear free energy relationships in heterogeneous catalysis. I. Dealkylation of alkylbenzenes on cracking catalysts. J. Catal. 1967, 7, 386−392. (28) Neurock, M.; Klein, M. T. Linear free energy relationships in kinetic analyses: Applications of quantum chemistry. Polycyclic Aromat. Compd. 1993, 3, 231−246. (29) 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, 4191−4210. (30) Sie, S. T.; Senden, M. M. G.; van Wechem, H. M. H. Conversion of natural gas to transportation fuels via the Shell Middle Distillate Synthesis Processes (SMDS). Catal. Today 1991, 8, 371− 394. (31) Denayer, J. F.; Baron, G. V.; Souverijns, W.; Martens, J. A.; Jacobs, P. A. Hydrocracking of n-alkane mixtures on Pt/H-Y zeolite: Chain length dependence of the adsorption and the kinetic constants. Ind. Eng. Chem. Res. 1997, 36, 3242−3247. (32) Sie, S. T. Acid-catalyzed cracking of paraffinic hydrocarbons. 3. Evidence for the protonated cyclopropane mechanism from hydrocracking/hydroisomerization experiments. Ind. Eng. Chem. Res. 1993, 32, 403−408. (33) Möller, K.; le Grange, P.; Acolla, C. A two-phase reactor model for the hydrocracking of Fischer−Tropsch-derived wax. Ind. Eng. Chem. Res. 2009, 48, 3791−3801. (34) Benson, S. W. Thermochemical Kinetics; John Wiley & Sons: New York, 1976. (35) Laxmi Narasimhan, C. S.; Thybaut, J. W.; Martens, J. A.; Jacobs, P. A.; Denayer, J. F.; Marin, G. B. A unified single-event microkinetic model for alkane hydroconversion in different aggregation states on Pt/H-USY-Zeolites. J. Phys. Chem. B 2006, 110, 6750−6758. (36) Li, B.; Calemma, V.; Gambaro, C.; Baron, G. V.; Denayer, J. F. Competitive adsorption of C20−C36 linear paraffins on the amorphous microporous silica−alumina ERS-8 in vapor phase and liquid phase. Ind. Eng. Chem. Res. 2010, 49, 7541−7549. (37) Schweitzer, J. M.; Galtier, P.; Schweich, D. A single event kinetic model for the hydrocracking of paraffins in a three-phase reactor. Chem. Eng. Sci. 1999, 54, 2441−2452. (38) Rossetti, I.; Gambaro, C.; Calemma, V. Hydrocracking of long chain linear paraffins. Chem. Eng. J. 2009, 154, 295−301. (39) Chen, J.; Jiang, W.; Yang, H.; Hawkins, R.; Ring, Z. Comparative Study of Vapor-Liquid Equilibrium during Hydroprocessing of Different Petroleum Feedstocks. Proceedings of 2006 AIChE National Spring Meeting, April 23−27, 2006, Orlando, FL, USA.

(40) Chen, J.; Wang, N.; Mederos, F.; Ancheyta, J. Vapor-Liquid Equilibrium Study in Trickle-Bed Reactors. Ind. Eng. Chem. Res. 2009, 48, 1096−1106. (41) Munteanu, M. C.; Chen, J. Vapor-Liquid Equilibrium (VLE)Based Modeling for the Prediction of Operating Regimes in a Heavy Gas Oil Hydrotreater. Energy Fuels 2012, 26, 1230−1236. (42) Chen, J.; Mulgundmath, V.; Wang, N. Accounting for VaporLiquid Equilibrium in the Modeling and Simulation of a Commercial Hydrotreating Reactor. Ind. Eng. Chem. Res. 2011, 50, 1571−1579. (43) Alvarez-Majmutov, A.; Chen, J. Modeling and simulation of a multibed hydrotreater with vapor-liquid equilibrium. Ind. Eng. Chem. Res. 2014, 53, 10566−10575.

148

DOI: 10.1021/acs.energyfuels.5b02084 Energy Fuels 2016, 30, 138−148