OPTIMIZATION OF KINETIC LUMPING MODEL PARAMETERS TO IMPROVE PRODUCTS QUALITY IN THE HYDROCRACKING PROCESS

Abstract

In this study a typical continuous lumping model with five parameters has been used for kinetic modeling of thermal and catalytic hydrocracking. Model parameters have been optimized according to experimental product distributions using a particle swarm optimization (PSO) algorithm. Experimental data from the hydrocracker setup have been employed to validate the proposed model. In this setup hydrogen and vacuum gasoil feed were introduced from the top of a vertical reactor and, after passing through a catalyst bed, the liquid and gas products were separated and analyzed. Temperature of the reactor was adjusted in the range of 440-470°C for thermal hydrocracking, and 410-430°C for catalytic hydrocracking. Liquid hourly space velocities (LHSV) were in the range of 0.5-1.5 feed flow rate per catalyst volume in both sets of experiments. Results of optimization showed that the parameters were only temperature dependent. The comparison between model results and experimental data indicates that the model is capable of predicting product yield with maximum errors of 0.986 and 0.041 for RMSE and AARE values, respectively.

Keywords:
Optimization; Hydrocracking Process; Model Parameters; Continuous Lumping Model; Tubular Reactor System

INTRODUCTION

Hydrocracking is a catalytic conversion process for producing lighter petroleum fractions such as diesel, kerosene, naphtha and gases from high-boiling constituent hydrocarbons in petroleum such as vacuum gasoil (VGO) feedstock, in the presence of hydrogen-rich gas. Normally, the hydrocracking process is carried out in a high pressure trickle bed reactor (Kumar and Balasubramanian, 2009).

Various reviews have been reported on hydrocracking technology in the literature. Different aspects of hydrocracking have been discussed, such as types of hydrocracking, process optimization, catalysis, catalyst acidity, pore diffusion, catalyst poisons and reactor type (Ramakanta et al., 2015; Angeles et al., 2014; Ahmed et al., 2013; Valavarasu et al., 2003; Qader and Hill, 1969; Sadighi and Ahmad, 2013; Matos and Guirardello 2002). These authors also emphasized that the available literature on the kinetics of hydrocracking is scarce, and that kinetic aspects and structural modeling needs to be studied in detail.

In the hydrocracking process, a kinetic model which can predict the performance of the process within the operating severity is very important. Moreover, this model could be extended beyond the experimental domain, for investigating the performance behavior (Moghadassi et al., 2011; Rana et al., 2007; Alhajree et al., 2011, Zhou et al., 2011).

Various kinetic models have been proposed by several researchers (Froment, 2005; Botchway et al., 2004; Browning et al., 2016; Ancheyta et al., 2005; Elizalde et al., 2010; Puron et al., 2014) to study the product distribution in hydrocracking. A commonly used approach is to consider the mixture in terms of selected lumps, which can be specified in terms of structural characteristics, such as boiling point. This approach is attractive for the kinetic modeling of complex mixtures because of its simplicity (Ancheyta et al., 1999). Depending on the mixture characteristic function form; there are two approaches for lumping models, which are continuous and discrete approaches (Canan and Arkun, 2012). In the literature, different lumping models have been reported in order to analyze the kinetics of hydrocracking reactions (Govindhakannan and Riggs, 2007; Khorasheh and Zainali, 2001; Adam et al., 2012; Stangeland, 1974; Chou and Ho, 1988; Stangeland and Kittrell, 1972; Becker et al., 2015; Elizalde et al., 2016). Becker et al. (2016) compared a continuous lumping model with a single events microkinetic model. They concluded that the continuous lumping model provides better results for the macroscopic effluent characteristics. Also, the single events model is far more complex and computationally expensive than the continuous lumping model. Furthermore, Elizalde et al. (2009) reported that the continuous kinetic lumping model versus discrete lumped model is able to predict the whole distillation curve of hydrocracked products. Moreover, the prediction of product composition can be performed with different boiling point ranges. In this work, a typical continuous lumping model with five parameters was employed. The main advantages of the lumping technique are its easy computational implementation and the small amount of data required for parameter estimation.

One of the major challenges in continuous lumping modeling is optimization of parameters (Elizalde and Ancheyta, 2011). In this study, the model parameters have been optimized according to experimental product distributions by a PSO algorithm. All of the experimental results have been obtained under limited operating conditions in the hydrocracker setup.

