Axial dispersion modeling of industrial hydrocracking unit and its multiobjective optimization

Accepted Manuscript Title: Axial Dispersion Modeling of Industrial Hydrocracking Unit and its Multiobjective Optimization Authors: Hima Harode, Manojkumar Ramteke PII: DOI: Reference:

S0263-8762(17)30129-6 http://dx.doi.org/doi:10.1016/j.cherd.2017.02.033 CHERD 2596

To appear in: Received date: Revised date: Accepted date:

27-11-2016 13-2-2017 28-2-2017

Please cite this article as: Harode, Hima, Ramteke, Manojkumar, Axial Dispersion Modeling of Industrial Hydrocracking Unit and its Multiobjective Optimization.Chemical Engineering Research and Design http://dx.doi.org/10.1016/j.cherd.2017.02.033 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Axial Dispersion Modeling of Industrial Hydrocracking Unit and its Multiobjective Optimization

Hima Harode1, Manojkumar Ramteke1,* 1

Department of Chemical Engineering,

Indian Institute of Technology Delhi, India *E-mail: [email protected], [email protected] Tel No.: 91-1126591026 (Office); Fax No.: 91-1126581120

1

Graphical abstract

Highlights

1) Industrial hydrocracker is simulated using a more realistic axial dispersion model 2) 25-lump kinetics is used with two separate lumps for sulfur and nitrogen compounds 3) A validated model is optimized multiobjectively using genetic algorithm 4) Results show considerable decrease in total hydrogen used in the process cycle

Abstract: Hydrocracking is one of the most important refining operations used to crack heavy crude oil components into valuable lighter products. It processes less expensive, difficult to crack heavy feedstock to produce low sulfur and low aromatics content hydrocarbon fuels. In this study, a two-stage industrial hydrocracker is modeled using axial dispersion model to incorporate the mixing effects of feed flow. A 25-lump kinetics is used to simulate and validate the plant data. This validated model is then used for three-objective optimization in which the yields of aviation turbine fuel (ATF) and diesel are maximized simultaneously with minimization of the total amount of hydrogen required in a process cycle using a real-coded elitist non-dominated sorting genetic algorithm (RNSGA-II). The Pareto optimal results

2

obtained show the 26 % reduction in total amount of hydrogen required in a process cycle, with nearly same diesel and marginally low ATF production compared to plant data.

Keywords: Axial Dispersion Model; Multiobjective Optimization; Hydrocracking; Genetic Algorithm; Industrial Optimization Keywords Hydrocracking; Axial Dispersion; Multiobjective Optimization; Genetic Algorithm; Plant Simulation 1. Introduction: The continuous growth in demand for transportation fuel stimulated the surge of interest in the conversion of less expensive but difficult to crack feedstock into refined clean products. One of such popular upgrading process is hydrocracking. It converts less expensive heavy crude comprising a high percentage of impurities into required and valuable transportation fuel by catalytic cracking in the presence of hydrogen gas. Further, strict regulations for fuel specifications featuring ultra-low-sulfur, high cetane number, and low polyaromatics content are continuously being imposed (“Worldwide Fuel Charter, 5th edition,” 2013). Hydrocracking is one of the most suitable processes to meet these constraints, as along with cracking it involves saturation of cyclic aromatic content and removal of impurities like sulfur and nitrogen from fuel. It also increases the octane number of fuel and prevents the buildup of coke and tar formation. Hydrocracking process is accomplished in a two-stage operation of two trickle-bed reactors in series with multiple catalyst beds in each (see Fig. 1). The feed [vacuum gas oil (VGO) with true boiling point (TBP) cut range of 653 - 838 K] to the first reactor having four catalyst beds undergoes both hydrotreating (removal of sulfur and nitrogen compounds) and hydrocracking (cracking the heavier crude into lighter fractions) in presence of hydrogen gas whereas the second reactor having three catalyst beds is fed with bottom fractions [boiling point above 653 K] from the fractionator prominently involves hydrocracking (Bhutani et al., 2006; Mohanty et al., 1991). Since the hydrogenation reactions are exothermic in nature, the reactors are operated in an adiabatic mode where hydrogen itself is used for intermediate quenching (Bhutani et al., 2006; Mohanty et al., 1991) as shown in the Fig. 1. Recycle of hydrogen is required to maintain hydrogen partial pressure and utilize unused hydrogen gas. Both reactors use bi-functional Silica-Alumina catalyst (Bhutani et al., 2006; Mohanty et al., 1991). The products from both reactors go to high-pressure separators to recover hydrogen gas and then to the fractionator where desired components like kerosene, naphtha, diesel and gases like butane and hydrogen sulfide, etc. are separated.

3

The hydrocracking process is modeled by assuming the plug flow behavior in the reactor by several researchers (Alhajree et al., 2011; Bhutani et al., 2006; Mohanty et al., 1991; Pacheco and Dassori, 2002). However, such idealized flow approximation of fluid is hardly realized in actual trickle bed reactors used for the hydrocracking process in which a counter-current flow between feed to catalyst bed and hydrogen is present. There often exists a mixing in the flow direction i.e. a deviation from plug flow. To account for such nonidealities, an axial dispersion model (Levenspiel et al., 1999) is preferred in which the axial dispersion (mixing) is incorporated through dispersion number. Though such axial dispersion models are used by several researchers (Krishna and Saxena, 1989; Kumar and Sinha, 2012; Matos and Guirardello, 2000; Sadighi, 2013) to simulate the hydrocracking process, these are applied on simpler pilot plant operations involving a limited number of lumps (5 – 9). The present study aims at fulfilling this gap by extending the use of axial dispersion model with the existing kinetics of 23 lumps (Mohanty et al., 1991) to simulate a more complex industrial VGO hydrocracking unit. Further, the model includes hydrotreating reactions such as hydrodesulphurization and hydrodenitogenation in addition to usual hydrocracking reactions using two additional lumps and also incorporates the thermal effects. Therefore, the proposed model represents a more realistic behavior of actual operation compared to lumped-kinetic models involving plug-flow behavior. Further, the model is validated on industrial data. This validated model is then used for multiobjective optimization involving maximization of ATF, maximization of diesel and minimization of total hydrogen used in the plant, simultaneously. For this, a popular evolutionary algorithm, namely, realcoded elitist non-dominated sorting genetic algorithm (RNSGA-II) is used. The rest of the paper is arranged as literature review in section 2, model development in section 3, simulation in section 4, multiobjective optimization formulation in section 5, results and discussion in section 6 followed by conclusion in section 7. 2. Literature Review: One of the prominent study on reaction kinetics of hydrocracking in a pilot plant operation is the one reported by Qadar and Hill (Qader and Hill, 1969). In this, the vacuum gas oil (VGO) was categorized into several lumps based on the boiling point ranges for which the first order kinetics was established and the rate constants were evaluated experimentally. Several such studies (Martínez and Ancheyta, 2012; Sanchez et al., 2005) were reported on the kinetics of hydrocracking which established the suitability of the first order lumped kinetics for hydrocracking. Similarly, several studies (Chu and Wang, 1982; Mosby et al., 1986; Tsamatsoulis and Papayannakos, 1998) were reported in literature on hydrotreating which established the suitability of pseudo-second and pseudo-first order kinetics for

4

