Dynamic Modeling and Simulation of a Slurry-Phase Reactor for

Article pubs.acs.org/EF

Dynamic Modeling and Simulation of a Slurry-Phase Reactor for Hydrotreating of Oil Fractions Cristian J. Calderón† and Jorge Ancheyta*,‡ †

Facultad de Química, UNAM, Ciudad Universitaria, México Distrito Federal 04510, México Instituto Mexicano del Petróleo, Eje Central Lázaro Cárdenas Norte 152, Col. San Bartolo Atepehuacan, Mexico City 07730, México

‡

ABSTRACT: The dynamic modeling and simulation of a slurry-phase reactor for catalytic hydrotreating of oil fractions obtained from heavy crude oil at 613 K and 5.3 MPa are reported. Compared with previous models, the proposed model incorporates the calculation of Henry constants for the volatile components through a state equation and the variation in superficial gas velocity regarding the G−L mass transfer. The mathematical model is based on the axial dispersion of gas and slurry phases and, therefore, allows approximation of the behavior of the slurry-phase system in a practical way. In order to model the hydrotreating process, different reactions are included such as hydrodesulfurization, hydrodenitrogenation, hydrodearomatization, olefins hydrogenation, and hydrocracking. Kinetic parameters and correlations for hydrodynamic and thermodynamic properties were taken from the literature. The effect of catalyst concentration, temperature, and the variation in gas phase velocity on impurities conversion is analyzed and discussed. Through the simulations, it is pointed out that under some considerations the full axial dispersion model can be simplified through ideal flow models, known as the plug-flow reactor model for the gas phase and stirred tank reactor model for the slurry phase. The need for generating experimental information in order to validate the correlations and the mathematical model is recognized.

1. INTRODUCTION Slurry-phase reactors (SPR) have gained great relevance in recent years due to a series of advantages compared with other types of reactors such as simple construction, low pressure drop, nearly isothermal operation even for highly exothermic reactions, large capacity and high conversion rates, operational flexibility, low catalyst concentration, and online catalyst addition.1,2 However, these reactors have also disadvantages such as the high back-mixing of the slurry phase which results in a distribution of residence times similar to those of a CSTR, whereby the conversion can be affected. Also, if the particle size is big enough, it would be necessary to consider abrasion and the effectiveness factor of the catalyst; if on the other hand the particles are small, this phenomenon is not considered. However, the separation of solid particles from the liquid phase after the reaction would become complicated. In some reactive systems such as FT synthesis, secondary products may be present due to the high L−S ratio in the reactor, but the main disadvantage is the lack of understanding of the hydrodynamics and consequently the scaling-up of this type of reactor.1 In the case of the petroleum industry, SPRs are mostly used for hydrocracking of heavy oils; however they can also be used for hydrotreating of lighter fractions. For example, Curtis et al.,3 Chai et al.,4,5 Deng et al.,6 and Xiang and Wang7 reported various experimental studies of slurry catalyst properties on hydrodesulfurization and hydrotreating of light oil fractions, in which general high conversions were obtained compared with those obtained in fixed-bed reactors. In the search of approaches to fulfill the stricter environmental regulations of fuel specifications, i.e., the so-called ultra low sulfur fuels, apart from developing a highly active catalyst, optimizing reaction conditions, and designing more efficient © 2017 American Chemical Society

reactor internals, among others, the utilization of other technological alternatives have been proposed for achieving such a purpose. Better catalyst effectiveness and ease of design are the main aspects that motivate the use of SPR as an option for hydrotreating of straight-run distillates. Despite slurry reactors being quite frequent on the industrial scale, mathematical models that describe their operation are scarce. This is due to the complexity in modeling of SPR, which involves many phenomena like mass and energy transfer between phases, chemical reactions, and the type of flow regime; even in some cases if the particle size is big enough, internal diffusion could be present. Various SPR models have been published in the literature in order to describe the dynamic behavior of different reacting systems such as Fischer− Tropsch (FT) synthesis, methanol synthesis, and dimethyl ether (DME) synthesis,2 which are represented by a global reaction where the gas and liquid feed to the reactor are pure components. Evidently the more components that are involved in the reacting system, the more complex is the modeling; such is the case for the hydrotreating of oil fractions. Although slurry-phase reactors are quite frequent on the industrial scale and an alternative to the hydrotreating of oil fractions, dynamical models capable of describing their operation are not reported in the literature; most of the modeling works are focused on fixed-bed or ebullated-bed reactors at steady state. A few authors have modeled slurry reactors for hydrocracking (HDC) and hydrotreating (HDT) of oil fractions with CFD techniques.8−11 Even though these models have shown Received: February 14, 2017 Revised: April 6, 2017 Published: April 10, 2017 5691

DOI: 10.1021/acs.energyfuels.7b00450 Energy Fuels 2017, 31, 5691−5700

Article

Energy & Fuels good results, they are focused on fluid dynamics and on the behavior of hydrodynamic variables; therefore despite being formulated in steady state, they are not capable of simulating and predicting detailed changes in product yields, product properties, mass and energy balances, etc. On the other hand, axial dispersion models (ADM) are the most common ones used in the literature dealing with SPR.12 ADM’s have been found to be more versatile than ideal flow models in which the gas phase is typically assumed to move as plug-flow and the slurry phase to move either as plug-flow or as a continuous stirred tank mode, if the Peclet number or dispersion coefficients correlations are accurate.13 Recently Khadem-Hamedani et al.14 presented an ADM for hydrodesulfurization of gas oil in a slurry reactor qualitatively validated with the experimental data of Deng et al.6 The reactor model takes into account the HDS of dibenzothiophene (DBT) as the main reaction. Therefore, the mass and energy balance equations are developed for DBT derivatives. However, this does not occur in reality, in which there are different sulfur components and a great number of reactions are taking place at the same time. Given that hydrotreating of middle distillates involves not only the hydrodesulfurization reaction but also other reactions such as hydrodenitrogenation, hydrodearomatization, and hydrocracking, it is then of high importance to include the kinetics of all of these reactions in the modeling in order to evaluate the slurry-phase reactor performance. This is the reason why the objective of this investigation is to propose a dynamic SPR model based on the ADM and on intrinsic kinetic reaction rates of HDT of gas oil.

