Development of a General Modelling Methodology for Vacuum Residue Hydroconversion.

— Development of a General Modelling Methodology for Vacuum Residue Hydroconversion — This work concerns the development of a methodology for kinetic modelling of reﬁning processes, and more speciﬁcally for vacuum residue conversion. The proposed approach allows to overcome the lack of molecular detail of the petroleum fractions and to simulate the transformation of the feedstock molecules into efﬂuent molecules by means of a two-step procedure. In the ﬁrst step, a synthetic mixture of molecules representing the feedstock for the process is generated via a molecular reconstruction method, termed SR-REM molecular reconstruction. In the second step, a kinetic Monte-Carlo method (kMC) is used to simulate the conversion reactions on this mixture of molecules. The molecular reconstruction was applied to several petroleum residues and is illustrated for an the conversion reactions with the same accuracy as the deterministic approach. The full-scale stochastic simulation approach using molecular-level reaction pathways provideshigh amounts of detail onthe efﬂuent composition and is brieﬂy illustrated for Athabasca VR hydrocracking.

the conversion reactions with the same accuracy as the deterministic approach. The full-scale stochastic simulation approach using molecular-level reaction pathways provides high amounts of detail on the effluent composition and is briefly illustrated for Athabasca VR hydrocracking. Normalised probability for the occurrence of the reaction v r v Probability that the reaction v will occur in the reaction volume V R j Deterministic production rate of lump j t Reaction time V Reaction volume Dt Reaction time step

INTRODUCTION
Nowadays, the world demand for high quality products such as gasoline and diesel continues to increase, while the available crude oil is becoming heavier. In this context, refining processes that convert heavy petroleum cuts into valuable products become increasingly more vital for the refining industry. Heavy petroleum fractions are complex mixtures of hydrocarbons that contain a large number of different chemical species. These hydrocarbons are mainly composed of carbon and hydrogen, but also contain heteroatoms such as sulphur, nitrogen and oxygen. Even metals such as nickel and vanadium are present [1]. Heavy oil conversion processes, such as residue hydrocracking or catalytic cracking, are based on the degradation of the largest molecules by thermal and/or catalytic cracking reactions at high temperature.
In order to improve the performance of conversion processes, reliable and accurate kinetic models are needed. Classic kinetic models for transformation of complex hydrocarbon mixtures generally use a lumped kinetics strategy, in which molecular components are grouped into several chemical families, according to their global properties (boiling point, solubility, etc.). However, the lumping approach assumes that the properties of these families do not change during the reaction, which is not true. Moreover, for heavy hydrocarbon mixtures, the number of families and of reaction pathways turns out to be so vast that this lumping approach is no longer manageable.
The present work is dedicated to the development of a methodology for modelling processes that treat heavy petroleum fractions, and is applied to the hydrocracking of vacuum residues. The proposed methodology allows to overcome the lack of molecular detail of the petroleum fractions and to simulate the degradation reactions of conversion processes through a two-step procedure. In the first step, a mixture of molecules representing the feedstock is generated via a molecular reconstruction method. In the second step, a method called ''Kinetic Monte-Carlo'' is utilised to simulate the degradation reactions of this mixture of molecules.