hydrodesulphurization and hydrodenitogenation, respectively. Further, LangmuirHinshelwood kinetics was also reported for hydrotreating in literature (Ancheyta and Sotelo, 2002; Korsten and Hoffmann, 1996; Schuit and Gates, 1973). Recently, Browning et al. (Browning et al., 2016) reported a detailed kinetic model of vacuum gas oil with 217 pseudocomponents based on carbon number and hydrocarbon family and validated against the experimental data reported by Henry et al. (Henry et al., 2012). Also, more accurate single event microkinetic models have been used in literature (Dewachtere et al., 1999; Laxmi Narasimhan et al., 2006; Van Geem et al., 2008) for modeling the hydrocracking process. However, these require detailed molecular composition of feedstock. In addition to the kinetics, the product distribution during cracking is an important aspect. One of the widely used models for feed distribution is given by Stangeland (Stangeland, 1974). This is a three- parameter model in which the first parameter gives values of rate constant according to boiling point whereas the remaining two parameters give the fractional distribution of cracking of each component into other components. Mohanty et al. (Mohanty et al., 1991) adopted this model for the distribution of component and used reaction kinetics proposed by Qadar and Hill (Qader and Hill, 1969) for two stage VGO hydrocracking with 23 pseudo components based on their boiling point ranges. The mass and energy balance were simulated and results were well in agreement with plant data. However, they did not include the hydrotreating reactions and assumed plug flow behavior. Pacheco and Dassori (Pacheco and Dassori, 2002) have proposed an improved version of the previous model (Stangeland, 1974) with two extra parameters. The model was applied to the industrial problem and gave better mass balance closure for the first reactor. Basak et al. (Basak et al., 2004) used continuum theory of lumping for modeling the commercial hydrocracking operation in which skewed Gaussian type distribution of yield to describe selectivity of hydrocracking reactions was incorporated. Valavarasu et al. (Valavarasu et al., 2005) modeled the hydrocracking operation using four-lump kinetics and tuned the model parameters using the experimental data available in the literature (Ali et al., 2002). Sildir et al. (Sildir et al., 2012) developed a dynamic, non-isothermal reactor model incorporating the continuous lumping of the components on the basis of a true boiling point. Further, the catalyst deactivation is also been incorporated in modeling of hydrocracking process in literature (Barkhordari et al., 1997; Moghadassi et al., 2011). The flow behavior in the reactor affects the performance of hydrocracker considerably. The countercurrent flow of gas and liquid develops mixing effect in the axial direction. To incorporate this mixing effect present in the reactor, axial dispersion model was applied by several researchers. Krishna and Saxena (Krishna and Saxena, 1989) simulated

5

VGO hydrocracking using axial dispersion model. In this model, the progress of cracking modeled by the distribution of boiling point around mid-boiling temperature and its decay kinetics. The Peclet number was calculated on the basis of the paraffinic content of component without incorporating the thermal effects. This study is extended by Matos and Guirardello (Matos and Guirardello, 2000) for different values of liquid superficial velocity, reactor length and diameter. They have shown that the extent of reaction decreases with the decrease in length and increase in velocity, due to a decrease in residence time. Recently, Kumar and Sinha (Kumar and Sinha, 2012) solved the five lump kinetic model for a tubular reactor using axial dispersion and showed that the heavy oil conversion increases with increases in temperature. Similarly, Sadighi (Sadighi, 2013) compared the effect of axial dispersion with ideal plug flow for a five lumped kinetic model of the hydrocracking reactor and found that the former is more efficient than the latter. Another interesting aspect is the optimization of hydrocracking process. The multiobjective optimization of the industrial hydrocracking unit was reported first time by Bhutani et al. (Bhutani et al., 2006). They have optimized the earlier model of Mohanty et al. (Mohanty et al., 1991) using a genetic algorithm for maximization of product yield and minimization of make-up hydrogen. Similarly, Li et al. (Li et al., 2008) proposed a six lump hydrocracking reactor model and tuned its kinetic parameter using Nelder-Mead minimization algorithm by minimizing the difference between plant and simulated values. Artificial neural network was applied by Alhajree et al. (Alhajree et al., 2011) to model the industrial hydrocracking operation. The model is then used to maximize the volume percent of gas oil, kerosene, light naphtha and heavy naphtha. Zhou et al. (Zhou et al., 2011) modeled and simulated hydrocracking reactions in hydrotreater and hydrocracker simultaneously using a discrete lumped kinetics. The model parameters were tuned and validated with industrial data using a genetic algorithm (GA). This validated model is then used for maximization of kerosene and diesel. The literature review illustrates that the use of lumped kinetic model with sufficiently large number of lumps and axial dispersion in feed flow is most suitable for the hydrocracking process. This led us to develop the 25-lump axial dispersion model by incorporating a more realistic effect of axial dispersion in the existing lumped kinetic model (Mohanty et al., 1991) for the industrial hydrocracking unit. The detailed model description is given next. 3. Model Development: The present work incorporates an axial dispersion in reactor model. Several authors have used axial dispersion effect (Krishna and Saxena, 1989; Kumar and Sinha,

6

2012; Matos and Guirardello, 2000; Sadighi, 2013) in material balance for pilot plant studies on hydrocracking. However, it is been implemented and validated on an industrial hydrocracking system for the first time in the present study. Moreover, the model used in the present study includes both hydrotreating and hydrocracking reactions simultaneously for industrial problem unlike previous studies (Alhajree et al., 2011; Bhutani et al., 2006; Mohanty et al., 1991) where just hydrocracking reactions were used. The axial dispersion model is formulated with assumptions that all cracking reactions follow the first order kinetics (Qader and Hill, 1969), the reactor operating adiabatically in steady state with negligible heat loss and no polymerization and isomerization reaction is taking place inside the reactor. Further, feed and products are assumed to be in liquid state and pure hydrogen is used as make-up and recycled gas. The kinetic scheme of 23 pseudo-components for hydrocracking (Mohanty et al., 1991) and two additional pseudo-components for hydrotreating is assumed. The concentration profile [i.e. variation of Ci (mass fraction, mass of component i /total mass of the feed)] with dimensionless reactor bed length z (length at axial location/total length of catalyst bed) of ith component in hydrocracking is obtained using the mass balance of ith component based on its cracking and formation from a cracking of all higher order jth components is given as follows: N   w 1 d 2Ci dCi   k C   0; i  1, 2,..., 23; j  6, 7,..., 23  i i  k j Pij C j  2 Pe dz dz  j  F vh Pe  D

(1)

In this, ε is the bed porosity, η is catalyst effectiveness, w is weight of catalyst (kg), F is total mass flow rate (kg/hr), k is reaction rate constant [(kg of feed)/(kg of catalyst)(h)], Pe is the Peclet number for the axial dispersion, v is the feed velocity (m/s), h is the length of the reactor (m) and D is the dispersion coefficient (m2/s). Similar to hydrocracking, the mass balance is performed for hydrotreating (viz. hydrodenitrogenation and hydrodesulfurization) reactions. These reactions considered to be pseudo-first and -second order respectively (Chu and Wang, 1982; Mosby et al., 1986; Tsamatsoulis and Papayannakos, 1998) as hydrogen is in excess. The mass balance of their corresponding lumps which account nitrogen and sulphur compounds is given as follows:

1 d 2CN dCN  w  (-k N CN ) 0 2 Pe dz dz F 1 d 2CS dCS  w  (-kS CS2 ) 0 2 Pe dz dz F

(2)

7

