Modeling and Simulation of Reaction and Fractionation Systems for the Industrial Residue Hydrotreating Process

Abstract

The residue hydrotreating process plays a significant role in the petroleum refining industry. In this process, modeling and simulation have critical importance for process development, control, and optimization. However, there is a lack of relevant reports of plant scale due to complexity in characterizing feedstock and determining reaction mechanisms. In this paper, reaction and fractionation models are constructed and simulated for a real-life industrial residue hydrotreating process based on Aspen HYSYS/Refining. Considering the heavier and inferior residue, analytical characterization is carried out for feedstock characterization based on laboratory analysis data. Moreover, two reactor models with parallel structures are proposed to implement the intricate reaction network, namely, a hydrocracker reactor and a plug flow reactor. The former simulates lighter petroleum hydrotreating based on the built-in reaction network. The latter emulates the conversion of a peculiar, heavier resin and asphaltene, using a six-lump model, which expands the scope of the feedstock and improves the accuracy of the model. To obtain a realistic simulation of fractionation, the database-based delumping method is adopted to model it with proper pseudo-components. The simulation results, including temperature rise, hydrogen consumption, temperature distribution, product yield, product properties, indicate that the model is capable of reflecting the realistic process accurately.

1. Introduction

The global energy crisis has become an urgent issue in the 21st century [1]. The worldwide growth in energy consumption and carbon emissions is about 1.5% and 1.4% annually over the past 5 years [2]. However, as the main form of energy, the remaining oil is forecasted to be available for only about 50 years [2]. Furthermore, the oil refineries face various challenges: growing demand for high-quality and low-boiling point products, deterioration of crude oil, ever-increasing difficulty of petroleum processing, and stringent environmental requirements. Since vacuum residue accounts for about half of the crude oil, more attention has been focused on vacuum residue hydroconversion processes, which are capable of converting heavy, inferior petroleum into light, valuable products such as gasoline, jet fuel, and diesel [3,4]. A good process model is important for monitoring the plant-level operation status, optimizing operating conditions, verifying the dynamic adjustment scheme, and maximizing the profit [5,6,7].

The residue hydrotreating process mainly consists of a reactor and a fractionation unit. The core of the process model is the kinetic model, which describes the catalytic cracking reactions based on mechanism knowledge. For the kinetic model, the inputs are the feedstock properties and reaction conditions. The model outputs are product yield, product properties, reactor temperature rise, and hydrogen consumption at different operation conditions. In the past years, two main kinds of methods have been studied to establish the kinetic model, which are based on molecular composition and lumped composition, respectively.

Based on molecular composition for the kinetic model, structure-oriented lumping (SOL) proposed by Quann et al. [8] is one of the most widely used techniques. In this technique, the feedstock mixture was expressed as a matrix, whose rows indicated different atomic elements like C, H, S, N, O, and columns represented different structural groups. This technique was then applied to model the fluid catalytic cracking (FCC) of the naphtha hydrodesulfurization [9]. Additionally, a two-step reconstruction algorithm was proposed by Hudebine et al. [10]. This algorithm could rebuild petroleum at a molecular level from overall petroleum analyses. Moreover, Peng et al. [11] investigated another approach called as the molecular type homologous series (MTHS), in which homologous series, along with carbon number information were considered for kinetic modeling. This approach obtained excellent predictions for a coal tar hydrogenation process [12]. However, these methods are with a shortcoming of the computation expensiveness since extensive molecular compositions are generated in the kinetic model. Additionally, another shortcoming is that routine chemical analysis and ¹³C nuclear magnetic resonance (NMR) are required. However, these analyses are mostly performed with low frequency, which limits these molecular-level methods on applying in residue hydroconversion process simulations.

The lumped method, proposed by Weekman et al. [13], was originally utilized for a three lumped model of a catalytic cracking process. In this method, the mixture of hydrocarbons was characterized by several lumped compositions classified by their reaction characteristics, which can alleviate the difficulty of building the complicated kinetic model. Then, the kinetic model is further extended to four [14], five [15], seven [16], and even thirty-seven lumped models [17]. However, those models only reflect the productivity, whereas some detailed information, such as sulfur and nitrogen content, is not included. Besides, once the lumps of kinetic model are determined, the cutting scopes of gasoline and distillate oil are permanent. However, the cutting scopes generally vary from refinery to refinery with different optimal economic benefits, which results in poor adaptability of the simulation models.

