Stochastic Modeling and Simulation Approach for Industrial Fixed-Bed

Subscriber access provided by Binghamton University | Libraries

Article

Stochastic Modeling and Simulation Approach for Industrial Fixed-Bed Hydrocrackers Anton Alvarez-Majmutov, and Jinwen Chen Ind. Eng. Chem. Res., Just Accepted Manuscript • Publication Date (Web): 31 May 2017 Downloaded from http://pubs.acs.org on June 5, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Industrial & Engineering Chemistry Research

Stochastic Modeling and Simulation Approach for Industrial Fixed-Bed Hydrocrackers Anton Alvarez-Majmutov* and Jinwen Chen Natural Resources Canada, CanmetENERGY-Devon, One Oil Patch Drive, Devon, AB, T9G 1A8, Canada © Her Majesty the Queen in Right of Canada, as represented by the Minister of Natural Resources, 2017.

Keywords Stochastic simulation; vacuum gas oil; hydrocracking; adiabatic reactor

Abstract In this study, we developed a comprehensive model to simulate vacuum gas oil hydrocracking. The model incorporates knowledge of feedstock composition and reaction chemistry at the molecular level and is tuned with pilot plant data. The molecular makeup of the feedstock is derived from its analytical characterization via statistical modeling. Hydrocracking reactivity is described based on the molecular structure of the reactants. Progression of chemical reactions is simulated using a kinetic Monte Carlo algorithm. A stochastic heat balance equation is integrated into the model to represent adiabatic operation. Vapor-liquid equilibrium calculations are executed in parallel with an in-house flash program. Hydrocracker simulations were conducted covering the typical range of conversion levels observed in commercial operation. The model generates detailed information on product distribution as well as product quality. In addition, the model is able to track the evolution of hydrogen consumption, hydrocarbon vaporization, and temperature along the reactor.

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

1. Introduction

The configuration of petroleum refineries nowadays is largely driven by market requirements. Hydrocracking has become a major conversion technology for the production of high-quality fuels and lubricants. Commercial hydrocracking units are able to take feedstocks comprised of a variety of difficult streams, including heavy vacuum gas oil (VGO), light cycle oil, coker gas oil, and deasphalted oil. Depending on the refining objectives and product specifications, hydrocracker designs can have one or two reaction stages, with once-through or recycle flow modes.1 The exothermic nature of hydrocracking reactions makes heat management an important aspect of the process. Large amounts of heat are produced through the chemical addition of hydrogen during the saturation of aromatic hydrocarbons and stabilization of cracking products.2 This heat accumulates in the reacting stream, resulting in sharp evolution of temperature along the catalyst bed. Maintaining the temperature rise within specific boundaries is crucial for optimal operation, as too-high temperatures can lead to loss of product selectivity, higher hydrogen consumption, and coke build-up on the catalyst surface.3 In commercial practice, hydrocracking reactors are designed with multiple catalyst beds for the purpose of injecting quench gas in between, thereby reducing temperature rise.4 Hydrocracking reactor design is an elaborate undertaking requiring reliable models with the capability to predict not only product composition and properties, but also the behavior of operational parameters. Because of the complexity of hydrocracking chemistry and feedstock composition, process models often rely on simplified reaction schemes formulated in terms of broad chemical lumps defined by bulk property measurements.5 To simulate adiabatic reactor operation, it is also common to use a lumped heat of reaction in the heat balance equation.6 This term is essentially a model parameter obtained by regression of plant data7 and is considered an approximate average of the net heat release of all reactions. A few examples of adiabatic hydrocracker simulations using this approach are available in the literature.8–11 The global heat of reaction reported in these studies is about 10 kcal per mole of H2 consumed, which corresponds to the range of normal alkane hydrocracking.2 Since the total heat evolved is determined by the extent of hydrogenation 2 ACS Paragon Plus Environment

Page 2 of 47

Page 3 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