Here, CN and CS are the mass fraction of nitrogen and sulphur present in the stream; kN and kS are the reaction rate constants of hydrodenitrogenation (Qader and Hill, 1969) and hydrodesulfurization (Tsamatsoulis and Papayannakos, 1998) reactions [(kg of feed)/(kg of catalyst)(h)]. In the mass balance equations, Peclet number of the axial dispersion for the corresponding pseudo-component is calculated using the dispersion coefficient. The corresponding dispersion coefficient is calculated using following empirical equation (Salmi et al., 2010):

D 3 107 1.35   ud Re2.1 Re1/8

(3)

Here, Re is Reynold’s number, u is the velocity of fluid (m/s) and d is the diameter of reactor (m). Pe varies in the range of 80-100 for the given system and a constant value of 85 is taken throughout the calculations. These mass balance equations are then coupled with model given by Mohanty et al. (Mohanty et al., 1991) (see Table S1 from supplementary material). The details are discussed next. In the mass balance equation (1), the first term includes the dispersion effect; the second term includes the variation of concentration of each component across the bed length whereas the third term comprises two subparts which correspond to the reaction. In the third term, the first part represents the disappearance of ith pseudo component whereas the second term gives the rate of formation of the ith pseudo component by cracking of components higher than it. The components 1 to 23 are in increasing the order of boiling point (BP). Only components having BP above 400K will undergo cracking. In this problem, the components 1 to 5 have BP below 400K hence cracking starts with the 6th component. Feed to the first reactor consists of components 13th to 23rd whereas product includes 1st to 23rd components. The reaction rate constant of vacuum gas oil with an average boiling point of 638K is reported by Qader and Hill (Qader and Hill, 1969) [see equation (S1) in supplementary material]. This average value is then multiplied with the relative rate constant (Ki) to give rate constant (ki) [see equation (S2)] for individual component i. The relative rate constant is calculated as per variation of the carbon number of normal alkanes C5/C10/C15/C20 with increment in their cracking rates as 1/32/72/120 (Rapoport, 1963). Product distribution correlations (Stangeland, 1974) are used to estimate the distribution of cracking [see equations (S5 – S8) from Table S1 of supplementary material]. In these, the mass fraction (P1j) of butanes and a lighter component formed from heavier components represented by j is estimated [equation (S5)] using boiling point, tbj, of

8

component j (ºC) and a constant parameter C which depends on characteristics of feed and catalyst. To calculate the yield of the compound heavier than butanes, first their boiling points are normalized [equation (S6)] to give normalized temperature (yij) for ith product formed from jth pseudo component. The cumulative yield till the ith pseudo component from the jth pseudo component (Pij) calculated using yij [equation (S7)] and constant parameter B. The subtraction of Pi-1, j from Pij [equation (S8)] gives the actual yield of ith component from component j. The parameter B depends on feed and catalyst. It gives the linearity of breaking. A zero value of B shows a linear distribution of breaking, whereas the value equals to -2 gives breaking into two equal halves. The constant C mainly defines the distribution of butane (Stangeland, 1974). Thermodynamic properties such as heat capacity and enthalpy are required to calculate the heat of reaction to obtain temperature distribution across the reactor length while solving energy balance. Empirical correlations for the molecular weight (MW), critical temperature (Tc), critical pressure (Pc), critical volume (Vc), critical temperature of mixture (Tcm), critical volume of mixture (Vcm) and compressibility factor (Zc), are calculated as a function of normal boiling point and density whereas acentric factor (ɷ) is being calculated as a function of the reduced boiling point [equation (S10-S17)] (Kesler et al., 1979; Lee and Kesler, 1975). Similarly, kinematic viscosity [equation (S18)] (Dutt, 1990) and ideal gas enthalpy for pure components (Hºi) and that of a mixture (Hºm) are calculated on the basis of the boiling points and specific gravities of components [equation (S27-S30)] (Weir and Eaton, 1932). Fugacity and excess enthalpy of components and mixture are calculated using Peng-Robinson’s equations of state [equations (S19-S27)]. The heat of reactions is calculated on the basis of hydrogen consumed in each reaction. It is assumed that feed consists mainly of paraffin and naphthenes, hence at standard condition the amount of heat released during any reaction is taken as 42 MJ/kmol of hydrogen consumed. For hydrodesulfurization reaction, the released amount of heat is taken as 62 MJ/kmol, whereas for hydrodenitrogenation an amount of 251 MJ/kmol is released (Mederos and Ancheyta, 2007). To estimate the hydrogen consumption, the value of C/H ratio plays a key role. This C/H value was obtained by taking an average of values from a monograph of API Technical Data Handbook,1977 and tabulated data of C/H values from Nelson on the basis of UOP characterization factor (Mohanty et al., 1991). The C/H ratio (Ri) and Pi,j gives product i formed from pseudo-component j. Using these values the total carbon content and hydrogen content in the products formed from component j is calculated [equations (S35-S36)] which further gives hydrogen content of pseudo-components [equation (S37)]. These values are then used to calculate hydrogen consumption per unit mass of product [equations (S38-S40)]. The enthalpy (Hi´) of reactant and product at reaction conditions is calculated using Peng-

9

Robinson’s equation of state (EOS) and the enthalpy of the liquid hydrocarbon (Hiº) at standard condition is calculated using the correlation [equation (S42)]. Butanes and lighter hydrocarbons combined to form pseudo-component 1. These are in a vapor phase at reference conditions, hence their enthalpy is calculated using heat capacity of butanes (CP1) [equation (S43)] (Reid and Prausnitz, 1977). Finally, the heat of reaction at reaction condition for hydrocracking is calculated [equation (S45)] using the enthalpy of hydrogen at standard conditions (HºH2) and reference conditions (HH2) calculated from heat capacity data (Vargaftik, 1975). Further details can be obtained from Mohanty et al. (Mohanty et al., 1991). 4. Simulation: The formulated reactor model is simulated in MATLAB 2013. The initial conditions such as plant operating conditions and feed characteristics are given in Table S2 and S3 (supplementary material) respectively. Thermodynamic properties, product distribution and critical properties are calculated using equations given in Table S1 (supplementary material). The mass balance [equations (1) – (3)] leads to second-order ordinary differential equations of boundary value type which are solved using finite difference method for each pseudocomponent. In contrast to this, the energy balance leads to a first-order ordinary differential equation of initial value type. This is solved using [equation (S9) from Table S1 of supplementary material] Runge –Kutta- Gill method to obtain temperature profiles across each catalyst bed. Since the energy and mass balances are coupled with each other, these are solved repeatedly by updating the concentration and temperature profiles into each other using successive substitution method [see Fig. 2]. Initially, a constant temperature is assumed throughout the reactor bed and concentration profile is calculated by solving mass balance equations. The mass balance simulation provides the concentration profiles across the catalyst beds. These concentration profiles are then used with interpolation to solve the energy balance for temperature profile across the bed length. The obtained temperature profile is then used in mass balance to calculate the concentration profiles across the bed. This cyclic procedure is repeated successively to get results with high precision. The first reactor consists of four catalyst beds. The temperature rise due to an exothermic reaction is controlled by hydrogen gas used for quenching after each bed. The inlet temperature of feed for the next bed is obtained by applying a heat balance between quenching hydrogen (354K) and mass flow from the previous bed. Using the procedure described above, all beds are simulated successively and final pseudo-component concentration and temperature are obtained at the exit of the first reactor. Also, the procedure is repeated for the second reactor with three catalyst beds.