Hence, it is necessary to construct a reasonable model based on proper material compositions. The compositions should contain theoretical, circumstantial information sufficient to reflect the process. Meanwhile, proper number of compositions are required to easily make the calculation and give good extrapolation results. Fortunately, another approach is utilized in the Aspen HYSYS/Refining hydrocracker model to simulate the reaction process at the molecular level, in which 97 lumps are selected to characterize petroleum. Also, these lumps take both the carbon number and molecular structure information into consideration, which are in accordance with the results in the reported literatures [18,19]. These lumps can be grouped into three classes of paraffin (P), naphthene (N), and aromatics (A). They can be tested using mass spectrometry and adopted to describe light petroleum (i.e., the PNA method). However, the heavy vacuum residue is usually characterized by the method four components of saturate, aromatic, resin, and asphaltene (SARA). The PNA characterization method is different from SARA, which may be improper for simulating the vacuum residue hydrogenation process. In addition, high-boiling point resin and asphaltene are not represented by the lumps in the hydrocracker [20,21], which results in the lower high-boiling points of simulated petroleum. Moreover, the colloidal structure of resin, and asphaltene (RA) have a great impact on the overall kinetics [22].

To alleviate the aforementioned problems, two parallel structure reactors model is proposed in this paper to simulate a real industrial residue hydrogenation process. The conversion process of the heavier petroleum regarded as an asphaltene lump is modeled in a plug flow reactor (PFR), which is in parallel with the original hydrocracker reactor (HCR) existing in Aspen HYSYS/Refining. In the PFR, the asphaltene conversion reaction network is established and its lumps are necessary to be characterized. To describe the asphaltene conversion, a six-lump reaction model is adopted. For the lump characterization, the property of asphaltene is calculated using the group contribution method, and the remaining lumps in the reaction model are represented utilizing the substitute mixtures of real components (SMRCs) method [23]. And the HCR simulates lighter petroleum hydrogenation based on the built-in reaction network and lumps.

As the outermost layer of residue hydrogenation model, the process model seeks to integrate the reactor model with the fractionation model. Due to the complexity of residues, most previous works paid attention to only one aspect of the reactor model or the fractionator model for the for the plant wide residue hydrogenation process [24]. To make contribution to this aspect, a residue hydroconversion process model is comprehensively developed by exploiting delumping, which is a committed step to concatenate the reactor with the fractionation model.

The paper is structured as follows: Section 2 describes the residue hydrogenation process in detail. Section 3 illustrates the kinetic model describing the hydrogenation process of built-in residue oil simulated in an HCR and asphaltene simulated in a PFR. Section 4 demonstrates the framework for building the whole process model, which includes characterization of the feedstock mixture, establishment of a reactor model of two parallel structure reactors, and delumping of the reactor model effluent. Section 5 presents the simulation results to clarify the model effectiveness. Section 6 provides the conclusion.

2. Process Description

The residue hydrogenation technology is one of the most effective ways to deeply hydro treat the heavy oil. The process is implemented in fixed beds for removal of the majority of impurities from oils, which can supply excellent raw materials to the subsequent catalytic cracking process. In the fixed bed reactors, both hydrotreating and hydrocracking reactions take place in appropriate temperature, pressure, hydrogen–oil ratio, and liquid hourly space velocity (LHSV) conditions. Generally, four reactors are adopted and connected in series, and the main function of each reactor is to remove some specific impurities.

A schematic diagram of the residue hydrogenation unit is displayed in Figure 1. The raw materials are blended with vacuum heavy wax oil, coke gas oil, tar wax oil and catalytic circulating oil. Before entering the reaction system, the blended oil is filtered in a buffer tank. Then the oil is boosted in the pump and exchanged heat with the reactor effluent. Afterwards, the oil goes through the heating furnace before entering the first reactor. Thus, the first reactor temperature is mainly affected by the heating furnace.

The temperatures of the rest ones are controlled by hydrogen quenching to attain the required operation conditions. In the hydrotreater (HT), sulfur, nitrogen, and oxygen compounds are mostly converted to hydrogen sulfide, ammonia, and water, while polycyclic aromatic hydrocarbons and some unsaturated hydrocarbons are hydro saturated.

The outflow from the last HT is then transferred to a high and low pressure separation system (HLPS) in which three-phase separation of the gas, oil and water takes place. The HLPS mainly consists of four separators, as illustrated in Figure 1. The gas of the separators then enters a desulfurization tower to desorb hydrogen sulfide. After this operation, the hydrogen purity of recycled gas is improved. Meanwhile, the off-gas discharges a part of the ${\mathrm{H}}_2$. The hydrogen is then compressed and blended with fresh hydrogen (${\mathrm{H}}_{2,\mathrm{makeup}}$).The blended gas recycles back to the four reactors, in which a part of the recycled gas (RG) is mixed with feed oil, and the other gas functioned as gas quenching ($Q_{HT}$) quenches the beds to ensure security and protect the catalyst. And the bottom stream of the HLPS is delivered to the downstream: the fractionation process, which includes a fractionator, a diesel stripper, a condenser, and a pump. The top product of the fractionator is naphtha, the middle one is diesel, and the bottom one is tail oil.