of the various hydrocarbon classes in the reacting mixture, models without detailed process chemistry have limited capacity to respond to variations in feedstock reactivity and operating conditions.12 Furthermore, other factors that affect hydrocracker performance but are usually ignored, for example vapor-liquid equilibrium (VLE), are needed to expand model predictive capabilities.13 Molecular-level modeling of refinery processes enhances simulation capability through more fundamental descriptions of the process. To overcome the difficulties of petroleum compositional analysis, composition modeling techniques have been developed and applied to derive consistent representations of refinery feedstocks based on thousands of molecules. These are subsequently used to model chemical reactions.14 The resulting reaction network in most cases is of considerable size, making its implementation through deterministic rate equations quite challenging. As an alternative, stochastic simulation algorithms have been proved to be capable of facilitating the task of dealing with complex chemical processes.15 Detailed composition-reaction models for various refinery hydroprocesses have recently been developed following this approach. 16–18

While the methodology has demonstrated sufficiently robust capabilities, the focus

so far has been limited to bench-scale isothermal reactors. In previous work18 we developed a molecular-level process model for vacuum gas oil hydrocracking using experimental data from a bench-scale reactor. The model accounts for acid- and metal-catalyzed reactions, chemical equilibrium, adsorption on catalytic sites, and VLE, but is restricted to isothermal operation. In the present study, we extend this model by introducing heat release kinetics to simulate adiabatic hydrocracking reactor operation with inter-bed gas quenching.

2. Modeling Approach

The developed modeling approach comprises two main components: feed composition modeling and process modeling. The first component consists of deriving a consistent molecular representation that mimics the chemical composition of the hydrocracker feedstock. The second component is for simulating chemical reactions at



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the molecular level, phase equilibrium, and thermal behavior of the hydrocracking reactor. The various elements of the hydrocracker model are described below.

2.1. Feed composition modeling The purpose of the composition model is to simulate the molecular composition of complex oil fractions by systematically generating virtual mixtures of representative molecules. The basic concept behind this approach is that the massive number of hydrocarbon molecules making up petroleum can be broken down into a limited number of structural fragments such as alkyl chains, benzene rings, cyclohexane rings, and heterocyclic rings, and that the frequency of occurrence of each of these fragments can be modeled using continuous statistical distributions.19 Thus, a wide spectrum of hydrocarbon molecules can be assembled from a limited number of fragments by computational sampling of their statistical distributions. The use of continuous distribution functions to represent oil composition has its origin in Boduszynski’s model, 20,21

in which petroleum composition is a continuum in molecular weight and structure. In the present work, compositional simulation of VGO feedstock is performed

using a methodology described in detail in our previous publications.22,23 Briefly, the VGO is represented by paraffins, cycloparaffins, aromatics, and their heterocyclic homologues (e.g. organosulfur and organonitrogen compounds). Specific building pathways and chemical rules are implemented to determine the molecular structure of each hydrocarbon class. For example, aromatics are constructed by first setting the number of benzene and cyclohexane rings and then establishing the number and length of the alkyl branches that are attached to the aromatic core structure. This process is executed through Monte Carlo sampling of a set of probability functions describing the distribution of each structural fragment or feature. The simulated feedstock is built on a large set (104) of computational molecules and optimized by applying simulated annealing and entropy maximization methods to match the available analytical data. The methodology was developed and validated by using middle distillate (204– 343°C) and vacuum distillate (343–550°C) samples having diverse geographic origins and processing histories. Input for the simulations was obtained from routine and advanced characterization techniques: liquid density, elemental analysis, simulated 4 ACS Paragon Plus Environment

Page 4 of 47

Page 5 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


distillation,

13

C NMR spectroscopy, SARA fractionation, and hydrocarbon type and

sulfur compound speciation by two-dimensional gas chromatography. Table 1 provides an example of the application of the composition model to one VGO feedstock. It is observed that the simulated VGO compares well with the real sample in terms of physical and chemical properties. The molecular representation contains sufficient information on the composition and chemical structure of the feed to enable detailed process modeling.

