Comprehensive Residual Oil Hydrocracking Reaction Kinetic

Article Cite This: Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

pubs.acs.org/IECR

Comprehensive Residual Oil Hydrocracking Reaction Kinetic Modeling Combined with Effective Phase Equilibrium Calculation Antti Pyhälahti,* Jaana Kanervo, and Susanna Kuitunen

Downloaded via KAOHSIUNG MEDICAL UNIV on August 9, 2018 at 00:47:41 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

Neste Engineering Solutions Oy, P.O. Box 310, FI-06101 Porvoo, Finland ABSTRACT: A model for a residual hydrocracking unit combining a comprehensive kinetic model with efficient phase equilibrium calculation has been created. For chemical reactions the model uses actual well-defined chemical components (base components) and kinetically limited stoichiometric reactions. For the phase equilibrium calculations base components are combined to lumps in order to accelerate the calculations. The physical properties of lumped component are based on those of the base components. The novelty involved in the model is that the physical properties of the lumps are continuously updated when the composition of the base components changes. The advantage of the model is that it avoids the tedious iteration of phase equilibrium with the full base component set, whereas allows full component set in kinetics calculation. The model is computationally very efficient and is demonstrated to work with a system of 19 444 base components and 486 079 reactions. For vapor−liquid equilibrium (VLE) calculations the components are lumped to 42 lump components, and for vapor−liquid−liquid equilibrium (VLLE) calculations they are lumped to 162 components.

■

INTRODUCTION Major challenges in the oil refining industry have been to maximize the distillate yield of crude oil while meeting more stringent product specifications. Residual hydrocracking (RHC) provides both simultaneous conversion of vacuum residuum into more valuable distillate products and the suppression of impurities in those. Understanding and optimization of residual hydrocracking call for experiments, analytics, and modeling, which have been challenged by the chemical and technological complexity of the system. Typical RHC reactors are ebullated bed reactors, containing vapor, liquid, and solid phases in continuously changing proportions and compositions, and the solution involves iteration. Thus, when simulating such a reactor, the reaction rates and phase equilibriums need to be solved in many points. Due to the complexity of the system, the efficiency of the algorithm is of crucial importance for an industrially useful system. This article presents a modeling methodology for a residue hydrocracking (RHC) reactor involving three phases (solid catalyst, liquid, and vapor) capable of handling a large number of chemical components and their reactions. The reaction kinetic model capturing conversion of well-defined reactants to well-defined products aims to mimic the true reactions as closely as possible, but lumped components, whose properties are derived from the well-defined components, are used in phase equilibrium and other calculations not involving chemical reactions. The difference from many lumped models available is that the real composition of each stream in the simulation is stored, and whenever physical operations are performed to a stream, the component lumping is performed and the physical properties of the lumps are upgraded. Obviously the question arises, is the benefit of lumping over © XXXX American Chemical Society

performing the vapor−liquid equilibrium (VLE) calculations with the complete component set lost along the computing necessary for performing the upgrading of the properties? However, this seemed not to be the case. This approach is by no means limited to the present system, but can be applied to modeling of many oil refining processes treating complex feeds and involving phase equilibrium calculation. Processing of biobased components involves often very complex mixtures, and this approach could be applied to them, as well. The approach has been selected because there are a huge number of individual chemical components in heavy oil mixtures. With a Fourier transformation ion cyclotron resonance mass spectrometer (FT-ICR-MS), several thousand gross formulas have been detected. Considering the size and complexity of the molecules, the number of possible structural isomers is easily up to hundreds of thousands. Because most of those components may undergo several reactions under RHC conditions, the number of reactions involved is likely to be many millions. Despite the large number of components, distillation of such mixtures has been successfully simulated using the so-called lumped component approach. The idea is to replace the actual components of the mixture with a limited number of lumped components, each of which represents a relatively narrow distillation fraction of the real mixture. The concentrations and properties of the lumped components are adjusted to so that Received: Revised: Accepted: Published: A

April 27, 2018 July 21, 2018 July 24, 2018 July 24, 2018 DOI: 10.1021/acs.iecr.8b01834 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