3. Kinetic Model of the Residue Hydrotreating Unit

In this paper, the reaction process is modeled with a two-parallel reactors model: an HCR and a user-defined PFR. For the HCR, the lighter petroleum hydrogenation process is simulated using the built-in residue oil hydrogenation reaction network with 97 lumps. For the PFR, the heavy asphaltene hydrogenation process is modeled with the six-lump kinetic model. Figure 2 macroscopically illustrates the reaction mechanism utilized in this paper.

3.1. Kinetic Model of the Built-In Residue Oil Hydrotreating Process

In the HCR model, the reaction network is extremely complicated with 177 reactions. According to reference [24], the reactions can be grouped into six types: (1) paraffin hydrocracking; (2) ring-opening; (3) dealkylation of aromatics, naphthenes, nitrogen lumps, and sulfur lumps; (4) saturation of aromatics, non-basic nitrogen lumps, and hindered sulfur lumps; (5) hydrodesulfurization (HDS) of unhindered sulfur lumps; and (6) hydrodenitrogenation (HDN) of nitrogen lumps. Each reaction rate equation is based on the classical theory of Langmuir–Hinshelwood–Hougen–Watson (LHHW) [25], which demonstrates that the heterogeneous catalytic reactions undergo three basic stages: adsorption, surface reaction, and desorption. Figure 3 illustrates the part of the reaction network employed in the HCR model of Aspen HYSYS for the reaction of nitrogen lumps.

According to the LHHW theory, reversible and irreversible reactions can be expressed by Equations (1) and (2).

$$R_r=R_{whole}w\frac{(AD{S}_i{C}_{i}AD{S}_{H_2}{(P_{H_2})}^m/eq-AD{S}_j{C}_j)}{ADS}$$

(1)

$$R_{ir}=R_{whole}w\frac{AD{S}_i{C}_{i}AD{S}_{H_2}{(P_{H_2})}^m}{ADS}$$

(2)

where $R_r$ is the reversible reaction rate and $R_{ir}$ is the irreversible reaction rate; $R_{whole}$ is the whole reaction rate; $w$ represents the inherent rate constant ascertained by basic study; $AD{S}_i$ and $AD{S}_j$ indicate the adsorption coefficients of lumps $i$ and $j$, respectively; lump $i$ and hydrogen $H_2$ are reactants; lump $j$ is the resultant; $C_i$ and $C_j$ are the concentration of lumps $i$ and $j$, respectively; $P_{H_2}$ is the partial pressure of hydrogen; $m$ is the reaction order of hydrogen; $eq$ is the equilibrium constant of the reaction; and $ADS$ indicates the whole LHHW adsorption coefficient, which expresses the competitive adsorption with diverse restrainers.

3.2. Kinetic Model of the Asphaltene Hydrogenation Process

The hydrogenation of resin and asphaltene with the high boiling point, and density also involves involute and heterogeneous catalytic reactions. Hence, by virtue of the complexity of high boiling point petroleum and its main products, a simple lump strategy is employed [26]. For simplification, resin and asphaltene are together regarded as the asphaltene lump. In the hydrogenation, three phase products are decomposed from asphaltene. The liquid mixtures are grouped into three lumps: the light oil (initial boiling point (IBP) < 350 °C) lump, the middle oil (350−540 °C) lump, and the heavy oil (>540 °C) lump. Gaseous gas and solid coke were taken as the gas lump and coke lump, respectively. Also, it is notable that the active hydrogen reaction is not considered in this paper since the tetralin (its donor), as shown in [26], does not exist in the real refinery production process. Figure 4 illustrates the reaction network of the asphaltene conversion process. The components substituting for lumps of the network are placed in the right and described in detail in Section 4.1.

Since the constants of LHHW in Equations (1) and (2) ascertained by basic research are scarcely found, general kinetic reactions are applied to describe the reaction network of petroleum based on the Arrhenius equation. And researchers have different opinions about the order of the reactions [27,28]. In this paper, it is presumed that the whole reaction pathway consists of first-order reactions. In addition, due to the low gas yield of the entire reaction system, the volume of the reaction system is assumed to be immutable.

The reaction rate $R$ can also be expressed in the form of the power law model:

$$R=-k{C}_i^{m,i}{C}_j^{m,j}$$

(3)

