Simulating vacuum residue hydroconversion by means of Monte-Carlo techniques

Catalysis Today 220–222 (2014) 208–220

Contents lists available at ScienceDirect

Catalysis Today journal homepage: www.elsevier.com/locate/cattod

Simulating vacuum residue hydroconversion by means of Monte-Carlo techniques Luís Pereira de Oliveira a,b , Jan J. Verstraete a,∗ , Max Kolb b a b

IFP Energies nouvelles, Rond-point de l’échangeur de Solaize, BP 3, 69360 Solaize, France Laboratoire de Chimie, École Normale Supérieure, 69364 Lyon Cedex 07, France

a r t i c l e

i n f o

Article history: Received 4 April 2013 Received in revised form 12 August 2013 Accepted 20 August 2013 Available online 23 October 2013 Keywords: Composition modeling Kinetic modeling Information entropy maximization Vacuum residue hydrocracking Stochastic Simulation Algorithm Kinetic Monte-Carlo

a b s t r a c t The present work focuses on the development of a novel methodology for the kinetic modeling of heavy oil conversion processes. The methodology models both the feedstock composition and the process reactions at a molecular level. The composition modeling consists of generating a set of molecules whose properties are close to those of the process feedstock analyses. This synthetic mixture of molecules is generated by a two-step molecular reconstruction algorithm. In its first step, an equimolar set of molecules is built by assembling structural blocks in a stochastic manner. In the second step, the mole fractions of the molecules are adjusted by maximizing an information entropy criterion. Once the composition of the feedstock is represented, the conversion process is simulated by applying, event by event, its main reactions to the set of molecules by means of a kinetic Monte Carlo (kMC) method. The methodology has been applied to hydroconversion of Ural vacuum residue and both the feed and the predicted effluents were favorably compared to the experimental yield pattern. © 2013 Elsevier B.V. All rights reserved.

1. Introduction Over the last decade, the world demand for high-quality lowboiling products such as gasoline, jet fuel and diesel has continually increased, while at the same time the available crude oils have become increasingly heavier. Both trends boost the importance of refining processes that are able to convert heavy petroleum fractions, such as vacuum residues, into lighter and more valuable clean products [1–3]. Petroleum residue conversion processes, such as residue hydrocracking or residue fluid catalytic cracking (RFCC), are based on the degradation of the largest molecules by thermal and/or catalytic cracking reactions at high temperature. To accurately predict the process performances, reliable

Abbreviations: AEBP, Atmospheric Equivalent Boiling Point; BT, benzothiophene; CME, chemical master equation; DBT, dibenzothiophene; GC, Gas Chromatography; GO, gas oil; HDA, hydrogenation of aromatic rings (hydrodearomatization); HDM, hydrodemetallization; HDN, hydrodenitrogenation; HDS, hydrodesulfurization; kMC, kinetic Monte Carlo; LAD, least absolute deviations; LCO, Light Cycle Oil; LFER, linear free energy relationships; LSQ, least squares; MS, mass spectrometry; MD, molecular discretization; NMR, nuclear magnetic resonance; PDF, probability distribution function; QS/RC, Quantitative Structure/Reactivity Correlation; REM, Reconstruction by Entropy Maximization; SARA, saturates–aromatics–resins–asphaltenes analysis; SR, stochastic reconstruction; SR-REM, coupled SR-REM method; SSA, Stochastic Simulation Algorithm; T, thiophene; VGO, vacuum gas oil; VR, vacuum residue. ∗ Corresponding author. Tel.: +33 04 37 70 25 07; fax: +33 04 37 70 20 08. E-mail address: [email protected] (J.J. Verstraete). 0920-5861/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cattod.2013.08.011

kinetic models are needed [4–6]. Kinetic models for heavy oil conversion processes have classically been based on a lumping strategy [7], in which molecular components are grouped into several chemical families, according to their global properties (boiling point, solubility, etc.). To improve the predictions of these lumped models, the number of lumps is generally increased, as exemplified by literature kinetic models for residue hydroprocessing [8–20]. Besides their feed dependence and lack of robustness, such models also seem to reach their limits due to very high number of reaction parameters (rate constants, adsorption parameters) that need to be identified when applied to the heavy oil conversion. The development of more detailed kinetic models containing molecule-based reaction pathways is therefore required [5,21] and has been ongoing over the last 2 to 3 decades for a wide variety of chemical processes, such as pyrolysis, steam cracking, catalytic cracking, hydrocracking, hydrotreating, catalytic reforming, alkylation, and many more [22]. Such models expect a molecular description of the feedstock, however. Unfortunately, even though the most advanced analytical techniques allow to identify a large number of compounds and classes of chemical families, the complete and quantitative molecular detail of heavy feedstocks still remains unknown [23–27]. Our work presents a novel two-step kinetic modeling methodology for heavy oil conversion processes that retains the molecular detail throughout the method (Fig. 1). In the proposed approach, the lack of molecular detail of the petroleum fractions is overcome by using a carefully selected synthetic mixture of representative molecules, the

L.P. de Oliveira et al. / Catalysis Today 220–222 (2014) 208–220

209

2. Description of the methodology a␯ Probability of reaction ␯ (s−1 ) c␯ Stochastic rate parameter of reaction ␯ cDehydro Stochastic rate parameter of a dehydrogenation reaction (s−1 ) Stochastic rate parameter of a hydrogenation reaccHydro tion (s−1 atm−n ) DR Normalized cumulative probability distribution for all possible reactions E(xi ) Shannon entropy criterion Ea Activation energy (J mol−1 ) Fj () Mixing rule for analytical property j Frep Replication factor h␯ Distinct combination of the reactant molecules of the reaction ␯ J Total number of analytical constraints in the REM method Equilibrium constant (–) Keq k␯ Deterministic rate parameter of the reaction ␯ Number of possible reactions M MW Molecular weight (g mol−1 ) N Number of molecules present in the predefined set of molecules Number of analyses present in the objective funcNA tion NM,i Number of measurements of analysis i NAR Number of aromatic rings Number of saturated rings NSR NTR Number of thiophenes fused to the aromatic ring nH Reaction stoichiometry with respect to molecular hydrogen Objective function OF pH2 Partial hydrogen pressure (atm) Value of analytical property j Pj P␯ Normalized probability of reaction ␯ Ideal gas constant (J K−1 mol−1 ) R RNi Random number i Vector of mole fractions xi x xi Mole fraction of molecule i (–) exp Xi,j Experimental value of measure j of analysis i calc Xi,j

ıi

HR0 HR

Calculated value of measure j of the analysis i Deviation between the calculated and experimental values of analysis i Heat of reaction at 25 ◦ C (J mol−1 ) Heat of reaction at reaction temperature (J mol−1 )

Modeling objective Feed

Product

Molecular Reconstruction

Property Calculations

Feed molecules

Molecule-based reaction and reactor model

Product molecules

Fig. 1. Molecule-based process modeling methodology (after [6]).

properties of which correspond closely to the available feedstock analyses. The conversion process is then simulated by transforming each molecule of this set by means of a kinetic Monte Carlo (kMC) method.

As mentioned in the introduction, the proposed methodology comprises two main steps. In the first step of the methodology, the feedstock composition is modeled by means of a set of molecules whose mixture properties are close to the process feedstock analyses. The set of molecules is generated using the SR-REM molecular reconstruction algorithm [28,29]. This approach results from the coupling of two methods: Stochastic Reconstruction (SR) [28,30–33] and Reconstruction by Entropy Maximization (REM) [28,34–38]. The SR algorithm generates an equimolar set of molecules by assembling structural blocks in a stochastic manner, while the REM algorithm adjusts the mole fractions of the set of molecules by maximizing a Shannon information entropy criterion. During previous work, this algorithm was first developed and validated for Light Cycle Oils (LCO) [28,33,39] and for vacuum gas oils (VGO) [40,41]. The SR-REM algorithm was also applied to various petroleum vacuum residues by Verstraete et al. [42] and de Oliveira et al. [29]. In the present work, the reconstruction of an Ural vacuum residue (VR) will be illustrated. In the second step of the methodology, the kinetics of the conversion process is modeled using a kinetic Monte Carlo (kMC) method, proposed by Gillespie [43] and termed Stochastic Simulation Algorithm (SSA). In this approach, the SSA identifies all possible reactions at each instant in time. For each reaction, its rate coefficient determines its reaction probability. Based on these probabilistic considerations, the algorithm then selects a reaction event, transforms the reactant molecule into its corresponding product, and increments the simulation time. The simulation proceeds event by event (reaction by reaction) until the final simulation time is reached, leading to the set of product molecules. In the present work, the hydroconversion of an Ural VR in a batch reactor will be simulated and compared to the experimental data. 2.1. Composition modeling 2.1.1. Stochastic reconstruction (SR) method The goal of the SR method is to create a representative ensemble of probability distribution functions (PDF) for molecular structural attributes (e.g. type of molecule, number of cores, number of aromatic rings, number of side chains, etc.) from available analytical information of the given petroleum fraction. Once the PDFs have been created, they are sampled to generate an equimolar set of molecules whose mixture properties are close to those of the given petroleum fraction. The transformation of the analytical information into PDFs of molecular attributes is carried out via an iterative procedure shown in Fig. 2. First, appropriate molecular attributes are chosen from the available analytical information (elemental analysis, 13 C NMR, GC–MS, etc.) of the feedstock and from expert knowledge. A PDF is selected to represent each molecular attribute. The PDFs employed may either be discrete, such as histograms, or continuous, such as normal distribution functions, gamma distribution functions, or exponential distribution functions. The choice of the type of PDF is based on its form, its flexibility and its number of parameters. A flexible PDF is desirable in order to easily pick up the correct form of the experimental observations, as well as a small number of parameters to characterize it. After the appropriate PDFs are found for each molecular attribute, they are sampled via a Monte Carlo procedure to determine the functional and structural attributes of a molecule. The sequence of the PDF sampling steps is defined through a building diagram. The selected molecular attributes are randomly assembled, yet respecting chemical rules. These rules avoid the creation of impossible and improbable molecules on chemical, thermodynamic or probabilistic grounds.