MOLECULAR REPRESENTATION OF VACUUM RESIDUES
The first step of the proposed methodology aims at determining an accurate representation of the process feedstock from partial analytical data. For this purpose, a molecular reconstruction algorithm, termed SR-REM, has been developed during previous works [21][22][23][24][25]. This algorithm results from the coupling of two methods, Stochastic Reconstruction (SR) [21,26,27] and Reconstruction by Entropy Maximisation (REM) [21,[28][29][30].
During previous works, the SR-REM algorithm was applied to various petroleum vacuum residues [23,24]. In the present work, the reconstruction of an Athabasca Vacuum Residue (VR) will be illustrated.

Brief Description of the SR-REM Algorithm
By using the SR and REM algorithms in series, the SR-REM molecular reconstruction approach generates a set of molecules with the same properties as the petroleum fraction to be represented.
The idea behind the SR step is to generate an equimolar set of molecules from a set of Probability Distribution Functions (PDF) for functional and structural attributes (e.g. molecule type, number of rings, number of side chains, length of the chains, etc.) of molecules. Each petroleum fraction is characterised by its own molecular attributes. By applying a Monte-Carlo procedure, the PDFs' are sampled, according to a sequence defined by a building diagram, to select a set of molecular attributes which are then assembled to obtain the structure of a molecule. The construction of a molecule is repeated N times to obtain a mixture of molecules. For each molecule, its properties are calculated either by direct inspection of its structure or by group contribution methods. Upon creating such a mixture of N molecules, its average properties are calculated from the pure component properties by means of mixing rules and compared to the available analyses through an objective function. The value of the objective function is then minimised by an elitist genetic algorithm that modifies the parameters of the distributions of structural attributes.
Starting from the representative mixture generated by the SR step, the REM step only needs to slightly modify the molar fractions of the molecules to represent the analytical data in an optimal way. Adjusting the molar fractions is done by maximising the Shannon's information entropy [31]. Entropy maximisation ensures that no molecule is preferred over any other, when no constraint (i.e. analytical information) is imposed. On the other hand, the introduction of constraints distorts the uniform distribution of the set of molecules until the entropic criterion is met. As opposed to the stochastic reconstruction technique, the entropy maximisation method uses classical optimisation techniques and its computational effort is much smaller.
For more details concerning to the SR-REM algorithm, the reader is referred to Verstraete et al. [23] and de Oliveira et al. [24].

Application to the Athabasca 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 characterisation cannot be achieved with the current analytical techniques. For this reason, VR streams are generally characterised by their bulk properties obtained from the elemental analysis (C, H, S, N, O), specific gravity, simulated distillation, 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 sulphur, nitrogen and oxygen. Even metals such as nickel and vanadium are present in the form of the porphyrin structures. VR fractions contain both acyclic (paraffins) and cyclic structures, the latter being combined into cores. These cores (poly-cyclic structures) may be composed of naphthenic rings, benzene rings and/or heterocycles, such as thiophenes, pyridines, pyrroles or furans. According to the 13 C NMR data, the alkyl chains are also present in these molecules. The boiling point of the molecules is typically higher than 350°C at atmospheric pressure, which corresponds to molecules more than 25 carbon atoms [24]. 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 [1,20,[32][33][34] have succeeded in obtaining molecular information by combining different analytical techniques, such elementary analysis, distillation, gas and liquid chromatography, mass spectrometry. In our study, 16 molecular attributes listed in Table 1 were selected to create a molecular representation of the VR fractions. 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 characterised 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 characterised 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 by alkyl chains. Their alkyl chains are described by a chain length and a degree of substitution by heteroatoms (S, N or O).
Three types of PDFs' were used as optimisation variables: histograms, exponential and gamma functions. Histograms are applied for molecular attributes containing a narrow range of possible values (less than three), while the exponential and gamma functions describe attributes with a range of more than three possible values. If large values of an attribute are improbable, the molecular attribute is characterised by an exponential function; otherwise this is described by a gamma function. A gamma distribution has two free parameters: a shape parameter (a) and a scale parameter (b). However, in order to reduce the number of parameters, it was assumed that the scale parameter is twice the shape parameter. For a more detailed description of the structural attributes and building diagram used in the present work, the reader is referred to Verstraete et al. [23] and de Oliveira et al. [24].
Concerning the pure component properties, chemical formula, molecular weight and aromatic carbon content ( 13 C NMR spectrum) are obtained by direct inspection of the structure. The normal boiling point is calculated through a group contribution initially developed by Hudebine and Verstraete [28] and extended by de Oliveira et al. [24]. The molecules are classified into the SARA fractions according to its molecular weight and hydrogen content by using a SARA diagram [24]. The SARA diagram is based on the solvent-resid phase diagram proposed by Wiehe [35]. The Athabasca VR was represented by a set of 5 000 molecules in order to ensure an optimal balance between the required CPU time and the accuracy of the feedstock representation. Elemental analyses (carbon, hydrogen and sulphur, 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 [36]. The comparison between analytical and calculated properties of the Athabasca VR is shown in Table 2.
As can be seen in the second column of Table 2 (SR output), a good agreement is observed between the properties of the equimolar mixture obtained after the SR step and experimental data. The elemental analysis, SARA fraction and 13 C NMR are well represented. However, average molecular weight and the simulated distillation show deviations between experimental and calculated values. These deviations may be caused by the group contribution method used to estimate the boiling point. This group contribution method has been developed on the basis of a large database of boiling points, but most of the data available pertain to molecules with less than 42 carbon atoms [24,28]. Therefore, extrapolation to very large compounds may increase the uncertainty of the method. Nevertheless, the application of the REM step leads to a clear improvement in the mixture properties, as illustrated in the third column of Table 2 (SR+REM output). The elemental analysis, SARA fractions separation and 13 C NMR were predicted perfectly. 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.
In conclusion, the results shown above illustrate that the SR-REM algorithm is able to generate a correct molecular representation of the process feedstock to be used in kinetic models.

SIMULATION OF THE DEGRADATION REACTIONS
The second step of the methodology simulates the effect of conversion reactions over the set of molecules generated in the previous step. The VR molecules, especially the asphaltenes, have a complex and large structure. Therefore, each molecule can undergo a large number of reactions and the kinetic network becomes too large to be simulated by traditional deterministic approaches. For this reason, the reactions of VR conversion are simulated using a kMC approach. In the present work, the classical kMC algorithm proposed by Gillespie [37,38], termed Stochastic Simulation Algorithm (SSA), was utilised. The kMC step was validated by simulating the kinetic network of an existing model for the VR hydroconversion.

Stochastic Simulation Algorithm
Instead of tracking the concentrations of the chemical species, the kMC method tracks the reactions of each molecule in a probabilistic way. 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) [38,39]. The CME is given in (1): For most reaction systems, the CME is extremely complicated to solve, both analytically and numerically. This fact motivated Gillespie [37] to develop the SSA. The SSA is a numerical procedure to determine the temporal trajectory of the molecular population in exact accordance with the CME. The theorems that prove the mathematical validity of the SSA to reproduce the solution of the CME have been demonstrated by Gillespie [37][38][39]. The various steps of the SSA are illustrated in Figure 1.
First, all potential reactions are identified from the reaction rules and the structure of each reactant molecule. The normalised probability of each reaction is determined based on the probability that the reaction v will occur in the reaction volume V (r v ), divided by the sum of all probabilities for all possible reactions, as shown in (2): This stochastic rate constant is defined as the product of the number of distinguishable combinations of the reactant molecules (h v ) and the stochastic rate parameter (c v ), as illustrated in (3): For monomolecular reactions, the number of distinguishable combinations of the reactant molecules (h v ) equals 1.
The stochastic rate parameters can be determined from the deterministic rate parameters (k v ), which are based on the average molecular concentration. As shown by Gillespie [37], for monomolecular reactions, k v and c v are equal; for bimolecular reactions, c v equals k v divided by the reaction volume V; for tri-molecular reactions, c v equals k v divided by V 2 .
In the SSA, a first Random Number (RN 1 ) is drawn to determine the reaction time step (Dt) between the current time t and the time at which the next, yet undefined, reaction will occur. This reaction time step is calculated using (4), as proposed by Gillespie [37]: The cumulative probability distribution D R at time t contains all reactions that can occur in the mixture at that reaction time. The probability of each reaction is given by their rate constant and the number of occurrences of this reaction. In the Gillespie's SSA, this cumulative probability distribution D R , is then randomly sampled in order to select the next reaction (l) that will be executed in the next reaction time step. The reaction selection procedure is computed by drawing a second Random Number (RN 2 ) between zero and one. The value of this random number selects the next reaction from the cumulative probability distribution D R . This is shown in (5).
Once the reaction time step Dt has been determined and the reaction has been selected, the reaction system is updated by executing the selected reaction and incrementing the simulation time by Dt. The simulation proceeds event-by-event (reaction-by-reaction) until the final simulation time is reached. The molecule mixture at the end of the simulation now represents the effluent of the process.

Stochastic Simulation of a Lumped Model for VR Conversion
As discussed in Section 1.2, the molecular composition of the heavy petroleum fractions remains unknown. This is why the kinetic modelling of the VR conversion is generally based on a lumping strategy. In lumped models, the molecules are grouped into several chemical families, termed lumps, according to their global properties (boiling point, solubility, etc.). These lumps are then connected by global reaction pathways in order to obtain the kinetic network. The reactions pathways are often described by a power law model (a rate constant and a reaction order) or by Langmuir-Hinshelwood kinetics.
To validate the kMC simulation procedure, a lumped kinetic model for VR hydrocracking [40] was simulated in a stochastic way (kMC simulation) and in a "classical" deterministic way, which corresponds to the numerical resolution of the mass balances. In our example, the kinetic network (Fig. 2) developed by Schweitzer and Kressmann [40] was used to simulate the reactions. Their  kinetic network is composed of 4 lumps that correspond to 4 distinct petroleum cuts: vacuum residue, vacuum gas oil, atmospheric gas oil and naphtha. This kinetic model was used to simulate a perfectly mixed batch reactor by solving the set of Ordinary Differential Equations (ODEs) given by the continuity equations for each lump: In our study, a set of 1 000 molecules representing the Athabasca VR has been generated by the SR-REM algorithm. The molecules were then classified a posteriori according to their boiling point into the lumps of the kinetic model, as illustrated in Table 3. The rate parameters for our stochastic simulation were obtained by using Equations (2), (3) and (4) and the values of the deterministic rate constants given in [40].
Both the deterministic and stochastic simulations were performed up to a conversion of 70% of the vacuum residue (lump A), which corresponds to the value that is typically obtained at the outlet of the VR hydrocracking unit. The comparison between the deterministic and the stochastic simulations is illustrated in Figure 3.
In Figure 3, the blue lines correspond to the vacuum residue (A), the green lines correspond to the vacuum gas oil (B), the red/orange lines represent the atmospheric gas oil (C) and the grey lines represent the evolution of the naphtha fraction (D).
For a single stochastic simulation, slight deviations are observed between the kMC method and the deterministic approach. This deviation is due to the stochastic nature of the kMC method. This result shows that the number of reaction sampling steps provided by a single stochastic simulation is not large enough to obtain an accurate representation of the VR conversion. However, by increasing the number of simulations and thus increasing the number of reaction sampling steps, a clear improvement is observed in the average value of the stochastic simulations with respect to the deterministic simulation. This is shown in the graphs "Average of 10 simulations" and "Average of 100 simulations" of Figure 3. These results illustrate that the kMC method is able to simulate the conversion reactions with the same accuracy as that of deterministic simulations. It needs to be stressed that the SSA reproduces the exact solution of the CME (Eq. 1), while the deterministic model solves a set of continuity equations (ODEs) given by Equation (6). Despite the fact that both approaches are based on very different hypotheses, they yield exactly the same solution.
As any solution technique, the SSA has advantages and limitations. The SSA is exact and mathematically rigorous. Unlike standard methods to numerically solve systems of ordinary differential equations, the stochastic simulation algorithm never approximates infinitesimal time increments dt by finite time steps Dt. This is especially advantageous when dealing with systems in which the molecular composition can change suddenly and sharply with time. The SSA is also very easy to code for any specified set of chemical reactions, and the corresponding computer code will normally require very little computer memory. When ensemble averages, such as means and variances, cannot be readily calculated from the CME, the stochastic simulation approach offers a universally straightforward way of numerically estimating them.
On the other hand, the SSA also has a number of limitations. For a single simulation, the computing time is directly proportional to the number of samplings of the cumulative probability distribution D R . Indeed, the denominator of Equation (4) depends on the number of reactions M and the reaction probabilities r i . For a given reaction system, the computing time is therefore  Figure 2 Kinetic network of VR hydroconversion [40].
roughly proportional to the total number of reactions, and hence approximately to the total number of molecules. Systems with very high reactivities and/or a very high number of molecules, each of which can undergo a large number of reactions, will not be simulated very efficiently using the SSA method. As with any stochastic procedure, the SSA requires quite some CPU time, mainly due to the fact that a relatively important number (100 to 1 000) of simulations need to be run before obtaining a stable ensemble average. This is somewhat offset by the fact that the SSA typically occupies only a small amount of computer memory, and that the various simulations can easily be run independently. Hence, the SSA computer code is therefore highly parallelisable.
On the numerical side, it needs to be stressed that the SSA requires a reliable uniform random number generator. Fortunately, many of the generators available today are quite good, as they run fast, have very large periods, exhibit low correlation between successive values, and are equidistributed even in high dimensions.
Given the molecular representation of a VR fraction generated by the SR-REM approach, one can also utilise the chemical structure of the molecules. Indeed, in the example above, the 1 000 generated molecules were simply classified according to their boiling point. The latter therefore defines their reaction pathways and their reactivity, ignoring the actual structure of each molecule. Instead of using global reaction pathways between lumps and average reactivities, the information available through the molecular representation can now be used to exactly identify the various molecular-level reactions of each molecule, and to attribute a molecule-dependent reactivity to each pathway. Such molecular-level simulations of hydrotreating and hydroconversion reactions are presented in detail elsewhere. The hydrogenation of benzene rings, dehydrogenation of naphthenic rings and hydrodesulphurization of thiophenes were added to the kMC algorithm, and the methodology was successfully applied to the hydrotreating of LCO gas oil [41]. Cracking reactions have also been included to be able Comparison between the deterministic simulation and the stochastic simulation.
to simulate VR hydrocracking [42][43]. Applying the full stochastic simulation methodology allows to utilise the molecular detail provided by the stochastic reconstruction algorithm and to generate much more details on the evolution of the composition. This is illustrated in Figure 4 and in Figure 5 for Athabasca vacuum residue hydrocracking. More detailed information on the implementation of the various reaction types can be found in [42]. The proposed overall approach has several advantages: first of all, the complex feedstock is represented by means of a set of molecules. The molecular description of the reaction system is retained throughout the entire reactor simulation. Finally, the simulation of the reac-tions can be performed without a pre-defined kinetic network. In the SSA method, the kinetic network is generated "on-the-fly" as the reactions proceed.

CONCLUSIONS
A novel two-step kinetic modelling strategy for heavy oil conversion processes has been presented in this paper.
The first step transforms the available analytical information into a set of molecules whose properties are the same as those of the process feedstock. This transformation is carried out by using the SR-REM molecular reconstruction algorithm, which results from the coupling two independent methods, Stochastic Reconstruction (SR) and Reconstruction by Entropy Maximisation (REM). The SR method is first applied to generate an initial equimolar set of molecules whose properties are close to the analytical data. The molar fractions of the generated molecules are then adjusted to obtain a mixture whose properties are closely aligned to those of the feedstock.
The molecular reconstruction step was applied to an Athabasca Vacuum Residue. As shown above, the generated set of molecules has the same properties as the feedstock. The results prove that the SR-REM algorithm is able to generate a correct representation of the VR feedstock.
The second step of the approach consists in simulating the effect of the conversion reactions on the set of molecules generated in the previous step by using a kinetic Monte-Carlo (kMC) method. To this aim, an algorithm based on the stochastic simulation algorithm (SSA) developed by Gillespie has been utilised.    Results of the full stochastic simulation at the molecular level: evolution of the molecular weight distribution of the hydrocarbon effluent during Athabasca vacuum residue hydroconversion.
In this work, the kMC simulation step was validated by simulating a lumped kinetic model for VR hydrocracking. The simulation results of the stochastic approach were compared to the results of the deterministic simulation, illustrating that the kMC method simulates the conversion reactions with the same accuracy as the deterministic approach, despite the fact that both approaches are based on very different hypotheses.
As the SR-REM algorithm provides a molecular representation of a VR fraction, one can also utilise the chemical structure of the molecules to simulate the conversion reactions by means of molecular-level reactions of each molecule and by assigning a moleculedependent reactivity to each pathway. This approach provides high amounts of detail and was illustrated for Athabasca VR hydrocracking.
The proposed overall approach has several advantages: first of all, the complex feedstock is represented by means of a set of molecules. This representation retains a molecular description of the reaction system throughout the entire reactor simulation. Moreover, the simulation of the reactions can be performed without a pre-defined kinetic network. In the kMC method, the kinetic network is generated "on-the-fly" as the reactions proceed.