where $k$ denotes the reaction rate constant; $C_i$ and $C_j$ are the molar concentrations of lump $i$ and $j$, respectively; $m,i$ and $m,j$ denote the $m_{th}$ reaction orders of lump $i$ and $j$, respectively. In this paper, all constants for $m,i$ are set to 1.

In fact, the data obtained from the factory is based on mass flow rate. It is difficult to acquire mole fractions of the lumps. Instead, the mass concentration is adopted to calculate the kinetic parameters. Accordingly, the rate equation can be rewritten as:

$$R_l=-k_l{y}_i^{m,i}{y}_j^{m,j}$$

(4)

where $y_i$ and $y_j$ indicate the mass concentrations of lump $i$ and $j$, respectively; $k_l$ is the rate constant of the $l_{th}$ reaction for lump $i$ and $j$, Thus, based on the reaction network of asphaltene conversion, the mathematical equations of the kinetic model can be rewritten in the following forms:

$$\frac{d{y}_1}{dt}=-(k_1{y}_1+k_2{y}_1+k_3{y}_1+k_4{y}_1)$$

(5)

$$\frac{d{y}_2}{dt}=k_1{y}_1-(k_5{y}_2+k_6{y}_2+k_7{y}_2)$$

(6)

$$\frac{d{y}_3}{dt}=k_2{y}_1+k_5{y}_2-k_8{y}_3$$

(7)

$$\frac{d{y}_4}{dt}=k_6{y}_2+k_8{y}_3$$

(8)

$$\frac{d{y}_5}{dt}=k_3{y}_1+k_7{y}_2$$

(9)

$$\frac{d{y}_6}{dt}=k_4{y}_1$$

(10)

4. Workflow of Developing an Integrated Residue Hydrogenation Process Model

The workflow for developing an integrated residue hydrogenation model is shown in Figure 5 with the Aspen HYSYS/Refining simulation software. A similar workflow for the hydrocracking process was proposed in [24]. However, it is different for the heavy petroleum refining process.

First of all, residue tends to become inferior and heavier since the crude oil gets heavier. However, the feed library for residue is predetermined in the simulation software. Furthermore, the asphaltene with higher boiling point is not considered in the feed library. In addition, owing to the absence of asphaltene, its product lumps are lacking in the feed library. Thus, it is required to make some adjustments for the feed in the HCR library and characterize the added six lumps. Moreover, due to the different feed properties and reaction conditions, it is necessary to determine the kinetic parameters of the reactor model. In the HCR model, the reaction activity variables are to be estimated for the built-in residue oil hydrogenation. Correspondingly, in the PFR model, the kinetic parameters are determined for the additional asphaltene hydrogenation. Finally, the reactor model is established with kinetic lumps grouped by reaction characteristics, whereas the fractionation model is built based on the pseudo-components according to fractionation characteristic (i.e., the boiling point). Therefore, the adopted kinetic lumps of reactor model may be inappropriate for a fractionator modeling of residue hydrogenation process. Consequently, it is prerequisite to transform kinetic lumps into suitable pseudo-components through a process called delumping. The detailed workflow of the residue hydrogenation model can be described as follows:

(1)

Obtain data from factories and select the required data for modeling. In this work, the set data (e.g., the feed property, flow) is the input of the model and the target data (e.g., the product yield, property) are the object value to be attained.

(2)

Characterize the feedstock based on laboratory testing data.

(3)

Develop the reactor model according to the real process data.

(4)

Delump the effluent from the reactor to build the fractionator.

(5)

Test the effectiveness of the model by comparing the model results with actual plant data.

4.1. Feedstock Mixture Characterization

Feed quality has a significant influence on residue hydrogenation process. Thus, it is fundamental to characterize the feedstock mixture according to its bulk properties. In this work, the “residue” fingerprint type is selected in the basic feed library of HCR for the residue hydrogenation process. And the feed data section inputs are the experimental bulk properties of feed, such as its density, boiling point, refractive index, sulfur content, and nitrogen content. The corresponding outputs are the simulated bulk properties and composition contents. Generally, the experimental bulk properties diverge from the simulated. Figure 6a show the residue properties of a dataset in the HCR. It can be found that the simulated initial boiling point (IBP) and final boiling point (FBP) are lower than their experimental values. Meanwhile, the simulated total light oil percentage is slightly higher in

Table 1. Simulated residue feed composition contents with basic and improved feed library in HCR.

	Naphtha	Distillate	Gas Oil	Total Light Oil
Basic library	0%	1.72%	12.04%	13.76%
Improved library	0%	2.79%	8.72%	11.01%

, which may have impact on the simulated reactor performance. It is because that the conversion of actual weight residue hydrogenation is only 15–20%, which is defined as the sum of the weight proportions of gas and light oil after reaction. Therefore, two problems demand to be solved.