210

L.P. de Oliveira et al. / Catalysis Today 220–222 (2014) 208–220

Analytical data

Building diagram

Set of PDFs Modified PDF parameters

N times

Sampling of PDFs and molecule assembly Calculation of pure component properties

Minimization using a robust global optimizer

Calculation of the mixture properties

No

Objective function

Converged?

Yes

Database of molecules

Fig. 2. Flow diagram of the stochastic reconstruction (SR) algorithm.

Molecule construction is repeated N times so as to obtain a mixture of molecules. The pure component properties of each molecule are calculated, either directly, by inspection of the structure (e.g. for chemical formula, molecular weight, elemental analysis, NMR analysis, etc.), or numerically, by group contribution methods or correlations (e.g. for density, boiling point distribution, SARA fractions, CCR content, etc.). The mixture properties are calculated from the pure component properties using linear mixing rules and are compared to the available analyses through an objective function. To compare the properties of this synthetic mixture to those of the actual feedstock, one first calculates the difference between the two. In our method, this can either be done by using absolute deviation errors (Eq. (1)), relative errors (Eq. (2)), squared errors (Eq. (3)), or squared relative errors (Eq. (4)), as shown below: NM,i





1   exp calc  ıi = Xj,i − Xj,i  NM,i

(1)

j=1





M,i  calc  1   Xj,i − Xj,i  ıi =   exp NM,i  Xj,i 

N

exp

(2)

j=1

NM,i

ıi =



1  exp calc Xj,i − Xj,i NM,i

2 (3)

j=1

NM,i

1  ıi = NM,i j=1



exp

calc Xj,i − Xj,i

2 (4)

exp

Xj,i

exp

where ıi is the deviation of analysis i, Xi,j is the experimental value calc Xi,j

of measure j of analysis i; is the calculated value of measure j of analysis i, and NM,i is the number of measurements of analysis i (e.g. each point of the simulated distillation or each atom content in the elemental analysis). Once the deviations ıi have been calculated for each analysis i, the objective function (OF) can be built. Various expressions are available in the literature to construct objective functions, such as the least squares function (LSQ) or L2 norm, the least absolute deviations (LAD) or L1 norm, Huber’s method, Andrew’s sine, Tukey’s

biweight, etc. [44]. In our method, the objective function is calculated as follows: 1  ıi NA NA

OF =

(5)

i=1

where NA is the number of analyses present in the objective function and ıi is the deviation of analysis i, which is obtained by one of the equations above (Eqs. (1)–(4)). Depending on the choice of the latter, minimizations can therefore be performed using either a least squares based objective function or a least absolute deviations based objective function. The value of the objective function is then iteratively minimized via an optimizer algorithm that modifies the parameters of the PDFs for the structural attributes. Although a large spectrum of optimizers exists [44–51], the stochastic nature of the Monte Carlo sampling process requires a not only a derivative-free optimizer, but also an optimizer that withstands noise on the objective function value. For these cases, robust global optimizers, such as simulated annealing [48], genetic algorithms [49,50] or particle swarm optimization [51], are well adapted and have therefore been implemented. Parameter estimation is continued until the molecular representation accurately matches the properties of the feedstock. 2.1.2. Reconstruction by entropy maximization (REM) method Now that the actual process feedstock has been replaced by a synthetic mixture of molecules generated by the stochastic reconstruction algorithm, the representativeness of this synthetic mixture can still be improved. This is done by adjusting the mole fractions of the molecules in order to achieve a close correspondence between the mixture properties of the set of molecules and those of the petroleum fraction. It should be noted that, given that the number of molecules is much larger than the number of analytical data, the problem is ill-posed and has therefore an infinite number of solutions, i.e. several combinations of the mole fractions allow adjusting the mixture properties to the petroleum analyses. Hence, adjusting mixture mole fractions is generally performed using a cost function, which can be based on Shannon entropy [28,35,38], Gibbs free energy [52], or other criteria [53–65]. In this way, the number of possible solutions is reduced to a single one. All molecular reconstruction techniques that modify mole fractions of

L.P. de Oliveira et al. / Catalysis Today 220–222 (2014) 208–220

a set of molecules therefore simply use a criterion to select a single molecular composition out of all possible molecular compositions that comply with the available analytical data. Reconstruction by Entropy Maximization (REM) [28,34–38] is a molecular reconstruction technique that adjusts the mole fractions of a pre-defined set of molecules by maximizing the Shannon entropy criterion E(xi ), originally formulated in the information theory developed by Shannon [66] and given by Eq. (6), subject to a mass balance (Eq. (7)) and to several analytical constraints (Eq. (8)). N 

E(xi ) = −

xi ln (xi )

(6)

i=1

with N 

xi = 1

(7)

i=1

and Fj (x) = Pj

j = 1, J

(8)

where N is the number of molecules present in the predefined set of molecules, xi represents the mole fraction of molecule i, x is the vector of mole fractions xi , Pj is the value of analytical property j, Fj () is the mixing rule for analytical property j, and J is the number of analytical properties or constraints. In our case, REM is applied to adjust the mole fractions of the molecules of the synthetic mixture so as to achieve a closer correspondence between the mixture properties of the set of molecules and those of the petroleum fraction. The maximization of this criterion ensures that in the absence of analytical information, no molecule is favored over the others, which leads to a uniform distribution of the mole fractions of the molecules. However, when constraints (analytical data) are added, the entropy maximization distorts the molar distribution of molecules so as to achieve an exact correspondence between the properties of the mixture and the user-provided analytical data. As opposed to stochastic reconstruction, REM uses classical optimization techniques, in our case a conjugate gradient method, and its computational effort is much smaller. For a more detailed description of the REM method used in the present work, the reader is referred to Hudebine et al. [28,35,38] and de Oliveira et al. [29]. 2.2. Reaction modeling In the second main step of the approach, the conversion process is modeled using a kinetic Monte Carlo (kMC) method, which is an alternative approach to the traditional deterministic kinetic modeling. Indeed, reaction systems are traditionally described in terms of continuous average compositions (molar concentrations) of the chemical species. For this purpose, a set of N coupled algebraic or differential equations (ODE) are developed and solved to obtain a temporal evolution of the N species in the system. Such an approach requires a pre-defined exhaustive reaction network, i.e. a list of all reactions that describe the conversion process throughout the simulation. As shown by several authors [22,67–72], the number of reactions increases significantly with the size of the molecules. For light or intermediate petroleum fractions, the molecules are small enough to generate a representative reaction network of the refining processes treating these mixtures. However, heavy petroleum fractions contain very large molecules for which the generation of the entire reaction network becomes very difficult or even impossible.

211