the heavy hydrocarbon mixtures containing polyaromatics and sulfur, nitrogen and oxygen heteroatoms is rather scarce. In the literature there are numerous models for hydrocracking of various heavy oil fractions. Some reviews have been published.2,3 The kinetic models of the reactions can be divided according to their fundamental approach:3 • Lumped models are where the feed is divided into a number of lumps which are then converted to product fractions. Usually the number of lumps is fairly small, varying from one to a few dozen. • Mechanistic models are where the rates of conversion are based on the actual reaction mechanisms including the reaction intermediates. • Molecular models are where the reaction mixture is described by means of molecules and the model includes kinetics for the conversion of reactant molecules to the product molecules. The lump models are based on dividing the mixture into a limited number of cuts and the conversion of the feed cuts to the product cuts is based directly on laboratory data. The molecular composition of the mixture is not considered in the model. The lump model is simple, and if sufficient laboratory data is available, its parameters can be identified in a relatively straightforward manner. A substantial advantage is that the analysis requirements are relatively easy to meet even if the necessary pilot runs may be time-consuming and expensive. Everyday standard characterization methods applied at oil refineries are adequate. Because the molecular composition of the lumps is not taken into account, it is unlikely that complex analysis technologies would bring much additional value. On the other hand, the predictive power of such models is minimal. If reaction conditions are outside the region covered by the experimental data and especially if the feed composition changes, new laboratory tests need to be performed and new model parameters fitted. Several lump models for hydrocracking have been published, but many of them are for heavy vacuum oil and lighter mixtures. For heavy component hydrocracking some kinetic lump models have been also presented.4−7 The mechanistic models take into account the detailed reaction mechanisms for thermal and catalytic reactions. The initial reactants and final products are molecules, as in the molecular model, but such a model includes the detailed reaction mechanistic steps, for instance production and consumption of the reaction intermediates such as free radicals. Such an approach is very comprehensive, but tedious. When the number of components increases, the complexity of the model increases substantially. Mechanistic models have been applied to systems involving light hydrocarbons, but as for residual oil, considering the size and complexity of the components, this approach seems infeasible with current computers. The molecular models are based on overall reactions where the reactant molecule or molecules convert to product molecules, but intermediate steps belonging to the fundamental reaction mechanism are not included. The reaction kinetics is based on the net reaction equation where reactants are converted to product molecules according to appropriate overall kinetic equations. The principal advantage of such models is that they take into account the structure and reactivity and most likely reaction paths of the molecules. Consequently, the predictive power of such models is superior

the mixture has properties similar to the actual mixture. The most usual properties used as basis of lumping are the normal boiling point (NBP) and density, but other properties can be used as well if measured data is available. Using this approach, the simulation of the distillation columns can be performed efficiently. Whereas physical bulk properties and separation processes can be well described by the lumped component approach, it has serious limitations when compositional changes via reactions need to be taken into account. Chemically very different components may have volatilities near each other and consequently behave in the same way in operations dependent on vapor−liquid equilibrium (VLE) properties such as distillation. However, in a reactive system their behaviors may differ from each other drastically. Thus, a simple lumped component model cannot describe properly a complex reactive system. A common method in modeling reactors treating heavy fractions in oil refinery processes is to implement them as black-box models, which do not contain any information about reaction mechanisms. The properties of these streams representing reaction products are based on laboratory or pilot results. Only the phase equilibrium operations are actually calculated. There appears to be a gap in existing modeling methodologies to take advantage of a proper description of the chemistry of the system without sacrificing the efficiency of the lumped components approach in physical operations. A serious obstacle for developing and applying more detailed models for such mixtures has been that it is difficult to get a detailed analysis of a mixture comprising a huge number of components. Due to the complexity of the mixture, the traditional approach has been to give up on analyzing the actual components, except for a limited number of light end components, and use various characterization methods for the bulk of the mixture. Those characterization methods include distillation curves, elemental composition, solubility classes, and various, technically useful bulk properties of the mixture, such as flash point or coke formation. However, they do not provide exact information about the actual chemical composition of the mixture. Recently much more powerful methods for the analysis of chemical composition have become available, such as the aforementioned FT-ICR-MS, NMR, and HPLC methods. Recently even real portraits of asphaltene molecules have been produced by atom force microscopy.1 However, if there is no connection from these analyses to the operation of the refinery, there is not much incentive to apply them in industry. The result is some kind deadlock: refined analysis methods are of very limited use, if their results cannot be applied in practice to improve plant operations via quantitatively predictive process simulation models and without analysis results more fundamental calculation methods cannot be established and applied. This work attempts to bridge the gap by showing that process simulations combining efficient physical property calculations along with detailed reaction chemistry taking advantage of modern analytical resources are possible in practice and to develop tools for that.

■

REVIEW OF PREVIOUS WORK The distillation models applying lumped components are wellestablished, and thus they are not reviewed further. Contrary to that, the knowledge concerning the reactions taking place in B

DOI: 10.1021/acs.iecr.8b01834 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