The first problem is that the lumped compositions whose normal boiling point is lower the real IBP value should be removed. Therefore, their proportion should be set to zero, which is based on the separation efficiency factor proposed by García et al. [29]. In this way, the proportions of remaining lumps inevitably increase. To handle this problem, it is necessary to adjust the proportions, which has been presented in the previous work [30]. Figure 6b and Table 1 show the results of the experimental value, simulated value with basic and improved library of the residue feed. Since the distillation range of naphtha is lower than that of residue, the contribution of naphtha is 0% in Table 1.

The other problem is that the simulated FBP of the feed is lower than its experimental value. Therefore, it is indispensable to add the composition representing the higher boiling point petroleum lump, and other lumps in the reaction network, which is displayed in Figure 4. To represent the added petroleum, the composition of Tahe residue heavy resin is adopted from [31], since its properties are closer to the real petroleum cut. Its properties of molecular weight and density are also given in reference [31]. Although the software has a powerful ability to calculate physical property, at least three properties are required to characterize the lumps more precisely. Fortunately, the normal boiling point of petroleum cut can be estimated using the group contribution method [32]. In this method, the composition (considered as a molecule here) is composed of 77 structural groups. Each group has its own contribution and can be referred to in reference [33]. The molecule property is calculated by summing the products over each group numbers and its corresponding contribution. Meanwhile, the method of first-order group contributions is applied to predict the normal boiling point, which is shown in the following equation [34]:

$$\exp (\frac{T_b}{T_{b,0}})={\displaystyle \sum _i(n_i\mathrm{\Delta }{T}_{b,i})+C{I}_{T,b}}$$

(11)

where $\mathrm{\Delta }{T}_{b,i}$ represents the group contribution of group $i$; $n_i$ represents the number of group $i$ in the molecule; $T_{b,0}$ represents a regulatable parameter ($T_{b,0}=307.63K$), and $C{I}_{T,b}$ represents a boiling point corrector, which can be computed as follows:

$$C{I}_{T,b}=0.2285N_R^2+0.4687N_R$$

(12)

where $N_R$ represents the number of rings in a molecule.

According to Equations (11) and (12), the normal boiling point of asphaltene is estimated, and other properties are calculated by the software. The results are shown in

Table 2. Lumps and their properties utilized in plug flow reactor (PFR).

Lump	Normal Boiling Point (°C)	Density (kg·m−3)	MW (g·mol−1)	Critical Temperatures (°C)	Critical Pressures (kpa)	Critical Volumes (m3·kg−1·mol−1)
ASP	889	1052	1057	989.6	378.7	2.78
C47	563	825	661	643	375	3.1217
n-C30	449.719	811.898	422.780	589.85	868	1.72353
n-C14	253.508	762.913	198.380	420.85	1620.18	0.82999
C2	−88.600	355.683	30.070	32.27801	4883.85	0.148
C	/	1642.06	/	/	/	/
ASP, asphaltene; C47: C47paraffins; n-C30: n-Triacontane; n-C14: 1,2-Tetradecanediol; C2: Ethane; C, Coke.

.

Then, the remaining five lumps need to be characterized in the reaction network. Here, the method of substitute mixtures of real components (SMRCs) is used for lump characterization. This method enhances comprehension of the mechanism of the reaction since the structure of a substituted component is included. In this method, each lump is substituted by one component whose normal boiling point approximates the average distillation range of the lump. Since multiple components may be available for a specified boiling point, it is important to select the proper components. The main principles are as follows: (1) The components exist in the system and can be detected by analysis; (2) The low-boiling-point components, the low-carbon components and their isomers can be tested by analytical techniques. (3) For the high-carbon components, it is preferable to choose paraffin and aromatics. All the selected lumps and their properties are shown in Table 2.

4.2. Reactor Model Construction with Two Parallel Structures

After feedstock mixture characterization, the next step is to develop reactor model. First, the process data (e.g., catalyst loading, feed rate, feedstock analysis, reactor inlet temperature, and reactor pressure) is adopted to synchronously identify the parameters of the rate equations and reactor design equations. Then, the deviation of model prediction from plant data is minimized by adjusting the reaction activity variable in the HCR model and the user-defined PFR. For the reactor model, modifying the activity factor is essential, since the feed property, reactor configuration, catalyst activity, and operating conditions differ greatly in various refineries. The procedure of adjusting the activity factor to lessen the disparity is called “calibration”. For the user-defined PFR, the kinetic parameters are identified similarly. The two parallel structure reactor model is demonstrated in Figure 7.

