Subscriber access provided by University of South Dakota
Article
Modeling of catalytic reactor for hydrotreating of straight-run gas oil blended with FCC naphtha and light cycle oil: The influence of vapor-liquid equilibrium Ivana M Mijatovic, Sandra B. Glisic, and Aleksandar Orlovic Ind. Eng. Chem. Res., Just Accepted Manuscript • DOI: 10.1021/ie503188p • Publication Date (Web): 17 Nov 2014 Downloaded from http://pubs.acs.org on November 22, 2014
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
Industrial & Engineering Chemistry Research is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
Modeling of catalytic reactor for hydrotreating of straight-run gas oil blended with FCC naphtha and light cycle oil: The influence of vapor-liquid equilibrium Ivana M. Mijatović, Sandra B. Glisic*, Aleksandar M. Orlović Faculty of Technology and Metallurgy, University of Belgrade, Karnegijeva 4, 11120 Belgrade, Serbia
ACS Paragon Plus Environment
1
Industrial & Engineering Chemistry Research
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 2 of 35
ABSTRACT Reducing the sulfur contents of diesel fuels may require adjusted operation of low pressure hydrotreater units. Mathematical model for co-hydrotreating of straight run gas oil blended with fluid catalytic cracking naphtha and light cycle oil was developed using axial distribution of phase equilibrium and effective wetting in the catalytic reactor. Model assumes that hydrodesulfurization (HDS) and hydrodearomatization (HDA) reactions occur on the catalyst surface which is in contact with vapor or liquid phase. Kinetic equations of HougenWatson type were used to describe HDS reactions for different classes of sulfur compounds. Model results were validated using the industrial test run data and very good predictions of overall sulfur conversion and reactor temperature were obtained. Simulation of reactor operation at different pressures, temperatures and hydrogen purities, confirms that reaction pressures of around 100 bar and high purity hydrogen streams are required for almost complete removal of sulfur compounds.
KEYWORDS Hydrotreating, Catalytic reactor, Middle distillates, Mathematical modeling, Light cycle oil
ACS Paragon Plus Environment
2
Page 3 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
1. INTRODUCTION Current environmental standards and regulations have become increasingly stringent leading to constantly decreasing concentration limits of sulfur in commercial diesel fuels. These global trends in standardization, followed by continuous advances in compression ignition engine technology, have resulted in increasing severity of hydrotreating process, or deep desulfurization of diesel fuels.1 The straight run gas oil (SRGO) which is obtained by crude oil distillation, still remains the dominant feedstock for diesel fuel production. Other possible streams like light cycle oil (LCO), coker gas oil (CGO)1 and visbreaker naphtha2 are also interesting options for refiners. However, these feedstocks are used to much lower extent since LCO and CGO properties like low cetane number as well as the high sulfur, nitrogen and aromatics content, make them difficult to treat and use as diesel blending components.3-5 Improved understanding of hydrotreating has been acquired over the past decades through systematic study of different aspects of this complex process.6,7 Mathematical models of the process have been successfully applied to study, simulate and optimize hydrotreater operation for processing of middle distillates, including LCO.8-13 However, much of the published data has been gathered through studies of synthetic fractions and model compounds in laboratory scale reactors and several pilot plant scale studies.14-18 Due to the complexity of hydrotreating it would be beneficial to develop further understanding of the process through mathematical models which are based on the industrial scale data and real feedstocks. Majority of kinetic models for hydrodesulfurization process are based on overall sulfur conversion.
Among
the
tested
models
the
best
predictions
were
obtained
using
pseudohomogeneous one dimensional reactor model incorporating kinetic model based on two first order parallel reactions.19-21 Anchyeta and coworkers have concluded that more detailed models allow more information to be gathered about the performance of HDS reactor, although some uncertainness may be associated to the kinetics and related model parameters.19,20 Recent advances in the study of hydrotreating reaction pathways and kinetics have resulted in kinetic equations which describe surface catalytic reactions of different classes of sulfur compounds normally found in middle distillates.22-24 Majority of mathematical models applied in the study of hydrotreating process have treated the organic phase, or middle distillate stream, as non-volatile phase.8-20 Using this assumption, the chemical reactions were limited to liquid phase since sulfur containing reactants and hydrogen were considered to be present only in
ACS Paragon Plus Environment
3
Industrial & Engineering Chemistry Research
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 4 of 35
liquid phase.25 Recently published study25 has shown that LCO and other petroleum fractions within similar boiling range can undergo phase change when contacted with hydrogen under elevated pressures and temperatures, which are normally encountered in hydrotreating process. This can be expected to occur even more frequently in industrial operation of hydrotreater reactor, due to the high hydrogen/oil ratios, presence of methane in hydrogen stream, presence of lower boiling point oil fractions, low LHSV values, and also due to the reactor temperature increase which might be accompanied by relatively low operating pressure. Under such circumstances catalytic hydrogenation reactions with reactants present in volatile liquid phase can undergo wet to dry transition when carried out in a trickle-bed reactor. This type of behavior leading to hysteresis loop, was first observed in cyclohexene hydrogenation.26,27 Elevated pressures and temperatures, along with the high hydrogen to hydrocarbon flow ratio, can lead to shift in vapor – liquid equilibrium and dry out locations in the catalyst bed. In those locations of the catalyst bed reaction rates are considerably higher28,29 resulting in a hysteresis profile of reactor performance, with one branch corresponding to dry reactor operation and high reactant conversion, and another branch corresponding to wetted particles throughout the reactor and a low reactant conversion.26-29 This study is dedicated to modeling and simulation of trickle-bed reactor operation for hydrotreating of SRGO blended with fluid catalytic cracking naphtha and light cycle oil (FCC NLCO). A one-dimensional pseudo-homogeneous model was developed using phase distribution obtained by vapor-liquid equilibrium calculations for the system containing hydrogen, methane and SRGO blended with FCC N-LCO. Catalyst bed was divided into dried out section and wetted section of the catalyst bed, as governed by the effective catalyst bed wetting under investigated conditions. Hougen-Watson type kinetic equations for different classes of sulfur compounds were applied in separate material balance equations for vapor and liquid phase. Mass transfer effects were taken into account using overall catalyst effectiveness factor for each of the reactions and overall energy balance equation was included in the model.
2. EXPERIMENTAL SECTION 2.1. Industrial test run data The industrial test run used for model development and validation was performed in industrial hydrotreater with conventional Co-Mo/γ-Al2O3 catalyst at Refinery Pancevo.
ACS Paragon Plus Environment
4
Page 5 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
Commercial Co-Mo/γ-Al2O3 hydrotreating catalyst used during the test run was in the form of trilobe particles (average diameter 0.0015 m and average length 0.007 m), with porosity of 0.52 cm3/g and other properties which can be found in the literature.30 The initial feed into the catalytic reactor was straight run gas oil (SRGO) while FCC naphtha and light cycle oil (FCC NLCO) stream was added to the feed gradually over the testing period. Total test run recorded time was 108 h and the 6 hours periods prior to the run and upon the run completion were used to stabilize the unit. Reactor pressure during the test run was 40 bar and the content of FCC N-LCO in the feed was gradually increased to 20 %vol. along with the inlet temperature increase from 327 to 334 oC (shown in Figure 1).
Figure 1. The data collected during industrial test run (analytical results and values recorded on HDT unit SCADA system) and six characteristic points used for validation of mathematical model. Hydrogen content of the hydrogen stream was 76 – 78 %mol and the remaining part was methane. Industrial catalytic reactor was adiabatic tubular reactor, with two catalyst beds (internal diameter of 2.135 m, the total length of 5.31 m, first bed length of 1.71 m and the
ACS Paragon Plus Environment
5
Industrial & Engineering Chemistry Research
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 6 of 35
second bed length of 3.6 m) separated by quench section (0.3% vol. of the reactor inlet hydrogen stream) and flow re-distributor. Additional description of the industrial test run can be found elsewhere.30 Very high sulfur conversions were achieved and temperature difference in the reactor increased from 4 to 12 oC during the test run (as shown in Figure 1).30 For most samples taken during the test run outlet sulfur concentrations were below 50 ppm. Six characteristic points of the test run were chosen for validation of the developed mathematical model, as shown in Figure 1 (M1 to M6). The GC-MS analysis of input streams30 was performed and concentration distribution of different classes of sulfur compounds is shown in Table 1. Table 1. Properties of input and output streams used in the test run relevant for the model development and validation Property Distillation, oC IBP 5 %vol. 10 %vol. 20 %vol. 50 %vol. 70 %vol. 90 %vol. 95 %vol. FBP Sulfur content, % wt. Sulfur content, mg/kg Thiol/mercaptan sulfur, %wt. Paraffin and naphtenic content, % vol. Olefins content, % vol. Aromatics content, % vol. Mono-aromatic content, % wt. Di-aromatic content, % wt. Tri+aromatic content, % wt. Sulfur distribution, wt% C1 and C2-BT C3 and C4-BT DBT and NT C1-DBT and C1-NT C2-DBT and C2-NT C3-DBT
Hydrotreated gas oil
SRGO
FCC N-LCO
208 229 236 245 265 281 304 315 326 / 98% 0.7305
33 63 101 170 242 274 333
0.0162 68.7 5.2 26.1 19.4 11.5 1.41
0.0271 32.1 16.9 51.0 Not available Not available Not available
37 0 73.0 0 27.0 27.4 4.6 0.70
12.73 60.68 8.39 8.24 9.56 0.39
24.92 45.73 4.66 11.04 10.79 2.85
1.30 5.40 0.50 14.33 76.40 2.07
335 / 91% 1.1400
ACS Paragon Plus Environment
196 223 231 241 259 275 294 302 321 / 98%
6
Page 7 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
The data concerning sulfur distribution is crucial for the model development and validation since the material balance equations were developed using specific rate equations for hydrogenation and hydrogenolysis of different classes of sulfur compounds. Classes of sulfur compounds used in the model were: substituted benzothiophenes (C1 to C4 BTs), dibenzothiophene and naphtothiophene (DBT and NT) and substituted DBTs and NTs (C1 to C3 DBT and NT).30 Other relevant properties of the inlet and outlet streams are shown in Table 1: distillation curve, sulfur contents and contents of aromatic compounds. These properties were measured on stream corresponding to blend M6 of the test run, shown in Figure 1. Properties of input streams for other characteristic points (M1-M5), like aromatic content, distillation range and sulfur distribution, were then calculated using data shown in Figure 1 and Table 1.
3. REACTOR MODEL
Catalytic hydrotreater was modeled using deterministic one-dimensional model of gasliquid-solid catalytic reactor under steady operation, or trickle-bed reactor (TBR). The following assumptions were used: • Trickle flow of reacting fluids; • The main reactions occurring in the reactor are: hydrodesulfurisation (HDS) and hydrodearomatisation (HDA); • Axial dispersion of the fluid flow is negligible31,32; • The reactor operates under adiabatic and isobaric conditions (as confirmed by the measured pressure drop); • There is no catalyst deactivation during the test run; • Catalyst particles are of uniform size and shape. Models of hydrotreating reactors are usually developed using the assumption that mass transfer of hydrogen to the non-volatile liquid phase of middle distillates takes place in the reactor and that reactions occur on the catalyst surface which is in contact with liquid phase only, i.e. the catalyst bed is completely wetted. This assumption is valid in some cases, but such type of hydrotreater model makes no account of the hydrocarbons and sulfur compounds which might be transferred to the vapor phase.
ACS Paragon Plus Environment
7
Industrial & Engineering Chemistry Research
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 8 of 35
Hydrogenation of liquids which contain volatile components was investigated in the past and interesting reactor behavior was reported.26-29 Following those experimental findings mathematical models which take into consideration evaporation of liquid phase were developed. This type of approach treats the fixed bed of catalyst as reaction domain in contact with reactants present in both phases. Such models were able to predict hysteresis type of behavior for tricklebed reactors quite successfully.28-29 The investigated reactions were hydrogenation of cyclohexene and hydrogenation of benzothiophene28-29, which are obviously susceptible to the evaporation of reactants and their migration to the gas phase under reaction conditions. Hydrotreating of middle distillates is conducted at different pressure and temperature conditions than hydrogenation of light organic compounds. Temperatures are usually around 300 o
C and higher, and pressures are normally 40 bar and above. Under these conditions it is hardly
expected to encounter considerable evaporation of hydrocarbons and sulfur compounds. However, hydrogen to middle distillates flow ratio can be very high and also hydrogen stream can contain relatively large fraction of methane and other lower hydrocarbons, especially in industrial units. These conditions can lead to quite different vapor-liquid equilibrium in the reactor than the expected one. This phenomenon can result in a situation where sulfur compounds react on the catalyst surface which is partially wetted and therefore in contact with reactants in both vapor and liquid phase, similar to the situation observed in hydrogenation of cyclohexene. Such behavior of reaction mixture was assumed and considered in this hydrotreater model and the following additional assumptions and procedures were introduced: Chemical reactions are occurring in the liquid and vapor phase. The reaction rate of any chemical species in vapor or liquid phase corresponds to concentration of that particular chemical compound in vapor or liquid phase; Reaction mixture vapor-liquid equilibrium (VLE) is approximated by the phase equilibrium of a four-, five-, or seven- component system, meaning that middle distillates stream is approximated with: n-dodecane and n-octadecane (two pseudo-components); n-dodecane, nhexadecane and n-octadecane (three pseudo-components); and n-dodecane, n-hexadecane, noctadecane, tetraline and phenantrene (five pseudo-components including aromatic compounds which were included in the model due to the presence of highly aromatic LCO). Two additional components used for VLE calculations for the overall mixture are methane and hydrogen;
ACS Paragon Plus Environment
8
Page 9 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
Catalyst bed wetting determines which part of the catalyst is in contact with reactants in vapor or liquid phase and bed wetting is calculated using equation for effective bed wetting in TBRs. The equation uses pressure drop across the catalyst bed which was measured during the industrial test run and applied in the model calculations; For every segment of the catalyst bed a new phase equilibrium is calculated using flash calculation based on Peng-Robinson Equation of State (PR EOS) for the above specified pseudocomponents, and Mass transfer resistances are taken into account by the calculation of overall effectiveness factor for all reactions in both phases. Thermodynamic properties and components distribution were calculated using UniSim software. Simulations were carried out in a user friendly program developed in MATLAB® & Simulink® Release 2010b. The system of ordinary differential equations was solved using a fourth order Runge – Kutta routine with variable step size. As a result of simulations, molar flows (or concentrations) of individual components in liquid and vapor phase were obtained, as well as reactor temperature, as a function of catalyst weight along the length of catalytic reactor.
3.1. Thermodynamic properties and VLE calculations Reaction mixture entering hydrotreater (HDT) reactor had four components: hydrogen (H2), methane (CH4), straight run gas oil (SRGO) and FCC naphtha and light cycle oil stream (FCC N–LCO). The distillation curves for both SRGO and FCC N-LCO are shown in Figure 2a. Since middle distillate stream entering the reactor during the test run was a blend of SRGO and FCC N-LCO with varying proportions of these two streams, six characteristic input points corresponding to mixtures M1-M6 in Fig.1 were analyzed and their calculated distillation curves are shown in Figure 2b.
ACS Paragon Plus Environment
9
Industrial & Engineering Chemistry Research
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 10 of 35
Figure 2. The distillation curves for a) neat SRGO and FCC N-LCO streams and b) six characteristic input mixtures M1-M6: M1- only SRGO, M2 – 91 vol% SRGO and 9 vol% FCC N-LCO, M3 – 91 vol% SRGO and 9 vol% FCC N-LCO, M4 – 87 vol% SRGO and 13 vol% FCC N-LCO, M5 and M6 – 80 vol% SRGO and 20 vol% FCC N-LCO
The number of components used to approximate SRGO and FCC N-LCO blends can have significant impact on simulation results. As can be seen from Figure 2a, SRGO and FCC N-LCO can be approximated with C12 to C18 hydrocarbons (n-dodecane, n-hexadecane and noctadecane). Another possibility which accounts for high aromatic content of LCO and which should be more accurate in describing distillation curve of different blends M1-M6 is shown in Figure 2b. The influence of aromatic compounds on the quality of simulation results is well known and it was previously reported in the literature.25 The approximation which uses aromatic pseudo-components, tetraline and phenantrene, was verified by the simulation of VLE for blends previously investigated in the literature.25 The difference between calculated and experimental weight percentage of evaporated LCO and white oil blend was less than 2%. Table 2 summarizes all compositions during the test run for each of the input mixtures M1-M6. Blends M1-M6 were approximated with two pseudo-components (n-dodecane and n-octadecane) as shown in the upper part of Table 2, as well as using three pseudo-components (n-dodecane, nhexadecane and n-octadecane) and five pseudo-components (n-dodecane, n-hexadecane and noctadecane, tetraline and phenantrene) as shown in the following part of Table 2. VLE and distribution of different classes of sulfur compounds corresponding to the selected pseudocomponents are also shown in Table 2 (all BTs are summarized as BT).
ACS Paragon Plus Environment
10
Page 11 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
Table 2. The characteristic inlet blends and components used for their approximation Input streams Tinput Toutput Sulfur conversion Input stream composition
Gas + M1 Gas + M2 Gas + M3 Gas + M4 Gas + M5 Gas + M6 K K mol% approx. compound
600.0 604.1 99.27
601.0 606.1 99.23
602.0 609.3 99.37
602.0 609.7 99.41
605.0 616.5 99.26
606.0 617.5 99.49
Mole fraction
Approximation of inlet diesel stream with two pseudo-components Hydrogen
H2
0.7215
0.7155
0.7120
0.7101
0.7096
0.7096
Methane
CH4
0.2155
0.2137
0.2126
0.2120
0.2119
0.2119
BT
C12
0.0046
0.0052
0.0060
0.0057
0.0072
0.0064
DBT1
C18
0.0004
0.0005
0.0006
0.0006
0.0007
0.0006
DBT2
C18
0.0006
0.0006
0.0007
0.0007
0.0009
0.0008
DBT3 non-sulfur light fraction non-sulfur heavy fraction Diaromatics
C18
0.0005
0.0006
0.0007
0.0007
0.0009
0.0008
0.0080
0.0113
0.0128
0.0142
0.0132
0.0137
0.0366
0.0368
0.0291
0.0295
0.0287
0.0292
C18
0.0019
0.0020
0.0074
0.0076
0.0075
0.0075
Monoaromatics
C12
0.0104
0.0138
0.0181
0.0189
0.0196
0.0196
C12 C18
Approximation of inlet diesel stream with three (five*) pseudo-components Hydrogen
H2
0.7215
0.7155
0.7120
0.7101
0.7096
0.7096
Methane
CH4
0.2155
0.2137
0.2126
0.2120
0.2119
0.2119
BT
C12
0.0046
0.0052
0.0060
0.0057
0.0072
0.0064
DBT1
C18
0.0004
0.0005
0.0006
0.0006
0.0007
0.0006
DBT2
C18
0.0006
0.0006
0.0007
0.0007
0.0009
0.0008
DBT3 non-sulfur light fraction non-sulfur mid fraction non-sulfur heavy fraction Diaromatics A1
C18
0.0005
0.0006
0.0007
0.0007
0.0009
0.0008
0.0080
0.0113
0.0128
0.0142
0.0132
0.0137
0.0337
0.0327
0.0291
0.0295
0.0287
0.0292
0.0029
0.0041
0.0000
0.0000
0.0000
0.0000
0.0019
0.0020
0.0074
0.0076
0.0075
0.0075
Monoaromatics A2 C16 (mono-A*) 0.0104 0.0138 0.0181 *-in the case when mixture is approximated with five compounds
0.0189
0.0196
0.0196
C12 C16 C18 C18 (di-A*)
The simulation of phase equilibrium was performed using UniSim software, PengRobinson equation of state (PR EOS) and physical and thermodynamic parameters based on DIPPR data. A slightly better performance of PR EOS than Soave-Redlich-Kwong EOS around critical point makes the PR EOS somewhat better suited to systems investigated in this study. Also, for a system consisting of a mixture of hydrogen and hydrocarbons, the description of
ACS Paragon Plus Environment
11
Industrial & Engineering Chemistry Research
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 12 of 35
vapor-liquid equilibrium can be approximated well using PR EOS and quadratic mixing rule.33 Peng-Robinson binary interaction parameters were obtained from DECHEMA Peng-Robinson Parameters - Interaction Parameters Data (IPD) using ChemSep – modeling separation processes software. Binary interactions parameters for all thermodynamic components are listed in Table 3. Table 3. Binary interaction parameters for all components used in VLE calculations lji=lij
H2
CH4
n-C12H26
n- C16H32
n-C18H38
tetraline
pyridine
H2
-
0.0026
0.0875
0.0895
0.0917
0.0032
0.0022
CH4
0.0026
-
0.0598
0.0625
0.0657
0.0745
0.0648
n-C12H26
0.0875
0.0598
-
0.1534
0.1705
0.0032
0.0223
n- C16H32
0.0895
0.0625
0.1534
-
0.18526
0.00123
0.00098
n-C18H38
0.0917
0.0657
0.1705
0.18526
-
0.00027
0.00013
tetraline
0.0032
0.0745
0.0032
0.00123
0.00027
-
0.0053
phenantrene
0.0022
0.0648
0.0023
0.00098
0.00013
0.0053
-
Distribution of components in the gas and liquid phase along the reactor length is calculated using equilibrium constants and Flash calculation. Ki =
yi xi
(1)
Equilibrium constants and liquid to vapor flow ratios were generated as polynomial functions dependent on the temperature (obtained using Flash calculation and PR EOS). These polynomials were used in the reactor model to calculate vapor-liquid distribution along the reactor length taking into account axial temperature gradient. Total length of the reactor is divided into ten equal segments and distribution of components in both phases along the reactor is calculated for every segment of the reactor where each segment represents 1800 kg of the catalyst. Ki = A + B ⋅T + C ⋅T 2 + D ⋅T 3
(2)
FL = a + b ⋅T + c ⋅T 2 V F
(3)
Densities of the vapor and liquid phase were calculated using following parametric equations as temperature dependent properties (Eq.4) and heat capacities were also estimated as functions of temperature (Eq.5).
ACS Paragon Plus Environment
12
Page 13 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
ρ i = Ai + Bi ⋅ T ; i = V , L
(4)
C p = α + β ⋅T + γ ⋅T 2 + δ ⋅T 3
(5)
3.2. Catalyst wetting efficiency Catalyst wetting efficiency (CWE) in trickle bed reactor can be predicted as a function of operating conditions. In trickle flow regime the catalyst pellets can be incompletely wetted and different approaches and correlations to estimate CWE are available17. The correlation shown in Eq.6 was developed using tracer technique34 in non-reactive system and was used in this TBR model: 1 + [(∆p / ∆z ) / ρ L g ] f = 1.104 ⋅ Re ⋅ Ga L
1/ 9
1/ 3 L
(6)
Catalyst wetting efficiency depends on the reactor design, catalyst shape and size, but mostly and essentially on the liquid flow rate and velocity in a catalyst bed. Re and Ga numbers for liquid phase were calculated using data from the test run and calculated flows of vapor and liquid phase resulting from the VLE calculations. The measured pressure drop across the catalyst bed obtained during the test run was also used in calculations. Other correlations and relationships17 for determination of CWE were also tested but the overall predictions of sulfur conversion and temperature difference in the TBR using the developed model were considerably worse than when using Eq. 6.
3.3. Material and energy balance equations Material balance equations for all reacting components involved in the process of hydrodesulfurization (BT, DBT1, DBT2 and DBT3), and in the process of hydrodearomatization (diaromatic and monoaromatic), as well as for H2 and H2S, were developed for both phases (vapor and liquid). Material balance equations for component i, where i=BT, DBT1, DBT2 and DBT3, in the vapor phase (V) are given by: −
dFiV = (1 − f ) ⋅ ΩVi ⋅ ri , HDS dW
(7)
Material balance equations for component i, where i=BT, DBT1, DBT2 and DBT3, in the liquid phase (L) are given by:
ACS Paragon Plus Environment
13
Industrial & Engineering Chemistry Research
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 14 of 35
L
dF L − i = f ⋅ Ωi ⋅ ri ,HDS dW
(8)
Material balance equations for aromatic compounds A1 (diaromatic) and A2 (monoaromatic) in vapor and liquid phase are given by: V
V (1 − f ) ⋅ Ω A1 ⋅ rA1, ARM dFA1 = dW ρ p ⋅ (1 − ε p )
(9)
L f ⋅ Ω A1 ⋅ rA1, ARM dFA1 = ρ p ⋅ (1 − ε p ) dW L
(10)
for reversible reaction of diaromatics hydrogenation. Balance equations for hydrogen and hydrogen sulfide are the sum of balance equations for HDS and HDA reactions multiplied by the corresponding stoichiometric coefficients. Since hydrogenation of diaromatics to monoaromatics is known to be the slowest reaction sequence in hydrogenation of condensed aromatic compounds and monoaromatics, this reaction was used to represent the overall HDA reaction. Consequently, averaged or “lumped” stoichiometric coefficient for the overall HDA reaction has been used in balance equation for hydrogen. This is an important aspect of hydrotreating process when processing cracked unsaturated streams like LCO since HDA reaction can substantially increase hydrogen consumption. Stoichiometric coefficients for hydrogen were averaged in all individual reactions in accordance with literature data35 and the values used in the model were: 3.6 for HDS of BT, DBT1, DBT2, and DBT3, while value of 6.3 was used for overall reaction of HDA. The overall energy balance shown in Eq.13 uses enthalpy changes of components involved in HDS and HDA reactions (Eqs.11 and 12). Heats of vaporization of organic compounds are not included in the energy balance since pressure conditions in the reactor are above critical values for these compounds, meaning that the actual values of latent heats are approaching zero. Reaction enthalpies of exothermic HDS and HDA reactions used in the model were -40 kJ/kmol and -85 kJ/kmol, respectively. i
i
1
1
i
i
1
1
∆E HDS = ∑ FiV ⋅ (1 − f ) ⋅ ΩVi ⋅ riV, HDS ⋅ ∆H HDS + ∑ Fi L ⋅ f ⋅ Ω iL ⋅ ri ,LHDS ⋅ ∆H HDS ∆E HDA = ∑ F jV ⋅ (1 − f ) ⋅ ΩVj ⋅ r jV, HDA ⋅ ∆H HDA + ∑ F jL ⋅ f ⋅ Ω Lj ⋅ r jL, HDA ⋅ ∆H HDA
ACS Paragon Plus Environment
(11) (12)
14
Page 15 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
dT = dW
∆E HDS + ∆E HDA ∑ Fi ⋅ Cpi,V + ∑ Fi L ⋅ Cpi ,L V
i
(13)
i
3.4. Kinetic equations for HDS and HDA reactions Hydrodesulfurization (HDS) is an exothermic and irreversible reaction under the investigated conditions and a number of kinetic studies dedicated to HDS can be found in the literature.22-24,
36-39
Sulfur compounds present in middle distillates are mainly thiophenes,
benzothiophenes, dibenzothiophenes and their alkyl derivatives and HDS reactivity depends strongly on the structure of sulfur compounds. Under typical hydrotreating conditions, thiophenes and benzothiophenes are highly reactive when compared with the DBTs species.36 The HDS rates of thiophene compounds, especially DBT compounds, are also affected by the presence of substituent groups on the positions that are adjacent to the S atom.39 The DBT compounds become less reactive when the number of substituent groups increases. The highly substituted DBTs, such as 4,6-dimethyldibenzothiophene (4,6DMDBT), are the least reactive group, while the non-substituted DBTs exhibit the highest reactivity.36-39 Two alternative reaction pathways are: thiophenic compounds desulfurized either by saturating the thiophene ring, as well the fused benzene rings, followed by removal of the sulfur atom (hydrogenation pathway proposed for highly substituted DBTs), and direct removal of the sulfur atom without ring hydrogenation (hydrogenolysis pathway constitutes the major route for the HDS of Ts, BTs and less substituted DBTs).37 The hydrogenation pathway is affected by thermodynamics because the hydrogenation reactions are thermodynamically limited under hydrotreating conditions.35-38 The 4-methyldibenzothiophenes (4-Me-DBT) are transformed by both pathways and no pathway is preferred over the other.37-39 Kinetic equations and parameters for HDS and HDA reactions were taken from the literature.22-24 Rate equations of the Hougen-Watson type for all reactions in the HDS reaction network were developed on a commercial CoMo/Al2O3 catalyst under operating conditions significant to industrial applications.22 Rate equations were also specifically developed for the HDS of 4-Me-DBT and 4,6-DiMeDBT for the same catalyst.23 The reaction rate equations for hydrodesulfurization of benzothiophene (BT), 4,6dimethyldibenzothiophene (DBT2) and trimethylbenzothiophenes (DBT3) in the vapor and in the liquid phase can be represented by the following equation (Eq.14):
ACS Paragon Plus Environment
15
Industrial & Engineering Chemistry Research
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
ri
L ,V
= Ci
L ,V
⋅ CH 2
L ,V
k K ⋅ K k K ⋅K ⋅ i ,σ i ,σ L,VH ,σ + i ,τ i ,τ L,VH ,τ DENτ DENσ
Page 16 of 35
(14)
The reaction rate equation for hydrodesulfurization of 4-dimethyldibenzothiophene (DBT1) is given by: rDBT1
L ,V
= C DBT1
L ,V
⋅ CH2
L ,V
k DBT1 ,σ K DBT1 ,σ ⋅ K H ,σ k DBT1 ,τ K DBT1 ,τ ⋅ (k1DBT1 ,τ + k 2 DBT2 ,τ ) ⋅ + L ,V L ,V DENσ DENτ
(15)
Temperature dependencies of the kinetic and adsorption coefficients can be found in the literature.
22-24
Additionally, adsorption parameter DENσ is influenced by the overall sulfur
conversion since the quantity of S containing molecules decreases with the extent of reaction and their products are less likely to be strongly adsorbed on the catalyst surface.12 This dependence was included in the model as proposed in the literature.12 Another important reaction during hydrotreating of gas oils is hydrogenation of benzene rings in aromatic molecules. Hydrogenation of aromatic rings is an exothermic and reversible reaction which takes place on the active sites of the catalyst and this reaction is thermodynamically limited under typical hydrotreating conditions, so, the reverse reaction (dehydrogenation of cyclohexane rings) has to be considered. Numerous studies on the HDA reaction have been published6,7,39-42 and the following conclusions can be summarized: no partially hydrogenated ring compounds were observed (hydrogenation occurs in a ring-by-ring manner) and rate of reaction increases with the number of aromatic rings. For the tri-aromatic and poly-aromatic species, hydrogenation of fused aromatic rings is preferred to hydrogenation of the external rings (central ring of anthracenic molecule is more reactive than that of phenanthrenic molecules) and the presence of alkyl chains and of cyclohexane rings increases the hydrogenation rate.37-40 Since the HDA is a monomolecular reaction, de Oliveira and coworkers37 used Quantitative Structure/Reactivity Correlations (QS/RCs) to calculate rate constants of the HDA for different types of molecules, observing that number of aromatic rings and saturations do not influence significantly the rate of the reaction. Rate equation of aromatics hydrogenation, used in this paper, was taken from the literature8 in the form of the Langmuir-Hinshelwood kinetics for single-site mechanism:
− rA =
k A K A K H pH2 C A
(16)
1 + K AC A + K H 2S p H 2S
ACS Paragon Plus Environment
16
Page 17 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
Values of apparent activation energy, heats of adsorption, expressions for kinetic parameters and their temperature dependences can be found in previously published work.8
3.5. Calculation of the overall catalyst effectiveness for vapor and liquid phase The influence of mass transfer resistances on the overall reaction rate for these catalytic reactions was incorporated into the model through overall catalyst effectiveness.43 The effectiveness factor of independent reactions is defined as the ratio of the volumetric average of the reaction rate into the particle to the reaction rate at the surface of the particle as proposed by Thiele (1939) and Zeldovich (1939).44,45 Analytical solutions for those equations are possible only for single reactions and for zero - and first - order rate expressions. Various correlations used to estimate the catalyst effectiveness factor for isothermal and irreversible reactions are available in the literature. For kinetic models other than the power - law expressions, such as Langmuir – Hinshelwood – Hougen – Watson (LHHW) - type kinetic expressions, there is no analytical solution of equation for effectiveness factor. An alternative method can be used to avoid the numerical integration, the Bischoff generalized modulus approach46, which enables an analytical solution for any type of rate equation and single reaction. Namely, the effectiveness factor in commercial HDS catalysts has been reported to be in the range 0.4 to 0.8.47-57 Bischoff (1965) proposed a general modulus to predict the effectiveness factor for any reaction type within a relatively narrow region.46 If reactions of order less than one - half are excluded, the spread between all different curves is around 15%. The mean deviation of values for η calculated from the empirical correlation proposed by Papayannakos and Georgiou (1988)55 is less than 2.4%, from those predicted with the normalized modulus for simple order reactions proposed by Froment and Bischoff (1990)56 in the range 0.05 < η < 0.99. Regarding the analyzed reaction system the following assumptions were applied: all reactions are occurring with large excess of H2, above stoichiometric, in liquid and vapour phase; and the amount of liquid phase is very small, 0.4 mol% at the reactor inlet and decreasing along the reactor length. So, the overall effectiveness factor for multiple reactions in HDS reactor can be calculated using following equation, as proposed by Aris, Eq.17 57-59
∑X Φ η = ∑X Φ f
Ω
f
f2 i
i
f
i
f i
(17)
f2 i
ACS Paragon Plus Environment
17
Industrial & Engineering Chemistry Research
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 18 of 35
where X
f
i
=
Defff ,i C i f
∑D
f eff ,i
Cif
,
(18)
with i=BT, DBT1, DBT2, DBT3 and A1, f=phase (liquid or vapor). The internal effectiveness factor η can be calculated by Eq.19:
ηif =
tanh(Φ if ) , Φ if
(19)
using generalized Thiele modulus for pseudo first order developed by Aris, Eq.2057-59
V Φ = p Sp f i
f kiI, ,APP Ci f ρ Defff , i
(20)
ki,APP uses kinetic expressions given by Eqs. 14-16. The overall and internal effectiveness factors calculation, for the liquid phase, follows essentially the same procedure as for the gas phase. Effective diffusivity for component i is calculated using following expression: V Deff ,i =
1
(21)
1 1 + V V DAB ,i DK ,i
where molecular diffusivity DAB,iV is calculated by the well known Fuller, Schettler and Giddings correlation60 which accounts for the molecular structure of diffusing molecule, pressure and temperature, while DK,iV represents Knudsen diffusivity given by:
DKV ,i = 2.056 ⋅10 −7 ⋅
T Mi
(22)
4. RESULTS AND DISCUSSION
4.1. Model validation The developed model was used to simulate the industrial trickle bed reactor operation during co-hydrogenation of SRGO and FCC N-LCO. The industrial test run data shown in
Figure 1 were used to validate model results for six characteristic points during reactor operation (M1-M6). The simulation results obtained with the proposed model were compared to the results of industrial test run (Figure 3) together with the influence of the number of pseudo-components
ACS Paragon Plus Environment
18
Page 19 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
on the simulation results. As can be seen the model results for overall sulfur conversion (Figure
3a) and reactor outlet temperature (Figure 3b) correlate well with the experimental data.
Figure 3. Model simulation results versus experimental data for a) overall sulfur conversion and b) temperature increase in reactor.
Total sulfur conversion deviates from the experimental data by up to 0.3%, while for most simulations the difference is around 0.1% (Table 4). Reactor outlet temperatures are predicted with high accuracy by the model using five pseudo-components while predictions with the models using two and three pseudo-components are less accurate, especially in case of two pseudo-components (Table 4).
Table 4. Deviation from the experimental values Time, h
Deviation from the experimental values, %
6
Sulfur conversion, % Two Three component component 0.16 0.13
Five component 0.15
∆T in reactor, K Two Three component component 14.63 14.63
Five component 17.07
18
0.32
0.17
0.17
0.00
3.92
11.76
30
0.15
0.01
0.17
21.92
9.59
8.22
38
0.06
0.02
0.11
12.99
1.30
1.30
ACS Paragon Plus Environment
19
Industrial & Engineering Chemistry Research
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 20 of 35
78
0.14
0.14
0.11
12.38
4.76
0.00
114
0.06
0.03
0.07
12.17
4.35
0.87
Reaction mechanism and kinetics of HDS reactions are very complex12,22-24 and different parameters can affect the overall reaction rate. The results shown in Figure 3 were obtained using detailed kinetic equations for individual classes of sulfur compounds13,22,23 and this type of approach to kinetics of HDS reactions requires the use of numerous rate parameters. Some of these parameters are difficult to generalize and quantify for broad range of possible applications as they can be influenced by the specific properties of the catalyst or reacting fluids. Dedicated investigations of these complex surface catalytic reactions showed that feed aromaticity61 and overall sulfur conversion12 can exhibit quantitative influence on the reaction rate due to the presence of numerous adsorbing species and their concentrations within the reactor. Having this in mind, it can be concluded that the agreement between simulation results and experimental data is very good. As shown in Figure 3, the best simulation results for overall sulfur conversion were obtained with the model which uses five pseudo-components. Both reactions taking place in the hydrotreater reactor are exothermic and the temperature increase is the result of conversion of sulfur containing and aromatic compounds. Individual conversions of different classes of sulfur compounds and aromatic compounds for M6 data set and the model with five pseudo-components are shown in Figure 4.
ACS Paragon Plus Environment
20
Page 21 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
Figure 4. Model simulated reactor profiles for the model with five pseudo-components and conditions M6: Conversions of sulfur compounds (BT, DBT1, DBT2, DBT3), aromatics (A1) and overall sulfur.
These conversion profiles indicate that the highest conversion of total sulfur (and consequently released reaction heat) is located in the entry section of the catalyst bed, which is typical for adiabatic TBR reactor and exothermic reaction. Conversion of DBT2, which is known to be the most refractory sulfur compound, steadily increases along the catalyst bed (final conversion reaches 93.6%) followed by the conversion of aromatics (conversion of aromatics reaches 71.7% while averaged experimental conversion of aromatics was around 67%). Hydrogen consumption calculated by the model for M6 was 291.4 std ft3/bbl, which after addition of hydrogen dissolved in the product35 (up to 6% of total hydrogen consumption) yields 308.9 std ft3/bbl and is in quite good agreement with previous studies35.
ACS Paragon Plus Environment
21
Industrial & Engineering Chemistry Research
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 22 of 35
Figure 5. Catalyst wetting efficiency along the catalyst bed: a) simulation of test run for M4; b) simulation of test run for M6; c) twofold and threefold increase of LHSV for conditions M4; and d) twofold and threefold increase of LHSV for conditions M6.
The data shown in Figure 5 a) and b) provide an indication of catalyst bed wetting efficiency in the reactor for test points M4 and M6, as calculated by the model (governing equations are PR EOS and Eq.6 for CWE). This wetting observation is in accordance with the data published by Henry and Gilbert (1973)62 and Sie and Krishna (1998)63-65 which state that commercial - scale reactors, operating at Reynolds number significantly less than 10 will produce poor wetting. The simulation which uses thermodynamic model with two pseudocomponents results in the existence of considerable quantity of liquid phase and therefore substantial fractional wetting of the catalyst bed across the entire reactor length. The thermodynamic model with three pseudo-components will always result in a completely dry bed of catalyst and reactions occurring only in vapor phase. Model with five pseudo-components predicts partially wetted bed and disappearance of the liquid phase within the reactor. Model results of CWE may appear unrealistically low, especially for the case of thermodynamic representation of middle distillates stream using three pseudo-components. Since the industrial test run which is used to validate the model results was performed at low LHSV (1.05 to 1.32 h1
), at high overall gas (hydrogen and methane) to liquid ratio (around 4300 std ft3/bbl) and at
high hydrogen to hydrocarbons ratio (around 3350 std ft3/bbl), the results of additional model simulations with increased LHSV are shown in Figure 5 c) and d) (simulation results were obtained only for model using five pseudo-components for VLE calculation). These results show CWE values for reactor operation which is more likely to be encountered in a typical industrial and laboratory scale operation. As can be seen, the increase of LHSV from 1.24 to 2.48 h-1 (overall inlet gas/liquid reduces to 2150 std ft3/bbl and inlet hydrogen/liquid reduces to 1675 std ft3/bbl) and to 3.72 (overall inlet gas/liquid reduces to 1230 std ft3/bbl and inlet hydrogen/liquid reduces to 1110 std ft3/bbl), will lead to high (twofold increase of LHSV) and almost complete (threefold increase of LHSV) catalyst wetting for other conditions unchanged and corresponding to M4 (Figure 5 c) and M6 (Figure 5 d). Simulation results for higher LHSV values indicate that typical industrial operation will be performed under the conditions of high catalyst wetting
ACS Paragon Plus Environment
22
Page 23 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
efficiency which is an assumption often used when modeling hydrotreater reactors8-20. Also, a shift from wet to dry reactor operation predicted by the model and Eq.6 should occur in case of mutually coupled low LHSV values, high inlet gas to liquid ratios and high inlet temperatures.
Values of η, internal effectiveness factor calculated by Eq.19 along with generalized Thiele modulus (Eq.20), are in agreement with previously published data.66 The impact of reaction mixture phase distribution on conversion is further influenced by the mass transfer calculated by the overall effectiveness factors (shown in Figure 6). Overall effectiveness factor has low values in liquid phase (approaching zero - Figures 6a), while in the gas phase it decreases along the reactor due to the temperature increase, the increased reaction rate and consequently growing influence of mass transfer resistance (Figures 6b).
Figure 6. Overall effectiveness factor in five component system for: a) liquid phase, b) vapor phase
As shown above, hydrotreating of SRGO and FCC N-LCO blends can be described better using VLE calculations with more pseudo-components. Model results indicate that five pseudocomponents can describe the system quite well, with good prediction of sulfur and aromatics conversion as well as reactor temperature rise. The prediction of most important parameters of the hydrotreating process, using model which accounts for the phase change of middle distillate stream, is very good. In order to be able to predict accurately the conversion of overall sulfur it is
ACS Paragon Plus Environment
23
Industrial & Engineering Chemistry Research
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 24 of 35
necessary to include distribution of sulfur compounds into different classes which are characterized by different reactivity. This type of model can be used to simulate hydrotreating process across a wide range of operating conditions.
4.2. Analysis of the influence of process parameters using model simulation The influence of relevant process parameters: temperature, pressure, hydrogen purity and different inlet concentrations of sulfur, was simulated using model with five thermodynamic pseudo-components. Increasing inlet temperature, at 40 bar reaction pressure and using 77 % mol. hydrogen in the hydrogen stream, will result in higher reactor temperature difference, higher aromatics conversion and practically the same level of sulfur conversion (Figure 7). Although much of the inlet feedstock is transferred into gas phase where reaction rates are higher significant gains in conversion of overall sulfur cannot be achieved by increasing temperature at 40 bar.
Figure 7. The influence of inlet temperature on the conversion of sulfur compounds, conversion of aromatics and outlet reactor temperature
Another strategy which can be used to bring outlet sulfur to very low concentrations, below 10 or 15 ppm, is to increase pressure and/or to increase hydrogen stream purity. Model simulations using high purity hydrogen (100% mol.) and higher pressures of 60 and 100 bar are shown in Table 5 and Figures 8. Increasing hydrogen purity and reactor pressure can result in increased sulfur conversion and outlet concentrations within 10 ppm range. Even when inlet
ACS Paragon Plus Environment
24
Page 25 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
sulfur concentrations are very high, the resulting outlet sulfur will be below or slightly above 10 ppm provided that reaction pressure is around 100 bar and that hydrogen stream contains low concentrations of lower hydrocarbons.1,55,67,68
Table 5. Model results for outlet sulfur concentration using different inlet temperatures, reaction pressures, hydrogen stream purities and inlet sulfur concentrations Tin, K
P, bar
Hydrogen purity 77%
Hydrogen purity 100%
Sinlet=7290 ppm
Sinlet=14580 ppm Outlet Sulfur concentration, ppm
Sinlet=7290 ppm
Sinlet=14580 ppm
606
40
43.06
161.89
30.65
65.61
60
34.30
125.39
19.68
59.81
100
5.84
14.58
1.09
10.21
40
39.50
164.87
25.51
58.32
60
35.72
123.97
26.97
56.93
100
4.37
16.04
5.1
11.66
40
37.18
182.25
26.97
61.24
60
43.01
137.05
27.07
58.32
100
5.10
13.12
7.29
8.75
616
626
Vapor and liquid phase distribution in the reactor is considerably influenced by the pressure, temperature and hydrogen purity. As shown in Figure 8 pure hydrogen stream will increase quantity of liquid phase for a given temperature and pressure, which will result in increased catalyst wetting. Also, increased reaction pressure will expectedly decrease evaporation of middle distillates stream and thereby result in increased catalyst wetting due to increased liquid flow in the catalyst bed, as discussed previously and as shown in Figure 5. Although reaction rates are higher in the gas phase (influence of mass transfer and higher hydrogen concentrations), the influence of higher hydrogen partial pressure on the hydrogen presence in liquid phase and reaction kinetics, is sufficient to almost completely convert even DBT2, as the most refractory sulfur compound. This influence of reaction pressure and hydrogen purity was previously reported in the literature1,55,65,66 and is currently used in modern hydrotreating technology to remove sulfur from diesel fuels and produce ultra low sulfur diesel (ULSD).
ACS Paragon Plus Environment
25
Industrial & Engineering Chemistry Research
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 26 of 35
Figure 8. Catalyst wetting efficiency for inlet temperature of 606 K: a) hydrogen stream contains 77% mol hydrogen and b) hydrogen stream contains 100% mol hydrogen.
Modeling of a hydrotreating reactor which uses deterministic approach through introduction of associated phenomena like phase equilibrium, mass transfer and surface catalytic reaction, can result in very good prediction of hydrotreater reactor operation. Additional improvement of hydrotreater modeling and simulation could be achieved by the inclusion of detailed kinetic expressions for complex and consecutive HDA catalytic reaction.
5. CONCLUSION
Mathematical model of catalytic reactor for hydrotreating of straight run gas oil (SRGO) containing FCC naphtha and light cycle oil fraction (FCC N-LCO) was developed. The model assumes phase change for middle distillates stream under low pressure operating conditions of the reactor, as predicted by the Peng-Robinson Equation of State (PR EOS). Phase equilibrium calculations covering typical operating conditions of the hydrotreater reactor were used to describe vapor-liquid equilibrium of the reaction mixture within the reactor catalyst bed. Benzothiophenes, dibenzothiophenes and their alkyl derivatives, as well as aromatic compounds, were assumed to react on the catalyst surface which is in contact with the liquid or with the vapor phase. Effective wetting of the catalyst bed was applied in material balance equations for each of
ACS Paragon Plus Environment
26
Page 27 of 35
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Industrial & Engineering Chemistry Research
the specific sulfur and aromatic compounds reacting in hydrodesulfurization (HDS) and hydrodearomatization (HDA) reactions, in order to account for reactions taking place in the wet or in the dry fraction of the catalyst bed. Similar approach was used previously to model hydrogenation of cyclohexene to cyclohexane. Mass transfer limitations were taken into account through the concept of overall catalyst effectiveness for both phases, while energy balance equation has been used to describe the temperature change within the adiabatic reactor. The results of industrial test run were used to validate model results. The agreement between test run data and model predictions was very good, especially when thermodynamic model with five pseudo-components was used to describe the middle distillate stream. The error between experimental and simulated results was around 0.1% for overall sulfur conversion, and around 6% for temperature difference in reactor and aromatics conversion. Model simulations predict the important influence of reaction pressure and hydrogen purity on conversion of the overall sulfur, as well as the fact that increasing reaction temperature (up to 626 K) at low to moderate pressure (40-60 bar) along with pure hydrogen will fail to decrease overall sulfur content to Ultra Low Sulfur Diesel (ULSD) levels, when processing blends containing Light Cycle Oil. In order to achieve very low sulfur concentrations (