base component composition changes. A typical property used as the basis for lumping is the NBP, which is useful for operations involving VLE calculations. However, the user may use any relevant property as the basis of lumping. Thus, for system where both vapor−liquid and liquid−liquid equilibriums were essential, the base components were regrouped to lumps based both on NBP and solubility factor. Thus, in spite of large number of components the model is conceptually simple, which is considered an advantage The reaction equations are normal stoichiometric equations with normal kinetics, and lumping means mixing of the components of the lump. Naturally the feed composition, component properties, and kinetics need to come from somewhere, but composing of the feed and building the kinetics is not part of the actual simulation model, although they are taken as given when solving material and energy balances. Nevertheless, a brief outline of methods used for setting up the feed is given below. Base Component Set. The base component set was generated based on the available information about the typical molecules present in the residual oil. The base component set is based on the core-side-chain model called the modified Yen model.18−21 Also, the concept of “naphthalene zigzag” reaction path22 has affected the selection of the basic polyaromatic hydrocarbons (PAHs). The resulting series of polyaromatic structures are of pericondensed type. The stability of pericondensed PAH molecules is found to be better under the hydrodemetallization conditions than that of katacondensed PAH molecules,23 which gave additional support to selecting this series as the basic structure of PAH cores. However, a number of kata-condensed PAH cores were included, as well. Direct evidence about the structure of asphaltene molecules is provided by atomic force microscopy.1 The base components consist of a core part to which is connected a certain number of carbons in alkyl side chains. The most common structure in the component set is the socalled island type, where there is only one aromatic ring structure in the molecule. Some structures are of archipelago type having two aromatic ring structures in the molecule. The cores consist of 0−11 rings containing 0−34 carbon atoms, which are either aromatic or naphthenic. In total the core contains one or two sulfur atoms and zero to one nitrogen atom. The total number of different cores is 417. The minimum number of carbons in the alkyl side chains is zero. The maximum number depends on the size of the core, because the maximum number of carbons in a single molecule is limited to 70. Consequently, the maximum acceptable number of carbons in the side chains is 70 minus the number of carbon atoms in the core. Each side chain carbon number covers all isomers with the same carbon number and molecule core. The basic component set resulting from these cores and side chains comprises 19 444 components. For example, a molecule may contain a hexahydropyrene core to which is attached alkyl side chains having 23 carbons. An example of molecules belonging to this category is heptyl-dioctyl-hexahydropyrene (Figure 1). The number and length of the side chains has an influence on the size distribution of the reaction products of the side chain cracking reactions. It is assumed that the number and length of the side chains attached to the cores are statistically distributed and product distribution is derived from that. The basic thermodynamic properties (critical properties, ideal gas heat capacity function, and enthalpy of formation) of

over the lump models, providing that the analytical, kinetic, and physical property data are sufficiently comprehensive and accurate. When moving from lumped models to molecular models, the necessary information concerning the composition of the mixtures, properties of the components, and kinetics of the principal reactions increase tremendously. Relevant data especially for the heavy components is very scarce. Monte Carlo method has been applied by a number of authors to develop component sets having proper properties and reaction networks with usable kinetics.8−11 Nevertheless, the positive aspect is that the data concerning real components’ properties and their reactions are not bound to a single system as is the case with lumped components, but once collected and stored to a databank, they can be utilized for handling new cases. Some models accounting for molecular structure have been published, e.g., for gasoil hydrotreating,12 for fluid catalytic cracker,13,14 and for hydrotreating of light cycle oil from fluid catalytic cracking and hydrocracking of vacuum residues.3 However, simultaneous VLE calculation was not included. A typical problem in residue hydrocracking is formation of coke to the equipment and sediment into the product. The most problematic parts of the process are not in the reactors but in the consequent separation stages where hydrogen pressure is reduced substantially.15 Studies about the reaction mechanisms in vacuum residue hydrocracking have been published.16,17 Availability of information on chemical species identified as potential sediment precursors, and on their reaction chemistry, holds a promise for better kinetic modeling of adverse sedimentation reactions. In practical engineering work the lump models seem to be most common. Considering the nature of the standard characterization methods applied to heavy oil analysis, the popularity of this approach is easy to understand. Typically the available data consists of distillation curves, elemental analysis, and solubility data. Moreover, the heavier and more complex the cuts become, the less detailed information is available. Thus, there is no sufficient data for creating any detailed model without adopting more powerful analysis methods.

■

DESCRIPTION OF THE MODEL The Model Treats Components at Two Levels: Base Components and Lumped Components. A rigorous but computationally effective model is obtained using component sets at two levels. The primary level is a set of the conventional components involved in reactions. These base components (as called in the present paper) have a well-defined chemical formula and physical properties. A model with almost 20 000 base components has been successfully tested, and this number seems not to be the upper limit. Even if the number of base components is so high, each base component still represents several components in the residual oil. The secondary level is set of lumped components of the base components. Lumping is done according to a suitable physical property or properties. In a typical case there are a few dozen lumped components. In the simulator input each base component is attached to a certain lumped component. The concentrations of the base components vary during the simulation because of the reactions and mixing taking place in the mixture. Consequently, the base component compositions of the lumped components change, as well. Thus, the properties of the lumped components must be updated using base components’ concentrations and properties whenever the C

DOI: 10.1021/acs.iecr.8b01834 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

components need to be updated. The same applies whenever reactions occur. One of the positive observations was that the upgrading of the lumped component properties did not slow down the simulation significantly. Reactions of Base Components. In this model the base components are well-defined chemical compounds such as hydrogen, butane, benzene, and butylbenzene. Consequently reactions between them have a proper stoichiometry. For example, the mentioned components may react according to the equation Figure 1. Heptyl-dioctyl-hexahydropyrene.

hydrogen + butylbenzene → butane + benzene