For the HCR model, the calibration procedure is already described in [24]. Here, there are some additional steps in the calibration procedure for the residue hydrogenation process. In the tuning step, it is discovered that the reaction intensity of the catalyst bed competes against each other as shown in Figure 8. The reaction intensity is reflected by the temperature rise since the reactions are overall exothermal. That is, when other parameters are unaltered, the global activity factor of the first catalyst bed will increase, which might result in a decline of temperature rise for the remaining catalyst bed, as demonstrated in Figure 8. Furthermore, it is recommended that the simulated temperature rises are overall lower than the real temperature rise. It may be because the exothermal hydrodemetallization reactions are excluded in the HCR model, which occur in practical production.

Since a PFR substitutes for four reactors of the real process, the arithmetic averages of operating parameters are adopted, such as temperature, pressure, volume, and so on. For the identification of kinetic parameters, reference [26] adopted genetic algorithm method for minimizing the deviation of the prediction of seven-lump kinetic model from the experimental data. In the established kinetic model, the activation energies were from 106.07 to 237.5 kJ. mol ⁻¹. However, due to the deviations of the feed, operation condition between the reference and the actual process, adjustments of kinetic parameters are required to match the real process. In addition, it is impossible to directly determine the kinetic parameters for the conversion of asphaltene because the petroleum is part of the inlet residue oil in a real industrial process. To handle this problem, the consistent rule and the whole asphaltene conversion are employed. And considering the existence of a catalyst in the real process, only the activation energy is changed. The consistent rule is that the adjusted distribution of reaction activation energy should be similar to that published in [26]. The kinetic parameters applied in this work are shown in

Table 3. Apparent activation energies with their frequency factors of the six-lump model in [26] and this paper. Reproduced with permission from copyright American Chemical Society, [2017].

			From [26]		This Paper
Reaction	Frequency Factor (h−1)		Activation Energy (kJ·mol−1)		Activation Energy (kJ·mol−1)
							Asp → HO	A1	6.70 × 108	E1	106.07	E1	82.2
Asp → MO	A2	7.78 × 108	E2	109.06	E2	83.05
Asp → G	A3	8.13 × 109	E3	175.56	E3	126
Asp → C	A4	5.06 × 109	E4	169.17	E4	125
HO → MO	A5	9.92 × 108	E5	130.78	E5	86.5
HO → MO	A6	9.92 × 109	E6	155.43	E6	112
HO → G	A7	9.91 × 108	E7	237.50	E7	156.75
MO → LO	A8	2.24 × 1011	E8	153.63	E8	110
Asp, asphaltene; HO, C47; MO, N-C30; LO, N-C14; G, C2; C, coke

. It can be observed that the adjusted value is almost 0.66–0.76 times its original value.

4.3. Delumping of the Reactor Model Effluent and Fractionator Model Development

Delumping the outflow of the reactor model works as a fundamental step combining the reactor model to the fractionator model. Firstly, lumps in the reactor model are more concerned with the carbon number and structure, whereas the lumps in the fractionator model pay more attention to the accuracy of thermodynamic behavior. Secondly, the lumps used in the kinetic model belong to two different reactor models and are not integrated in this paper.

Numerous researchers have contributed to the development of this aspect. Haynes et al. [35] applied the Gauss–Legendre quadrature to calculate the vapor–liquid equilibrium (VLE) of the petroleum whose properties are expressed by critical properties and acentric factors. Yang et al. [36] determined appropriate fractionation lumps and their physical properties by interpolation, cutting, and normalization of the distillation curve of reactor effluent for the residue hydrogenation process. Chang et al. [24] utilized the Gauss–Legendre quadrature to delump the overflow of the reactor model into 20 pseudo-components for the medium-pressure hydrocracking (MP HPR) and high-pressure hydrocracking (HP HPR) units.

However, a common disadvantage in these methods is that some of the related correlations are controversial. Additionally, the residue feed covers a wider range of distillation and contains more complex compositions, making the delumping tougher. Moreover, for the high and heavy fractions, property correlation may generate questionable results because the formulas are usually summarized by low boiling-point fractions. In this paper, a data-based method is applied to delump the reactor model outflow and develop the fractionator model. The component list of the method is named “Assay Components Celsius to 850 C” in ASPEN/HYSYS software. The software has collected and checked the data from the world’s most respected sources [37]. Thus, its data source has high credibility.

Due to page limitations,

Table 4. Part of boiling point-based pseudo-components and their properties.