10

The heavier components (boiling above 643 K), a bottom product, obtained from the distillation of effluent mixture from the first and second reactor, is sent to the second reactor for further cracking. The pseudo-components obtained from fractionation, are distributed into products as per their boiling point range. Pseudo-component 1 gives light-ends like butane; components 2 – 6 form Naphtha and components 7 – 10 and 2/5th of the component 11 form kerosene, components 12 – 15, 3/5th of 11th and 1/5th of the 16th component form diesel. The components 17 – 23 fall in the range of recycling components. The successive substitution method is used to calculate the mass flow rate of recycle stream. The hydrogen consumption is given by the difference between hydrogen content in product and reactant after cracking. The sulfur in the feed is assumed to be in the form of benzothiophene which requires 3 moles of hydrogen per mole of sulfur consumed. Similarly nitrogen compound requires 4 moles of hydrogen per mole of nitrogen consumed. Hence the difference between the initial and final value of sulfur and nitrogen multiplied with their respective hydrogen requirement number add up to give overall hydrogen consumption. This simulated model of industrial hydrocracking process is then used for multiobjective optimization and the problem formulation is described next. 5. Multiobjective Optimization Formulation: Multiobjective optimizations of industrial systems are widely reported (Chaudhari and Gupta, 2012; Gujarathi and Babu, 2016; Guria et al., 2014; Halim et al., 2015; Kasat and Gupta, 2003; Mogilicharla et al., 2014; Nandasana et al., 2003; Patel and Padhiyar, 2016; Rangaiah, 2010, 2008; Rangaiah and Bonilla-Petriciolet, 2013; Sharma and Rangaiah, 2016; Sultana et al., 2016) using Genetic Algorithm (GA) and its variants (Coello et al., 2002; Deb, 2001; Holland, 1975). One of such popular variant is real-coded elitist non-dominated sorting genetic algorithm (RNSGA-II) (Deb, 2001). RNSGA-II starts with a random population of chromosomes constituted by decision variables of a problem. This random population undergoes the genetic operations such as selection, simulated binary crossover and polynomial mutation to produce the offspring population. This offspring population is then combined with the parent population and the best solutions from these are allowed to survive. These best elite solutions are then used as a parent population in the next generation and the process is repeated till the optimal results are obtained or the user defined a maximum number of generations are executed. The details of the algorithm can be obtained from Deb (Deb, 2001) and Coello Coello et al. (Coello et al., 2002). In the present study, RNSGA-II is used for optimizing the industrial hydrocracking system. The hydrocracking process produces high speed diesel, aviation turbine fuel, naphtha, butane and lighter gases as a product. As per market value, the increase in production of the selective component with minimum utilization of resources can make more profit. Since multiple objectives are naturally present, a multiobjective optimization is required to gain

11

optimal conditions to meet demands of particular products. The aviation turbine fuel (ATF) or kerosene makes good market demand and higher cost for its use in aircraft and domestic purpose. Therefore, the first objective of the present study is to maximize the production of ATF. Similarly, high speed diesel being widely used as a transportation fuel has very high market demand. Therefore, the second objective is to maximize the production of diesel. Since most of the two-stage hydrocracking plants operate with the conversion range of 98 – 99 %, there is a little scope for improving the yield of particular product. Often this is at the expense of other products. However, the operating expense of H2 is the real concerns in the industry. Hydrogen is used for both quenching and hydrogenation purpose. The total hydrogen used in the process cycle comprises the total of five quench and two reactant streams. The hydrogen utilization for hydrocracking is around 1% of total hydrogen used and rest of it is used for quenching (Bhutani et al., 2006). Hence, hydrogen is recycled continuously that increases the compressing cost. Thus, any minimization in aggregate hydrogen in the process cycle can earn value to the plant. Therefore, minimization of the total hydrogen used in a process cycle is an important objective. Based on the above description, a three-objective optimization problem consistent with the industrial requirement is formulated for the first time as follows: Problem 1: Max F1 = ATF (Aviation turbine fuel)

(4)

Max F2 = Diesel

(5)

Min F3 = H2 in process cycle

(6)

The decision variable selected are the feed input temperatures, [T1, T2 (K)] in reactor 1(R1) and 2 (R2), inlet hydrogen flow rate (HPM1, HPM2 m3H2/m3feed hr-1) in both the reactors, flow rate of quenching [Q1 (i), Q2 (i), m3/hr] hydrogen in all seven catalyst beds. The upper bounds and lower bounds for these variables are taken from literature (Bhutani, 2007), which are as follows:

650  T1  680 625  T2  655 600  HPM1  900 600  HPM1  900

35, 000  Q1 (i)  50, 000;i=1,..3 7000  Q2 (i)  9000;i=1,..2

(7)

The outlet temperature (t1, t2, K) of catalyst beds are kept below their maximum operating values. Hence following constraints are added to the model

t1 (i)  714;i=1,..,4 t2 (i)  700;i=1,..,3

(8)

12

6. Validation and Results and Discussion: The developed model is first tuned using the plant data reported by (Mohanty et al., 1991) and then extended to multiobjective optimization. Such tuning is also been reported in the literature (Bhutani et al., 2006; Mohanty et al., 1991). Similar to these, the rate constants [equation (S2)] are tuned by multiplying the constant factors CF1 and CF2 for first and second reactors, respectively as these are obtained from the best fit polynomial equation [equation (S2)] based on Rapaport principle (Rapoport, 1963). In addition to these, parameters B and C of product distribution equations [(S5 and S7)] are also tuned to match the obtained product distribution as suggested by (Stangeland, 1974). In order to select the best values of these, CF1, CF2, B and C are varied between the ranges: 0.9 – 1.0, 0.7 – 1.0, 0.7 – 1.0, and 0.25 – 0.4, respectively (Mohanty et al., 1991). After the rigorous tuning, the best values obtained for these are 0.9790, 0.860, 0.9900 and 0.2980, respectively. However, the best values reported by (Mohanty et al., 1991) were 0.9906, 0.7844, 0.7 and 0.37, respectively for plug flow model in which hydrotreating, catalyst effectiveness and bed porosity are not been included, unlike the developed axial dispersion model. With the tuned model, the obtained results are compared with Plant data and the results obtained by Mohanty et al. (Mohanty et al., 1991) using pure plug flow model without hydrotreating as shown in Table 1. The results are well in agreement with the plant data except for the error of 9.3 % in hydrogen consumption in the second reactor. Also, the obtained results of the proposed model are better than that reported in (Mohanty et al., 1991). Further, the obtained results show the conversion of ~ 99 % which is in good agreement with industrial practice. For further validation of axial dispersion model with reactor operation, temperature profiles for catalyst beds were plotted in Fig. 3 and are shown by black colored lines. In this figure, the temperature profiles of pure plug flow model (Mohanty et al., 1991) are also shown by red lines for comparison. Fig. 3 (a) shows the temperature profile along the catalyst bed length for the first reactor. Due to exothermic nature of the hydrocracking reactions, the steep increase in temperature was expected and is captured by both models. The quenching with hydrogen between the catalyst beds brings down the temperature in the range of 680 K. The temperature variation in the second reactor [see Fig. 3(b)] is quite different from that of the first reactor. The temperature drop along the bed of reactor is quite less between the beds due to less amount of hydrogen available for quenching. These graphs [Fig. 3 (a –b)] show the same pattern as expected for exothermic reactions along with quenching. However, the temperature prediction of axial dispersion model is higher than that of the plug flow model. Since, the CF1 used in the present study is slightly lower than that used by Mohanty et al. (1991), the rate of reactions are expected to be lower and therefore the temperature profiles are also expected to be lower. However, the contrasting trends in temperature profiles of first stage are obtained since the plug flow model [Mohanty et al.

13

(1991)] used for comparison did not accounted the hydrotreating reactions unlike the developed axial dispersion model. So the increase in temperature is primarily due to additional contribution from hydrotreating reactions which predominantly occur in the first stage. However, the temperature profiles for second stage which does not involve hydrotreating are also higher than plug flow model. This increase is primarily due to use of higher CF2 value than that of Mohanty et al. (1991). This clearly illustrates that the temperature profiles obtained are along the expected lines. Similar to temperature profiles, the concentration profiles (mass fractions) for one component each representing butane (C1), naphtha (C4), kerosene (C9), diesel (C14), feed (C20) and higher components (C23) are plotted across the beds of first and second reactors as shown in [Fig. 4(a-b)] for axial dispersion and pure plug flow model without hydrotreating. These profiles show the slightly lower rate of formation of components corresponding to products (i.e. C1, C4, C9, C14) and consumption of components corresponding to feed (i.e. C20, C23) in axial dispersion model even though the higher temperature profiles compared to plug flow model compensates reduction in rate due to catalyst effectiveness factor and bed porosity. The additional concentration profiles (mass fractions) for sulfur and nitrogen of axial dispersion model are shown in Fig. 5. These profiles match with the industrial operation in which the nitrogen and sulfur mostly removed in first two catalyst beds (Mohanty, 1989). Thus, the obtained concentration profiles further validate the correctness of the developed model. Since the axial dispersion model depicts a more realistic flow behavior and also includes the hydrotreating in the kinetics, the results shown in Figs. 3 and 4 are more realistic than the pure plug flow model without hydrotreating used in the literature (Mohanty et al., 1991).