A proper stoichiometry enables the estimation of the heat of reaction from the properties of the base components participating in the reaction. Currently the selection of the reaction consists of the following types: • cracking of alkyl side chains (thermal and catalytic) • cracking of archipelago type cores (thermal and catalytic) • dearomatization/aromatization of ring structures (catalytic, equilibrium) • opening of naphthenic rings (thermal and catalytic) • hydrodesulfurization (catalytic) • hydrodenitrogenation (catalytic) • regrouping of alkyl radicals resulting from cracking reactions Most of the reactions were treated by neglecting the possible radical reaction intermediates. However, the main mechanism of alkyl cracking is thermal and it is likely that a substantial portion of the immediate products are radicals which then react further in a very complex way. The modeling approach will be explained below. Kinetic data for large polyaromatic hydrocarbon molecules is scarce and does not cover the temperature and pressure region relevant in RHC processing (around 700 K and 15 MPa hydrogen pressure). Thus, the existing data which is mainly measured for smaller molecules and at lower temperature regime must be extrapolated considerably. Some kinetic data concerning the hydrogenation of benzene, naphthalene, and phenanthrene was found28 in addition to data concerning the hydrogenation of pyrene.11 The rates of reaction depend inevitably on the activity of the catalyst. The catalyst activity is difficult to predict theoretically, especially that of deactivated catalyst. Moreover, it is not constant in the reactors; thus it needs to be adjusted based on experience. The kinetic expressions used were the simplest possible, i.e., the homogeneous power-law type. Aromatization and dearomatization were treated as kinetically limited equilibrium reactions. Hydrogen partial pressure was used as basis of the equilibrium. Many reactions are thermal and independent of catalyst. For such reactions the power-law kinetics is acceptable. For catalytic reactions such as hydrodearomatization more advanced models such as the Langmuir−Hinshelwood mechanism might be better. Due to a lack of data concerning most of the components, such a refinement was considered useless. The model was not satisfactory without allowing recombination of the alkyl radicals in the reaction mixture. After some tests a simple approach was selected. In each calculation step in all alkyl chain cracking reactions, the formed alkyl molecule was placed in a “radical pool”. When the rates of all nonradical

the base components are necessary for estimating the physical properties of the lumped components. Unfortunately, measured data about physical properties of the large organic molecules found in residual oil is very limited. Thus, most of the properties must be estimated using group contribution methods. The methods used were the Marrero−Gani method24,25 and the Joback method.26 For the gases and light hydrocarbons up to C16 components, data was available in the literature and it was used instead of estimation. The PCSAFT parameters were estimated based on the method of Gonzalez Rodriguez27 with modifications for heteroatoms. The solubility factors were estimated using Hildebrand solution theory and the PC-SAFT equation of state for predicting the heat of vaporization. The actual VLLE calculation was based on the PC-SAFT equation of state; the solubility parameters were used only as a basis of classification when performing lumping for the sedimentation simulation. However, the architecture of the model is not dependent on the method used for generation of the feed composition and properties. An alternative approach for feed composition generation would be the Monte Carlo methods presented in the literature.3,8−11 Lumped Components for Phase Equilibrium. The phase equilibrium calculations are performed with lumped components. Each base component is explicitly attached to a certain lumped component in the simulator input. Thus, the user can select any basis for the grouping considered relevant. It is assumed that in the phase equilibrium operations all basic components belonging to a certain lumped component distribute like that lumped component. Thereby the tedious iteration of equilibrium composition between thousands of basic components is avoided. For VLE operations the NBP is a natural basis for lumping. For operations involving both VLE and VLLE operations, the lumping can be based for example on the NBP and the solubility factor. The complication resulting from the two levels of components is that for all streams present in the simulation both the lumped component composition and base component composition must be stored. The benefit is that the lumping and delumping are lossless with respect to the base component composition. Thus, every individual stream has its own base component composition and corresponding lumps, which are different from the lumps of other streams with different base component composition and during the iteration those lumps change whenever the base composition changes. Moreover, physical properties of the lumped components must be updated when the base component composition of the stream changes. For example, when two streams with different base compositions are mixed, the lumped components of the resulting stream have different physical properties than the original streams. Thus, the properties of the lumped D

DOI: 10.1021/acs.iecr.8b01834 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

e.g., 10 slices, each slice contains exactly 10% of the total catalyst load of the reactor. In an ebullated reactor the catalyst is in fluidized state. Consequently, the heights of the slices depend on flow conditions and must be reevaluated in each iteration round. This implies that the liquid and vapor holdup of each slice varies, as well. This has an impact on the relative rates of thermal and catalytic reactions. The remaining empty space above the catalyst bed and up to the liquid surface is divided into slices, and thermal reactions may proceed in these sections as well. The reason to select this method of dividing the reactor volume is to avoid discontinuities near the top of the fluidized bed. There is a sharp change in reaction rates at the catalyst bed top, and when the limit between the sections is set at the top of the bed, this change occurs at the section limit. If the heights of the slices were fixed, there would be at least one slice where this change occurs within the slice. The effluent separation is done at the outlet flash and the liquid phase is divided between recycle and product streams in a fixed ratio. The ratio needs to be adjusted so that the catalyst bed is appropriately expanded. Reaction Stage Internal Structure. Each reactor slice consists internally of hydraulic, VLE, and reaction calculation steps as schematically presented in Figure 3. They are solved sequentially, not simultaneously, which simplifies the calculations substantially.