2.2. Reaction modeling The set of molecules describing the VGO feedstock forms the basis to simulate the hydrocracking process. Instead of using deterministic equations to simulate the progression of chemical reactions, a variable time step kinetic Monte Carlo algorithm is employed.16,24 By this method, reactions of the discrete population of VGO molecules are simulated one at a time, meaning that through each reaction step one reactant molecule is replaced by new product molecules. For example, the change of unimolecular reaction system “a → b” from the current reaction step s to the next one (s+1) is executed as follows:

= − 1

(1)

= + 1

(2)

where ma and mb are the number of reactant “a” molecules and product “b” molecules, respectively. The reaction event to be executed is selected by random sampling of a cumulative probability distribution, which is formed by accumulating the rates of all potential reactions. To this end, it is necessary to define the possible reactions (e.g., hydrogenation, cracking, and isomerization) for each molecule present in the population. After executing the selected reaction, a time increment is set to advance the reaction time. The time increment (∆ts) between the current reaction step s and the next one is defined as: ∆ = −

ln() ∑

(3)

where rn is a random number uniformly distributed between 0 and 1 and ∑ is the summation of all the potential reaction probabilities, or in this case reaction rates ri. Since

the molecular population evolves with each reaction step, a new distribution of potential 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 47

reaction probabilities is generated each time to define the corresponding reaction event and time increment. To simulate hydrocracking reactions in the axial reactor coordinate, an appropriate form of equation 3 is formulated. Assuming that the hydrocracking reactor behaves as a pseudo-homogeneous plug-flow reactor, the mass continuity equation for hydrocarbon component i is:

=

(4)

where CLi is the component concentration in the liquid phase, τ is space time, rik is the component rate to form hydrocracking product k, and ρb is the catalyst bulk density. Space time can be defined in terms of liquid phase velocity (uL) and axial reactor distance (z) as follows: =

(5)

Considering that the time coordinate used in equation 3 is in essence equivalent to space time in a plug-flow reactor, the stochastic axial position increment (∆zs) between the current reaction step s and the next one becomes: ∆ = −

ln() 1 ∑ ∑

(6)

The stochastic approach is well-suited for handling large-scale systems where numerous molecular species participate in many highly coupled chemical reactions. A particular advantage over deterministic modeling is that the hydrocracking reaction network is generated while reactions are moving forward.16 However, we must note that since the model is derived from a pseudo-homogeneous plug-flow reactor equation, the effects of external mass transfer limitations are excluded by default. This means that this aspect of the process is implicit in the model parameters. On the other hand, given that the pilot plant unit used in this study operates in up-flow mode (section 3), external mass transfer is expected to improve, particularly at the liquid-solid interphase.25 Another related assumption in the model is full catalyst wetting. This is justified by the fact that up-flow operation provides complete wetting of the catalyst bed. Therefore, it can be considered that correction by wetting efficiency to represent commercial-scale reactor operation under full catalyst wetting conditions is not necessary. Likewise, it is implicitly 6 ACS Paragon Plus Environment

Page 7 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


assumed that diffusional limitations inside the catalyst particle are minimal owing to the wide pore structure of the zeolite catalyst matrix.26,27 It should be noted that pilot plant experiments were conducted using commercial size catalyst particles, meaning that the measured kinetics reflect catalyst performance in its working state including any possible diffusional effects.