This validated model is then used for multiobjective optimization using RNSGA II. The GA parameters used for optimization are selected based on the values reported in the literature (Deb, 2001) and are listed in Table 2. Since the simulation is computationally intensive, the population size of 50 and number of generations equals to 100 are used. The computational time taken for optimization is ~60 hrs on a desktop computer with a configuration of Intel Xeon E3 – 1225 [email protected] GHz processor, 8 GB RAM and Windows 7 operating system. Fig. 6 represents the Pareto optimal solutions for minimization of total hydrogen in the process with maximization of HSD and ATF simultaneously, normalized with their plant value. These normalized values are shown by asterisk *. The optimal results corresponding to the robust operating point (i.e. a point with inflation is shown by a red circle) show a decrease of 26 % in total H2 in process cycle compared to industrial practice for a nearly same yield of HSD and marginal decrease (6 %) in ATF yield. The results corresponding to a robust operating point are analyzed in detail.

14

As the optimum result shows a significant decrease in total hydrogen in process cycle compare to the increase in production of valuable products, all decision variables are plotted (Fig. 7 and 8) against the total hydrogen used in a process cycle normalized with its plant value. With the increase in inlet temperature in both reactors, reaction rates increases which increase the production and therefore increases inlet hydrogen for reaction (HPM1 and HPM2). This also increases the amount of hydrogen required for quenching in both reactors. Therefore with the increase in inlet temperature in both reactors, the total hydrogen flow increases as shown in Figs. 7a and 7b. In these figures, the optimum inlet temperature corresponding to the robust operating point is also shown by a red circle. This point also shows the tradeoff characteristics as like Pareto optimal plot is shown in Fig. 6. Here one can see that the robust temperature value is very close the inflation point on the curve. Similarly, the inlet hydrogen for reaction in both reactors (HPM1 and HPM2) is plotted against total hydrogen in the process (see Figs. 7c and 7d). In this, the increase in former increases the latter since HPM1, HPM2 and quenching hydrogen requirement gives total hydrogen in the process. Further, the HPM1 and HPM2 values corresponding to robust optimization solution are shown by a red circle in Figs. 7c and 7d. In a similar manner, the effect quenching hydrogen requirement on production is shown in Fig. 8. Since reactor 1 has three catalyst beds and reactor 2 has two catalyst beds, their individual requirement for quenching hydrogen with respect to aggregate hydrogen in the process are shown in Figs. 8a and 8b. Also, the total quenching hydrogen requirement in both the reactors against total hydrogen in the process is shown Figs. 8c and 8d. Similar to the HPM1 and HPM2, increase in quenching hydrogen increases the aggregate hydrogen in the process. Also, the values corresponding to the robust optimum point (shown by a red circle) are close to the lower bounds.

Further, the temperature and concentration profiles of selected components corresponding to the robust optimum point are compared to the profiles corresponding industrial simulation as shown in Figs. 9 (a-b) and 10 (a-b). Fig. 9 (a-b) illustrates that the temperature to be maintained in the first reactor should be lower in all 4 beds whereas it should be slightly higher in first two beds and lower in third bed in the second reactor compared to that of industrial simulation. Therefore it requires lower quenching hydrogen which reduces the total hydrogen requirement in the process. Also, the concentration profiles of selected components corresponding to different products and feed are compared with the industrial simulation in Fig. 10(a-b). These clearly illustrate the slightly lower production in first and nearly same production in the second reactor which is quite expected as per the temperature profiles [Fig. 9(a-b)]. However, it significantly reduces the total amount of H2 in a process cycle with the same yield of HSD and the marginal decrease in ATF yield. Also, the

15

conversion corresponding to this operating point is 59 and 98.8 % in reactor 1 and reactor 2 which is acceptable in the industrial practice. The concentration profiles for the sulfur and nitrogen in the first reactor are also shown for industrial simulation and those corresponding to the robust operating point in Fig. 5. Since hydrotreating occurs only in the first reactor, these profiles are shown only for the first reactor. From the profile, it can be seen that these occur primarily in first two catalyst beds. Also, there optimum concentration profiles for sulfur and nitrogen are found to be coinciding with the industrial simulation profiles.

Conclusion: In the present study, a more realistic axial dispersion model with 25-lump kinetics is used to model, simulate and optimize the industrial hydrocracking unit. The developed model includes both hydrotreating and hydrocracking reactions of hydrocracker. In this, 23 lumps are used to represent paraffin and naphthenes contents whereas 24th and 25th lumps are used to represent nitrogen and sulfur contents respectively. The model provides values of product formation, consumed hydrogen and temperature across reactor beds and the prediction shows a close match with the plant data. It also shows better results compared to earlier plug flow model. The validated model is then used for multiobjective optimization using real coded NSGA-II. A realistic three-objective optimization problem is formulated based on the industrial practice. An operating point in the obtained Pareto optimal front is identified which shows 26 % decrease in total hydrogen in a process cycle (H2) compared to industrial practice with an almost same value of HSD and a marginal decrease in ATF.

Acknowledgement: The authors gratefully acknowledge the partial financial support from Science and Engineering Research Board, Department of Science and Technology, Government of India, New Delhi [through grant No. SB/FTP/ ETA-0125/2013, dated June 5, 2014].

16