reactions had been determined, a given fraction of carbon atoms of these molecules in the radical pool were counted and redistributed to a new set of alkanes. The chain lengths of the new molecules were determined by a given distribution (e.g., Rayleigh distribution) and given parameters for that distribution. The rest of the molecules in the radical pool were assumed to remain as they were. Naturally the total amount of hydrogen bound into these molecules depends on the size of the molecules. Thus, after solving the final molecule sizes the amount of hydrogen consumed by the reactions was adjusted. Reactor Model. The modeling methodology was applied to an ebullated bed reactor train treating residual oil. A schematic cross section of such a reactor is presented in Figure 2. A schematic figure showing the how the reactor is divided

Figure 2. Schematic presentation of RHC reactor and simulation scheme.

into computational sections is also shown. The feed oil and hydrogen are fed to the bottom chamber below the catalyst bed. There they mix with the recycle stream and flow through the distributor grid, above which there is the fluidized catalyst bed. In the upper part of the reactor is empty space above the catalyst bed where catalyst disengages from the fluids. Furthermore, at the top there is a vapor−liquid separation section where vapor and liquid are disengaged. Most of the liquid is recycled to the bottom section via the recycle pipe to the suction of the recycle pump. The flow rate through the bed determines the bed expansion. The part of the liquid which is not recycled and the vapor phase are recovered as two-phase flow via the reactor product line. In the simulation the reactor is treated as series of adiabatic reactive flash stages. The catalyst bed is divided into slices, which are considered fully mixed. The division is made so that the first slice is the bottom chamber of the reactor below the catalyst support plate where the feed mixture and liquid recycle are mixed together and thermal reactions may proceed. The catalyst containing part of the reactor is divided so that in each slice there is the same amount of catalyst. Thus, if there are,

Figure 3. Blocks of reaction/flash stage.

Hydraulic Calculation Step. The hydraulic model is solved first, because bed expansion is needed for reaction rate calculations. The hydraulic model is approximate, and its principal goal is the calculation of the local bed expansion and the estimation of the liquid holdup. A fully satisfactory model for the bed expansion was not found from the literature. Here bed expansion is calculated from correlation with one of them29 even if this case is outside of the regime for which it was developed. Flash Stage. The VLE calculations are performed using the lumped components as explained earlier. Both ideal thermodynamics and the Soave−Redlich−Kwong equation of E

DOI: 10.1021/acs.iecr.8b01834 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

and quantitative description of the sedimentation for the used tuning data. The Algorithm. The reaction section is solved using the forward Euler method. Its great advantage is avoiding the evaluation of Jacobian matrix for the huge equation system. The Euler method is naturally rather inaccurate, but on the other hand, the flow pattern through the fluidized bed is probably not ideal plug flow. The real drawback of the Euler method is its limited range of stability. However, no substantial problems were encountered with the present kinetic model. Convergence of the recycle stream composition was achieved using the Wegstein method. Implementation of the Model. The reactor model was implemented at an in-house simulator of Neste Engineering Solutions called Flowbat. Flowbat is a sequential modular process simulator. For this model many parts of the simulator code had to be modified to accept the exceptional amount of components (20 000 base components) as well as the two levels of components (base and lumped) as explained. The lumped components were given the status of normal simulation components, and special treatment was implemented for the base components. Handling of the physical properties had to be revised thoroughly. The properties of the lumped components change, whenever chemical reactions modify the composition of the base component set. Thus, routines for upgrading the physical properties of the lumped components were written. The stoichiometry and kinetic data were placed in a separate reaction database file. The separator section following the reactors was placed in a separate simulation stage which could be solved with a commercial simulation program because the reactions were not taken into account. Examples of the Results. Results of simulating a hypothetical RHC reactor train with two hypothetical vacuum residue feeds at two temperature levels are presented. Thus, the equipment was the same in each case, and the variables were feed and temperature. The sample case involves 19 444 base components and 486 079 reactions. The base components are grouped to 42 lumped components for VLE calculations. The basis of lumping has been the NBP. The first lumps in the light end contain only a single base component each, but when the molecular masses increase, there are more components having boiling points near each other. Thus, in the highest boiling lumps each of the groups contains several hundred components. The volatility of the last groups is so low that they are almost absent from the vapor phase. The algorithm fulfilled the expectations considering efficiency of computation. For example, after changing the VLE model from ideal to Soave−Redlich−Kwong, achieving a converged solution required 244 s of CPU time with an 8 year old desktop computer based on an AMD Phenom II X955 processor. Thus, the requirements for computing hardware are moderate. The sedimentation calculation was performed in a separate simulation run. For that the components were reordered to different lumps based on NBP and solubility factor. Some key properties of the feeds are presented in Table 1. Distillation curves of the feeds and products are presented in Figure 4. Each of the curves consists actually of 19 444 points plotted on cumulative weight−NBP coordinates. The main difference between the feeds is that the first one contains more aromatic carbons and sulfur, about 32.4% aromatic carbons of all carbon atoms and 3.5% sulfur. The second feed contained

