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

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.


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 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.

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. 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.

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 H 2 . The hydrogen is then compressed and blended with fresh hydrogen (H 2,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.

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. 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 2 H . The hydrogen is then compressed and blended with fresh hydrogen ( 2, makeup H ).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 ( HT Q ) 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.

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.

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) ringopening; (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.

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).
where r R is the reversible reaction rate and ir R is the irreversible reaction rate; whole R is the whole reaction rate; w represents the inherent rate constant ascertained by basic study; 2 H P is the partial pressure of hydrogen; m is the reaction order of hydrogen; eq is the equilibrium constant of the reaction; and AD S indicates the whole LHHW adsorption coefficient, which expresses the competitive adsorption with diverse restrainers.

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 reaction network of nitrogen lumps adopted in the HCR model [24]. According to the LHHW theory, reversible and irreversible reactions can be expressed by Equations (1) and (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; ADS i and ADS 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.

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.
Processes 2020, 8, x FOR PEER REVIEW 6 of 20 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: where k denotes the reaction rate constant; i C and j C are the molar concentrations of lump i and j , respectively; , m i and , m j denote the th m 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: where i y and j y indicate the mass concentrations of lump i and j , respectively; l k is the rate constant of the th l 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:  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: 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: 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:

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.

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: 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.

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, 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. 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.

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, 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  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]: where ∆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 CI T,b represents a boiling point corrector, which can be computed as follows: 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. 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.

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.

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.  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.

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 −1 . 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. It can be observed that the adjusted value is almost 0.66-0.76 times its original value.

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 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.

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).

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 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.

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. Figures 10 and 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.

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. Figures 10 and 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.  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  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 (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.

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 2,makeup H 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.

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 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 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.  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.   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.   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.     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.   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.   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.

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 databased 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 2,makeup H flow rate with the ARD of 2.27%, the temperature distribution of the column with Pearson coefficients over

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 H 2,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%.

Conflicts of Interest:
The authors declare no conflict of interest.