The kMC method models the chemical reaction system at the molecular level by tracking the transformation of a discrete population of molecules. The evolution of the system does not occur continuously in the time, but by discrete events. Each event describes the transition from one state to another resulting from a chemical reaction. For each transition (or reaction), a probability function is associated, which depends on the molecules involved in the reaction and the rate constant. As opposed to the deterministic approach, the kMC method does not need the exhaustive reaction network, but a list of reaction rules that describe the possible transformations for each molecule. In this way, the network can be generated “on-the-fly”. Thus, the kMC method becomes better suited for the simulation of processes treating heavy fractions than the deterministic method. In the kMC method, the temporal evolution of the species population is analytically described by a single differential–difference equation for a probability function, termed the chemical master equation (CME) [43,73]. This CME determines the states through which the system moves over time. The CME can rarely be solved, either analytically or numerically, due to its high complexity. To overcome this drawback, two distinct Monte Carlo procedures have been proposed in the literature. The first, developed by Gillespie [43], corresponds to an event-space procedure, while the second technique was formulated by McDermott [74] and is a fixed timestep procedure. In the event-space procedure, the evolution of the system is followed event by event, i.e. reaction by reaction, resulting in disparate time intervals between two subsequent reactions. Indeed, the time intervals depend on the number of possible reactions and on their associated rate constants. On the other hand, McDermott’s procedure follows the states of the system at fixed time intervals and calculates which reactions occur in each time interval. In the present work, the complete molecular pathways were initially not available for each molecule, and it is hence difficult to implement the fixed time-step procedure. For this reason, the simulations of the reactions are carried out using Gillespie’s algorithm based on the event-space procedure, termed Stochastic Simulation Algorithm (SSA). The various steps of the implemented algorithm are shown in Fig. 3. The first step of the algorithm is the “molecular discretization” of the feedstock that was generated previously by molecular reconstruction. After SR-REM, the process feedstock is represented by a set of molecules and their mole fractions. However, the kMC method does not follow the mole fractions of the molecules, but the transformation of the discrete molecules. Therefore, the molecular representation of the feedstock must be discretized. The “molecular discretization” is based on the same principle as the quantitative approach used by Hudebine and Verstraete to create a representative database of molecules for gasoline fractions [38]. This approach consists of replicating the molecules according to their mole fractions. The number of replications is calculated as the product of the mole fraction of a molecule with a replication factor (Frep ), and rounding off this product to its nearest integer value. The replication factor Frep defines the maximum number of molecules in the discrete mixture, as well as the lowest mole fraction that is retained. Indeed, all molecules with a mole fraction lower than 0.5/Frep will be discarded from the mixture. For example, in case the replication factor Frep is set to 1000, all molecules whose mole fraction is lower than 0.0005 (i.e., 0.5 × 1/1000 or half of the reciprocal value of Frep ) are removed from the mixture of molecules. The factor 0.5 orginates from the rounding off procedure to the nearest integer. Once the “molecular discretization” step has been performed, all potential reaction events need to be identified. This is performed using a set of reaction rules that generate all reactions by inspecting the structure of the reactant molecules. Once all potential reactions events have been listed, the normalized probability for each reaction event is determined based on the

212

L.P. de Oliveira et al. / Catalysis Today 220–222 (2014) 208–220

Discrete set of reactant molecules Reaction conditions (tinitial, tfinal, T, P, ...)

Reaction rules and rate parameters

Identification of all possible reaction events Calculation of reaction probabilities Calculation of the reaction probability distribution DR Generation of two Random Numbers RN1

Update reactions and reaction probabilities RN2

Determination of the reaction time step ∆ t

Selection of the reaction event by sampling DR No

Execute reaction event

tfinal ?

Yes

Product molecules

Fig. 3. Flow diagram of the SSA algorithm used in the present work.

probability that the reaction ␯ will occur (a␯ ) divided by the sum of all probabilities for all M possible reactions, as shown in Eq. (9). P =

a

M i

(9)

ai

the next, yet undefined, reaction will occur. This reaction time step is calculated using Eq. (12), as proposed by Gillespie [43]: t = −

ln(RN1 )

M i

The probability av that reaction ␯ will occur is defined as the product of the number of distinguishable combinations of the reactant molecules (h␯ ) and the stochastic rate parameter (c␯ ), as given by the following equation: a = h c

(10)

For monomolecular reactions, the number of distinguishable combinations of the reactant molecules (h␯ ) equals 1, while for bimolecular reactions it is equal to the number of available molecules of the second reactant. The stochastic rate parameters (c␯ ) can be determined from the deterministic rate parameters (k␯ ), which are based on the average molecular concentration. As shown by Gillespie [43], for monomolecular reactions, k␯ and c␯ are equal; for bimolecular reactions, c␯ equals k␯ divided by the reaction volume V; for tri-molecular reactions, c␯ equals k␯ divided by V2 . In the simulations, the effect of reaction temperature on the stochastic rate parameters is accounted for by means of an Arrhenius expression:

 E 1 a

c = c (Tref ) exp −

R

T

−

1 Tref



ai

(12)

The normalized cumulative probability distribution DR at time t contains all M reactions that can occur in the mixture at that reaction time and is the sum of the normalized reaction probabilities Pv for each reaction event v. This normalized cumulative probability distribution DR is then randomly sampled in order to select the next reaction (␮) that will be executed in the next reaction time step. This reaction selection step is performed by drawing a second random number (RN2 ) between zero and one. The value of this random number selects the next reaction from the cumulative probability distribution DR , as shown in Eq. (13). DR ( − 1) < RN2 ≤ DR ()

(13)

Once the reaction time step t has been determined and the reaction has been selected, the reaction system is updated by executing the selected reaction, i.e. replacing the reactant(s) by the product(s), and incrementing the simulation time by t. The simulation proceeds event by event (reaction by reaction) until the final simulation time is reached, leading to the set of product molecules. Because of the stochastic nature of the method, several runs are performed and averaged to obtain an accurate representation of the evolution of the reaction system.

(11)

where c␯ is the stochastic rate constant of the reaction type ␯ at the reaction temperature, c␯ (Tref ) is the stochastic rate constant of the reaction type ␯ at the reference temperature Tref , Ea is the activation energy, R is the ideal gas constant, T is the temperature of the system and Tref is the reference temperature. A first random number (RN1 ) is drawn to determine the reaction time step (t) between the current time t and the time at which

3. Simulation of the hydroconversion of Ural vacuum residue Residue hydroconversion, also called residue hydrocracking, is a refining process that converts heavy petroleum fractions into valuable products, such as gasoline, kerosene and gas oil in the presence of a catalyst and of hydrogen at high temperature (410–450 ◦ C) and high pressure (110–200 bar). The catalysts utilized usually

L.P. de Oliveira et al. / Catalysis Today 220–222 (2014) 208–220

contain oxides of Mo, Co, Ni, and/or W, which are supported on alumina, silica, or silica/alumina, and which are used in their sulfided form. The main reactions involved in the residue hydroconversion process are cracking, hydrodesulfurization, hydrodenitrogenation, hydrodemetallization and hydrogenation of aromatic rings. During previous work, the proposed methodology was applied to the Athabasca vacuum residue (VR) hydroconversion [75]. In this paper, the methodology is applied to the hydroconversion of the vacuum residue from a different crude oil, an Ural crude oil. As already discussed in previous work [29], the Ural VR has very different characteristics than the Athabasca VR. The Athabasca VR is a heavy residue that contains high amounts of heteroatoms, such sulfur, nitrogen, oxygen and metals, and very low hydrogen content. In contrast, the Ural VR exhibits a somewhat lower molecular weight, high hydrogen content and low amounts of heteroatoms. The Ural VR hydroconversion experiments were carried out in a batch reactor using a 300 ml autoclave, according to the procedures described elsewhere [76]. The experiments with Ural VR were performed at two temperatures (395 ◦ C and 410 ◦ C) and one reaction time (2 h). The reactor was charged with 50 g of feedstock and 0.5 g of crushed equilibrium industrial sulfided NiMo catalyst, pressurized with hydrogen in excess and quickly heated to reaction temperature (less than 15 min from 100 ◦ C to 420 ◦ C). At the end of the reaction (after a few hours), the reactor was quickly cooled down with an air vortex system (less than 15 min from 420 ◦ C to 100 ◦ C). 3.1. Composition modeling of the Ural vacuum residue A VR is a petroleum fraction that is drawn from the bottom of a vacuum distillation column. This petroleum cut is an extremely complex mixture of hydrocarbons that contains several thousands to millions of different species. As mentioned above, a complete molecular characterization cannot be achieved with the current analytical techniques. VR streams are generally characterized by their bulk properties obtained from the elemental analysis (C, H, S, N, O), specific gravity, average molecular weight, 13 C NMR, and SARA analysis. The SARA analysis classifies molecules according to their polarity into Saturates (non-polar), Aromatics (slightly polar), Resins (polar) or Asphaltenes (highly polar). The analytical data show that VR molecules are mainly composed of carbon and hydrogen, but also contain heteroatoms such as sulfur, nitrogen and oxygen. Even metals such as nickel and vanadium are present, generally in the form of the porphyrin structures.

213