References: Alhajree, I., Zahedi, G., Manan, Z.A., Zadeh, S.M., 2011. Modeling and optimization of an industrial hydrocracker plant. J. Pet. Sci. Eng. 78, 627–636. doi:10.1016/j.petrol.2011.07.019 Ali, M.A., Tatsumi, T., Masuda, T., 2002. Development of heavy oil hydrocracking catalysts using amorphous silica-alumina and zeolites as catalyst supports. Appl. Catal. A Gen. 233, 77–90. doi:10.1016/S0926-860X(02)00121-7 Ancheyta, J., Sotelo, R., 2002. Kinetic modeling of vacuum gas oil catalytic cracking. J. Mex. Chem. Soc. 46, 38–42. Barkhordari, A., Fatemi, S., Daneshpayeh, M., Zamani, H., 1997. Development Of A Discrete Kinetic Model For Modeling Of Industrial Hydrocracking Reaction Accompanied With Catalyst Deactivation 1–6. Basak, K., Sau, M., Manna, U., Verma, R.P., 2004. Industrial hydrocracker model based on novel continuum lumping approach for optimization in petroleum refinery. Catal. Today 98, 253–264. doi:10.1016/j.cattod.2004.07.056 Bhutani, N., 2007. Modelling, simulation and multi-objective optimization of industrial hydrocrackers. National University of Singapore. Bhutani, N., Rangaiah, G.P., Ray, A.K., 2006. First-principles, data-based, and hybrid modeling and optimization of an industrial hydrocracking unit. Ind. Eng. Chem. Res. 45, 7807–7816. doi:10.1021/ie060247q Browning, B., Afanasiev, P., Pitault, I., Couenne, F., Tayakout-fayolle, M., 2016. Detailed kinetic modelling of vacuum gas oil hydrocracking using bifunctional catalyst : A distribution approach. Chem. Eng. J. 284, 270–284. doi:10.1016/j.cej.2015.08.126 Chaudhari, P., Gupta, S.K., 2012. Multiobjective Optimization of a Fixed Bed Maleic Anhydride Reactor Using an Improved Biomimetic Adaptation of NSGA-II. Chu, C.I., Wang, I., 1982. Kinetic study on hydrotreating. Ind. Eng. Chem. Process Des. Dev. 21, 338–344. Coello, C.A.C., Van Veldhuizen, D.A., Lamont, G.B., 2002. Evolutionary algorithms for solving multi-objective problems. Springer. Deb, K., 2001. Multi-objective optimization using evolutionary algorithms. John Wiley & Sons. Dewachtere, N. V., Santaella, F., Froment, G.F., 1999. Application of a single-event kinetic model in the simulation of an industrial riser reactor for the catalytic cracking of vacuum gas oil. Chem. Eng. Sci. 54, 3653–3660. doi:10.1016/S0009-2509(98)00518-1 Dutt, N.V.K., 1990. A Simple Method of Estimating the Viscosity of Petroleum Crude- Oil and Fractions. Chem. Eng. J. Biochem. Eng. J. 45, 83–86. Gujarathi, A.M., Babu, B. V, 2016. Evolutionary Computation: Techniques and Applications. Guria, C., Goli, K.K., Pathak, A.K., 2014. Multi-objective optimization of oil well drilling using elitist non-dominated sorting genetic algorithm. Pet. Sci. 11, 97–110. doi:10.1007/s12182-014-0321-x

17

Halim, I., Adhitya, A., Srinivasan, R., 2015. A novel application of genetic algorithm for synthesizing optimal water reuse network with multiple objectives. Chem. Eng. Res. Des. 100, 39–56. doi:http://dx.doi.org/10.1016/j.cherd.2015.05.015 Henry, R., Tayakout-Fayolle, M., Afanasiev, P., Couenne, F., Lapisardi, G., Simon, L.J., Souchon, V., 2012. Methodology for the Study of Vacuum Gas Oil Hydrocracking Catalysts in a Batch Reactor – Coupling of GC-2D Data with Vapour-Liquid Equilibrium, in: Material Sciences and Technology, Advanced Materials Research. Trans Tech Publications, pp. 207–213. doi:10.4028/www.scientific.net/AMR.560561.207 Holland, J.H., 1975. Adaptation in natural and artificial systems: an introductory analysis with applications to biology, control, and artificial intelligence. U Michigan Press. Kasat, R.B., Gupta, S.K., 2003. Multi-objective optimization of an industrial fluidized-bed catalytic cracking unit (FCCU) using genetic algorithm (GA) with the jumping genes operator. Comput. Chem. Eng. 27, 1785–1800. doi:http://dx.doi.org/10.1016/S00981354(03)00153-4 Kesler, M.G., Lee, B.I., Sandler, S.I., 1979. A Third Parameter for Use in Generalized Thermodynamic Correlations. Ind. Eng. Chem. Fundam. 18, 49–54. doi:10.1021/i160069a012 Korsten, H., Hoffmann, U., 1996. Three‐ phase reactor model for hydrotreating in pilot trickle‐ bed reactors. AIChE J. 42, 1350–1360. doi:10.1002/aic.690420515 Krishna, R., Saxena, A.K., 1989. Use of an axial-dispersion model for kinetic description of hydrocracking. Chem. Eng. Sci. 44, 703–712. doi:10.1016/0009-2509(89)85045-6 Kumar, A., Sinha, S., 2012. Steady state modeling and simulation of hydrocracking reactor. Pet. Coal 54, 59–64. Laxmi Narasimhan, C.S., Thybaut, J.W., Martens, J.A., Jacobs, P.A., Denayer, J.F., Marin, G.B., 2006. A unified single-event microkinetic model for alkane hydroconversion in different aggregation states on Pt/H-USY-zeolites. J. Phys. Chem. B 110, 6750–6758. doi:10.1021/jp054981j Lee, B.I., Kesler, M.G., 1975. A generalized thermodynamic correlation based on threeparameter corresponding states. AIChE J. 21, 510–527. doi:10.1002/aic.690210313 Levenspiel, O., Wiley, J., Hepburn, K., Levenspiel, O., 1999. Chemical Reaction Engineering - third edition, 3th edition. doi:10.1016/0009-2509(64)85017-X Li, Q., Jiang, Q., Cao, Z., Ai, L., Zhang, Y., 2008. Modeling and simulation for the hydrocracking reactor. Proc. 27th Chinese Control Conf. CCC 204–208. doi:10.1109/CHICC.2008.4604960 Martínez, J., Ancheyta, J., 2012. Kinetic model for hydrocracking of heavy oil in a {CSTR} involving short term catalyst deactivation. Fuel 100, 193–199. doi:http://dx.doi.org/10.1016/j.fuel.2012.05.032 Matos, E.M., Guirardello, R., 2000. Modelling and simulation of the hydrocracking of heavy oil fractions. Brazilian J. Chem. Eng. 17, 79–90. doi:10.1590/S010466322000000100007 Mederos, F.S., Ancheyta, J., 2007. Mathematical modeling and simulation of hydrotreating reactors: Cocurrent versus countercurrent operations. Appl. Catal. A Gen. 332, 8–21. doi:10.1016/j.apcata.2007.07.028

18