state were applied to VLE calculations. The results with these models were clearly different, but the impact on overall time of the solution was minor. The results in tables are from simulations with the ideal model. In order to accelerate the calculations further, the components are ordered in the simulation input so that the base components belonging to a certain lumped component are stored in a continuous block in computer memory, which improves performance when properties are updated. Reaction Stage. The reactions rates are evaluated last. There are both catalytic and thermal reactions which are assumed to proceed in the liquid phase; thus both the amount of the catalyst and the liquid volume of the section affect the reaction rates. Because vapor occupies part of the volume of the section and thereby reduces the liquid volume, it has an impact on the reaction rate. The reactions involved are not fast; thus separating VLE calculation and reaction calculation from each other is not expected to cause prohibitive errors. At least inside the catalyst pores, it is not even self-evident that vaporization and reaction actually proceed simultaneously. Sedimentation Model. A specific challenge with RHC operations is the sediment formation especially during the physical separations after the actual hydrocracking.30 The exact reasons and chemical precursors for sediment formation remain unclear, but the phenomenon typically limits the RHC conversion target attainable for a given feed type. Presumably asphaltenes form nanoaggregates via van der Waals interactions that subsequently undergo chemical transformations that lead to species that can be classified as coke. Nanoaggregate coalescence is believed to lead to liquid−liquid separation of a mesophase that is insoluble in the main bulk of liquid. A model that could predict sedimentation tendency or even quantitatively predict the extent of sedimentation would considerably improve plant operations. For detecting the possibility of separation of sediment from the reactor effluent, a process model based on lumped components was developed. The sedimentation kinetic model was incorporated into a process simulation model of a separation section (Aspen Plus) consisting of several physical separation steps leading to a variety of hydrocarbon products. The input stream for this model was the simulated RHC reactor effluent, grouped into 162 lumped components according to solubility factor and NBP. The number of lumped components was kept this low to allow liquid−liquid equilibrium calculations to be performed efficiently enough. The physical properties within the lumped components were derived similarly from those of the base components as in the RHC simulations. Unlike in the RHC reactor part, the stream composition changes in a separation section involve mainly physical processes such as flashes and distillations instead of chemical reactions. The sedimentation was described as a kinetic process generating a new component from the existing ones, which represents formation of nanoaggregates from asphaltenes. The functional form of the kinetic model was established based on data analysis of the behavior of components during the separations. The parameters of the thermodynamic model (PC-SAFT31) were tuned so that the generated sediment components separated from the bulk liquid as a second liquid phase under the chosen conditions. Kinetic parameters of the sedimentation model were tuned to match modeled sediment amount with the amounts found in plant operations. Altogether the model provided a valid qualitative F

DOI: 10.1021/acs.iecr.8b01834 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research

components. They are based on mass fractions and NBPs of all base components. Table 2 presents some key numbers about the changes taking place in various cases. The reduction of the components boiling above 783.15 K is larger in the low aromatic feed than in the high aromatic feed. Likewise the cracking rate at high temperature (HT) is higher than at low temperature (LT). The higher cracking rate of low aromatic feed can be explained by the larger proportion of aliphatic side chains which are easier to break in thermal cracking. An interesting result is the higher conversions of aromatic carbons in the low aromatic feed. The reason is obviously the character of the aromatic cores. The proportion of more stable large pericondensed types is higher in the high aromatic feed. The sediment formation is clearly higher with high aromatic feed than with low aromatic feed, because sediment is formed from large polyaromatic molecule cores. At both temperatures the conversions of aromatic carbon types are lower for the high aromatic case than for the low aromatic feed. This results from the larger average size of the polyaromatic cores which makes them less reactive. On the other hand, the relative cracking rate of the aliphatic side chains is higher. The sedimentation model was applied with the simulated RHC effluent as feed (four cases as above). The conditions when entering the separation section were 688 K and 16.3 MPa. First a high-pressure flash was simulated for the feed, and thereafter a midpressure flash (688 K and 2.3 MPa) was evaluated. Sediment formed in the system according to the defined first-order kinetics assuming coalescence of eight asphaltene entities from one lumped component group. The separation of mesophase took place only in the stream composition corresponding to the midpressure flash, where the key changes in the composition were the escape of dissolved hydrogen and other light components from the liquid stream. The amounts of formed sediment per total feed (ppm) are provided in Table 2. With the present model formulation and parameter tuning, the amounts of sediment vary considerably with the feed type entering the separation section. This example serves as a qualitative demonstration of the usability of the presented modeling methodology as applied to the description of the sedimentation.

Table 1. Properties of Feeds of the Simulated Cases Case

Aromatic carbons of all carbons in oil feeda (%)

av MW of feed

783 K+ fraction in oil feeda (%)

Sulfur content of oil feeda (%)

32.40

570.9

89.60

3.51

17.10

620.2

84.80

1.01

high aromatic (HA) low aromatic (LA) a

Weight basis.

Figure 4. Distillation curves of total feeds and products of demonstration cases on weight percent basis, with vertical axis temperature in kelvin. HA = 32.4% aromatic carbons, LA = 17.1% aromatic carbons, LT = 673.15 K, and HT = 698.15 K.

■