VR fractions contain both acyclic (paraffins) and cyclic structures, the latter being combined into cores. These cores (or poly-cyclic structures) are composed of any combination of naphthenic rings, benzene rings and/or heterocycles, such as thiophenes, pyridines, pyrroles or furans. According to Mass Spectrometry data and 13 C NMR data, alkyl chains are also present in these molecules. The Atmospheric Equivalent Boiling Point (AEBP) of the molecules is typically higher than 350 ◦ C at atmospheric pressure, which corresponds to molecules with more than 25 carbon atoms. As the VR contains non-volatile molecules, the final boiling point cannot be determined and consequently the maximum number of carbon atoms remains unknown. Despite the complexity of VR fractions, several authors [23–27] have succeeded in obtaining molecular information by combining different analytical techniques, such elementary analysis, distillation, gas and liquid chromatography, mass spectrometry. In the present work, it is assumed that a vacuum residue can be described by 16 structural attributes shown in Table 1. Metals were not accounted for during the molecular reconstruction due to their low abundance. VR molecules may have 0 cores (paraffins), one core (naphthenes and aromatic mono-core molecules) or multiple interconnected cores (archipelago structures). Each molecule is therefore initially defined according to its type, which may vary between Paraffins (type 0), Naphthenes (type 1), Aromatic Mono-Core (type 2) or Aromatic Multi-Core (type 3). Given the low concentration of heteroatoms in the paraffins and naphthenes, these elements were not taken into account for this type of molecules. The aliphatic structures (paraffins and alkyl chains) are considered to be linear. Paraffins are therefore described by a chain length, while naphthenes are characterized by a number of cycles, a degree of substitution by alkyl chains and a length for each of their alkyl chains. For aromatic molecules, each core is characterized by a number of naphthenic rings, a number of benzene rings, a number of heterocycles (thiophene, pyridine, pyrrole or furan), and a degree of substitution with alkyl chains. Their alkyl chains are described by a chain length and a degree of substitution by heteroatoms (S, N or O). Three kinds of PDFs were used as optimization variables: histograms, exponential functions and gamma functions. The histograms were used for structural attributes containing a narrow range of possible values (less than four), while the exponential and gamma functions describe attributes with a wide range of possible values. The criterion to choose between the two latter PDFs was the shape of the distribution. For a more detailed description of

Table 1 Structural attributes for vacuum residue molecules. Structural attribute 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 * ** *** ****

*

Type of molecule Number of cores Type of heterocycle** Number of benzene rings per core Total number of rings per core Number of thiophenes per core Number of pyridines per core Number of pyrroles per core Number of furans per core Acceptance probability for a peripheral carbon Length of the paraffinic chains Length of an alkyl chain (lateral and inter-core) Probability of sulfur substitution for aliphatic CH3 or CH2 Substitution probability of a carbon atom by a heteroatom Type of heteroatom substitution*** Type of oxygen group****

Type of molecule: 0, paraffin; 1, naphthene; 2, aromatic mono-core; 3, aromatic multi-core. Type of heterocycle: 0, thiophene; 1, pyridine; 2, pyrrole; 3, furan. Type of heteroatom: 0, nitrogen; 1, oxygen. Type of oxygen group: 0, ether function; 1, carbonyl function.

Values

Distribution

Parameters

0, 1, 2 or 3 >1 0, 1, 2 or 3 >0 >1 0, 1 or 2 0, 1 or 2 0, 1 or 2 0, 1 or 2 0 or 1 >1 >1 0 or 1 0 or 1 0 or 1 0 or 1

Histogram Exponential Histogram Exponential Gamma Histogram Histogram Histogram Histogram Histogram Gamma Gamma Histogram Histogram Histogram Histogram

0, 1 and 2 3 4, 5 and 6 7 8 9 and 10 11 and 12 13 and 14 15 and 16 17 18 19 20 21 22 23

214

L.P. de Oliveira et al. / Catalysis Today 220–222 (2014) 208–220

Table 2 Comparison between experimental and calculated properties after the stochastic reconstruction (SR) step, the entropy maximization (REM) step and molecular discretization (MD). Experimental data

Simulation (SR output)

Simulation (SR/REM output)

Simulation (MD output)

Elemental analysis Carbon Hydrogen Sulfur Nitrogen Oxygen

wt% wt% wt% wt% wt%

85.5 10.6 2.7 0.6 0.6

85.4 10.6 2.7 0.6 0.6

85.5 10.6 2.7 0.6 0.6

85.5 10.6 2.8 0.6 0.6

SARA Saturates Aromatics Resins Asphaltenes

wt% wt% wt% wt%

11.7 46.1 37.6 4.6

11.4 46.5 38.3 3.8

11.7 46.1 37.6 4.6

11.5 46.2 37.7 4.6

at% at%

72.8 27.2

72.9 27.1

72.8 27.2

72.7 27.3

13

C NMR Saturated carbon Aromatic carbon

Simulated distillation 0% wt 5% wt 10% wt 20% wt 30% wt 40% wt Other Average molecular weight

◦

C C ◦ C ◦ C ◦ C ◦ C

384 496 520 550 574 598

379 449 498 581 650 715

379 496 520 550 574 601

379 493 520 549 574 598

g mol−1

706.8

724.8

713.2

709.9

◦

the structural attributes and building diagram used in the present work, the reader is referred to Verstraete et al. [42] and de Oliveira et al. [29,75]. All pure component properties (chemical formula, aromatic carbon content, molecular weight, etc.) are calculated by inspection of the structure of the molecule, except for the normal boiling point, or AEBP, and the SARA analysis. The AEBP is calculated by means of a group contribution method initially developed by Hudebine et al. [38] and extended by de Oliveira et al. [29]. The classification of molecules into SARA fractions is based on their molecular weight, their hydrogen content and their nitrogen content by means of a SARA diagram. The SARA diagram was derived from the solventresid phase diagram proposed by Wiehe [77] and is described elsewhere in more detail [29]. Once all pure component properties are available, the average properties of the mixture are calculated using linear mixing rules. These mixture properties are compared to the analytical data through an objective function (Eq. (5)) that contains the sum of the absolute values of the relative errors (Eq. (2)). By combining those two expressions, the objective function used in this application is therefore a least absolute deviations objective function. This objective function was minimized in order to determine the optimal values of the PDF parameters using an elitist genetic algorithm [78,79] evolving over 100 generations. The Ural VR was represented by a set of 5000 molecules in order to ensure an optimal balance between the required CPU time and the accuracy of the feedstock representation, as illustrated elsewhere [75]. Elemental analyses (carbon, hydrogen, sulfur, nitrogen and oxygen), 13 C NMR, SARA analysis, simulated distillation and average molecular weight were used as input to the SR-REM algorithm. All experimental data were obtained at IFP Energies nouvelles. The “experimental” value for the average molecular weight was obtained from an API correlation that is based on specific gravity and simulated distillation [80]. The comparison between analytical and calculated properties of the Ural VR is shown in Table 2. As can be seen in Table 2, the properties of the mixture obtained after the SR step (SR output) show a good agreement with most of the analytical data. The elemental composition of the set of

molecules is very close to the experimental data. The SARA fraction and 13 C NMR are well represented. Some differences are, however, observed in the average molecular weight and partial simulated distillation curve. As discussed in previous works [29,75], these deviations are probably caused by various hypotheses used during the calculation of the simulated distillation curve. First of all, the AEBP of each molecule is obtained by means of a group contribution method and adjusted for an additional retention from interactions with the not completely apolar GC column. Then, the simulated distillation curve is calculated assuming that the molecules are perfectly separated by increasing adjusted boiling points, which is not strictly true, especially in the region near the initial and final boiling point [28,81]. In addition, the group contribution method was developed on the basis of a large database of boiling points, but containing mainly molecules with less than 42 carbon atoms [29,38]. For very large molecules, the boiling point is therefore estimated by extrapolation, which may increase the imprecision of the method. Nevertheless, the properties of the synthetic mixture of molecules are well improved when the REM step is applied, as illustrated in the third column of Table 2. The elemental analysis, SARA fractions separation and 13 C NMR are nearly perfectly adjusted to the experimental data. In the case of the molecular weight and simulated distillation, the observed deviations were significantly reduced except for the initial boiling point. This is due to the fact that the initial boiling point is determined by a single molecule, the lightest one. The results shown above illustrate that the SR-REM algorithm generates a set of molecules whose properties are identical to those of the Ural vacuum residue analyses. This mixture can therefore be used as an input for the kinetic model.

3.2. Simulation of residue hydroconversion reactions To simulate vacuum residue hydroconversion by means of a stochastic simulation method, the molecular representation of the feedstock must be discretized, a set of reaction rules must be specified, together with the rate constants of the various reactions.

L.P. de Oliveira et al. / Catalysis Today 220–222 (2014) 208–220

215

Table 3 Reaction types and their reaction rules for hydroconversion. Reaction type

Reaction rules

Ring dealkylation (saturated and aromatic rings)

+

+

H2

Cleavage of the C C bond near the ring when the side chain length is higher than 2 Ring opening

+

H2