Table 1. Dimensionless Numbers and Variables

εL

∂τ

=−

2 ) εG ∂ yi 1 ∂(yu i G + − Sti(yi − xi) PeG ∂z 2 uG0 ∂z

volatile liquid phase concentration

xi =

nonvolatile liquid phase concentration

xi =

reactor length

z=

h L

dimensionless time

τ=

tuG0 L

Peclet number for liquid phase

PeL =

Peclet number for gas phase

PeG =

Stanton number

Sti =

CG0 CiHi CG0 Ci C L0

uG0L DL uG0L DG

(KLa)i L HiCG0uG0

ρ εPεLL ∂xi u ∂x ε ∂ 2xi = − L0 i + L + Sti(yi − xi) + P 0 2 PeL ∂z ∂τ uG ∂z HiCGuG

NR

∑ νirj ij

(2) Nonvolatile liquid phase εL

ρ εPεLL ∂xi u ∂x ε ∂ 2xi = − L0 i + L + P0 PeL ∂z 2 ∂τ uG ∂z C LuG

NR

∑ νirj ij

(3)

with the following set of initial and boundary conditions:

2.1. Model Assumptions. In order to formulate the mathematical model, it is necessary to consider the main phenomena taking place in the reactor. However, developing a model in as much detail as possible could be a problem, especially for the HDT in a SPR, in which various reactions occur simultaneously.15 Nonetheless, some considerations have been proposed by several authors in order to simplify the modeling of SPR1,14,16−19 for HDT of oil fractions,20,21 which are as follows: • The effects of radial dispersion are not considered. • There is isothermal operation. • There is no evaporation of the liquid phase; thus the total liquid volume remains constant. This assumption is justified through a simulation in Aspen Hysys in which the hydrogen and gas oil feed were mixed under the operating conditions (5.3 MPa and 613 K). A 0.19 wt % of gas oil evaporation was observed; then it could be neglected in the modeling. • Gas side mass transfer resistance is negligible. • Due to the small catalyst particle present in the slurry phase, intraparticle concentration profiles are neglected. • All of the solid surface is in contact with the reactive phase. • There is good homogeneity between solid and liquid phases. • The solid particles follow the liquid motion, and the effects of terminal and settling velocities are negligible, which generate constant distribution of particles. • Given the range of operating conditions in HDT of oil fractions, homogeneous flow of the gas phase is assumed. • Due to the small velocity in the liquid phase, no variation in liquid velocity along the reactor is assumed. 2.2. Model Equations. On the basis of the previous considerations and taking into account the dimensionless variables given in Table 1, the dynamic equations of the ADM are as follows: Gas phase

∂yi

yi =

Volatile liquid phase

2. MODEL FORMULATION

εG

Ci

gas phase concentration

Gas phase: τ = 0, 0 ≤ z ≤ 1, yi = yi0

(4)

Liquid phase: τ = 0, 0 ≤ z ≤ 1, xi = xi0

(5)

Gas phase: 0 ≤ τ ≤ τf , z = 0,

εG d yi = (yi − yi0 ) PeG dz

Liquid phase: 0 ≤ τ ≤ τf , z = 0,

εL dxi = (xi − xi0) PeL dz

Gas and liquid phases: 0 ≤ τ ≤ τf , z = 1,

dy dxi = i =0 dz dz

(6)

(7) (8)

Equations 1−3 are the general mass balance equations. The subscript i represents each component in the system: the volatile components (H2, H2S, NH3, CH4) and the nonvolatile components (S, N, O, Poly, Di, Mono, Naph, GO, Naphtha) in the gas and liquid phases. Equations 4 and 5 establish the initial concentration along the reactor length. Equations 6 and 7 are the Dackwerts type border conditions at the inlet of the reactor, and eq 8 is the outlet boundary condition in which there are not mass fluxes. The first term of the right-hand side of eq 1 describes the convective flux for each of the components in the gas phase. Although most ADMs consider a constant superficial gas velocity1,19 there could be variations that lead to an increase or decrease in mass transfer and dispersion of both phases inside the reactor. This situation affects directly the reaction rate; therefore the effects associated with the gas phase conversion and molar flow change are considered. Then, the necessity of an extra equation that describes this phenomenon becomes evident. Khadem-Hamedani et al.14 presented a modeling and validation work based on experimental data for an HDS slurry reactor. In their model, the variation of superficial gas velocity is considered by a correlation that involves the previous knowledge of a contraction factor associated with volume change in the gas phase.

(1) 5692


Article

Energy & Fuels Table 2. Correlations Used in the Proposed Modela parameter

correlation

gas phase holdup16

⎛u ⎞ εG = 0.053⎜ G0 ⎟ ⎝ uG ⎠

gas axial dispersion coefficient25

⎛ u ⎞3 DG = 5 × 10−4⎜ G ⎟ DR 1.5 ⎝ εG ⎠

liquid axial dispersion coefficient16

DL = 3.676uG 0.32DR 1.34

gas−liquid mass transfer coefficient19

⎛ DM, i ⎞0.5 ⎟⎟ (KLa)i = (KLa)ref ⎜⎜ ⎝ DM,ref ⎠

molecular diffusivity coefficient26

⎛ v 0.267 ⎞⎛ T ⎞ DM, i = 8.93 × 10−8⎜⎜ L0.433 ⎟⎟⎜⎜ ⎟⎟ ⎝ vi ⎠⎝ μL ⎠

1.1

ρL (P , T ) = ρ0 + ΔρP − ΔρT ΔρP = (0.167 + 16.181 × 10−0.0425ρ0) × liquid density at operating conditions27

⎛ P ⎞ ⎜ ⎟ ⎝ 1000 ⎠

− 0.01 × (0.299 + 263 × 10−0.0603ρ0) ×

⎛ P ⎞2 ⎜ ⎟ ⎝ 1000 ⎠

ΔρT = (0.0133 + 152.4 × (ρ0 + ΔρP )−2.45 ) × (T − 520) − (8.1 × 10−6 − 0.0622) × 10−0.764(ρ0 +ΔρP) × (T − 520) μL = 3.141 × 1010(T − 460)−3.444 [log10 API]a a = 10.313[log10(T − 460) − 36.447]

liquid viscosity at operating conditions27

API = a

141.5 − 131.5 sg

For liquid and viscosity correlations, T is in R and P in psia.

Thus, a factor reported in the literature was used for the specific case of a Fischer−Tropsch system. However, for the HDT system, this value could not be accurate due to differences between both systems, such as the composition and molar variation of the gas phase. In order to address this situation, the equation for gas velocity is obtained by taking into consideration the total mass balance for the gas phase so that the variation is ∂uG = ∂z

NVC

∑ KLai(yi − xi)L i

RT P

reported in the literature, obtaining good results. Although these conditions remain below the operating conditions of the hydrotreating system, they have been used in different models obtaining acceptable results. For instance, in the ADM reported by Khadem-Hamedani et al.14 for hydrodesulfurization of gas oil in a slurry reactor, the most common correlations used in slurry reactors were taken into account for parameters that depend on the type of reactor configuration, such as axial dispersion coefficient and phase holdups, and parameters that depend on the reacting system, such as gas−liquid mass transfer coefficients; vapor−liquid equilibrium calculation; and liquid, gas, and slurry properties. This model was successfully validated with the experimental data of Deng et al.6 On the basis of this, the same set of correlations is implemented in the slurry reactor model for HDT. All correlations are presented in Table 2. On the other hand, if the variation in gas flow rate is going to be considered, it is necessary to implement an equation of state (EoS) to take into account all the volatile components involved in the reactions and for a better approximation of the vapor−liquid equilibrium calculation. In this work, the Peng−Robinson state equation was used. 2.4. Reaction Kinetics. In HDT of gas oil, various reactions take place at the same time, such as hydrodesulfurization (HDS), hydrodenitrogenation (HDN), hydrodearomatization (HDA), olefins hydrogenation (HGO), and hydrocracking (HDC).21 These reactions are represented as follows:

(9) 22

This approximation had only been used by Jansen et al. in a CSTR for hydrocracking involving G−L mass transfer and for Parulekar and Shah23 in an ebullated-bed reactor without mass transfer as a function of gas phase conversion. The introduction of eq 9 in the modeling allows evaluation in a more precise way of the effect of the variation of the superficial gas phase velocity on the conversion rate. 2.3. Model parameters. Despite the considerations that lead to simplifying the reactor model, it is still necessary to evaluate a series of parameters and chemical properties of the system. To do so, there are a great number of correlations proposed in the literature for different reactive systems; however, there is a lack of correlations when dealing with more complex systems operating under elevated conditions.15 Due to the absence of such experimental information for HDT systems, a good approximation is to select those correlations that have been developed for similar fluids under similar conditions compared to the system under study. Due to the complexity of performing experiments with systems operating at high pressure and temperature, it is usually preferred to use data under ambient conditions for simple air−water systems in order to generate correlations for the different parameters involved in the modeling; however, for the particular case of gas holdup and liquid axial dispersion coefficient, which are two of the most important parameters in the model, those were formulated from experimental data for an F−T system operating above 250 °C. These correlations were compared with different experimental data

kHDS

A r−S + 2H 2 ⎯⎯⎯⎯→ A r−H + H 2S

(10)

kHDN

A r−N + 3H 2 ⎯⎯⎯⎯→ A r−H + NH3

(11)

kPoly + , kPoly −

Polyaromatics (Poly) + H 2 ←⎯⎯⎯⎯⎯⎯⎯⎯⎯→ Diaromatics (Di) kDi + , kDi −

Diaromatics (Di) + 2H 2 ←⎯⎯⎯⎯⎯⎯→ Monoaromatics (Mono) 5693

(12) (13)


Article

Energy & Fuels Table 3. Kinetic Parameters Reported by Mederos et al.20,28 reaction j

k0,j

Eaj

2.639287 × 1017

150.1

kinetic model

′ jCSm1C Hm22 k in,

K0,j

ΔHAd

ΔHr,j

HDS

r′j =

HDN

rj = k inC N

1.553296 × 1012

172.28

HDA Poly

r′j = k in, j +C Poly − k in, j −C Di

2.413720 × 10

10

156.09

9.856780

58.61

HDA Di

r′j = k in, j +C Di − k in, j −CMono

2.855135 × 10

8

131.91

66.570899

58.61

HDA Mono

r′j = k in, j +CMono − k in, j −C Naph

54.31

438.006631

58.61

HGO

r′j = k1CO

(1 + K H2SC H2S)2

11.143786 9.463471 × 10

8

HCR

rj = k1CGO − k 2CGO

1.384870 × 10

HCR

rj = − k1CGO + k 3C Naphtha

2.391649 × 10

6

139.46

HCR

rj = − k 2CGO − k 3C Naphtha

2.212240 × 10

16

273.78

(14) kHGO

k1

Gas oil (GO) → Naphtha k2

Gas oil (GO) → Light gases k3

Naphtha → Light gases

(15)

3. RESULTS AND DISCUSSION 3.1. Validation of the Model with HDS Reaction. Deng et al.6 conducted an experimental study of the hydrodesulfurization of gas oil in a slurry-phase reactor with particles of supported NiMoS/γ-Al2O3 catalyst of 80 μm diameter. This particle size allows elimination of the external mass transfer phenomena between the liquid phase and the catalyst particles; therefore the HDS reaction occurs in the liquid phase. The experiments were performed in a continuous stirred tank reactor (CSTR) perfectly mixed of 1 L with a stirring speed high enough to eliminate the G−L mass transfer. Thus, the liquid phase is saturated with hydrogen; it also promotes a homogeneous catalyst distribution within the reactor. The process diagram is presented in Figure 1. Although their work focuses on

(16) (17) (18)

Equations 10 and 11 correspond to HDS and HDN reactions in which the terms Ar−S and Ar−N represent all the sulfur and nitrogen compounds in the feed, respectively. Eqs 12−14 stand for the reversible reactions of HDA of polyaromatics, diaromatics, and monoaromatics. Equation 15 represents the olefins’ hydrogenation, and eqs 16−18 deal with the HDC reactions. The corresponding kinetic expressions, reaction rate coefficients, and other parameters were taken from Mederos et al.,20 which were successfully applied to describe the dynamical behavior of a trickled-bed reactor. The experimental conditions under which all of the parameters were obtained ensure that they are free of internal and external diffusion limitations; therefore the kinetic parameters are considered to be intrinsic. All reactions and kinetic parameters are presented in Table 3. The intrinsic reaction rate coefficients and adsorption and equilibrium constants for HDS and HDA reactions are calculated with the following equations: k in, j = k 0, j e(−Eaj / RT )

(19)

K H2S = K 0,H2S e(−ΔHAd / RT )

(20)

k in, j +

= K 0, j e[−ΔHr,j / R(1/ T0− 1/ T )]

56.41

volatile and nonvolatile components involved in the reacting system. After the discretization of the spatial coordinates, the ODE system is solved numerically as an initial value problem through MATLAB software, particularly the ode23 subroutine.

Monoaromatics (Mono) + 3H 2 ←⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯→ Saturates (Sat) R−CHCH−R′ + H 2 ⎯⎯⎯⎯→ R−CH 2−CH 2−R′

34.02

138.15

−1

kMono + , kMono −

5.1668888

(21)

Figure 1. Schematic representation for the experimental reactor of Deng et al.6

2.5. Model Solution. The dynamic model of the slurry reactor is defined by eqs 1−3 along with the initial and boundary conditions given by eqs 4−8. This set of equations forms a system of hyperbolic partial differential equations (PDE) in which a numerical solution is needed. Therefore, the system of PDE needs to be reduced into an ordinary differential equation (ODE) system in order to find the numerical solution. The method of orthogonal collocation, which is based on the concept of interpolation of unequally spaced points, is used. The orthogonal collocation consists of choosing an orthogonal polynomial function (typically the Legendre polynomials) to approximate the solution of a differential equation for a given interval. The dimension of the resulting matrix system of ODE depends on the grade of the polynomial to be chosen and on the number of state variables. In this work, there are 18 state variables that stand for the

evaluating the effects of catalyst concentration (4.8−23.1 wt %), temperature (320−360 °C), and pressure (3−5 MPa) on conversion rate, their results, particularly the conversion achieved at steady state, can be used as a reference point to compare with a reactor model. This procedure was used by Khadem et al.14 in order to validate their reactor model under the same operating conditions, particularly at 350 °C and 5 MPa. Under these conditions and with a WHSV of 1.67 h−1, conversions of 98.3% are expected in the HDS reaction on a laboratory scale. Hence, by maintaining the volume and space velocity relationship in the industrial scale model, the system would be expected to reach a similar steady state. The outline of this process is

Kj =

k in, j −

5694


Article

Energy & Fuels presented in Figure 2. One of the assumptions made by the authors in the reactor modeling is that the hydrodesulfurization

Figure 2. Schematic representation of an industrial slurry-phase HDT reactor.

reaction is controlled by the hydrodesulfurization of DBT compounds; therefore the kinetic model of Vanrysselberghe and Froment24 presented in eqs 22−25 as well as the reaction rate expressions in eqs 26−29 are used in their dynamic model. kDBT, σ

DBT + 2H 2 ⎯⎯⎯⎯⎯→ BPH + H 2S (σ site) kDBT, τ

DBT + 5H 2 ⎯⎯⎯⎯⎯→ CHB + H 2S (τ site) kBPH, τ

BPH + 3H 2 ⎯⎯⎯⎯⎯→ CHB (τ site) k CHB, τ

CHB + 3H 2 ⎯⎯⎯⎯⎯→ BCH (τ site) rDBT, σ =

Figure 3. Comparison of dynamic simulation of HDS of gas oil at z = L. (···) Model proposed by Khadem-Hamedani et al.14 (−) Model proposed in this work. (○) Experimental value.

exit of the reactor is represented by a void circle. As can be seen, for both simulations, the dynamic behavior is similar. The main difference is in the steady-state prediction, at which the proposed model tends to slightly improve the prediction of the experimental value. It should be addressed that the numeric method of solution could have an influence on the simulation results: orthogonal collocation in this work versus finite differences in the work of Khadem-Hamedani et al.14 It is then concluded that the proposed reactor model properly represents the literature results; however, the simulations did not match the experimental value with good accuracy. This is because the kinetic model used to describe the hydrodesulfurization reaction is based on one component, in this case the DBT. In order to have a better approximation of the HDS reaction, a model based on the total sulfur composition of the feed should be considered, or kinetic models based on pseudocomponents with different types of sulfur compounds used. In addition, in a hydrotreating process, different reactions occur simultaneously, not only the HDS reaction as discussed above; therefore the inclusion of more kinetic expressions to improve the predictions would be necessary. 3.2. Simulation for All HDT Reactions. Once the reactor model is validated, it was used to simulate the most important reactions taking place in the hydrotreating of gas oil. The kinetic models used were those previously reported by Mederos et al.20 for gas oil obtained from heavy crude oil, which were obtained at 613 K and 5.3 MPa under the kinetic regime, so that they can be implemented directly in the slurry-phase reactor model considering the same pressure and temperature. It was assumed that the catalyst used for the simulations was the same as the one used by Mederos et al.20 in their experimental work. It was assumed that a trilobular cobalt molybdenum supported on alumina catalyst (equivalent particle diameter of 0.254 cm, total pore volume of 0.5 cm3/g, bulk density of 0.92 g/cm3) was used, with the exception that for the slurry-phase reactor, this catalyst needs to be crushed in order to originate a small particle diameter, about 0.008 cm, in order to eliminate the effects of internal

(22) (23) (24) (25)

kDBT, σK H2, σKDBT, σC DBTC H2 (1 + KDBT, σ + (K H2, σC H2)1/2 + KBPH, σC BPH + K H2SC H2S)3

(26) rDBT, τ =

rBPH, τ =

rCHB, τ =

kDBT, τK H2, τKDBT, τC DBTC H2 (1 + KDBT, τ + (K H2, τC H2)1/2 + KBPH, τC BPH)3

(27)

kBPH, τK H2, τKBPH, τC DBTC H2 (1 + KDBT, τ + (K H2, τC H2)1/2 + KBPH, τC BPH)3

(28)

k CHB, τK H2, τK CHB, τCCHBC H2 (1 + KDBT, τ + (K H2, τC H2)1/2 + KBPH, τC BPH)3

(29)

with DBT, BPH, CHB, and BCH being dibenzothiophene, biphenyl, cyclohexyl benzene, and bicyclohexyl, respectively. In order to corroborate the performance of the dynamic model proposed by Khadem-Hamedani et al.,14 their simulations were reproduced with the exact same data of operating conditions and feed properties. After that, the proposed model given by eqs 1−9 is used for validation. All the simulations were carried out for a feed with the following properties: an API gravity of 34 and sulfur content of 2200 ppm at a temperature of 623 K and pressure of 5 MPa. Figure 3 shows a comparison between the simulated results at the exit of the reactor with the model reported by KhademHamedani et al.14 and that proposed in this work for the hydrodesulfurization of gas oil considering the kinetic model of Vanrysselberghe and Froment.24 The experimental value at the 5695


Article

Energy & Fuels

HDA reactions (Figure 5), the behavior is different because they are reversible. While the concentration of polyaromatics

mass transfer and to consider the pseudohomogeneous (L−S) phase. Feed properties and operating conditions used in the simulations for the slurry-phase reactor are presented in Table 4. Table 4. Feed Properties and Operating Conditions parameter liquid phase (gas oil) density (g/cm3) molecular mass (g/mol) viscosity (cp) inlet superficial velocity (m/s) sulfur concentration (mol/cm3) nitrogen concentration (mol/cm3) olefins concentration (mol/cm3) polyaromatics concentration (mol/cm3) diaromatics concentration (mol/cm3) monoaromatics concentration (mol/cm3) naphthenes concentration (mol/cm3) gas phase (hydrogen) density (g/cm3) inlet superficial velocity (m/s) operating conditions temperature (K) pressure (MPa)

value 0.873 248.7 7.2554 0.005 4.9 × 10−4 1.7 × 10−5 1.4 × 10−4 1.9 × 10−4 2.3 × 10−4 6.4 × 10−4 1.35 × 10−3 0.0013 0.1 613 5.3

Figure 5. Dynamic profiles of aromatic components at z = L. (−) εP = 0.25, (···) εP = 0.5. (+) Poly, (×) Di, (▽) Mono, (Δ) Sat.

Figures 4 shows the results of the dynamic simulation for the liquid phase components (sulfur, nitrogen, olefins) and the

and diaromatics reduces, those of monoaromatics and saturates increase. In all cases, the dynamic behavior reaches the steady state after 2000 s. The dynamic profiles of liquid concentration of gas oil, naphtha, and light gases at the exit of the reactor are shown in Figure 6. Evidently, the change of concentration produced by the cracking of gas oil to naphtha is bigger than the concentration of light gases dissolved in the liquid phase due to the high volatility of these components. In addition, it can be seen that the light gases have faster dynamics and reach the steady state before naphtha and gas oil do; this means that the gas oil

Figure 4. Dynamic profiles of liquid phase concentrations of sulfur, olefins, and nitrogen at z = L. (−) εP = 0.25, (···) εP = 0.5. (○) Sulfur, (∗) olefins, (□) nitrogen.

effect of catalyst concentration on conversion at the exit of the reactor. As expected, an increase in catalyst concentration favors the conversion of all reactions. However, it should be noted that the influence of particle concentration does not affect in the same way all the reactions. HDS reaction is the least affected due to the inhibition effect originated by H2S formation in the liquid phase, while HGO and HDN reactions are greatly influenced by catalyst concentration. In the case of

Figure 6. Dynamic profiles of hydrocracked products at z = L. 5696


Article

Energy & Fuels is saturated with the light gases and the production of naphtha continues. The variation of partial pressure of gas phase components is presented in Figure 7. Hydrogen partial pressure decreases

Figure 8. Axial concentration profiles of liquid and gas phase components at steady state. (−) Non volatile components. (○) Sulfur, (∗) olefins, (□) nitrogen, (+) Poly, (×) Di, (▽) Mono, (Δ) Sat. (- -) volatile components. (○) H2, (∗) H2S, (+) NH3, (×) light gases.

Figure 7. Dynamic profiles of partial pressure of gas phase components at z = L.

rapidly around 200 s due to the high absorption and hydrogen consumption in the liquid phase caused by the reactions taking place, mainly HDS. Once the HDS reaction reaches the steady state, the liquid starts to saturate with hydrogen and the partial pressure of the gas phase increases until it reaches the steady state. This effect can also be seen in Figure 4. Inversely to H2 behavior, at the beginning of the reaction, the partial pressure of H2S increases considerably until the steady state is reached and the production of H2S is constant. Moreover, the rates of HDN and HCR reactions are slower than that of the HDS reaction; that is why the productions of NH3 and light gases are considerably lower. Figure 8 shows the axial concentration profiles of liquid phase components. It can be seen that the concentration of liquid components remains practically constant along the reactor length, which agrees with literature reports for other slurry-phase systems in which, due to the hydrodynamic conditions inside the reactor, it is assumed that the liquid phase can be modeled as an ideal dispersion. Therefore, a CSTR model can be used instead of a full ADM. This would simplify in a great manner the model equations and the numeric method of solution. Figure 9 shows the profile of partial pressure of gas phase components; as can be noted, those that are present in high amounts (H2 and H2S) exhibit a major concentration profile along the reactor length. As expected, the partial pressure of H2 shows a decrease with respect to the reactor axial position due to its absorption and consumption in the reacting system. Partial pressure of H2S, NH3, and the light gases showed an evident increase as the reactions occurred. The observed profiles within the reactor indicate that the convective terms in the general mass balance can be assumed to be dominant over the dispersion effects; therefore an obvious simplification could be the use of a PFR model instead of the full ADM.

Figure 9. Axial partial pressure profiles of gas phase components at steady state.

The variation in gas phase velocity along the reactor is presented in Figure 10. As mentioned before, at the beginning of the reactions there is a high gradient in hydrogen concentration (Figure 7) which generates a sudden change in gas velocity. Then, the system starts to stabilize, generating a constant superficial gas velocity profile in each point of the reactor. While there is a change in this velocity, it does not affect to a great extent the conversion of the liquid phase, as can be seen in Figure 11, in which the difference between the simulations with constant and variable velocity for gas and liquid components is plotted. Figure 12 shows the influence of inlet gas velocity in the HDS reaction as well as in dissolved hydrogen and H2S in the liquid phase. The enhanced sulfur conversion at higher inlet 5697


Article

Energy & Fuels

Figure 10. Dynamic profile of gas phase velocity at different reactor lengths. (- -) z = 0.5 m, (- ·) z = 2.3 m, (···) z = 4.7 m, (−) z = 6.5 m. Figure 12. Effect of inlet superficial gas velocity on sulfur conversion and H2 and H2S compositions in the liquid phase.

Figure 11. Difference between simulations with constant velocity (−) and variable velocity (- -) at z = L for sulfur and gas components.

superficial gas phase velocity is because an increase in this velocity improves the dispersion of the liquid phase and provokes better mixing between gas and liquid phases, which facilitates the mass transfer of H2 and H2S to the low concentration phase, that being the liquid phase for H2 and the gas phase for H2S. Having a greater amount of H2 in the liquid phase and a smaller amount of H2S benefits the conversion of sulfur. The temperature effect on HDT reactions can be seen in Figures 13 and 14. As expected, increasing the temperature favors the HDS, HDN, and HGO reaction rates. However, special attention should be noted for the HDA reactions. On the one hand, polyaromatic and diaromatic HDA reactions are favored by the increase in temperature, although not to a great extent. On the other hand, for the monoaromatic HDA reaction, a negative effect in conversion is observed even at high temperatures due to hydrogenation of polyaromatics, and diaromatics is faster than monoaromatics hydrogenation. This has been previously observed and explained in the literature by Chowdhury et al.29 and Mederos et al.20 in which the

Figure 13. Effect of temperature on HDS, HGO, and HDN reactions at steady state. (○) Sulfur, (∗) olefins, (□) nitrogen.

thermodynamic equilibrium is found at 360 °C (633 K), which is not the case for the slurry phase reactor.

■

CONCLUSIONS The following conclusions can be pointed out: • It is possible to model a slurry-phase reactor in which multiple reactions occur simultaneously using an axial dispersion model with intrinsic reaction rate expressions and appropriate hydrodynamic correlations. • Introducing an equation of state improves the robustness of the model, especially in multiphase and multicomponent systems where volatile and nonvolatile components are involved. • Assuming that only the effects of gas−liquid mass transfer affect the superficial gas velocity increases considerably the complexity of the reactor model; however, considering that the superficial gas velocity remains constant yields similar results. 5698


Article

Energy & Fuels

DBT = Dibenzothiophene DM,i = Molecular diffusion coefficient of component i, cm3f / cms·s DGi = Axial dispersion coefficient of gas phase, cm2r /s DLi = Axial dispersion coefficient of liquid phase, cm2r /s Eaj = Activation energy for j reaction, J/mol HCR = Hydrocracking reaction HDA = Hydrodearomatization reaction HDN = Hydrodenitrogenation reaction HDS = Hydrodesulfurization reaction HGO = Olefins hydrogenation reaction Hi = Henry’s law coefficient of component i, MPa·cm3L/mol i ki,w = Rate coefficient of component i, on w site for HDS of DBT reactions, kmol/kgcat·h kj = Intrinsic reaction rate constant for j reaction, cm3L/(mol· gs·s) k0,j = Frequency factor for j reaction, mol/(cm3L·s) K0,j = Equilibrium constant for j reaction, dimensionless KH2S = Adsorption-equilibrium constant for H2S, cm3L/mol Ki,w = Adsorption coefficient of component i on w site for HDS of DBT reactions, m3L/kmol KL = Overall gas−liquid mass transfer coefficient, cm3L/(cm3I · s) L = Reactor longitude, cm MONO = Monoaromatic compounds NVC = Number of volatile components NR = Number of reactions P = Pressure, MPa P0 = Inlet pressure of gas phase, MPa PeL = Peclet number of liquid phase, dimensionless PeG = Peclet number of gas phase, dimensionless Poly = Polyaromatic compounds R = Gas law constant, J/mol·K rj′ = Rate of reaction j per unit of catalyst mass in the liquid phase, mol i/(gS·s) rj = Rate of reaction j per unit of volume in the liquid phase, mol i/(gS·s) sg = Specific gravity of i component, dimensionless t = Time, s T = Temperature, K T0 = Reference temperature, K u0G = Inlet superficial gas velocity, cm/s uG = Superficial gas velocity, cm/s u0L = Inlet superficial liquid velocity, cm/s uL = Superficial liquid velocity, cm/s V = Reactor volume, cm3r xi = Mole fraction of i component in the liquid phase, moli/ molL x0i = Inlet mole fraction of i component in the liquid phase, moli/molL yi = Mole fraction of i component in the gas phase, moli/ molG y0i = Inlet mole fraction of i component in the gas phase, moli/molG z = Axial position, cm

Figure 14. Effect of temperature on HDA reactions at steady state. (+) Poly, (×) Di, (▽) Mono, (Δ) Sat.

• The reactor model can be simplified similarly to other systems in the slurry phase wherein all the components in the liquid phase are considered as a completely mixed reactor (CSTR) and the gas phase components as an ideal plug-flow reactor (PFR). • There is an evident lack of experimental information that allows for better model validation and for the formulation of correlations that can be applied to different oil fraction systems.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Jorge Ancheyta: 0000-0001-9626-637X Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank the Mexican Institute of Petroleum for financial support. Cristian Calderón also thanks Consejo Nacional de Ciencia y Tecnologiá (CONACYT) for his Ph.D. scholarship.

■

NOMENCLATURE

Symbols

a = Interfacial area per unit reactor volume, cm2/cm2r API = API gravity ArN = Nitrogen compounds ArS = Sulfur compounds BPH = Biphenyl BCH = Bicyclohexyl C0G = Inlet gas phase concentration, mol/cm3 CGi = Concentration of i component in the gas phase, mol/ cm3 C0L = Inlet liquid phase concentration, mol/cm3 CLi = Concentration of i component in the liquid phase, mol/ cm3 CHB = Cyclohexyl benzene Di = Diaromatic compounds

Subscripts

G = Referring to gas phase i = Referring to i component L = Referring to liquid phase P = Referring to solid particles S = Referring to solid phase 5699


Article

Energy & Fuels Superscripts

(17) Maretto, C.; Krishna, R. Modelling of a Bubble Column Slurry Reactor for Fischer−Tropsch Synthesis. Catal. Today 1999, 52, 279− 289. (18) Coselli Vasco de Toledo, E.; Leite de Santana, P.; Regina Wolf Maciel, M.; Maciel Filho, R. Dynamic Modelling of a Three-Phase Catalytic Slurry Reactor. Chem. Eng. Sci. 2001, 56, 6055−6061. (19) de Swart, J.; Krishna, R. Simulation of the Transient and SteadyState Behaviour of a Bubble Column Slurry Reactor for Fischer− Tropsch Synthesis. Chem. Eng. Process. 2002, 41, 35−47. (20) Mederos, F. S.; Ancheyta, J.; Elizalde, I. Dynamic Modeling and Simulation of Hydrotreating of Gas Oil Obtained from Heavy Crude Oil. Appl. Catal., A 2012, 425, 13−27. (21) Ancheyta, J. Modeling of Processes and Reactors for Upgrading of Heavy Petroleum; CRC Press Taylor and Francis Group: Boca Raton, FL, 2013. (22) Jansen, T.; Guerry, D.; Leclerc, E.; Ropars, M.; Lacroix, M.; Geantet, C.; Tayakout- Fayolle, M. Simulation of Petroleum Residue Hydroconversion in a Continuous Pilot Unit Using Batch Reactor Experiments and a Cold Mock-Up. Ind. Eng. Chem. Res. 2014, 53, 15852−15861. (23) Parulekar, S. J.; Shah, Y. T. Steady-State Behavior of Gas-LiquidSolid Fluidized-Bed Reactors. Chem. Eng. J. 1980, 20, 21−33. (24) Vanrysselberghe, V.; Froment, C. Hydrodesulfurization of Dibenzothiophene on a CoMo/Al2O3 Catalyst: Reaction Network and Kinetics. Ind. Eng. Chem. Res. 1996, 35, 3311−3318. (25) Wilkinson, P. M.; Spek, A. P.; van Dierendonck, L. L. Design Parameters Estimation for Scale-up of High Pressure Bubble Columns. AIChE J. 1992, 38, 544−554. (26) Poling, B. E.; Prausnitz, J. M.; O’Connell, J. P. The Properties of Gases and Liquids; 5th ed.; McGraw-Hill: New York, 2004. (27) Ahmed, T. H. Hydrocarbon Phase Behavior; Gulf Publishing Company: Houston, TX, 1989; Vol 7. (28) Jaffe, S. B. Kinetics of Heat Release in Petroleum Hydrogenation. Ind. Eng. Chem. Process Des. Dev. 1974, 13, 34−39. (29) Chowdhury, R.; Pedernera, E.; Reimert, R. Trickle-Bed Reactor Model for Desulfurization and Dearomatization of Diesel. AIChE J. 2002, 48, 126−135.

0 = Inlet conditions Greek letters

α = Gas phase contraction factor, dimensionless ΔHAd = Adsorption enthalpy of H2S, J/mol μL = Liquid phase viscosity cP νfi,j = Stoichiometric coefficient of component i in reaction j in the f phase, dimensionless ρP = Solid phase density, g/cm3P ρL = Liquid phase density, g/cm3L τ = Dimensionless time τf = Dimensionless final time of simulation εG = Holdup of gas phase, cm3G/cm3r εL = Holdup of liquid phase, cm3L/cm3r εP = Solid fraction in the liquid phase, cm3S/cm3L

■

REFERENCES

(1) Rados, N.; Al-Dahhan, M.; Dudukovic, M. P. Dynamic Modeling of Slurry Bubble Column Reactors. Ind. Eng. Chem. Res. 2005, 44, 6086−6094. (2) Wang, T.; Wang, J.; Jin, Y. Slurry Reactors for Gas-to-Liquid Processes: A Review. Ind. Eng. Chem. Res. 2007, 46, 5824−5847. (3) Curtis, C. W.; Chen, J. H.; Tang, Y. Hydrodesulfurization of Model Systems Using Slurry-Phase Catalysts. Energy Fuels 1995, 9, 195−203. (4) Chai, Y.; Xiang, C.; Kong, H.; Liu, Y.; Liu, C. A Study on SlurryBed Hydrotreating Process of Distillate Oil: I Catalyst preparation. J. Fuel. Chem. Technol. 2008, 6, 014. (5) Chai, Y.; Xiang, C.; Kong, H.; Liu, Y.; Liu, C. A Study on Process for Slurry-Bed Hydrotreating of Distillate Oil: II Processing Conditions and Kinetics Analysis of FCC Diesel Oil. J. Fuel. Chem. Technol. 2009, 1, 010. (6) Deng, Z.; Wang, T.; Wang, Z. Hydrodesulfurization of Diesel in a Slurry Reactor. Chem. Eng. Sci. 2010, 65, 480−486. (7) Xiang, H.; Wang, T. Kinetic Study of Hydrodesulfurization of Coker Gas Oil in a Slurry Reactor. Front. Chem. Sci. Eng. 2013, 7, 139− 144. (8) Carbonell, M. M.; Guirardello, R. Modelling of a Slurry Bubble Column Reactor Applied to the Hydroconversion of Heavy Oils. Chem. Eng. Sci. 1997, 52, 4179−4185. (9) Matos, E.; Guirardello, R. Modeling and Simulation of a PseudoTwo-Phase Gas-Liquid Column Reactor for Thermal Hydrocracking of Petroleum Heavy Fractions. Braz. J. Chem. Eng. 2002, 19, 319−334. (10) Matos, E.; Nunhez, J. The Effect of Different Feed Flow Patterns on the Conversion of Bubble Column Reactors. Chem. Eng. J. 2006, 116, 163−172. (11) Matos, E.; Guirardello, R.; Mori, M.; Nunhez, J. Modeling and Simulation of a Pseudo-Three-Phase Slurry Bubble Column Reactor Applied to the Process of Petroleum Hydrodesulfurization. Comput. Chem. Eng. 2009, 33, 1115−1122. (12) Duduković, M. P.; Larachi, F.; Mills, P. L. Multiphase Catalytic Reactors: A Perspective on Current Knowledge and Future Trends. Catal. Rev.: Sci. Eng. 2002, 44, 123−246. (13) Rados, N.; Al-Dahhan, M.; Dudukovic, M. P. Modeling of the Fischer−Tropsch Synthesis in Slurry Bubble Column Reactors. Catal. Today 2003, 79−80, 211−218. (14) Khadem-Hamedani, B.; Yaghmaei, B. K.-H.; Fattahi, M.; Mashayekhan, S.; Hosseini-Ardali, S. M. Mathematical Modeling of a Slurry Bubble Column Reactor for Hydrocracking of Diesel Fuel: Single and Two Bubble Configurations. Chem. Eng. Res. Des. 2015, 100, 362−376. (15) Calderón, C. J.; Ancheyta, J. Modeling of Slurry-Phase Reactors for Hydrocracking of Heavy Oils. Energy Fuels 2016, 30, 2525−2543. (16) Deckwer, W. D.; Serpemen, Y.; Ralek, M.; Schmidt, B. Modeling the Fischer−Tropsch Synthesis in the Slurry-Phase. Ind. Eng. Ind. Eng. Chem. Process Des. Dev. 1982, 21, 231−241. 5700


Dynamic Modeling and Simulation of a Slurry-Phase Reactor for

Recommend Documents