CONCLUSIONS The goal of this work has been to develop a model for a multiphase RHC process able to handle the vast amount of reactive components and reactions and to perform the VLE and VLLE calculations efficiently. The solution is to use a lumped model for treating the physical phenomena and a detailed component model for modeling the chemical reactions. The essential part of the model is that the complete

about 17.1% aromatic carbons of all carbons and its sulfur content was 1.01%. The temperatures were defined as outlet temperatures of the reactors: T1 = 673.15 K and T2 = 698.15 K. Figure 4 presents the true boiling point distillation curves for total feed and total reactor effluent for all cases. The graphs cover the whole feed and product streams including light

Table 2. Simulated Conversions of Classes of Carbon Atoms in Different Cases, Reduction of Components Having NBPs Higher than 783.15 K and Formation of Sediments Aromatic carbons in feeda 32.4 32.4 17.1 17.1

T (K)

Secondary aromatic carbonsb

External tertiary aromatic carbonsb

Internal tertiary aromatic carbonsb

Naphthenic carbonsb

Open chain aliphatic carbonsb

783.15 K+ fraction (wt %)

Sediment per residue feed (ppm)

673.15 698.15 673.15 698.15

8.3 4.5 17.4 14.2

19.4 19.7 21.1 22.8

7.5 5.4 8.5 8.3

20.7 37.8 21.8 39.0

−26.1 −37.2 −20.3 −30.5

42.7 52.8 51.8 59.1

273 535 43 90

a

Weight percent. bAll conversions of carbon classes are percent of the amount present in the feed. G

DOI: 10.1021/acs.iecr.8b01834 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Industrial & Engineering Chemistry Research

■

composition of each stream is stored and lump properties are upgraded accordingly for each stream individually when necessary. The motivation of creating such a model is to take advantage of the advanced analysis methods such as FT-ICR mass spectrometry, NMR, and HPLC. Performing phase equilibrium calculations with thousands of components is very tedious. Use of lumped components for phase equilibrium calculations solves this problem. The advantage of the complex reaction model is its predictive potential. Use of actual chemical components as the basis of the model makes it possible to solve heat and material balances and to predict the physical properties using thermodynamic methods. The model gives detailed compositions and properties of products. Moreover, their generation in the process as well as consumption of hydrogen can be traced. The methodology presented here was successfully applied to a system involving vapor−liquid−solid catalyst mixture with almost 20 000 components and 490 000 reactions. As a demonstration, results for two hypothetical feed compositions at two temperature levels are presented. For the VLE calculations for RHC reactors the base components were lumped to 42 pseudocomponents. With a standard laptop computer the heat and material balances of a reactor train with three reactors can be solved in a few minutes. For VLLE calculations for the separation train following the RCH reactors, the reactor effluent was divided into 162 lumped components based on their NBPs and solubility factors. The simulated RHC conversions and mesophase amounts show distinct results for each feed type and reaction temperature. The modeling methodology can be applied to many systems comprising large number of reactive components and involving phase equilibrium calculations such as diesel hydrotreating or upgrading of biomaterials such as lignocellulosic feeds, which are highly complex mixtures.

■

Article

REFERENCES