Moghadassi, A., Amini, N., Fadavi, O., Bahmani, M., 2011. Hydrocracking Lumped Kinetic Model with Catalyst Deactivation in Arak Refinery Hydrocracker Unit. J. Pet. Sci. Technol. 1, 31–37. Mogilicharla, A., Mitra, K., Majumdar, S., 2014. Modeling of propylene polymerization with long chain branching. Chem. Eng. J. 246, 175–183. doi:10.1016/j.cej.2014.02.052 Mohanty, S., 1989. Modelling simulation of hydrocracking and hydrotreating processes. PhD Thesis, Indian Inst. Technol. Kanpur, India. Mohanty, S., Saraf, D.N., Kunzru, D., 1991. Modeling of a hydrocracking reactor. Fuel Process. Technol. Elsevier. Mosby, J.F., Buttke, R.D., Cox, J.A., Nikolaides, C., 1986. Process characterization of expanded-bed reactors in series. Chem. Eng. Sci. 41, 989–995. doi:10.1016/00092509(86)87184-6 Nandasana, A.D., Ray, A.K., Gupta, S.K., 2003. Applications of the Non-Dominated Sorting Genetic Algorithm ( NSGA ) in Chemical Reaction Engineering Applications of the Non-Dominated Sorting Genetic Algorithm ( NSGA ) in Chemical Reaction Engineering. Pacheco, M. a., Dassori, C.G., 2002. Hydrocracking: An improved Kinetic Model and Reactor Modeling. Chem. Eng. Commun. 189, 1684–1704. doi:10.1080/00986440214584 Patel, N., Padhiyar, N., 2016. Multi-Objective Optimal Control Study of fed-batch bioreactor. IFAC-PapersOnLine 49, 97–102. doi:10.1016/j.ifacol.2016.07.223 Qader, S.A., Hill, G.R., 1969. Hydrocracking of gas oil. Ind. Eng. Chem. Process Des. Dev. 8, 98–105. Rangaiah, G.P., 2010. Stochastic global optimization: techniques and applications in chemical engineering. World Scientific. Rangaiah, G.P., 2008. Techniques and Applications in Chemical Engineering (With CDROM). Rangaiah, G.P., Bonilla-Petriciolet, A., 2013. Multi-objective optimization in chemical engineering: developments and applications. John Wiley & Sons. Rapoport, I.B., 1963. Chemistry and technology of synthetic liquid fuels. Reid, Rc., Prausnitz, J.M., 1977. TK Sherwood The properties of gases and liquids. MacGraw-Hill, New York. Sadighi, S., 2013. Modeling a vacuum gas oil hydrocracking reactor using axialdispersion lumped Kinetics. Pet. Coal 55, 156–168. Salmi, T.O., Mikkola, J.-P., Warna, J.P., 2010. Chemical reaction engineering and reactor technology. CRC Press. Sanchez, S., Rodriguez, M.A., Ancheyta, J., 2005. Kinetic mode for moderate hydrocracking of heavyoils. Ind. Eng. Chem. Res. 44, 9409–9413. Schuit, G.C. a, Gates, B.C., 1973. Chemistry and Engineering of Catalytic H yd rodesulf uriza tion. AIChE J. 19, 417–438.

19

Sharma, S., Rangaiah, G.P., 2016. Designing , Retro fi tting , and Revamping Water Networks in Petroleum Re fi neries Using Multiobjective Optimization. doi:10.1021/acs.iecr.5b02840 Sildir, H., Arkun, Y., Cakal, B., Gokce, D., Kuzu, E., 2012. A dynamic non-isothermal model for a hydrocracking reactor: Model development by the method of continuous lumping and application to an industrial unit. J. Process Control 22, 1956–1965. doi:10.1016/j.jprocont.2012.08.016 Stangeland, B., 1974. A kinetic model for the prediction of hydrocracker yields. … Eng. Chem. Process Des. … 13, 71–76. doi:10.1021/i260049a013 Sultana, N., Hajra, B., Saxena, V.K., Guria, C., 2016. Optimal Synthesis of Sal ( Shorea robusta ) Oil Biodiesel Using Recycled Bentonite Nanoclay at High Temperature. doi:10.1021/acs.energyfuels.5b02227 Tsamatsoulis, D., Papayannakos, N., 1998. Investigation of intrinsic hydrodesulphurization kinetics of a VGO in a trickle bed reactor with backmixing effects. Chem. Eng. Sci. 53, 3449–3458. doi:10.1016/S0009-2509(98)00144-4 Valavarasu, G., Bhaskar, M., Sairam, B., Balaraman, K., Balu, K., 2005. A Four Lump Kinetic Model for the Simulation of the Hydrocracking Process. Pet. Sci. Technol. 23, 1323–1332. doi:10.1081/LFT-200038172 Van Geem, K., Reyniers, M.-F., Marin, G., 2008. Challenges of modeling steam cracking of heavy feedstocks. OIL GAS Sci. Technol. L Inst. Fr. DU Pet. 63, 79–94. Vargaftik, N.B., 1975. Handbook of physical properties of liquids and gases-pure substances and mixtures. Weir, H.M., Eaton, G.L., 1932. Heat content of petroleum-oil fractions at elevated temperatures. Ind. Eng. Chem. 24, 211–218. Worldwide Fuel Charter, 5th edition, 2013. Zhou, H., Lu, J., Cao, Z., Shi, J., Pan, M., Li, W., Jiang, Q., 2011. Modeling and optimization of an industrial hydrocracking unit to improve the yield of diesel or kerosene. Fuel 90, 3521–3530. doi:10.1016/j.fuel.2011.02.043

20

Fig. 1: (a) A schematic diagram for two stage VGO hydrocracking unit (Mohanty et al., 1991) (b) the fourth catalyst bed of first-stage is enlarged to highlight non-ideal flow behavior in the reactor

21

Temperature and component concentrations in the feed, i=1

j=1

Constant temperature profile, T0(z) assumed across the ith catalyst bed

Solution of mass balance using finite difference method to obtain concentration profiles, C(j)(z) of ith bed

Solution of energy balance using Runge –Kutta- Gill method to obtain temperature profiles, T(j)(z) of ith bed

j = j+1

No

T(j-1)(z) ≈ T(j)(z) or j ≥ 25

Yes

Final C(j)(z) and T(j)(z) profiles of ith bed i = i+1

No Heat balance between outlet temperature of ith bed and quench hydrogen

i = Total number of beds

Yes

Complete Concentration and Temperature profiles across the reactor Fig. 2: Flow chart of the simulation procedure [superscripts j denote the cycle number and T0(z) is same as inlet temperature, obtained by heat balance between outlet temperature, t(i), of previous bed and hydrogen quenching]

22

Axial Dispersion Model

705

(a)

700 695 690 685 680 675 Plug Flow Model without Hydrotreating 670 0

1

2

4

3

Temperature across second reactor T2, K

Temperature across first reactor T1, K

710

690

(b)(b) Axial Dispersion Model

680

670

660

650 Plug Flow Model without Hydrotreating 640 0

1

Catalyst bed

3

2

Catalyst bed

Fig. 3: Temperature profiles across catalyst beds a) first reactor and b) second reactor [profiles corresponding to tuned axial dispersion model are shown by black lines whereas those corresponding to plug flow model without hydrotreating (Mohanty et al., 1991) are

Concentration across second reactor (mass fraction)

shown by red lines]. Concentration across first reactor (mass fraction)

0.16

(a)

0.14

C(20)

0.12

C(23)

0.10 0.08

C(14)

0.06 C(9)

0.04

C(4)

0.02

C(1)

0.00 0

1

2

Catalyst bed

3

4

0.14

(b)

0.12

C(20)

0.10 0.08 0.06 C(14)

0.04

C(23) C(9)