Pseudo-Components	TBP (°C)	ILD (kg·m−3)	MW (g·mol−1)	Pseudo-Components	TBP (°C)	ILD (kg·m−3)	MW (g·mol−1)
340–350 C *	350	877	299	520–540 C *	540	958	538
350–360 C *	360	882	309	540–560 C *	560	967	565
360–370 C *	370	887	320	560–580 C *	580	975	593
370–380 C *	380	892	330	580–600 C *	600	984	622
380–390 C *	390	894	341	600–625 C *	625	990	668
390–400 C *	400	898	353	625–650 C *	650	997	718
400–410 C *	410	903	367	650–675 C *	675	1006	770
410–420 C *	420	908	381	675–700 C *	700	1015	828
420–430 C *	430	912	394	700–725 C *	725	1019	911
430–440 C *	440	917	408	725–750 C *	750	1023	1002
440–450 C *	450	921	421	750–775 C *	775	1030	1093
450–460 C *	460	926	434	775–800 C *	800	1036	1190
460–480 C *	480	932	461	800–825 C *	825	1043	1292
480–500 C *	500	941	487	825–850 C *	850	1049	1398
500–520 C *	520	950	512	850+C *	770	1062	998

shows only part of the pseudo-components from 340–850+C * and its property of true boiling point (TBP), ideal liquid density(ILD), molecular weight(MW). More detailed information about the complete components from 36–850+C * can be found in the Supplementary Materials Table S1. For the distillation range from 36 to 460 °C, each temperature interval with 10 °C is represented as a lump; for the distillation boiling point range from 460 to 600 °C, approximately an interval of 20 °C was regarded as a lump; for the range from 600 to 850 °C, each interval with 25 °C is regarded as a lump; the range over 850 °C is considered as a lump. The stream cutter is adopted to connect the two models belonging to different physical property packages. The purpose of the stream cutter is to recalculate stream enthalpy. In the calculation, component A (TBP) in physical property pack 1 is mapped to pseudo-component A* in physical property pack 2 based on the boiling point.

5. Simulation Results

In this section, five operation points are used to demonstrate the efficacy of the established model. Three parts are displayed for the simulation: comparison of a single HCR and the proposed two-reactor model, comparison of the default HCRSRK method in Aspen HYSYS/Refining and the proposed delumping method, and the prediction of the process model for the whole residue hydrotreating conversion. The HCRSRK is a built-in fluid package in the HCR model (with a special HCR.cml component list).

5.1. Comparison of a Single HCR Model with the Proposed Two Reactor Model of HCR and PFR

Five datasets are utilized to verify the effectiveness of the model with two parallel structure reactors. It is easily seen from Figure 9 that compared with a single HCR, the proposed model can provide more precise prediction at high boiling points for tail oil and slightly better results for the lighter product.

Table 5. Comparison of the prediction average relative deviations of boiling point for tail oil with the HCR model and HCR + PFR model.

	Dataset 1	Dataset 2	Dataset 3	Dataset 4	Dataset 5
HCR	12.97%	14.20%	13.50%	14.85%	14.54%
HCR + PFR	5.92%	6.22%	6.06%	7.00%	6.05%

presents the prediction average relative deviation of boiling points for tail oil with the HCR and the proposed model in detail. The proposed reactor model extends the simulated distillation range of the products through the added simulation of asphaltene conversion process. Hence, the distillation curve of the tail oil can be predicted better by the proposed reactor model. However, the predicted initial boiling point of naphtha is lower than the real value in Figure 9. In the simulated naphtha, the additional gas that does not exist in the actual stream contributes to the strange result.

5.2. Comparison of the HCRSRK Method and the Proposed Delumping Method

Sensitivity tests and predictions of properties of the products are performed to demonstrate the rationality and accuracy of results with the HCRSRK and the proposed delumping method. Sensitivity tests are essential to ensure that delumped lumps are capable of producing rational results [38]. A rational result is that the relationships of the side draw rate change to the corresponding distillation curve and temperature are consecutive rather than stepwise. Figure 10 and Figure 11 show the sensitivity test results of the side products naphtha and diesel for the HCRSRK method (a default method in HCR model) and the data-based delumping method. Apparently, the HCRSRK method generates stepwise relationships at most distillation points with the change of draw rates, especially in naphtha, as seen in Figure 10a. The data-based delumping method produces continuous relationships at most distillation points, except the acceptable ASTM D86 10% of naphtha and ASTM D86 70% of diesel. Actually, the ASTM D86 10% of naphtha is easily affected by the contents of its light components. Meanwhile, the actual value of ASTM D86 10% of naphtha is very small. Thus, the relative deviation may seem large even for a small absolute deviation. As can be seen from Figure 10d, the large relative deviation for ASTM D86 70% appears after the change of side draw is large than 6%. Therefore, the stepwise relationship is acceptable in this case. Figure 11 illustrates that the relationship between side draw rate and temperature is smooth with the two methods. In this way, the proposed delumping method is able to predict rational result. The results also show that it is necessary to delump the reactor model.