(1) Schuler, B.; Meyer, G.; Peña, D.; Mullins, O. C.; Gross, L. Unraveling the Molecular Structures of Asphaltenes by Atomic Force Microscopy. J. Am. Chem. Soc. 2015, 137 (31), 9870−9876. (2) Ancheyta, J.; Sánchez, S.; Rodríguez, M. A. Kinetic modeling of heavy oil fractions: A review. Catal. Today 2005, 109, 76−92. (3) Pereira de Oliveira, L. C. Développement d’une méthodologie de modélisation cinétique de procédés de raffinage traitant des charges lourdes. Ph.D. Dissertation, Ecole normale supérieure de lyonENS LYON, 2013. (4) Sanchez, S.; Rodriguez, M. A.; Ancheyta, J. Kinetic Model for moderate hydrocracking of heavy oils. Ind. Eng. Chem. Res. 2005, 44, 9409−9413. (5) Benito, A. M.; Callejas, M. A.; Martinez, M. T. Kinetics of asphaltene hydroconversion 2. Catalytic hydrocracking of a coal residue. Fuel 1997, 76, 907−911. (6) Callejas, M. A.; Martinez, M. T. Hydrocracking of a Maya residue. Kinetics and product yield distributions. Ind. Eng. Chem. Res. 1999, 38, 3285−3289. (7) Wiehe, I. Process Chemistry of Petroleum Macromolecules; CRC Press: 2008; 456pp. (8) Neurock, M.; Libanati, C.; Nigam, A.; Klein, M. T. Molecular structure and reactivity in modeling heavy oils. Chem. Eng. Sci. 1990, 45 (8), 2083−2088. (9) Pereira de Oliveira, L. P.; Verstraete, J. J.; Kolb, M. Simulating vacuum residue hydroconversion by means of Monte-Carlo techniques. Catal. Today 2014, 220−222, 208−220. (10) Pereira de Oliveira, L.; Verstraete, J. J.; Kolb, M. Development of a General Modeling Methodology for Vacuum Residue Hydroconversion. Oil Gas Sci. Technol. 2013, 68 (6), 1027−1038. (11) Korre, S.; Klein, M.; Quann, R. J. Polynuclear Aromatic Hydrocarbons Hydrogenation. 1. Experimental Reaction Pathways and Kinetics. Ind. Eng. Chem. Res. 1995, 34, 101−117. (12) López-García, C.; Hudebine, D.; Schweitzer, J.-M.; Verstraete, J. J.; Ferré, D. In-depth modeling of gas oil hydrotreating: From feedstock reconstruction to reactor stability analysis. Catal. Today 2010, 150 (3−4), 279−299. (13) Quann, R. J.; Jaffe, S. B. Structure-oriented lumping: describing the chemistry of complex hydrocarbon mixtures. Ind. Eng. Chem. Res. 1992, 31 (11), 2483−2497. (14) Quann, R. J.; Jaffe, S. B. Building useful models of complex reaction systems in petroleum refining. Chem. Eng. Sci. 1996, 51 (10), 1615−1635. (15) Kunnas, J.; Ovaskainen, O.; Respini, M. Mitigate fouling in ebullated bed hydrocrackers. Hydrocarbon Process. 2010, 89 (10), 59− 64. (16) Gray, M. R.; McCaffrey, W. C. Role of Chain Reactions and Olefin Formation in Cracking, Hydroconversion, and Coking of Petroleum and Bitumen Fractions. Energy Fuels 2002, 16, 756−766. (17) Bi, W.; McCaffrey, W. C.; Gray, M. R. Agglomeration and Deposition of Coke during Cracking of Petroleum Vacuum Residue. Energy Fuels 2007, 21, 1205−1211. (18) Dickie, J. P.; Yen, T. F. Macrostructures of asphaltic fractions by various instrumental methods. Anal. Chem. 1967, 39, 1847−1852. (19) Mullins, O. C.; Sheu, E. Y.; Hammami, A.; Marshall, A. G. Asphaltenes, Heavy Oils, and Petroleomics; Springer: 2007, 670pp. (20) Mullins, O. C. The Modified Yen Model. Energy Fuels 2010, 24, 2179−2207. (21) Mullins, O. C.; Sabbah, H.; Eyssautier, J.; Pomerantz, A. E.; Barré, L.; Zare, R. N.; et al. Advances in Asphaltene Science and the Yen−Mullins Model. Energy Fuels 2012, 26, 3986−4003. (22) Sullivan, R. F.; Boduszynski, M. M.; Fetzer, J. C. Molecular transformations in hydrotreating and hydrocracking. Energy Fuels 1989, 3, 603−612. (23) Seki, H.; Kumata, F. Structural change of petroleum asphaltenes and resins by hydrodemetallization. Energy Fuels 2000, 14, 980−985. (24) Marrero, J.; Gani, R. Group-contribution based estimation of pure component properties. Fluid Phase Equilib. 2001, 183−184, 183−208.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Author Contributions

The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes

The authors declare no competing financial interest.

■ ■

ACKNOWLEDGMENTS The authors thank Neste Oyj and Neste Engineering Solutions Oy for permit to publish this article. ABBREVIATIONS FT-ICR-MS = Fourier transform ion cyclotron resonance mass spectrometry NBP = normal boiling point NMR = nuclear magnetic resonance HPLC = high performance liquid chromatography PAH = polyaromatic hydrocarbon PC-SAFT = perturbed chain statistical associating fluid theory RHC = residue hydrocracking VLE = vapor−liquid equilibrium VLLE = vapor−liquid−liquid equilibrium H

DOI: 10.1021/acs.iecr.8b01834 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX

Article

Industrial & Engineering Chemistry Research (25) Hukkerikar, A. S. Development of pure component property models for chemical product-process design and analysis. Ph.D. Dissertation, DTU, Lyngby, Denmark, 2013. (26) Joback, K. G.; Reid, R. C. Estimation of pure component properties from group-contributions. Chem. Eng. Commun. 1987, 57, 233−243. (27) Gonzalez Rodriguez, D. L. Modeling of asphaltene precipitation and deposition tendency using the PC-SAFT equation of state. Ph.D. Dissertation, Rice University, Houston, TX, 2008, 238pp. (28) Ancheyta, J.; Speight, J. G. Hydroprocessing of Heavy Oils and Residua; CRC Press: Boca Raton, FL, USA, 2007; 345pp. (29) Nikov, I.; Grandjean, B. P. A.; Carreau, P. J.; Paris, J. Viscosity effects in cocurrent three-phase fluidization. AIChE J. 1990, 36, 1613−1616. (30) Stanislaus, A.; Hauser, A.; Marafi, M. Investigation of the mechanism of sediment formation in residual oil hydrocracking process through characterization of sediment deposits. Catal. Today 2005, 109, 167−177. (31) Gross, J.; Sadowski, G. Application of the perturbed-chain SAFT Equation of State to Associating Systems. Ind. Eng. Chem. Res. 2002, 41, 5510−5515.

I

DOI: 10.1021/acs.iecr.8b01834 Ind. Eng. Chem. Res. XXXX, XXX, XXX−XXX