2.2.1. Reaction kinetics Hydrocracking reactions are promoted in the presence of dual-functional catalysts. The metal phase catalyzes hydrogenation and dehydrogenation reactions, while the acid function of the catalyst governs cracking and isomerization. Figure 1 illustrates the possible reaction pathways for aromatics hydrocracking, with naphthalene as the representative compound.28–31 Naphthalene is first saturated on a terminal ring via naphthalenic-type hydrogenation to form tetralin. The saturated ring of the tetralin molecule is subsequently isomerized into a five-membered ring with a methyl group substituent, followed by ring opening into butylbenzene. The final step is the removal of the butyl branch attached to the benzene ring. The tetralin molecule from the initial step can also be saturated into decalin through a benzenic-type hydrogenation, which then follows a similar sequence of isomerization and cracking steps. In an analogous way, paraffins undergo various skeletal isomerizations to form mono- and multibranched structures before being cracked into smaller fragments via β-scissions (Types A, B1, and B2), as shown in Figure 2.32,33 Given the recursive pattern of the hydrocracking network, reactions can be efficiently organized into families according to the types of metal- and acid-catalyzed transformations34: aromatic ring hydrogenations (benzenic, naphthalenic, and phenanthrenic)28, paraffin isomerization, paraffin cracking, ring isomerization, ring opening, and ring dealkylation. The rate of hydrocracking reactions on metal and acid sites is formulated in terms of the Langmuir-Hinshelwood-Hougen-Watson (LHHW) formalism as follows:35 Metal sites:

=

! "# $%" & '() *# − $ # 1 + ∑ $%"

+,"#


(7)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Acid sites:

, =

! ./ $%. & − $ / +,./

Page 8 of 47

(8)

1 + ∑ $%.

where ksMij and ksAiq are the rate constants for metal and acid site reactions, respectively;

CLi, CLj, CLq are component concentrations in the liquid phase; '() is hydrogen partial

pressure and nij is the hydrogen stoichiometry; KeqMij and KeqAiq are chemical equilibrium constants for metal and acid site reactions, respectively; KadMi and KadAi are the metal and acid site adsorption equilibrium constants, respectively. Considering the size of the hydrocracking reaction network, applying Equations 7 and 8 involves an overwhelming number of model parameters (rate constants, adsorption and chemical equilibrium parameters). Strategic simplifications are needed to reduce the number of model parameters to a manageable level. A consistent way of doing this is to use quantitative structure/reactivity correlations (QSRC) whereby reactivity parameters are correlated with molecular descriptors pertaining to each reaction family.36,37 Korre et al.29,38 developed robust QSRCs for hydrocracking of aromatic hydrocarbons in a comprehensive study using model compounds. These correlations enable systematic descriptions of component reaction rate, adsorption, and equilibrium parameters for metal and acid sites. Table 2 lists the QSRCs for the hydrocracking reaction families. As observed, reactivity parameters are described by using molecular structure descriptors and thermodynamic quantities that can be estimated through semi-empirical quantum chemistry calculations (e.g. AMPAC software). The expressions used for paraffin reactions are given in Table 3. Paraffin reactivity is assumed to increase as a function of paraffin chain length.39,40 Structural proportionality factors (NPCi - aiq) indicating the number of reactive paraffinic carbons are introduced for approximating isomerization and cracking rate parameters.41 The rates of B1- and B2-type cracking are related to that of type A using factors derived from experimental data.42 Isomerization and cracking reactions are assumed to have equal activation energies.33,43 Adsorption on acid sites is considered to be invariant with respect to paraffin chain length, as suggested by experimental evidence44,45, and therefore is represented by a single fitting coefficient. Paraffin isomerization equilibrium coefficients are estimated using Benson’s group contribution method. Reference rate and adsorption 8 ACS Paragon Plus Environment

Page 9 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


parameters are obtained from model compound studies43 of representative long-chain paraffins.

2.2.2. Vapor-liquid equilibrium Hydrocarbon vaporization is a key factor influencing conversion of reactants, fluid dynamics, and the heat balance in a hydroprocessing reactor. Flash calculations are therefore necessary to quantify the extent of feed vaporization in the reactor as well as the changes in phase compositions and properties. Assuming that the vapor and liquid phases achieve equilibrium at any position (i.e. reaction step) inside the hydrocracking reactor, the distribution of hydrocarbon components and hydrogen between the vapor and liquid phases is given by: 1 = $34 2

(9)

where yi and xi are component mole fractions in the gas and liquid phases, respectively,

and $34 is the component VLE partition coefficient. Mole fractions in each phase are