The experiments of Figure 12 are performed with the feed of fractionator from a HCR since the added ashphaltene lump is out of the range of the HCRSRK. Figure 12 shows the simulation accuracy comparison of products properties, including the distillation curve and standard gravity with these two methods. Figure 12a,d show that the data-based method has better prediction for the distillation range of naphtha, diesel fuel and the standard density of the products. For the distillation curve of tail oil, the two methods have almost the same prediction. The detailed average relative deviations (ARD) of distillation with the HCRSRK and data-based method are 5.59% and 2.56% for naphtha, respectively; they are 6.48% and 3.61% for diesel fuel, respectively; and they are 16.41%, 13.19% for tail oil, respectively. The deficiency of the higher boiling point RA accounts for large ARD of tail oil. Meanwhile, the average relative deviations (ARD) of density of three products are 6.4% and 2.8% with the HCRSRK and data-based method. The lumps with the data-based method, which possesses more number and better distillation properties account for the higher accuracy.

5.3. The Predictions of the Whole Residue Hydrotreating Process Model

The industrial residue hydrotreating process is modeled with a reactor and fractionator model in this paper. The model predictions of four HT reactor temperature rises are displayed in Figure 13. In the reactor model, the catalyst reactor temperature rise is computed with the given inlet temperature. The average absolute deviations (AAD) of each reactor temperature rise is 1.22, 0.95, 2.10, and 1.80 °C for the residue hydrorefining reactor. Thus, it is concluded that the model generates good predictions for the temperature rise of the HT reactor, which plays an important role in evaluating product conversion. The temperature rises of reactor R3 and R4 are evidently higher than R2 and R1. This is because different catalysts are loaded in the reactor. Reactor R3 and R4 mainly removes sulfide and nitrogen while R2 and R1 mainly remove metals. Figure 14 shows the prediction of ${\mathrm{H}}_{2,makeup}$ flow rate for the residue hydrorefining process and the ARD is 2.27%. Obviously, the excellent hydrogen consumption prediction is achieved by the developed model.

To assess the performance of the simulated fractionator, temperature profile, product yields, distillation curves, and density of liquid products are important indicators for model evaluation. Temperature distribution is meaningful for evaluating energy consumption, which is also important to energy optimization and product cutting for engineers. Figure 15 shows the prediction of the temperature distribution of the column. The Pearson coefficients of the temperature distribution are 0.977, 0.985, 0.987, 0.983, 0.985 for dataset 1 to 5, respectively. Apparently, the model demonstrates good capability in predicting column temperature distribution.

Three products are produced in the residue hydrogenation process, namely, naphtha, diesel fuel, and tail oil. Figure 16 illustrates the model predictions on product yields: naphtha, diesel fuel, and tail oil. And

Table 6. Prediction relative deviations of the product yield.

	Naphtha	Diesel Fuel	Tail Oil
Dataset 1	−0.32%	4.79%	−2.11%
Dataset 2	−0.25%	−0.91%	−1.99%
Dataset 3	4.08%	1.52%	0.27%
Dataset 4	−0.01%	−2.52%	0.06%
Dataset 5	−2.62%	−0.07%	−0.57%

shows the results of five datasets in detail. It is observed that the model can accomplish accurate prediction on product yields with consideration of the relative deviation of each product, which can reveal more detailed information.

Figure 17 illustrates model predictions on the distillation curves of major products: naphtha, diesel fuel, and tail oil (in Dataset 1 and 2). Figure 18 shows the density predictions of three major products. The pinpoint predictions indicate the proposed method can well model the distillation curves and density predictions of major products.

6. Conclusions

In this paper, a comprehensive residue hydrotreating process model is established based on plant data with consideration of the product yields, qualities, and process operations. A novel two parallel structure is developed for reactor model. In this model, resin and asphaltene with high boiling point and density are added as a lump, which makes it more consistent with the actual feedstock and production conditions. The proposed reactor model can predict the boiling point of product more accurate. Especially, the average relative deviations of tail oil are only about at almost 6.25%. The accuracy is improved by 55% compared with the original HCR model. For the fractionator, a database-based delumping method is proposed for transforming kinetic lumps of HCR into data-based pseudo-components. The method not only provides a consecutive response to the changes of column specification, but also has a slightly higher prediction accuracy. Five datasets are used to verify the effectiveness of the established integrated model. The developed process model can predict temperature rise of the four catalyst beds with the AAD smaller than 2.10 °C, the ${\mathrm{H}}_{2,\mathrm{makeup}}$ flow rate with the ARD of 2.27%, the temperature distribution of the column with Pearson coefficients over 0.977, the product yield with the ARD smaller than 5%.