The C C ring bond cleavage takes place near the bridgehead atom. Ring opening occurs only when no alkyl chain is connected to ring

Hydrogenation of an aromatic ring

+ 2H2 Hydrogenation of a benzene ring into a cyclohexane ring (ring by ring) Dehydrogenation of a saturated ring

+ 2H2 Dehydrogenation of a cyclohexane ring into a benzene ring (ring by ring) Hydrodesulfurization of non-heterocyclic compounds

SH

+ H2

+ H2S Cleavage of the C S bond on a chain, removing the sulfur atom as H2 S

Hydrodesulfurization of heterocyclic compounds

+ 3H2

+ H2S

S

Removal of the sulfur atom from a thiophene ring as H2 S

The molecular discretization was carried out with a replication factor equal to 10,000 and resulted in a discrete population of 9984 molecules. The mixture properties of this discrete population of molecules are also shown in Table 2. It should be noted that the molecular discretization does not introduce any significant changes to the properties of the mixture. The reactions rules for vacuum residue hydroconversion used in the present work are summarized in Table 3. The present simulation includes reaction rules for the dealkylation of rings, ring opening, hydrogenation of aromatic rings, dehydrogenation of saturated rings, and hydrodesulfurization of the non-heterocyclic and heterocyclic compounds. A justification for the choice of these reaction rules is given elsewhere [75]. The last step concerns the definition of the rate parameters of the various reaction events. Concerning the hydrogenation reactions, the reaction probabilities av is calculated from Eq. (14): aHydro = cHydro pnHH

(14)

2

where cHydro is the stochastic rate parameter for hydrogenation, pH2 is the hydrogen partial pressure, and nH is the reaction stoichiometry with respect to hydrogen molecule. The stochastic rate parameter for the hydrogenation reactions is calculated using a Quantitative Structure/Reactivity Correlation (QS/RC) that was developed by de Oliveira et al. [82] and given by: ln(cHydro ) =

3.26 − 7.06n1.18 H

   HR0    + 0.616NAR − 1.08  R  1  E 1

+ 0.330NSR − 2.80NTR −

a

R

·

T

−

Tref

(15)

where nH is the reaction stoichiometry with respect to hydrogen molecule, HR0 is the heat of hydrogenation reaction at 25 ◦ C, NAR and is the number of aromatic (benzene) rings of the molecule, NSR

is the number of saturated (cyclohexane) rings of the molecule, NTR is the number of thiophene rings fused to the aromatic ring, Ea is the activation energy for hydrogenation, R is the ideal gas constant, T is the reaction temperature, and Tref is the reference temperature. For the dehydrogenation reactions, the reaction probabilities av and stochastic rate parameters cv are calculated to be thermodynamically consistent with the hydrogenation reactions. This therefore leads to the following expression for the reaction probability: aDehydro = cDehydro =

cHydro

(16)

Keq

where cDehydro is the stochastic rate parameter for dehydrogenation, and Keq is the equilibrium constant. For the equilibrium constant Keq , the Quantitative Structure/Reactivity Correlation (QS/RC) given by has been developed by de Oliveira et al. [82]:

   HR0   R 

ln(Keq ) = 2.952 − 13.215nH + 5.196 × 10−3  − 0.784NSR −

HR R

1 T

−

1 TRef

 (17)

where nH is the reaction stoichiometry in respect to hydrogen molecule, HR0 is the heat of hydrogenation reaction at 25 ◦ C, NSR is the number of saturated rings (cyclohexane) of the molecule, R is the ideal gas constant, T is the reaction temperature, and Tref is the reference temperature. For ring dealkylation, ring opening and HDS reactions are considered to be first order or pseudo-first order reactions with respect to the reacting molecule. Hence, for each molecule, the probability av that the reaction ␯ will occur is therefore obtained using Eq. (10), in which the number of distinguishable combinations of the reactant molecules (h␯ ) is equal to 1.

216

L.P. de Oliveira et al. / Catalysis Today 220–222 (2014) 208–220

Sulfur removal 80 410°C

60

HDS (%)

VR Conversion (%)

VR conversion 80

40 395°C 20

60 410°C 40 395°C

20

0

0 0

1

2

3

0

4

Time (h) Simulated

1

2

3

4

Time (h)

Experimental

Simulated

Experimental

Fig. 4. Comparison between experimental data and simulation results for residue conversion and sulfur removal of Ural VR hydroconversion at 410 ◦ C and 395 ◦ C.

Product distribution

Product Yields (%)

100 80

Gas (C1 – C4)

395°C

Naphtha (C5 – 220°C)

60

410°C

GO (220°C – 380°C)

40

VGO (380°C – 520°C) VR (520° C+)

20 0 0

20

40

60

80

VR Conversion (%) Fig. 5. Comparison between experimental data and simulation results for the products yields as function of VR conversion of Ural VR hydroconversion at 410 ◦ C and 395 ◦ C.

The stochastic rate parameters at the reference temperature c␯ (Tref ) are fitted to the experimental data on Ural vacuum residue hydroconversion at 410 ◦ C (Figs. 4 and 5) and their values are presented in Table 4, together with the reference temperature. To reduce the number of degrees of freedom during fitting, the number of adjustable parameters was reduced as much as possible. For the HDS of heterocyclic compounds, only an HDS activity parameter was fitted, and the rate parameters have been adjusted for reactivity/activity differences in such a way that the ratios between the reaction rates are equal to the ratios obtained in previous work [39,75,82]. It is well-known that the HDS reaction family can follow a hydrogenolysis pathway or a hydrogenation pathway [83,84]. The rate parameters given in Table 4 correspond to the hydrogenolysis pathway. For hydrogenation pathway, saturation of the benzene rings fused to the thiophene ring is assumed to follow a classical hydrogenation reaction and its rate constant is estimated by Eq. (15). The thiophene saturation followed by removal of the sulfur atom is assumed to occur in a single step and its rate parameter is equal to that of the hydrogenolysis of the thiophene. Concerning