defined in terms of the discrete population of reacting molecules as follows: 1 = 2 =

5,

∑ 5,

(10)

,

∑ ,

(11)

where mGi,s and mLi,s are the number of component i molecules in the gas and liquid phases, respectively, at current reaction step s. Component liquid phase concentrations defined in rate equations 7 and 8 are evaluated using the following expression: = 1

7

(12)

where ρL is the liquid phase density and ML is the liquid phase molecular weight. Calculation of VLE partition coefficients for hydrocarbons and hydrogen and phase compositions at reactor conditions is done by solving the Peng-Robinson equation of state (EoS). Computations are executed using an in-house flash calculation program designed specifically for hydrogen-petroleum systems. The critical properties of hydrocarbons and other related parameters used in the EoS are estimated from their 9 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 47

boiling point, density, and molecular weight. For hydrogen, such information is available the literature. Liquid and gas phase thermophysical properties, such as enthalpies and heat capacities, are estimated by means of various correlations. The key feature of the flash program is that the binary interaction coefficients between hydrogen and hydrocarbons of the EoS were determined from experimental VLE data sets for a variety of middle and heavy distillates over wide ranges of temperatures and pressures. A number of previous studies documenting the development of this program and its application to hydroprocessing reactor modeling are available in the literature.13,46–48

2.2.3. Heat balance Hydrocracking reactions are exothermic overall since the heat generated by the hydrogenation of aromatic rings and unsaturated cracking products exceeds the heat consumed by C-C bond scission during cracking reactions.2 The consequent temperature rise along the hydrocracking catalyst bed can be modeled with the following pseudohomogeneous energy balance equation for fixed-bed adiabatic reactors: 895 : 5 + 9 : ;

< @3 < + =3 >? − 8−ΔBC*D ; = 0 <

(13)

where FL and FG are the liquid- and gas-phase superficial mass velocities, respectively,

: and : 5 are the liquid- and gas-phase heat capacities, respectively, T is reactor temperature, λV is the heat of vaporization, G0 is the superficial mass velocity of the

hydrocarbon feed, φV is the fraction of vaporized hydrocarbons, and ΔBC*D is the heat

release of hydrocracking reaction ik, which comprises any possible metal and acid center reaction. The first term in equation 13 represents convective transport, the second term is the heat consumed by vaporization, and the third term is net heat generation. Heat capacities of each phase, heats of vaporization, and the vaporized fraction at process conditions are calculated by the VLE program described in section 2.2.2. Heats of reaction conventionally are estimated by thermodynamic calculations using the molecular structure of reactants and products. However, when taking into account the number of reactions involved, it is preferable to reduce the problem in a systematic manner. Jaffe2 identified from various sets of thermochemical data of model hydrocracking reactions49 10 ACS Paragon Plus Environment

Page 11 of 47

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


that the heat generated per mole of hydrogen consumed is practically constant throughout distinct reaction families, namely paraffin cracking, ring opening, aromatics hydrogenation, and olefin saturation. Also, isomerization reactions were determined to have a minor contribution, as there is no hydrogen addition involved. Such a pattern indicates that the heat released in the hydrocracking process can be efficiently described using representative values for each reaction family that consumes hydrogen. The heat release data for each reaction family considered in the model is listed in Table 4 (from Stull et al.49). Reported values are relative to the moles of hydrogen consumed in each reaction family, hence the hydrogen stoichiometric coefficients (nik) in the third term of equation 13. Also noted is that the heat release of paraffin cracking, ring opening, and dealkylation is an overall pathway value considering cracking and product saturation. Deterministic equation 13 is reformulated into the corresponding stochastic expression. The following approximation is proposed to calculate the temperature step increase (∆Ts) between reaction steps: ∆< = F∆ G

8−ΔBC*D ;

H95 : 5 + 9 : + =3 >?

@3 I

Stochastic Modeling and Simulation Approach for Industrial Fixed-Bed

Recommend Documents