The present work has used the continuous lumping concept for kinetic modeling of hydrocracking of a typical refinery vacuum gasoil. The considered five parameter model used TBP and reactivity (k) of the species as continuous variables for hydrocarbon mixture characterization. Next, Optimization of parameters was performed using a PSO algorithm as a novel optimization method. Finally, optimized parameters were used to predict product yields at various severities.

METHODOLOGY

Continuous Lumping Concept (Continuous Lumping Model Formulation)

The continuous lumping concept was first proposed by De Donder (1931). This concept is used when the mixtures include many different components and their difference is partially detectable. Under these conditions, a continuous variable, such as boiling point, molecular weight or other specifications can be used to describe mixture composition. In order to apply a continuous model, a continuous distribution function for the reactant mixture based on one of their properties needs to be defined. During the reaction, the distribution function changes. The aim of modeling is to investigate the distribution function variation and to obtain its final values. The continuous kinetic model is applicable for hydrocracking feed and products.

In this study, a typical proposed model based on continuous lumping for hydrocracking of petroleum fractions has been applied (Laxminarasimhan et al., 1996). The model considered the dimensionless boiling point (θ) and reactivity (k) as continuous variables for components identification. The dimensionless boiling point for component i is defined as follows:

(1) ${\mathit{\theta }}_i=\frac{{\mathrm{TBP}}_i-\mathrm{TBP}(l)}{\mathrm{TBP}(h)-\mathrm{TBP}(l)}$

where TBP(l) and TBP(h) indicate the lowest and highest possible boiling point in the reactant mixture, respectively. Usually experimental results of hydrocracking have been expressed in the form of cumulative weight percent (cum wt%) of components versus true boiling point(TBP). On the other hand, material balance and reactor performance equations have been written in the form of concentration (C) or non-cumulative weight percent versus reactivity (k) of components. Thus, for solving material balance equations and coordinate transformation, the continuous model considered a simple power law form relation between k_i and θ_i as below:

(2) $\frac{k_i}{k_{\max }}=\mathit{\theta }_i{}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathit{\alpha }$}\right.}$

Here k_max is the reactivity of the species with the highest boiling point, corresponding to θ =1, whereas α is a parameter of model. Both k_max and α are positive constants. Also this equation for θ = 0 yields the value of k equal to zero. In fact, the cracking rate constant for the lightest species is equal to zero.

Chou et al. (1988) described the coordinate transformation from the i axis (discrete mixture index) to the k axis (continuous mixture variable). They used a component type distribution function D(k) as the correction term for scale similarity along the axes.

(3) $\mathrm{di}=D(k\mathrm{).}\mathrm{dk}$

D(k) was considered as a Jacobian in coordinate transformation from i to k as follows:

(4) $D(k)=\frac{\mathrm{di}}{\mathrm{dk}}=\frac{\mathrm{di}}{d\mathit{\theta }}\frac{d\mathit{\theta }}{\mathrm{dk}}$

Assuming that the total number of components in the mixture is equal to N and (N→ ∞), then di/dθ can be replaced by N. Thus, the component type distribution function D(k) can be defined as follows:

(5) $D(k)=\frac{N\mathit{\alpha }}{k_{\max }^{\mathit{\alpha }}}{k}^{\mathit{\alpha }-1}$

The mass balance equation as a function of reactivity of species i (k_i ) and residence time t (or alternatively the inverse of space velocity, 1/LHSV), can be derived as:

(6) $\frac{\mathrm{dC}(k_i,t)}{\mathrm{dt}}=-k_{i}C(k_i,t)+{\int }_{\mathrm{ki}}^{k_{\max }}p(k_i,K)\mathrm{KC}(K,t)D(K)\mathrm{dK}(i=1\mathrm{...}N)$

In the above equation i=1 indicates the lightest component and i=N represents the heaviest component in the mixture. Furthermore, it is assumed that the model is one dimensional and homogeneous and neglects mass transfer resistances in boundary layer and pores of catalyst. The left side of the equation is the accumulation rate of species k_i. The first term on the right side represents the consumption reaction rate of k_i species during cracking to lower boiling point species, which assumed first order. The second term shows the generation rate of k_i species from heavier boiling point species with reactivity k_i < K < k_max. In the second term p(k_i,K) is the yield distribution function, that represents k_i species formation from K species cracking reactions. The function p (k_i, K) must have the following conditions:

1. The value of p (k_i, K) for k_i = K or k_i > K should be zero because the component with reactivity ki cannot be generated from itself or lighter components during hydrocracking.

2. The function p(k_i, K) must satisfy the following mass balance equation:

   (7) ${\int }_{\mathrm{ki}}^{k_{\max }}P(k,K\mathrm{).}D(k\mathrm{).}\mathrm{dk}=1$

3. Function p (k_i, K) should have a finite non-zero value when k_i = 0. Experimental findings show that when a component with reactivity K cracks, the lightest component in the mixture (k_i =0) is produced in traces.

4. The value of p (k_i, K) must always be positive.

5. A skewed Gaussian type distribution function such as Equation (8) can be considered to best depict the yield distributions and satisfy the above mentioned conditions.

(8) $P^{\left(k_i,K\right)}=\frac{1}{S_0·\sqrt{2\mathit{\pi }}}·\left[\exp \left(-{\left(\frac{{\left(k_i/K\right)}^{a_0}-0.5}{a_1}\right)}^2\right)-A+B\right]$

Terms A and B and S₀ of Equation (8) have been formulated as follows:

(9) $A=e^{-}{}^{{\left(\raisebox{1ex}{$0.5$}\!\left/ \!\raisebox{-1ex}{$a_1$}\right.\right)}^2}$

(10) $B=\mathit{\delta }\left(1-\left(\frac{k_i}{K}\right)\right)$

(11) $S_0=\underset{0}{\overset{K}{\int }}\left[\frac{1}{\sqrt{2\mathit{\pi }}}·\left(\exp \left(-{\left(\frac{{\left(k/K\right)}^{a_0}-0.5}{a_1}\right)}^2\right)-A+B\right)\right]·D\left(K\right)·\mathrm{dk}$

The parameters a₀, a₁ and δ are specific for each system and are used for model tuning.

Model Solution

Assuming that the model parameters including a₀, a₁, δ, k_max, α are known, then the solution method for Equation (6) is as follows:

For i=N or the component with maximum reactivity k_N (k_max), the integral term of Equation (6) will be equal to zero because it is not generated during cracking reactions. Therefore, this equation is converted to the following form:

(12) $\frac{\mathrm{dC}\left(k_N,t\right)}{\mathrm{dt}}=-k_{N}C\left(k_N,t\right)$

Thus the concentration of the heaviest component at time t, C(kN,t), is obtained from its value at time t-δt by the following equation:

(13) $C\left(k_N,t\right)=C\left(k_N,t-\mathit{\delta }t\right)e^{-k_N\mathit{\delta }t}$

For components i=1...N-1, the following implicit (backward Euler) finite difference equation was considered for numerical solution:

(14) $\frac{C\left(k_i,t\right)-C\left(k_i,t-\mathit{\delta }t\right)}{\mathit{\delta }t}=k_{i}C\left(k_i,t\right)+{\int }_{k_i}^{k_{\max }}p\left(k_i,K\right)\mathrm{KC}\left(K,t\right)D\left(K\right)\mathrm{dK}\left(i=1\mathrm{...}N-1\right)$

Considering an approximation, replacing concentration values at time t with the values at time t - δt in the integral term, the above equation is changed to Equation (15). This equation provides the possibility of calculating the concentration values of species i at time t from their values at time t - δt. Using small time steps δt in the calculation procedure could minimize this approximation error.

(15) $C\left(k_i,t\right)=\frac{C\left(k_i,t-\mathit{\delta }t\right)+\mathit{\delta }t·{\int }_{k_i}^{k_{\max }}p\left(k_i,K\right)\mathrm{KC}\left(K,t-\mathit{\delta }t\right)D\left(K\right)\mathrm{dK}}{1-k_i\mathit{\delta }t}\left(i=1\mathrm{...}N-1\right)$

The calculation procedure starts from residence time t=0 where the component concentrations are known and equal to the feed. Using Equation (2) C versus θ data are converted to C versus k form. After a δt time step, the component concentrations at new time t are calculated from Equation (13) for i=N and Equation (15) for i=1...N-1. Since the functionality of C(K,t) versus K in Equation (15) is unknown, therefore the integral term has no analytical solution. For node i=N-1 numerical integration from the trapezoidal rule was considered because the integration range included two nodes i=N-1 and i=N. For i=1...N-2 the more accurate cubic spline method for integration was considered because the integration range included three nodes or more. Cubic spline uses a third order polynomial for estimation of the internal function of the integral term.