0.02

C(4)

0.00

C(1) 0

1

2

3

Catalyst bed

Fig. 4: Concentration profile across reactor beds in a) first reactor and b) second reactor for selected components: C1, C4, C9, C14, C20 and C23 representing butane, naphtha, kerosene, diesel, feed and higher components respectively [profiles corresponding to tuned axial dispersion model are shown by black lines whereas those corresponding to plug flow model without hydrotreating (Mohanty et al., 1991) are shown by red lines].

23

0.014

Industrial Simulation Robust Operating Point

0.012 0.010 0.008 0.006 0.004 0.002 0.000 0

1

2

Catalyst bed

3

4

Concentration of nitrogen across first reactor (mass fraction)

Concentration of sulfur across first reactor (mass fraction)

0.016

6e-4 5e-4

Industrial Simulation Robust Operating Point

4e-4 3e-4 2e-4 1e-4 0

0

1

2

3

4

Catalyst bed

Fig. 5: Concentration profiles across the beds for a) sulfur b) nitrogen in the first reactor [profiles corresponding to plant operating point based on industrial simulation are shown by black lines whereas those corresponding to robust operating point of problem 1 are shown by red lines].

24

0.770 0.765

H2*

0.760 0.755 0.750 0.745

Robust Operating Point

HS D*

0.740

0.975 0.980 0.985 0.990

0.735 0.91

0.92

0.995 0.93

0.94

ATF*

1.000 0.95

0.96

0.97

1.005

Fig.6: Pareto optimal plot for three-objective optimization (minimization of H2*, maximization of ATF* and maximization of HSD*). The robust operating point is shown by the red circle is identified at inflation and is further analyzed in detail.

25

(a) 667

666

665 Robust Operating Point

664

663 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80

Inlet temperature in second reactor, T2 (K)

Inlet temperature in first reactor, T1 (K)

668

644 643

(b)

642 Robust Operating Point

641 640 639 638

637 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80

H2*

H2*

660

(c)

Inlet hydrogen in second reactor, HPM2 (m3/m3 feed hr-1)

Inlet hydrogen in first reactor, HPM1 (m3/m3 feed hr-1)

612

650

610

(d)

640

608

630

606

Robust Operating Point

620

604

610

602

600

600 Robust Operating Point

598 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80

H2*

590 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80

H2*

Fig. 7: Plot for decision variables a) inlet temperature of reactor 1, b) inlet temperature of reactor 2, c) hydrogen flow at the inlet of first reactor (HPM1), d) hydrogen flow at the inlet of the second reactor (HPM2) vs. total hydrogen in a process cycle (H2*). The values corresponding to the robust operating point as shown in Fig. 6 are shown by the red circle.

26

Quenching hydrogen Q2(i) after each catalyst bed in second reactor (m3/hr)

Quenching hydrogen Q1(i) after each catalyst bed in first reactor (m3/hr)

44000

(a) 42000

40000

38000

36000

34000 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80

9500

(b) 9000 8500 8000 7500 7000 6500 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80

H2*

H2*

Q1(1) Q2(2)

Q2(1) Q2(2)

Q3(3)

Optimal quench flow rate 1st bed Optimal quench flow rate 2nd bed

Optimal quench flow rate for 1st bed Optimal quench flow rate for 2nd bed

Total quenching hydrogen entering in first reactor (m3/hr)

17500 17000

1.22e+5

(c)

1.20e+5

16500

(d)

1.18e+5 1.16e+5

16000

1.14e+5

15500

1.12e+5

15000

1.10e+5

14500 14000

Total quenching hydrogen entering in second reactor (m3/hr)

Optimal quench flow rate for 3rd bed

Robust Operating Point

1.08e+5 1.06e+5

13500 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80

H2*

Robust Operating Point 1.04e+5 0.73 0.74 0.75 0.76 0.77 0.78 0.79 0.80

H2*

Fig. 8: Plot for decision variables a) quench hydrogen flow rate between each bed of reactor 1, b) quench hydrogen flow rate between each bed of reactor 2, c) total quench hydrogen flow rate in reactor 1, d) total quench hydrogen flow rate in reactor 2 vs. total hydrogen in a process cycle (H2*). The values corresponding to the robust operating point as shown in Fig. 6 are shown by the red circle.

27

(a)

Plant Operating Point

700

690

680

670

Robust Operating Point

660 0

1

2

Catalyst bed

3

4

Temperature across second reactor, T2(K)

Temperature across first reactor,T1 (K)

710

700

(b) 690 680 Robust Operating Point 670 660 650 Plant Operating Point

640 630 0

1

2

3

Catalyst bed

Fig. 9: Temperature profiles across catalyst beds a) first reactor and b) second reactor [profiles corresponding to plant operating point based on industrial simulation are shown by black lines whereas those corresponding to robust operating point of problem 1 are shown by red lines].

28

(a)

0.14

C(20)

0.12 0.10 C(23)

0.08

C(14)

0.06

C(9)

0.04

C(1) 0.02

C(4)

0.00 0

1

2

3

4

Concentration across second reactor (mass fraction)

Concentration across first reactor (mass fraction)

0.16

0.14

(b)

0.12

C(20)

0.10 0.08 0.06 C(14) 0.04

C(23) C(4)

0.02

C(9) C(1)

C(1) C(4)

0.00 0

Catalyst bed

1

2

3

Catalyst bed

Fig. 10: Concentration profile across reactor beds in a) first reactor and b) second reactor for selected components: C1, C4, C9, C14, C20 and C23 representing butane, naphtha, kerosene, diesel, feed and higher components respectively [profiles corresponding to plant operating point based on industrial simulation are shown by black lines whereas those corresponding to robust operating point of problem 1 are shown by red lines].

29

Table 1: Comparison of model predicted results with plant data (Mohanty et al., 1991) Axial Dispersion Model (with Hydrotreatin g)

Simulation Results using Plug Flow Model without Hydrotreatin g (Mohanty et al., 1991)

Plant Data

Error of Axial Dispersion Model (with Hydrotreatin g) with Plant Data (%)

Error of Plug Flow Model without Hydrotreatin g with Plant Data (%)

Total feed to second stage (kg/hr)

183761.4

183236

183385

+0.20

-0.08

Hydrogen consumptio n (kg/hr)

1st Stage

324 3

1st Stag e

2816

1st Stag e

3267

-0.73

-13.8

2nd Stage

123 5

2nd Stag e

1196

2nd Stag e

1363

-9.3

-12.2

Reactor outlet temperatur e (K)

1st Stage

704

1st Stag e

699. 3

1st Stag e

714(upp er limit)

within the given limits

2nd Stage

686

2nd Stag e

677

2nd Stag e

700(upp er limit)

High-speed diesel (wt. %)

48.79

48.79

50.5

-3.3

-3.3

Aviation turbine fuel (wt. %)

31.0

30.53

29.4

+5.4

+3.8

Naphtha (wt. %)

15.6

16.17

15.8

-1.2

+2.3

Butanes and lighter fraction

4.46

4.51

4.5

-0.88

+0.02

(wt. %)

30

Table 2: Values of NSGA parameters

Parameter

Value/Specification

Population size

50

Number of generation

100

Encoding technique

Real Coding

Distribution index for real-coded mutation

20

Crossover Probability

0.9

Distribution index of real-coded crossover

20

Selection strategy

Tournament selection

Mutation probability

1/number of decision variables

31