the hydrogenation and dehydrogenation reactions, only an overall hydrogenation activity parameter was adjusted (the constant in the QS/RC given in Eq. (15). Finally, the activation energies were not estimated, but directly taken from the literature, with the activation energy for ring dealkylation estimated by Danial-Fortain [85] on the same data set, while the activation energies for the HDS reactions and the hydrogenation reactions were proposed by López-García [86] and obtained for hydrotreating of gas oils [21,39,86]. The activation energies for the various reaction types are also listed in the Table 4. The Ural VR hydroconversion was simulated at 410 ◦ C and 395 ◦ C. For each temperature, the simulations are repeated 5 times and the mixture properties of the population of molecules are averaged. The comparison between the predicted and experimental data is shown in Figs. 4 and 5. The left hand side graph of Fig. 4 corresponds to the 520 ◦ C+ conversion of the VR, i.e. of the molecules boiling above 520 ◦ C, while the temporal evolution of the sulfur removal is represented in the right hand side graph. In both graphs, the lines represent the simulation results and the symbols correspond to the experimental data. For both temperatures, the prediction of the conversion of VR and the sulfur removal is quite accurate. These results are even more remarkable considering the fact that the rate parameters were adjusted at one single temperature (410 ◦ C) and that the performances at the second temperature (395 ◦ C) were predicted using values for the activation energies obtained from different literature sources. Fig. 5 represents the yields of the various petroleum products as function of the 520 ◦ C+ conversion of the VR. To this aim, the product molecules are grouped into 5 families according their boiling point: gas, naphtha, gas oil (GO), vacuum gas oil (VGO) and unconverted vacuum residue (VR). The gas is the lightest fraction containing the molecules with number of carbons between 1 and 4 which has a boiling point lower than 0 ◦ C. The naphtha, gasoil and vacuum gasoil corresponds to petroleum fractions that

Table 4 Reaction families and their rate parameters, activation energies and reference temperature. Reaction families

c␯ (h−1 )

Ea (J mol−1 )

Tref (◦ C)

Hydrogenation/dehydrogenation Ring dealkylation Ring opening Hydrodesulfurization of non-heterocyclic compounds

Eq. (15) + Eq. (17) 8.0 × 10−1 7.2 × 10−1 4.0 × 10−1

50,942 217,000 50,000 50,000

320 410 410 410

Hydrodesulfurization of heterocyclic compounds Hydrogenolysis of thiophene and their derivatives Hydrogenolysis of BT and their derivatives Hydrogenolysis of DBT and their derivatives Hydrogenolysis of 4-alkyl-DBT and their derivatives Hydrogenolysis of 4,6-dialkyl-DBT and their derivatives

6.4 × 10−2 4.0 × 10−2 6.4 × 10−3 4.0 × 10−5 4.0 × 10−6

50,000 70,000 84,991 113,340 120,416

320 320 320 320 320

L.P. de Oliveira et al. / Catalysis Today 220–222 (2014) 208–220

Fig. 6. Comparison between the simulated and experimental composition in terms of paraffins, naphthenes and aromatics for the naphtha fraction at 49 wt% conversion of the Ural VR.

Fig. 7. Comparison between the simulated and experimental composition in terms of paraffins, naphthenes, mono-, di- and poly-aromatics of the atmospheric gas oil fraction at a conversion level of 49 wt% of Ural VR.

contains the molecules with a boiling point lower than 220 ◦ C, 380 ◦ C and 520 ◦ C. For low VR conversion, the product yields are well predicted. At intermediate conversions for the VR, the gas yield and gas oil yield are well predicted, but the deviations are observed in the prediction of the naphtha and vacuum gas oil fractions. To understand this imbalance between naphtha and vacuum gas oil, the composition of the naphtha and atmospheric gas oil fractions was determined in terms of paraffins, naphthenes and aromatics. For the simulations, this composition is obtained by grouping the various product molecules according to their PONA type and their AEBP. The experimental data were obtained by gas chromatography analysis for the naphtha fraction and by mass spectrometry for the GO fraction. Fig. 6 shows the comparison between the simulated and the experimental composition of the naphtha fraction at a conversion level of 49 wt% of Ural VR, while Fig. 7 compares the simulated and experimental composition of the GO fraction for the same

condition. According to the experimental data, the naphtha fraction is composed by 50 wt% of non-cyclic compounds and 50 wt% of non-cyclic compounds, while atmospheric gas oil is mainly composed by cyclic compounds (>70 wt%). However, the population of molecules obtained by stochastic simulation contains a very high amount of paraffins for both petroleum fractions. These originate from the fact that the transformation of heavy fractions into lighter fractions is mainly done by dealkylation of the rings. Even though the formation of the lighter fractions can also be achieved by ring opening of VR and VGO molecules, the implemented reaction rules only open a ring when no alkyl chains are connected to the ring. This limitation imposed on the ring opening reactions means that the ring can be opened only after all chains connected to the ring have been removed. This then explains why the paraffins content of the naphtha and atmospheric gas oil cuts is much higher than the experimentally observed data. As illustrated above, this molecule-based kinetic modeling approach not only allows predicting the global information of the process, such as its yield structure and hydrotreating performances, but it also provides many molecular details throughout the reactor simulation. Indeed, one can easily “regroup” the molecules of the synthetic mixture so as to assess the evolution of the concentrations of various chemical families. Similarly, the composition and various properties (density, molecular weight, sulfur content, aromatics content, CCR content, etc.) of selected petroleum cuts can easily be determined, as illustrated in Figs. 6 and 7. All this information is available throughout the reactor simulation, and not only for the effluent, which allows following their evolution and optimizing the conversion process. The proposed methodology also provides the molecular structural information. For example, the temporal evolutions of the some molecular attributes are illustrated in Fig. 8. The left hand side graph corresponds to the sulfur-containing structures that are grouped into heterocyclic and non-heterocyclic structures. The non-heterocyclic structures are nearly completely removed from the VR, which is expected due to the high reactivity of these structures. The amount of heterocyclic structures is gradually reduced over the reaction time as a result of the hydrogenolysis of the thiophene, benzothiophene, and dibenzothiophene type compounds. The second graph in Fig. 8 illustrates the temporal evolution of the core structures in terms of the average number of cores per molecule and the average number of rings per core. The substantial decrease in the average number of cores implies the separation of the cores by cracking of the inter-core chains. The number of rings per core also decreases along the reaction time. Ring opening reaction, as well the HDS of the heterocyclic structures are responsible by this reduction.

Sulfur-containing structures

Structural attributes of cores

3

6

Average Number

Fraction (%wt)

217

2

1

5 4 3 2 1 0

0 0

1

2

3

4

0

1

Total

Heterocyclic

Non-heterocyclic

2

3

4

Time (h)

Time (h) Cores

Rings per core

Fig. 8. Average temporal values of the sulfur-containing structures and of the structural attributes for the cores for the Ural VR.

218

L.P. de Oliveira et al. / Catalysis Today 220–222 (2014) 208–220

4. Conclusions In this paper, a novel two-step kinetic modeling strategy for heavy oil hydroconversion processes has been described and illustrated. In the first step, a SR-REM molecular reconstruction algorithm is used to generate a synthetic mixture of molecules that represent the process feedstock. This algorithm couples two independent molecular reconstruction methods: Stochastic Reconstruction (SR) and Reconstruction by Entropy Maximization (REM). The SR method is employed first to generate an initial equimolar set of molecules that are representative of the process feedstock and whose mixture properties correspond well to the analytical data. The mole fractions of the generated molecules are then adjusted by the REM method to obtain a synthetic mixture whose properties closely match those of the process feedstock. The second step of the modeling strategy consists in simulating the effect of the conversion reactions on the generated set of molecules using a kinetic Monte Carlo (kMC) method, termed Gillespie’s Stochastic Simulation Algorithm (SSA). The SSA is a Markovian process that transforms the molecules, event by event, based on probability considerations. Given the discrete nature of this approach, the SSA is preceded by a additional step, termed “molecular discretization”, to transform the set of molecules generated by the SR-REM algorithm into a discrete population of molecules. The methodology has been successfully applied to the Ural vacuum residue (VR) hydroconversion. The SR-REM algorithm was used to generate a molecular representation of this VR starting from the available analytical information (CHNOS analysis, 13 C NMR analysis, SARA analysis, and simulated distillation). These results illustrate that the SR-REM algorithm is able to generate a very good molecular representation of the processed feedstock, but they also confirm the importance of using both the SR and the REM method. Indeed, the SR method is not able to generate a set of molecules whose properties are closely adjusted to some of the analytical data, such as the distillated simulation. By adding the REM method as second step to the SR-REM algorithm, the deviations between the experimental and simulated properties are almost completely eliminated. The effect of the conversion reactions was simulated over the synthetic mixture of molecules using our implementation of Gillespie’s SSA. In the present work, ring dealkylation, ring opening, hydrogenation of aromatics (HDA) and hydrodesulfurization (HDS) of non-heterocyclic and heterocyclic sulfur compounds were considered as the main reactions occurring during hydroconversion. The rate constants for HDA were calculated by using Quantitative Structure/Reactivity Correlations (QS/RCs), while the other reactions were assigned a rate parameter. The parameters of QS/RCs and the rate parameters were obtained in previous work by fitting to experimental data of the Athabasca VR hydroconversion, so that only two activities (HDS and dealkylation) had to be adjusted to simulate the hydroconversion of Ural VR. For the sulfur removal and overall VR conversion, the agreement between the simulation results and the experimental values was quite accurate. Concerning the yield structure, slight deviations are observed for some products at the highest conversion. These deviations may result from the absence of some reaction pathways in the reaction network. The effect of operating temperature was very well represented, especially given the fact that the rate parameters were adjusted at one single temperature (410 ◦ C) and that the performances at the second temperature (395 ◦ C) were predicted using values for the activation energies obtained from different literature sources. Given that the modeling methodology is molecule-based, a vast array of detailed information is available throughout the reactor

simulation. Indeed, the temporal evolution of the various structural attributes can be recalculated. Also, one can easily “regroup” the molecules of the synthetic mixture so as to assess the evolution of the concentrations of various chemical families. Similarly, the composition and various properties (density, molecular weight, elemental analysis, PIONA analysis, SARA analysis, NMR analyses, CCR content, etc.) of selected petroleum cuts can easily be determined. As this information is available throughout the reactor simulation, and not only for the effluent, one can then follow the evolution of various properties and hence optimize the process. In conclusion, the proposed molecule-based modeling approach has several advantages. First of all, the complex feedstock is represented by means of a synthetic mixture of molecules, which allows to retain the molecular detail throughout the entire reactor simulation. Secondly, only a limited number of reaction types need to be defined. This allows tracking the transformation of each molecule, event by event, thereby generating the reaction network “on the fly” as the reactions proceed. Finally, by following the modifications to the synthetic mixture throughout the reactor simulation, the evolution of all product yields and product properties can be assessed by evaluating the yields, composition and characteristics of sub-groups of molecules. Hence, by means of fundamental composition, kinetic and reactor modeling techniques, one can acquire a more fundamental understanding of very complex processes. References [1] M.S. Rana, V. Sámano, J. Ancheyta, J.A.I. Diaz, A review of recent advances on process technologies for upgrading of heavy oils and residua, Fuel 86 (2007) 1216–1231. [2] A.-Y. Huc (Ed.), Heavy Crude Oils: From Geology to Upgrading, Editions Technip, Paris, 2011, ISBN 978-2-7108-0890-9. [3] H. Toulhoat, P. Raybaud (Eds.), Catalysis by Transition Metal Sulphides: From Molecular Theory to Industrial Application, Editions Technip, Paris, 2013, ISBN 978-2-7108-0991-3. [4] R.J. Quann, S.B. Jaffe, Structure-oriented lumping: describing the chemistry of complex hydrocarbon mixtures, Industrial and Engineering Chemistry Research 31 (1992) 2483–2497. [5] R.J. Quann, S.B. Jaffe, Building useful models of complex reaction systems in petroleum refining, Chemical Engineering Science 51 (1996) 1615–1635. [6] M.T. Klein, G. Hou, R.J. Bertolacini, L.J. Broadbelt, A. Kumar, Molecular Modeling in Heavy Hydrocarbon Conversions, CRC Press Taylor & Francis Group, Boca Raton, 2006, ISBN 978-0-8247-5851-6, pp. 5851–5856. [7] V.W. Weekman, Lumps, models, and kinetics in practice, AIChE Monograph Series 75 (11) (1979) 1–29. [8] F. Mosby, R.D. Buttke, J.A. Cox, C. Nikolaides, Process characterization of expanded-reactors in series, Chemical Engineering Science 41 (1986) 989–995. [9] R. Krishna, A.K. Saxena, Use of an axial-dispersion model for kinetic description of hydrocracking, Chemical Engineering Science 44 (1989) 703–712. [10] M.R. Gray, Lumped kinetics of structural groups: hydrotreating of heavy distillate, Industrial and Engineering Chemistry Research 29 (4) (1990) 505–512. [11] M.A. Callejas, M.T. Martínez, Hydroprocessing of a Maya residue: intrinsic kinetics of sulfur-, nitrogen-, nickel- and vanadium-removal reactions, Energy and Fuels 13 (3) (1999) 629–636. [12] F.X. Haulle, Modélisation cinétique de l’hydrotraitement en lit fixe des résidus pétrolièrs: Etude de la réactivité des composés soufrés, Université Paris VI, Paris, 2002 (PhD thesis). [13] C. Botchwey, A.K. Dalai, J. Adjaye, Product selectivity during hydrotreating and mild hydrocracking of bitumen-derived gas oil, Energy and Fuels 17 (2003) 1372–1381. [14] L.O. Oyekunle, B.O. Kalejaiye, Kinetic modeling of hydrodesulphurization of residual oils. I. Power law model, Petroleum Science and Technology 21 (9–10) (2003) 1475–1488. [15] J.-M. Schweitzer, S. Kressmann, Ebullated bed reactor modeling for residue conversion, Chemical Engineering Science 59 (22–23) (2004) 5637–5645. [16] J. Ancheyta, M.A. Rodríguez, S. Sánchez, Kinetic modeling of hydrocracking of heavy oil fractions: a review, Catalysis Today 109 (2005) 76–92. [17] S. Sánchez, M.A. Rodríguez, J. Ancheyta, Kinetic model for moderate hydrocracking of heavy oils, Industrial and Engineering Chemistry Research 44 (2005) 9409–9413. [18] J.J. Verstraete, K. Le Lannic, I. Guibard, Modeling fixed-bed residue hydrotreating processes, Chemical Engineering Science 62 (18–20) (2007) 5402–5408. [19] A. Alvarez, J. Ancheyta, Modeling residue hydroprocessing in a multi-fixed-bed reactor system, Applied Catalysis A 351 (2) (2008) 148–158. [20] C. Ferreira, J. Marques, M. Tayakout-Fayolle, I. Guibard, F. Lemos, F. Ramôa Ribeiro, Modeling residue hydrotreating, Chemical Engineering Science 65 (1) (2010) 322–329.

L.P. de Oliveira et al. / Catalysis Today 220–222 (2014) 208–220 [21] C. López-García, D. Hudebine, J.-M. Schweitzer, J.J. Verstraete, D. Ferré, In-depth modeling of gas oil hydrotreating: from feedstock reconstruction to reactor stability analysis, Catalysis Today 150 (2010) 279–299. [22] H. Cochegrue, P. Gauthier, J.J. Verstraete, K. Surla, D. Guillaume, P. Galtier, J. Barbier, Reduction of single event kinetic models by rigorous relumping: application to catalytic reforming, Oil and Gas Science and Technology – Revue d’IFP Energies Nouvelles 66 (3) (2011) 367–397. [23] M.M. Boduszynski, Composition of heavy petroleums. 1. Molecular weight, hydrogen deficiency, and heteroatom concentration as a function of atmospheric equivalent boiling point up to 1400 ◦ F (760 ◦ C), Energy and Fuels 1 (1987) 2–11. [24] M.M. Boduszynski, Composition of heavy petroleums. 2. Molecular characterization, Energy and Fuels 2 (1988) 597–613. [25] E.Y. Sheu, Petroleum asphaltene properties, characterization, and issues, Energy and Fuels 16 (2002) 74–82. [26] I. Merdrignac, D. Espinat, Physicochemical characterization of petroleum fractions: the state of the art, Oil and Gas Science and Technology – Revue d’IFP Energies Nouvelles 62 (2007) 7–32. [27] A.M. McKenna, G.T. Blakney, F. Xian, P.B. Glaser, R.P. Rodgers, A.G. Marshall, Heavy petroleum composition. 2. Progression of the Boduszynski model to the limit of distillation by ultrahigh-resolution FT-ICR mass spectrometry, Energy and Fuels 24 (2010) 2939–2946. [28] D. Hudebine, J.J. Verstraete, Molecular reconstruction of LCO gas oils from overall petroleum analyses, Chemical Engineering Science 59 (2004) 4755–4763. [29] L. de Oliveira, A. Trujillo Vazquez, J.J. Verstraete, M. Kolb, Molecular reconstruction of petroleum fractions: application to various vacuum residues, Energy and Fuels 27 (7) (2013) 3622–3641. [30] M. Neurock, A. Nigam, D.M. Trauth, M.T. Klein, Molecular representation of complex hydrocarbon feedstocks through efficient characterization and stochastic algorithms, Chemical Engineering Science 49 (1994) 4153–4177. [31] D.M. Trauth, S.M. Stark, T.F. Petti, M. Neurock, M.T. Klein, Representation of the molecular structure of petroleum resid through characterization and Monte Carlo modeling, Energy and Fuels 8 (1994) 576–580. [32] F. Khorasheh, R. Khaledi, M.R. Gray, Computer generation of representative molecules for heavy hydrocarbon mixtures, Fuel 77 (1998) 247–253. [33] D. Hudebine, C. Vera, F. Wahl, J.J. Verstraete, Molecular representation of hydrocarbon mixtures from overall petroleum analyses, AIChE Spring Meeting (2002) 10–14. [34] S.B. Jaffe, H. Freund, W.N. Olmstead, Extension of structure-oriented lumping to vacuum residua, Industrial and Engineering Chemistry Research 44 (26) (2005) 9840–9852. [35] K.M. Van Geem, D. Hudebine, M.-F. Reyniers, F. Wahl, J.J. Verstraete, G.B. Marin, Molecular reconstruction of naphtha steam cracking feedstocks based on commercial indices, Computers and Chemical Engineering 31 (9) (2007) 1020–1034. [36] K.M. Van Geem, M.-F. Reyniers, G.B. Marin, Challenges of Modeling Steam Cracking of Heavy Feedstocks, Oil and Gas Science and Technology – Revue d’IFP Energies Nouvelles 63 (1) (2008) 79–94. [37] S.P. Pyl, K.M. Van Geem, M.-F. Reyniers, G.B. Marin, Molecular reconstruction of complex hydrocarbon mixtures: an application of principal component analysis, AIChE Journal 56 (12) (2010) 3174–3188. [38] D. Hudebine, J.J. Verstraete, Reconstruction of petroleum feedstocks by entropy maximization. application to FCC gasolines, Oil and Gas Science and Technology – Revue d’IFP Energies Nouvelles 66 (2011) 437–460. [39] L. de Oliveira, J.J. Verstraete, M. Kolb, A Monte Carlo modeling methodology for the simulation of hydrotreating processes, Chemical Engineering Journal 207–208 (2012) 94–102. [40] J.J. Verstraete, N. Revellin, H. Dulot, Molecular reconstruction of vacuum gas oils, Preprints of Papers – American Chemical Society, Division of Fuel Chemistry 49 (2004) 20–21. [41] N. Charon-Revellin, H. Dulot, C. López-García, J. Jose, Kinetic modeling of vacuum gas oil hydrotreatment using a molecular reconstruction approach, Oil and Gas Science and Technology – Revue d’IFP Energies Nouvelles 66 (3) (2011) 479–490. [42] J.J. Verstraete, P. Schnongs, H. Dulot, D. Hudebine, Molecular reconstruction of heavy petroleum residue fractions, Chemical Engineering Science 65 (2010) 304–312. [43] D.T. Gillespie, A general method for numerically simulating the Stochastic time evolution of coupled chemical reactions, Journal of Computational Physics 22 (1976) 403–434. [44] W.H. Press, S.A. Teukolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes: The Art of Scientific Computing, third ed., Cambridge University Press, Cambridge, UK, 2007, ISBN 978-05218-8068-8. [45] J. Bard, Nonlinear Parameter Estimation, Academic Press, New York, 1974, ISBN 978-0-1207-8250-5. [46] G.A.F. Seber, C.J. Wild, Non Linear Regression, Wiley-Interscience, New York, 1989, ISBN 978-0-4714-7135-6. [47] P. Englezos, N. Kalogerakis, Applied Parameter Estimation for Chemical Engineers, Marcel Dekker Inc, New York, 2000, ISBN 978-0-8247-9561-0. [48] S. Kirkpatrick, C.D. Gelatt, M.P. Vecchi, Optimization by simulated annealing, Science 220 (1983) 671–680. [49] D.E. Goldberg, Genetic Algorithms in Search, Optimization, and Machine Learning, Addison-Wesley Longman Publishing Co. Inc., Boston, MA, 1989, ISBN 978-0-2011-5767-3.

219

[50] M. Mitchell, An Introduction to Genetic Algorithms, MIT Press, Cambridge, MA, 1998, ISBN 978-0-2626-3185-3. [51] J. Kennedy, R.C. Eberhart, Swarm Intelligence, Morgan Kaufmann Publishers, San Francisco, CA, 2001, ISBN 978-1-5586-0595-4. [52] Z.Y. Ha, Z. Ring, S.J. Liu, Derivation of molecular representations of middle distillates, Energy and Fuels 19 (2005) 2378–2393. [53] D.T. Allen, D.K. Liguras, Structural models of catalytic cracking chemistry: a case study of a group contribution approach to lumped kinetic modeling, in: A.V. Sapre, F.J. Krambeck (Eds.), Chemical Reactions in Complex Mixtures: The Mobil Workshop, Van Nostrand Reinhold, New-York, 1991, pp. 101–125. [54] T.A. Albahri, Molecularly explicit characterization model (MECM) for light petroleum fractions, Industrial and Engineering Chemistry Research 44 (24) (2005) 9286–9298. [55] I.P. Androulakis, M.D. Weisel, C.S. Hsu, K.N. Qian, L.A. Green, J.T. Farrell, K. Nakakita, An integrated approach for creating model diesel fuels, Energy and Fuels 19 (2005) 111–119. [56] M.M.S. Aye, N. Zhang, A novel methodology in transforming bulk properties of refining streams into molecular information, Chemical Engineering Science 60 (2005) 6702–6717. [57] E. Eckert, T. Vanˇek, New approach to the characterisation of petroleum mixtures used in the modelling of separation processes, Computers and Chemical Engineering 30 (2005) 343–356. [58] E. Eckert, T. Vanˇek, Extended utilization of the characterization of petroleum mixtures based on real components, Chemical Papers 59 (6) (2005) 428–433. [59] T.A. Albahri, Enhanced method for predicting the properties of light petroleum fractions, Fuel 85 (2006) 748–754. [60] E. Eckert, T. Vanˇek, Mathematical modelling of selected characterisation procedures for oil fractions, Chemical Papers 62 (1) (2008) 26–33. [61] J. Gomez-Prado, N. Zhang, C. Theodoropoulos, Characterisation of heavy petroleum fractions using modified molecular-type homologous series (MTHS) representation, Energy 33 (2008) 974–987. [62] Y. Wu, N. Zhang, Molecular management of gasoline streams, Chemical Engineering Transactions 18 (2009) 749–754. [63] E. Eckert, T. Vanˇek, Improvements in the selection of real components forming a substitute mixture for petroleum fractions, Chemical Papers 63 (4) (2009) 399–405. [64] M.I. Ahmad, N. Zhang, M. Jobson, Molecular components-based representation of petroleum fractions, Chemical Engineering Research and Design 89 (2011) 410–420. [65] E. Eckert, T. Vanˇek, Z. Bˇelohlav, P. Zámostny, Effective characterization of petroleum C7+ fractions, Fuel 102 (2012) 545–553. [66] C.E. Shannon, A mathematical theory of communication, The Bell System Technical Journal 27 (379–423) (1948) 623–656. [67] M.A. Baltanas, G.F. Froment, Computer generation of reaction networks and calculation of product distributions in the hydroisomerization and hydrocracking of paraffins on Pt-containing bifunctional catalysts, Computers and Chemical Engineering 9 (1985) 71–81. [68] E. Vynckier, G.F. Froment, Modeling of the kinetics of complex processes upon elementary steps, in: G. Astarita, S.I. Sandler (Eds.), Kinetic and thermodynamic Lumping of Multicomponent Mixtures, Elsevier Science Publishers BV, Amsterdam, 1991, pp. 131–161. [69] D. Guillaume, Network generation of oligomerization reactions: principles, Industrial and Engineering Chemistry Research 45 (2006) 4554–4557. [70] J.R. Shahrouzi, D. Guillaume, P. Rouchon, P. Da Costa, Stochastic simulation and single events kinetic modeling: application to olefin oligomerization, Industrial and Engineering Chemistry Research 47 (2008) 4308–4316. [71] E. Valéry, D. Guillaume, K. Surla, P. Galtier, J.J. Verstraete, D. Schweich, Kinetic modeling of acid catalyzed hydrocracking of heavy molecules: application to squalane, Industrial and Engineering Chemistry Research 46 (14) (2007) 4755–4763. [72] D. Guillaume, E. Valéry, J.J. Verstraete, K. Surla, P. Galtier, D. Schweich, Single event kinetic modelling without explicit generation of large networks: application to hydrocracking of long paraffins, Oil and Gas Science and Technology – Revue d’IFP Energies Nouvelles 66 (3) (2011) 399–422. [73] D.T. Gillespie, A rigorous derivation of the chemical master equation, Physica A: Statistical Mechanics and Its Applications 188 (1992) 404–425. [74] J.B. McDermott, C. Libanati, C. LaMarca, M.T. Klein, Quantitative use of model compound information: Monte Carlo simulation of the reactions of complex macromolecules, Industrial and Engineering Chemistry Research 29 (1990) 22–29. [75] L. de Oliveira, J.J. Verstraete, M. Kolb, Molecule-based kinetic modeling by Monte Carlo methods for heavy petroleum conversion, Science China Chemistry (2013), http://dx.doi.org/10.1007/s11426-013-4989-3 (in press). [76] P. Danial-Fortain, T. Gauthier, I. Merdrignac, H. Budzinski, Reactivity study of Athabasca vacuum residue in hydroconversion conditions’, Catalysis Today 150 (2010) 255–263. [77] I.A. Wiehe, The pendant-core building block model of petroleum residua, Energy and Fuels 8 (1994) 536–544. [78] D.E. Goldberg, K. Deb, J.H. Clark, Genetic algorithms, noise, and the sizing of populations, Complex Systems 6 (1992) 333–362. [79] J.R. Parker, Genetic algorithms for continuous problems, Advances in Artificial Intelligence (2002) 176–184. [80] API 2B2.1, “API procedure 2B2.1 for estimating the molecular weight of a petroleum fraction”. API Technical Handbook, 1987.

220

L.P. de Oliveira et al. / Catalysis Today 220–222 (2014) 208–220

[81] D. Hudebine, J.J. Verstraete, T. Chapus, Statistical reconstruction of gas oil cuts, Oil and Gas Science and Technology – Revue d’IFP Energies Nouvelles 66 (3) (2011) 461–477. [82] L. de Oliveira, Développement d’une méthodologie de modélisation cinétique de procédés de raffinage traitant les charges lourdes, Ecole Normale Superieure de Lyon, Lyon, 2013 (PhD thesis). [83] M.J. Girgis, B.C. Gates, Reactivities, reaction networks, and kinetics in highpressure catalytic hydroprocessing, Industrial and Engineering Chemistry Research 30 (1991) 2021–2058.

[84] D.D. Whitehurst, T. Isoda, I. Mochida, Present state of the art and future challenges in the hydrodesulfurization of polyaromatic sulfur compounds, Advances in Catalysis 42 (1998) 345–471. [85] P. Danial-Fortain, Étude de la réactivité des résidus pétroliers en hydroconversion, Université Bordeaux 1, Bordeaux, 2010 (PhD thesis), http://tel.archives-ouvertes.fr/tel-00839871/ [86] C. López-García, Analyse de la réactivité des composés soufrés dans les coupes pétrolières: Cinétique et modélisation de l’hydrotraitement, Université Claude Bernard, Lyon, 2000 (PhD thesis).