The calculation procedure continued up to the reactor ultimate residence time with δt step times. Finally, all component concentrations in the reactor outlet will be determined. Due to the high calculation volume, a code was carried out using MATLAB 7.10 simulation software. A personal computer equipped with Intel (R) Core (TM) 2 CPU (3.0 GHz) and 3.25 GB of RAM was used.

In the above procedure for product concentration prediction, it was assumed that the model parameters are known. But usually existing data are in the form of feed and product concentrations and model parameters are unknown. Therefore, another program was developed for parameter optimization, so that the predicted product concentration had minimum difference with experimental values. Thus, the computer program was developed in two modes: the prediction mode and the optimization mode. The developed program was validated with data reported by Elizalde et al. (2009). Then it was employed for the parameter estimation of the experimental data.

In the present work, a PSO algorithm has been applied for the parameter optimization. This method has many advantages including: fast convergence, finding the global optimum in the presence of many local optima, and simple programming (Shokri et al., 2014; Santos et al., 2015). Figure 1 shows a flowchart for the optimization mode of the program.

Figure 1
Flowchart for the optimization mode of the program.

PSO method

The PSO method is a population based search optimization method that was originally proposed by Kennedy and Eberhart (1995). It is developed from swarm intelligence and is based on bird flocking and fish schooling to a promising position to achieve precise objectives in a multidimensional space. Particles fly around the parameter space with velocities which are dynamically adjusted according to their historical behavior and the particles have a tendency to fly towards the global solution over the course of the search process. In PSO, each potential solution of the problem is called a particle and the population of particles is called a swarm.

In the standard PSO, each particle i has a position X_i and velocity V_i that change according the following equality:

(16) $\left\{\begin{array}{c}{V}_i\left(m+1\right)=W\left(m\right)V_i\left(m\right)+C_1{R}_1\left(m\right)\left(P_{\mathrm{besti}}\left(m\right.\right)-X_1\left(\left.m\right)\right)+C_2{R}_2\left(m\right)\left(G_{\mathrm{besti}}\left(m\right.\right)-X_i\left(\left.m\right)\right)\hfill \\ X_i\left(m+1\right)=X_i\left(m\right)+V_i\left(m+1\right)\hfill \end{array}\right.$

In this equality, m is the iteration number; W represents an inertial weight used as a trade-off between global and local exploration capabilities of the swarm. Small values of W result in more rapid convergence usually on a suboptimal position, while a too large value may prevent divergence. P_besti and, G_besti represent the position with the 'best' objective value found so far by particle i and the entire population respectively. C₁ and C₂ represent two positive acceleration coefficients regulating the speed of the particle's flying with respect to the best global and local positions respectively, R₁(m) and R₂(m) represent random numbers uniformly distributed in the closed interval of [0, 1].

This algorithm terminates when no significant improvement in the objective function is observed over a number of iterations. There are different ways in which improvement can be measured. For example, if the average change in particle positions is small, the swarm can be considered to have converged. Alternatively, if the average particle velocity over a number of iterations is approximately zero, only small position updates are made, and the search can be terminated.

Experimental Setup

In this study, a hydrocracking setup was utilized to obtain experimental data sets. The schematic diagram of the system is depicted in Figure 2. The core of the set up is the isothermal tubular reactor having an internal diameter of 2.0 cm and total length of 2.0 m. The reactor material was stainless steel. Thermocouples for temperature indication were inserted in the stainless steel thermowell at the center of reactor. Reaction temperature was supplied by an electrical heating jacket in the four zones.

Figure 2
Schematic diagram of the hydrocracker setup.

A commercial NiMo supported on γ-alumina catalyst was used. In order to make the flow regime similar between the setup reactor and the industrial reactor, the ratio between the catalyst diameter and the reactor diameter should be approximately equal. For this purpose, the catalyst pellets were crushed to 0.8 mm.

According to this Figure for the process, the feed and high-pressure hydrogen are mixed before entering the reactor section. In the operating conditions of the reactor, a trickle bed flow regime is governing. Inside the reactor, heavier hydrocarbons are converted to lighter liquid and gas hydrocarbons in the presence of hydrogen and catalyst at high pressure. Moreover, organic sulfur and nitrogen are transformed into H₂S and NH₃ respectively. Reactor product is sent to the cooler. There it is cooled to ambient temperature. Thus, condensable hydrocarbon fractions are separated as liquid.

Then, a two-phase stream is sent into a high-pressure separator (HPS) for separation into hydrogen-rich gas and liquid hydrocarbon. Liquid hydrocarbon from the HPS is routed to the low pressure separator (LPS). The remaining dissolved gases are separated from the liquid in the LPS. Utilizing an Agilent gas chromatograph liquid products were analyzed by the simulated distillation (ASTM D2887) method. Analysis results are in form of true boiling point (TBP) versus distilled cumulative weight percent.

The properties of the vacuum gasoil feed are presented in Table 1.

Table 1
Characteristics of feed.

	Temperature(ºC)	Fraction(Cumulative wt %)
Simulated distillation results	260.2	IBP
	302.5	10
	351.7	50
	406.7	90
	466.5	FBP
Density at 15.6°C (g/cm3)	0.825
Kinematic viscosity (cSt)	9.86

The experiments were performed at different space velocities of feed and different temperatures for thermal and catalytic cracking. In thermal cracking experiments, for hydrodynamic similarity 0.8 mm quartz pellets were loaded inside the reactor.

Table 2 shows the operating conditions of the experiments. The volume ratio of hydrogen to liquid hydrocarbon feed at the reactor inlet was approximately 1200:1. Moreover, the reactor pressure was maintained at 110 bars in all experiments.

Table 2
Operating conditions of the experiments.

Experiment No.	Temperature (ºC)	LHSV (hr-1)
Catalytic cracking:
1	410	1.5
2	410	1.0
3	410	0.5
4	430	1.5
5	430	1.0
6	430	0.5
Thermal cracking:
1	440	1.5
2	440	1.0
3	440	0.5
4	470	1.5
5	470	1.0
6	470	0.5

Figures 3(a) and 3(b) give TBP versus cumulative weight percent for feed and product of catalytic and thermal cracking at 410°C and 440°C, typically. As can be seen in these Figures with increasing residence time or decreasing LHSV, the curves move down.

Figure 3
TBP versus cumulative weight percent: (a) catalytic cracking at 410°C; (b) thermal cracking at 440°C.

Cumulative weight percent (Cum wt%) versus θ are shown in Figures 4(a) and 4(b) for catalytic and thermal cracking at 410°C and 440°C, respectively.

Figure 4
Cumulative weight percent versus : (a) catalytic cracking at 410°C; (b) thermal cracking at 440°C.

RESULTS AND DISCUSSION

The continuous lumping model was applied to the results obtained from catalytic and thermal cracking experiments. Optimization of model parameters were performed for all operating conditions individually, using a PSO algorithm. Tables 3 and 4 show results of the model parameter optimization. As observed in these tables, the parameters are temperature dependent and the LHSV has no effect on them.

Table 3
Optimum values of parameters for catalytic cracking experiments.

RUN	1	2	3	4	5	6
T(C)	410	410	410	430	430	430
LHSV(hr-1)	1.5	1	0.5	1.5	1	0.5
α	1.391	1.41	1.398	1.181	1.183	1.189
a0	10.115	10.112	10.098	11.851	11.847	11.814
a1	29.21	29.19	29.14	38.361	38.334	38.231
kmax	0.243	0.242	0.24	0.793	0.791	0.788
δ ∙ 10e +5	2.314	2.318	2.32	2.744	2.748	2.754

Table 4
Optimum values of parameters for thermal cracking experiments.

RUN	1	2	3	4	5	6
T(C)	440	440	440	470	470	470
LHSV(hr-1)	1.5	1	0.5	1.5	1	0.5
α	1.358	1.362	1.354	1.205	1.207	1.21
a0	10.464	10.461	10.446	11.756	11.748	11.734
a1	30.147	30.138	30.128	38.09	37.96	37.81
kmax	0.251	0.250	0.248	0.781	0.780	0.776
δ ∙ 10e +5	2.33	2.331	2.34	2.723	2.726	2.731

For evaluation of the model performance, in any temperature condition parameters optimized according to experimental results at LHSV=1(hr^-1) were considered as the optimum set. Then by using this optimum set as input in the prediction mode of the computer program, product concentrations were calculated in other severities including LHSV=0.5 and LHSV=1.5. Typically, Figures 5(a) and 5(b) show predicted and experimental results for catalytic and thermal cracking at 430°C and 470°C, respectively.

Figure 5
Predicted and experimental results: (a) catalytic cracking at 430°C; (b) thermal cracking at 470°C.

In Figures 5(a) and 5(b), it is observed that the model gives acceptable product predictions at LHSV= 0.5 and 1.5, which indicates it suitable performance. Furthermore, the model results show good agreement with experimental results at LHSV=1. Consequently, the optimization algorithm has good performance.

In this study, the prediction performance of the model was evaluated based on two criteria: Root Mean Square error (RMSE) and Average Absolute Relative Error (AARE). These criteria are tabulated in Table 5. In this table y_i shows the experimental data and ŷ_i shows the model predicted data. A comparison of RMSE and AARE values for the model prediction is summarized in Table 6 for all experiments.

Table 5
Model performance criteria.

Performance criteria
Average Absolute Relative Error (AARE)	$\mathit{AARE}=\frac{1}{\mathit{n}}{\sum }_{\mathit{i}=1}^{\mathit{n}}\left|\frac{y_{\mathit{i}}-{\widehat{y}}_{\mathit{i}}}{y_{\mathit{i}}}\right|$
Root Mean Square Error (RMSE)	$\mathit{RMSE}=\sqrt{\frac{1}{\mathit{n}}{\sum }_{\mathit{i}=1}^{\mathit{n}}{\left(y_{\mathit{i}}-\widehat{y_{\mathit{i}}}\right)}^2}$

Table 6
Comparison of the performance of continuous lumping model.

Experiment Number	RMSE	AARE
Catalytic Cracking
1	0.902	0.039
2	0.845	0.029
3	0.986	0.041
4	0.892	0.035
5	0.747	0.024
6	0.861	0.026
Thermal cracking
1	0.921	0.036
2	0.842	0.029
3	0.896	0.032
4	0.969	0.037
5	0.762	0.022
6	0.941	0.024

The results reveal that the continuous lumping model indicates good prediction for LHSV= 1.5 and 0.5, using optimized parameters at LHSV=1. Maximum RMSE and AARE values were 0.986 and 0.041, in experiment Number 3.

Tables 7 and 8 compare product fraction yields (gasoline, diesel and residue) predicted by the model with those of the experimental data at different severities for catalytic and thermal cracking at 430°C and 470°C, respectively. The TBP range for gasoline is considered to be 30-150°C, and 150-370°C for diesel. The fraction of 370°C⁺ is classified as residue. Moreover, the normalized TBP (θ) range for these fractions is given in the tables. As observed, the model shows the best prediction for the diesel fraction.

Table 7.
Fraction yields prediction for catalytic cracking at 430°C.

Fraction	θ(range)	LHSV= 0.5(hr-1)			LHSV= 1(hr-1)			LHSV= 1.5(hr-1)
		Exp	Model	AARE	Exp	Model	AARE	Exp	Model	AARE
Gasoline (30-150°C)	0-0.22	5.89	6.23	0.059	7.64	7.83	0.024	15.40	16.13	0.047
Diesel (150-370°C)	0.22-0.77	69.76	70.73	0.014	71.28	72.08	0.011	73.52	71.74	0.024
Residue (370°C+)	0.77-1	24.35	23.04	0.054	21.08	20.09	0.047	11.08	12.13	0.095

• Exp=experimental

Table 8
Fraction yields prediction for thermal cracking at 470°C.

Fraction	θ(range)	LHSV= 0.5(hr-1)			LHSV= 1(hr-1)			LHSV= 1.5(hr-1)
		Exp	Model	AARE	Exp	Model	AARE	Exp	Model	AARE
Gasoline (30-150°C)	0-0.22	5.53	5.73	0.037	7.77	7.97	0.025	14.83	15.55	0.049
Diesel (150-370°C)	0.22-0.77	71.54	70.07	0.021	71.61	70.33	0.018	72.83	70.87	0.027
Residue (370°C+)	0.77-1	22.93	24.20	0.055	20.62	21.70	0.052	12.34	13.58	0.1

• Exp=experimental

CONCLUSIONS

A continuous lumping kinetic model was applied for the hydrocracking of vacuum gas oil feedstock in an experimental setup. The model parameters were optimized based on the experimental data for LHSV=1(hr^-1) utilizing a PSO method. Then, the model was used to predict the product distribution for other severities. In order to evaluate the performance of the proposed model, a wide range of experimental data was taken from a hydrocracker setup. Some statistical criteria (AARE, RMSE) were used to evaluate the prediction performance of models.

The comparison between the continuous model results and the experimental data shows that this model is capable of predicting the kinetic behavior well with an acceptable accuracy (maximum RMSE, 0.986 and maximum AARE, 0.041) for complex mixtures. Therefore, it was observed that the model can predict product yields with a highly acceptable degree of accuracy.

REFERENCES

• Adam, M., Arrighi, V., Ocone, R., Continuum lumping modelling for step growth polymerisation mechanism. Chemical Engineering Research and Design, 90, 2287-2292 (2012).
• Ahmed, H. S., Shaban, S. A., Menoufy, M. F., and El Kady, F. Y., Effect of catalyst deactivation on vacuum residue hydrocracking. Egyptian Journal of Petroleum, 22, 3, 367-372 (2013).
• Angeles, M. J., Leyva, C., Ancheyta, J. and Ramirez, S., A review of experimental procedures for heavy oil hydrocracking with dispersed catalyst. Catalysis Today, 220-222, 274-294 (2014).
• Ancheyta, J., Sanchez, S. and Rodriguez, M. A. Kinetic modeling of hydrocracking of heavy oil fractions: A review. Catalysis Today, 109, 76-92 (2005).
• Alhajree, I., Zahedi, G., Manan, Z. A. and Mohammad Zadeh, S., Modeling and optimization of an industrial hydrocracker plant. Journal of Petroleum Science and Engineering, 78, 627-636 (2011).
• Ancheyta, J., Lopez, F. and Aguilar, E., 5- lump kinetic model for gas oil catalytic cracking. Applied Catalysis A: General, 177, 227-235 (1999).
• Becker, P. J., Serrand, N., Celse, B., Guillaume, D. and Dulot, H., Comparing hydrocracking models: continuous lumping vs. single events. Fuel, 165, 306-315 (2016).
• Becker, P. J., Celse, B., Guillaume, D., Dulot, H., and Costa, V., Hydrotreatment modeling for a variety of VGO feedstocks: A continuous lumping approach, Fuel, 139, 133-143 (2015).
• Botchway, C., Dalai, A. K. and Adjaye, J., Kinetics of bitumen-derived gas oil upgrading using a commercial NIMO/Al2O3 catalyst. Canadian Journal of Chemical Engineering, 82, 478-487 (2004).
• Browning, B., Afanasieva, P., Pitaultb, I., Couenneb, F. and Tayakout-Fayollea, M., Detailed kinetic modelling of vacuum gas oil hydrocracking using bifunctional catalyst: A distribution approach. Chemical Engineering Journal, 15, 270-284 (2016).
• Canan, U. and Arkun, Y., Steady-state modeling of an industrial hydrocracking reactor by discrete lumping approach. Proceedings of the World Congress on Engineering and Computer Science (2012).
• Chou, M. Y., and Ho, T. C., Continuum theory for lumping nonlinear reactions. AIChE Journal, 34, 1519-1527 (1988).
• DeDonder, T., L'Affinite, Chapter 3, 2nd edition, Gauthier-Villars, Paris (1931).
• Elizalde, I., Rodriguez, M. and Ancheyt, J., Application of continuous kinetic lumping modeling to moderate hydrocracking of heavy oil. Applied Catalysis A: General, 365, 237-242 (2009).
• Elizalde, I. and Ancheyt, J., On the detailed solution and application of the continuous kinetic lumping modeling to hydrocracking of heavy oils. Fuel, 90, 3542-3550 (2011).
• Elizalde, I., Rodriguez, M. A. and Ancheyta, J., Modeling the effect of pressure and temperature on the hydrocracking of heavy crude oil by the continuous kinetic lumping approach. Applied Catalysis A: General, 382, 2, 15, 205-212 (2010).
• Elizalde, I., Mederos, F.S., and Mena-Cervantes, V.Y., Dynamic modeling of adiabatic reactor for hydrocracking of VGO by using of the continuous lumping approach. Fuel Processing Technology, 152, 200-206 (2016).
• Froment, G. F. Single event kinetic modeling of complex catalytic processes. Catalysis Reviews: Science and Engineering, 47, 83-124 (2005).
• Govindhakannan, J. and Riggs, J. B. On the construction of a continuous concentration-reactivity function for the continuum lumping approach. Industrial & Engineering Chemistry Research, 46(5), 1653-1656 (2007).
• Kumar, V. and Balasubramanian, P., Kinetic parameter estimation in hydrocracking using hybrid particle swarm optimization. Fuel, 88, 2171-2180 (2009).
• Khorasheh, F., Zainali, H., Chan, E. C. and Gray, M. R., Kinetic modeling of bitumen hydrocracking reactions. Petroleum & Coal, 43 (3-4), 208-18 (2001).
• Kennedy, J., Eberhart, R. C., Particle Swarm Optimization. Proceedings of IEEE Conference on Neural Network, 4, 1942-1948 (1995).
• Laxminarasimhan, C.S., Verma, R.P., Ramachandran, P.A., Continuous lumping model for simulation of hydrocracking, AIChE Journal. 42, 9, 2645-2653 (1996).
• Matos, E. M. and Guirardello, R., Modeling and simulation of a pseudotwo-phase gas-liquid column reactor for thermal hydrocracking of petroleum heavy fractions. Brazilian Journal of Chemical Engineering, 19, 3, 319-334 (2002).
• Moghadassi, A. Amini, N., Fadavi O. and Bahmani, M., Hydrocracking lumped kinetic model with catalyst deactivation in arak refinery hydrocracker unit. Journal of Petroleum Science and Technology, 1, 1, 31-37 (2011).
• Puron, H., Arcelus-Arrillaga, P., Chin, K.K., Pinilla, J.L., Fidalgo, B. and Millan, M., Kinetic analysis of vacuum residue hydrocracking in early reaction stages. Fuel, 117, 408-414 (2014).
• Qader, S. A. and Hill, G. R., Hydrocracking of gas oil. Industrial & Engineering Chemistry Process Design and Development, 8(1), 98 (1969).
• Rana M. S., Samano V., Ancheyta J. and Diaz J. a. I., A review of recent advances on process technologies for upgrading of heavy oils and residua. Fuel, 86, 12, 16-31 (2007).
• Sahu, R., Song, B.J., Im, J.S., Jeon, Y.P., Lee, C.W., A review of recent advances in catalytic hydrocracking of heavy residues. Journal of Industrial and Engineering Chemistry, 27, 12-24 (2015).
• Santos, L. S., Pagano, R. L., Calado, V. M. A. and Biscaia Jr., E. C., Optimization of a pultrusion process using finite difference and particle swarm algorithms. Brazilian Journal of Chemical Engineering, 32, 2, 543-553 (2015).
• Sadighi, S. and Ahmad, A., An optimization approach for increasing the profit of a commercial VGO hydrocracking process. The Canadian Journal of Chemical Engineering, 91, 1077-1091 (2013).
• Shokri, S., Sadeghi, M. T., Ahmadi Marvast, M., High reliability estimation of product quality using support vector regression and hybrid meta-heuristic algorithms. Journal of the Taiwan Institute of Chemical Engineers, 45, 2225-2232 (2014).
• Stangeland, B. E., Kinetic Model for prediction of hydrocaracker yields. Industrial & Engineering Chemistry Process Design and Development, 13, 72 (1974).
• Stangeland, B. E. and Kittrell, J. R., Jet fuel selectivity in hydrocracking. Industrial & Engineering Chemistry Process Design and Development, 11(1), 15-20 (1972).
• Valavarasu, G., Bhaskar, M. and Balarman, K. S., Petroleum Science and Technology, 21, 1185-1205 (2003).
• Zhou, H., Lu, J., Cao, Z., Shi, J., Pan, M., Li, W. and Jiang, Q., Modeling and optimization of an industrial hydrocracking unit to improve the yield of diesel or kerosene. Fuel, 90, 3521-3530 (2011).

Publication Dates

• Publication in this collection
  Apr-Jun 2018

History

• Received
  14 Aug 2016
• Reviewed
  15 Jan 2017
• Accepted
  03 Mar 2017