Detailed kinetic modelling of vacuum gas oil hydrocracking using bifunctional catalyst: A distribution approach

Chemical Engineering Journal 284 (2016) 270–284

Contents lists available at ScienceDirect

Chemical Engineering Journal journal homepage: www.elsevier.com/locate/cej

Detailed kinetic modelling of vacuum gas oil hydrocracking using bifunctional catalyst: A distribution approach B. Browning a,⇑, P. Afanasiev a, I. Pitault b, F. Couenne b, M. Tayakout-Fayolle a a b

Institut de Recherche sur la Catalyse et l’Environnement de Lyon, CNRS, Université de Lyon, 43 bd du 11 novembre 1918, 69616 Villeurbanne, France LAGEP, CNRS, Université de Lyon, 43 bd du 11 novembre 1918, 69616 Villeurbanne, France

h i g h l i g h t s
Hydrocracking of a real vacuum gas oil was modelled based on GC  GC analysis results.
The model includes 217 lumps, classified by carbon number and hydrocarbon family.
17 kinetic parameters were estimated and the model results match the data well.
The estimated parameters are coherent with the literature.
The effects of molecular structure are observed in the model results.

a r t i c l e

i n f o

Article history: Received 7 April 2015 Received in revised form 26 August 2015 Accepted 27 August 2015 Available online 3 September 2015 Keywords: Kinetic modelling Hydrocracking Lumped model GC  GC

a b s t r a c t A kinetic model for hydrocracking a real vacuum gas oil at 120 bar and 400 °C in a semi-batch reactor has been constructed based on analysis results from two dimensional gas chromatography (GC  GC). The model has 217 pseudo-component lumps classified by carbon number and hydrocarbon family. Hydrogenation and cracking reactions are included separately and product distributions are generated using a probabilistic approach. Mass transfer resistances and vapour–liquid equilibrium are also accounted for. 17 parameters were estimated and the resulting model is able to fit the experimental data quite well. It was found that it is important to consider the average molecular structure within each lump. To this end, the non-reactivity of components within the paraffin lumps is taken into account in the probability matrix. A parameter has been included for the number of branches removed from the aromatics and naphthenes and it was also necessary to separate the aromatics and naphthenes into mono-ringed and multi-ringed lumps. This many lumped model allows us to approach the results found in studies with model molecules. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction Hydrocracking is used in industry to upgrade vacuum gas oil fractions. This is of increasing importance because of the need to meet a growing demand for middle distillate (MD), mainly due to the rapid growth of developing Asian markets [1]. Whilst hydrocracking is more costly for upgrading oil fractions than the traditional method of catalytic cracking, it offers refiners more flexibility through improved product quality and the possibility of processing heavier feedstocks. A typical hydrocracking process is carried out over a bifunctional catalyst in a trickle bed reactor operating under high temperature and pressure, between 340 °C and 450 °C and 80–180 bar. ⇑ Corresponding author. Tel.: +33 4 72 44 53 26. E-mail address: [email protected] (B. Browning). http://dx.doi.org/10.1016/j.cej.2015.08.126 1385-8947/Ó 2015 Elsevier B.V. All rights reserved.

Because oil fractions are complex mixtures and hydrocracking is carried out under severe conditions, it is difficult to unravel the influences of composition, mass transfer, thermodynamics and intrinsic reaction kinetics on product composition. At the same time, it is important to the oil industry to maximise selectivity for MD. Indeed, much laboratory research has been given to understanding hydrocracking kinetics. Many key experimental studies [2,3] for hydrocracking reaction mechanisms over bifunctional catalysts have used model molecules. The strengths of this method are that, by isolating representative components, elementary hydrocracking steps can be studied and analysed in detail and precise reaction rate parameters can be estimated [4–6]. The drawbacks for extrapolating to real feedstocks are that the molecules chosen for study tend to be relatively small hydrocarbons and do not contain traces of sulphur or nitrogen containing compounds, which in real feeds may affect reaction kinetics, for example, by

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284

271

Nomenclature Aro Aro= API C C0 D DiB De d F IP IP= i j k kB kLa ksL ksz m MB n nC nH nDB NP NP= N N= Nimp Nps pi r r0 rh Re s S Sc Sh t

aromatic aromatic with double bond on branch API gravity (–) concentration (mol m3) concentration (molC m3) objective function final value diffusion coefficient of component i in reaction mixture (m2 s1) effective diffusion coefficient (m2 s1) diameter (m) flux (mol s1) iso-paraffin dehydrogenated iso-paraffin any hydrogenated lump dehydrogenated lump corresponding to i reaction constant (s1) Boltzmann constant (1.38  1023 m2 kg s2 K1) vapour–liquid mass transfer coefficient (s1) liquid–solid mass transfer coefficient (m s1) zeolite phase mass transfer coefficient (m s1) mass (g) molecular weight (g mol1) number of data points (–) carbon number (–) number of hydrogen atoms (–) number of double bonds (–) N-paraffin dehydrogenated n-paraffin naphthene naphthene with double bond on branch impeller frequency (s1) number of pseudocomponents (lumps) (–) parameter production rate (mol s1) carbon production rate (molC s1) hydrodynamic radius (m) Reynold’s number (–) surface tension (kg s2) surface area (m2) schmidt number (–) Sherwood number (–) time (s)

poisoning active catalyst sites [7]. Work with real feedstocks allows these aspects to be considered and process plant conditions to be replicated closely, potentially leading to useful data for design and development of new units and catalysts. The bifunctionality of the catalyst comes from the presence of both Brönsted acid sites which catalyse the cracking and isomerisation reactions, and sites for hydrogenation/dehydrogenation in the form of a metal sulphide or noble metal [3]. The accepted mechanism for hydrocracking over a bifunctional catalyst is reported by Weitkamp [8] and includes several elementary reaction steps. Firstly, the reacting molecule is adsorbed onto a metallic site and dehydrogenated to form a double bond. It then migrates to an acid site, the double bond is protonated and a carbenium ion forms. Isomerisation occurs and the molecule may also undergo b-scission, i.e. cracking. The product(s) then migrate to another metallic site where they are hydrogenated before desorbing and diffusing out of the catalyst. Iso-paraffins undergo b-scission more readily than normal paraffins [2]. Opening of aromatic and naphthenic rings can occur but requires strong hydrogenolytic

T V v We X Y

a

b

C K

e

h

l q ri s


u £

temperature (K) volume (m3) molar volume (m3 mol1) Weber number (–) conversion (%) selectivity (%) fraction of Aro= and N= producing NP= (–) fraction of Aro and N producing > 1 NP or IP (–) gamma function (–) number of parameters porosity (–) parameter for gamma distribution (–) viscosity (mPa s1) density (kg m3) standard deviation tortuosity (–) stoichiometric coefficient parameter for gamma distribution (–) association factor (–)

Subscripts cat catalyst c cracking hydrogen H2 hyd/dehyd hydrogenation/dehydrogenation reaction i any hydrogenated lump j dehydrogenated lump corresponding to i imp impeller in entering reactor l liquid phase MD middle distillate p pore r reactor s metal phase v vapour phase v volumetric VGO vacuum gas oil z zeolite phase 0 initial condition/composition 1 mono-ringed Aro or N 2 multi-ringed Aro or N

functionality, provided by noble metals [9]. However, at lower conversion, the branched chains of naphthenes can easily crack and the paring reaction leads to an isomerization of rings before cracking [8,10]. The homogeneity of the catalyst and the relative proportions of the different types of sites strongly influence the catalyst selectivity and performance [11]. Zeolite crystallites are used to provide the acidic sites, and mass transfer limitations are a potential issue here due to the very small pore sizes (1 nm) [12]. Martens and Marin found an efficiency of 0.8 for hydrocracking components of size C18+ over a bifunctional catalyst [13], which gives an idea of the expected degree of internal mass transfer resistance. Modelling of real feedstocks and effluents is advancing with computing power and improvements in analysis tools and techniques. Lumping techniques have been used historically to represent petroleum conversion processes such as Fluid Catalytic Cracking (FCC), hydrocracking and hydroconversion. In their 2005 review, Ancheyta et al. [14] classified the types of models types used in the literature. The simplest and most common type of

272

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284

model uses a traditional lumping method [15–20]. The reaction mixture composition is analysed and lumped by True Boiling Point (TBP). Each lump is equivalent to a pseudo-component defined only by its TBP range. The number of lumps varies but five are often used, including one to represent gaseous products. The temperature limits of these TBP ranges vary. Within this group of models, there exists a wide range of possible reaction schemes, with products being formed directly from the feedstock and/or via the other lumps. Classification and lumping of the reaction mixture by carbon number has also been used with this type of model, as that developed for hydrocracking Fischer–Tropsch paraffin products by Pellegrini et al. [21]. For vacuum gas oils (VGO), cracking reactions are always considered to be macroscopically irreversible and first order [15,22–24]. Mostly, these models consider only cracking reactions. However, some authors do take the hydrogenation reactions into account separately [25–27]. Use of the discrete lumping method with a more detailed analysis of the reaction mixture allows information about the molecular structures to be used in conjunction with the TBP data [28–30]. For example, Sun et al. [31] recently constructed and validated a 7-lump kinetic model for the hydrotreatment of coal tar, based on TBP and SARA (saturated hydrocarbon, aromatic hydrocarbon, resin, asphaltene) analysis. The benefits and drawbacks of the traditional and discrete lumping methods lie in their simplicity. With relatively few equations to be solved and few parameters to be estimated, computing times are reasonable as is the amount of experimental data required. On the other hand, information about the reaction mixture is lost. The number of lumps can be increased to better exploit the experimental data but this rapidly leads to a very large number of reaction rate parameters. To overcome this problem, methods were developed where the rate constants are described by a function. Stangeland [32] maintained the concept of a lumped boiling point distribution and introduced a correlation for hydrocracking rate as a function of TBP. Laxminarasimhan et al. [33] considered the boiling point range as a continuous function and represented it as a statistical distribution. Both of these methods, based only on TBP, have proved useful for modelling industrial units [34] and continue to be developed [35–40]. As with the traditional lumping method, a more detailed analysis of the feedstock and reaction mixture allows these models to be extended so that hydrocarbon structure is also taken into account [41,42]. The distribution approach that we take for our model has the advantages of continuous lumping, in that little information is lost and the number of parameters required is not excessive. It also has the benefit that, through considering many discrete lumped pseudocomponents, attention is drawn to the molecular structure. This should allow us to produce a model which is less empirical and can take account of the physical differences between the lumps. Also, since the experimental results for hydrocracking are often reported in the literature as distributions, it seems logical to work with the same format for modelling. As opposed to the empirical lumping models described above, Single Event MicroKinetic (SEMK) models approach the problem in the opposite way by taking the underlying elementary reactions as a starting point. The elementary reactions are then classified by type to reduce the number of kinetic parameters. Both Ancheyta et al. [14] and Froment et al. [43] review this method which continues to be developed [44,45]. Due to the substantial experimental effort made towards this method, SEMK allows for a detailed vision of complex reaction systems. It would be of great interest to compare the lumped distribution approach with SEMK. Another important aspect of hydrocracking modelling is the vapour–liquid equilibrium (VLE) inside the hydrocracking reactor, which affects the reaction kinetics. Chavez et al. have recently published a detailed review on the subject [46]. Typically, most of the

hydrogen is in the vapour phase with some dissolved in the liquid, and the heavy hydrocarbons are mainly in the liquid phase which is in contact with the catalyst. However, both experimental work [47–49] and modelling studies [50,51] have shown that there is sufficient vaporisation of the hydrocarbons to impact apparent conversion in hydrocracking reactors. Therefore taking into account VLE seems mandatory. Earlier, we carried out VGO hydrocracking experiments in a laboratory scale batch reactor under conditions close to those found in industrial units [52,53]. In that work, cooled reactor contents and liquid phase samples were analysed by two-dimensional gas chromatography (GC  GC) which quantifies components by hydrocarbon family and carbon number. The experimental data obtained in [52,53] allow therefore a model of the lumped type to be constructed with more detail on the physical structure of the pseudo-components than has hitherto been possible. In this work, we construct a kinetic model which includes each pseudo-component identifiable by GC  GC. There are 217 discrete hydrocarbon lumps in total, of which 105 are present in the initial feedstock. As well as the reaction kinetics, the model takes into account internal and external mass transfer and VLE effects. In the next section of this article, the experimental method and results will be briefly summarised. The kinetic model will then be presented followed by a discussion of the results with parameter estimation. The mass balances for the reactor and the correlations used to estimate mass transfer coefficients are presented separately in Appendices A and B. 2. Materials and methods 2.1. Experimental set-up The model developed in this article is based on the experimental work and set-up described in detail by Henry et al. [52]. The reactor was a 300 mL Parr 4842 bench scale pressure reactor equipped with a Robinson–Mahoney stationary catalyst basket. It incorporates an on-line sampling system which allowed composition data to be collected from the hot reactor whilst in operation. The reactor was used to contact 120 g of hydrotreated VGO with 10 g of alumina supported Ni–Mo–S and zeolite bifunctional hydrocracking catalyst supplied by IFP Energies nouvelles. The operating conditions were defined to give a selectivity and conversion similar to a particular industrial unit and were 400 °C and 120 bar with the pressure maintained by hydrogen injection. The VGO had a boiling range of 370–510 °C and density of 0.8589 g cm3 at 20 °C which corresponds to carbon numbers in the range 22–42. There are minimal impurities due to the feedstock having been previously hydrotreated. 1.4 mL of Dimethyldisulphide (DMDS) decomposing to methane and hydrogen sulphide and 0.35 mL of aniline decomposing to ammonia are added. The reaction time was measured from the moment the reactor reached its operating temperature of 400 °C, after 23 min of heating. The final mass balance was generally accurate to 97.5% with a small amount of material lost as volatiles on opening the reactor and as liquid remaining on the reactor cover and agitator. Vapour samples were analysed by gas chromatography for hydrogen, hydrogen sulphide and light hydrocarbons: methane, ethane, propane, normal and iso-paraffins C4–C6, methylcyclopentane and cyclohexane. GC  GC analysis [54] was used for the feedstock and liquid products. Based on this, six hydrocarbon families are considered: n-paraffins (NP), iso-paraffins (IP), mononaphthenes (N1), naphthenes with two fused rings or more (N2), mono-aromatics (Aro1) and aromatics with two fused rings or more (Aro2). Component lumps are defined by hydrocarbon family and carbon number. According to this method, each NP is included

273

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284 Table 1 Reaction times and values for conversion and MD selectivity calculated from experimental data. Reaction time (min)

Conversion (%)

Selectivity (%)

0 28 68 158

0 ± 4.00 20.72 ± 3.17 37.18 ± 2.51 54.75 ± 1.81

– 84.69 ± 21.0 78.40 ± 11.8 67.92 ± 8.11

where mVGO;0 and mMD;0 represent the masses of VGO and MD respectively at the start of the experiment at operating conditions. In terms of carbon number, VGO is defined as the lumps C23 or greater for the IP and C22 or greater for the other hydrocarbon families. MD is defined as the lumps C10–C22 for the IP and C9–C21 for the other hydrocarbons. The maximum errors DX and DY are determined from Eqs. (3) and (4):

DX ¼  as a specific molecule. For the other hydrocarbon families the same definition gives discrete lumps of multiple isomers.

X¼

mVGO;0  mVGO mVGO;0

ð1Þ

Y¼

mMD  mMD;0 mVGO;0  mVGO

ð2Þ

1 jmVGO;0  mVGO j

 
mMD  mMD;0

ðjDmVGO;0 j þ jDmVGO jÞ  jDmMD j þ jDmMD;0 j þ

mVGO;0  mVGO


Fig. 1 gives an example of the GC  GC liquid composition analysis. The results shown in Fig. 1(a) and (b) are given for an experiment where the reaction has been stopped immediately the reactor reached operating conditions. The initial feed contained only hydrocarbons greater than C21, and a small amount of cracking has occurred during the heating step. A more detailed analysis was carried out for the initial feedstock and the mass fractions of Aro and N with more than one (fused) ring were determined. These are 71% and 39% respectively. Components with rings of each type

(b)

6%

5%

Mass fracon (wt%)

ð4Þ

Aromacs

t0

4%

3%

2%

6%

5%

Naphthenes Mass fracon (wt%)

(a)

1%

4%

t0 Isoparaffins N-paraffins

3%

2%

1%

0%

0%

5

8

11

14

17

20

23

26

29

32

35

38

41

5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41

Carbon number

Carbon number

6%

(d)

5%

Mass fracon (wt%)

Aromacs

4%

3%

2%

1%

6%

5%

Naphthenes

t = 158 min

t = 158 min Mass fracon (wt%)

(c)

ð3Þ

DY ¼ 

2.2. Experimental results The repeatability of the experiments was verified and experiments of increasing duration were carried out to observe the progression of the hydrocracking reaction. Table 1 shows the reaction times used and the conversion, X, and selectivity, Y, for MD, calculated from the following expressions:




mVGO DmVGO;0


þ jDmVGO j

mVGO;0 mVGO;0 1

4%

Isoparaffins N-paraffins

3%

2%

1%

0%

0%

5

8

11

14

17

20

23

26

29

Carbon number

32

35

38

41

5 7 9 11 13 15 17 19 21 23 25 27 29 31 33 35 37 39 41

Carbon number

Fig. 1. Liquid phase hydrocarbon distributions based on analysis of cooled reactor contents by family and carbon number: (a and b) immediately after reaching reaction conditions and (c and d) at 54.75% conversion.

274

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284

are classified as Aro. Fig. 1(c) and (d) shows the results for an experiment of 158 min duration, with many more hydrocracking products present. In the range of the initial feedstock, above C21, the two graphs show a similar pattern with reduced quantities but little change to the distribution shape. The Aro and N mass fractions are reduced more than for the NP and IP suggesting they are more reactive. For carbon numbers of C21 and below there is a marked difference between the two graphs. For the Aro and N there is a gap. Aro and N of size C18–C20 are only formed in tiny quantities and from C7 to C17 the general shape of the product distribution resembles that of the feedstock. There are few N and Aro products lighter than C9, and none of size C6 or smaller, as the minimum size is limited by the aromatic or naphthenic rings. For the paraffins, the product distribution is much less constrained. Comparison of the GC  GC results for experiments of different duration also showed that the number of aromatic and naphthenic rings remained constant. This is in accordance with the literature as ring opening is reported to be a difficult reaction, occurring only in the presence of noble metals [8–10]. These observations will form the basis for the hydrocracking kinetics used to construct the model. Note that no olefins were measured in the GC  GC analysis so, under the experimental conditions, hydrogenation must be very rapid. However, as the hydrogenation/dehydrogenation reactions are important steps in hydrocracking they are included in the model as intermediate products. The vapour phase analysis provided useful data for the lightest components measured. However, for C5 and C6 components there were significant losses to the liquid phase and also to the condenser on the vapour outlet line, so these results were not used. The analysis confirmed that negligible amounts of methane and ethane were produced by the hydrocracking reaction.

3. Model construction Fig. 2 shows the concentration profile for a typical pseudocomponent in the model through the four phases present (vapour, liquid and two solids, mechanically mixed). The catalyst is formed from the two solids: zeolite crystallites, which contain the Brönsted acid sites and alumina supported Ni–Mo–S. In the model these are considered as two separate phases. The cracking reaction, which is represented by a distributed approach described in Sections 3.1.2 and 3.1.3, occurs only in the zeolite and hydrogenation/dehydrogenation is supposed to occur only in the metal phase. Furthermore, as the amount of zeolite is much smaller than that of alumina supported catalyst, the crystallites of zeolite are mostly embedded within the body of supported catalyst granules. Therefore we assume that to reach the zeolite any component must pass through the supported catalyst phase. Mass transfer resistances in the solid and at the solid boundaries are included. The liquid and vapour phases are also considered separately. Vapour (+H2 supply)

1. Isothermal reactor. 2. Perfect mixing in liquid and vapour phases. 3. VLE represented by Peng Robinson equation of state with Grayson Streed activity coefficient model [55]. 4. No contact between vapour and catalyst. 5. No contact between the zeolite and the external liquid phase. 6. Mass transfer resistances within metal and zeolite phases of the catalyst can be represented by a linear driving force model [56]. Secondly, assumptions pertaining to the kinetic model are as follows: 1. Hydrocarbons can be grouped into 6 families: NP, IP, N1, N2, Aro1, Aro2. 2. In order to have a reasonable number of parameters, adsorption/desorption processes are not considered in the model. 3. Hydrogenation/dehydrogenation reversible reactions only occur in the metal phase of the catalyst. 4. The same hydrogenation and dehydrogenation rate parameters can be used for all lumps. 5. Hydrogenation/dehydrogenation of aromatic and naphthenic rings is not considered because it is not observed under our operating conditions. 6. Cracking only occurs in the zeolite phase. 7. Only dehydrogenated hydrocarbons undergo cracking. 8. No dehydrogenated hydrocarbons are present in the vapour phase. 9. All reactions are pseudo-first order. 10. All stoichiometric coefficients are +1 for products or 1 for reactants. 11. No ring opening occurs. 3.1. Kinetic model 3.1.1. Hydrogenation/dehydrogenation reactions The hydrogenation/dehydrogenation is said to be a reversible pseudo-first order reaction, occurring only on the active sites in the metal phase of the bifunctional catalyst. Only hydrogenation/ dehydrogenation of chains or branches is considered here. The model does not include hydrogenation or dehydrogenation of aromatic or naphthenic rings (assumption 5 in kinetic model). Indeed, the feedstock considered here is already hydrotreated and the ratio of naphthenes to aromatics in it is not far from equilibrium under HCK reaction conditions. For any lump, i (Aro1, Aro2, N1, N2, NP or IP), a dehydrogenated form, j (Aro1=, Aro2= N1=, N2=, NP= or IP=), can be produced:

Catalyst (alumina support)

Liquid

pi

The list of assumptions for the model has been divided into two groups. Firstly, those regarding heat and mass transfer are:

Catalyst (zeolite phase)

•NiMoS acve sites •Hydrogenaon/ Dehydrogenaon C i,S

Ci,L

•Bronsted acid sites •Isomerisaon/Cracking

C i,Z kLa

ksL

Vapour/Liquid interface at equilibrium

Liquid/Alumina interface

kSZ Alumina/Zeolite interface

Fig. 2. Concentration profile for a typical pseudo-component in the model.

275

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284 N1 toIP

Table 2 Hydrocracking incidence matrix for the model by hydrocarbon family. Parents Aro1= Products

NP¼ toNP¼

carbon number. For the paraffins, the constants kc Aro2=

N1=

U U

U U

U U

Aro1= Aro2 = N1= N2

Aro2= to IP =

N2 = to NP =

=

Aro2

N2

=

Aro2= to NP =

k

Aro2 = N1=

NP=

N2+=

NP= to NP=

kc

NP= IP=

k

N2= to IP =

= kc

Aro1=

kc

N2= to NP = c

N1=

kc

Aro1= to NP=

k

U U

N1= k cN1= to IP=

Aro1=

kc

N1= to NP= c

U U

N1 = to NP=

Aro2= to NP =

= kc

U U U

=

Aro1= to IP =

kc

Aro2= to IP =

kc

N1= to IP =

kc k

N2= to IP = c

IP =

2

3 ¼7 r 0Arox c  6. 7 6. 7 ¼
Arox¼ kArox¼ to c 4. 5

IP= to IP =

kc

Products are shown in boxes and are always smaller than parent molecules

¼17 r 0Arox c

Fig. 3. Hydrocracking scheme by product family.

kdehyd

i ¢ j þ H2

ð5Þ

khyd

The overall reaction rate of lump, i, to or from its dehydrogenated form, j, is given by Eq. (6). The hydrogenation and dehydrogenation rate constants, khyd and kdehyd , are assumed to be independent of hydrocarbon structure. The global hydrogen consumption rate is given by Eq. (9)

ri ¼

þ


hyd r ihyd i

ð6Þ

with

r idehyd

¼ kdehyd C i;s

ð7Þ

r ihyd ¼ khyd C j;s C H2 ;s

r H2

Aro1¼ toNP¼

kc

NP¼

Nx¼ to

þ kc

Aro2¼ toNP¼

kc

,

Aro1¼ toIP¼

kc

,

Aro2¼ toIP¼

kc

,

N1¼ toNP¼

kc

,

N2¼ toNP¼

kc

,

32 Nx¼ 3 C 20 Probability  7 76 IP¼ 6 6 7 4 matrix 54 ... 5 Nx¼ to Nx¼ C Nx¼ 42

ð11Þ 2

¼3 r 0NP c

6 6. 6 .. 4

¼39 r 0NP c

3

2

7 X 7 Nx to
Nx¼ kc ¼ 7¼ 5 x¼1;2

NP¼ 6 6

Nx¼ to NP¼ 2

X

32 C Nx¼ 3 20 7 76 7 6 . 76 . 7 54 . 5

Probability

4 matrix

Arox¼ toNP¼ 6 6


Arox¼ kc

Probability

¼ C Nx 42 32

4 matrix

x¼1;2

3 ¼ C Arox 20 7 6 76 . 7 76 . 7 54 . 5

Arox¼ to NP¼ 2

,

ð10Þ

2

¼17 r 0Nx c

k¼1

3.1.2. Cracking reactions In the model, the dehydrogenated hydrocarbons are the species which undergo cracking and this occurs only in the zeolite phase. The possible cracking products are shown by family in Table 2. An Aro= or N= must come from a larger Aro= or N= whereas an NP= or an IP= can come from any parent molecule. The corresponding schemes are shown in Fig. 3. There are two groups of parameters. For the Aro and N families the parameters



Probability

3 ¼7 r 0Nx c  6. 7 6. 7 ¼
Nx¼ kNx¼ to c 4. 5

þ ð9Þ

Arox¼ to IP¼

þ kc

2

ð8Þ

Nps   X dehyd k k ¼
hyd r dehyd H2 r hyd þ
H2

NP¼

2 3 3 C Arox¼ 20 6 7 7 6 76 .  4 matrix 7 56 .. 4 5 Arox Arox¼ to Arox¼ C 42 ¼ 2


dehyd r idehyd i

,

NP= to IP =

kc

NP= IP=

IP= to NP= c

N2

NP¼ toIP¼

, kc

IP toIP kc ¼ ¼

relate to the product mass balances. For any reaction, an overall rate constant is applied and the variation in reaction rate for components within that family is described using a probability matrix. Eqs. (10)–(13) give the formation rates for each hydrocarbon family. The product formation rates and parent concentrations for each lump are written as vectors with dimension equal to the number of lumps produced or consumed. The terms in the probability matrices are the probabilities of forming any product for any given lump. If a reaction is impossible, its probability is set to zero. To represent the cracking reactions as simply as possible, the rate equations are based on moles of carbon as this quantity is conserved. The minimum and maximum carbon numbers consumed and produced for each hydrocarbon family are included in Eqs. (10)–(13). These are discussed in the next section. The model structure is not identical for all the hydrocarbon families because there is more information about the Aro= and N=. They are not formed in hydrocracking so their source must be an Aro= or N=, whereas an NP= or IP= can come from any parent. So, to meet the requirement of reaction stoichiometry, the probability matrices used for the Aro= and N= in Eqs. (12) and (13) are found from Eqs. (10) and (11) which are valid for both mono(x = 1) and multi-ringed species (x = 2).

kc

kc

Aro2

IP=

U

Aro1= to NP= to IP =

NP=

U

kc

Aro1= k cAro1=

N2=

IP toNP¼ kc ¼ ,

U

Aro1= Aro2= N1= N2= NP= IP=

N2 toIP

kc ¼ ¼ , kc ¼ ¼ are typical reaction constants and Aro and N form a smaller Aro or N plus one or two paraffins to conserve the total

NP¼ to NP¼ 6 6

þ
NP¼ kc

Probability

4 matrix

32 C NP¼ 3 6 7 76 7 6 . 76 . 7 54 . 5

¼ C Arox 42

NP¼ to NP¼

IP to

þ
IP¼ kc ¼

¼ C NP 42 2 3 2 3 C IP¼ Probability 6 7 76 7 NP¼ 6 . 6 matrix 76 . 7 4 56 . 4 5 IP¼ IP¼ NP¼ C

42

ð12Þ

276

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284

3

2

¼4 r 0IP c 7 X 6. Nx to 7¼ 6.
Nx¼ kc ¼ 5 4. ¼39 r 0IP c

x¼1;2

2 IP¼ 6

Probability

4 matrix Nx¼ to IP¼

32 Nx¼ 3 C 20 7 76 . 7 6 54 .. 5

¼ C Nx 42 32 Arox¼ 3 C 20 Probability X 7 76 Arox¼ to IP¼ 6 . 7 þ
Arox¼ kc 4 matrix 56 5 4 .. x¼1;2 Arox¼ Arox¼ to IP¼ C 42 2 32 IP¼ 3 C6 Probability 7 76 IP toIP 6 . 7 þ
IP¼ kc ¼ ¼ 4 matrix 56 4 .. 5 IP¼ IP¼ to IP¼ C 42 2 32 NP¼ 3 C6 Probability 7 76 NP to IP¼ 6 . 7 6 ð13Þ þ
NP¼ kc ¼ matrix 4 54 .. 5 NP¼ to IP¼ C NP¼

2

42

3.1.3. Probability matrices As discussed above, the experimental results show that the hydrocarbon families behave differently, with more constraints on the Aro= and N= product distributions. Eq. (14) shows the probability matrix for Arox= to Arox=. The matrix for Nx= to Nx= is similar. The first important point is that only components from C20 to C42 react and only products from C7 to C17 for mono-species and from C10 to C17 for multi-species are formed. Referring back to Fig. 1(c), both Aro= and N= form two distinct and separate distributions during hydrocracking with some potential products not formed at all and it has already been established that ring opening can be neglected. To represent this behaviour in the model, the parent molecule is said to react in a single step and to give a final product in the range C7–C17. This imposes a considerable limitation, the consequences of which will be discussed later, in the results section.

2

3

Probability 6 7 4 matrix 5 Arox¼ to Arox¼ 2 pðArox¼20 ! Arox¼7 Þ 6 pðArox¼20 ! Arox¼8 Þ 6 ¼6 6 .. 4.

parameter, b, which allows the NP= and IP= products from cracking Aro= and N= to vary between two limits as shown in Fig. 4. For b = 0, each cracking reaction forms a single paraffin and for b = 1, two smaller paraffins are formed. The smaller paraffins are half the size of the single paraffin, or as close as possible for odd carbon numbers. The probability matrix for Arox= to NP= is shown Eq. (16). Eq. (17) shows how the probability matrix for b > 0 is constructed from the matrix with b = 0. 2

3 Probability 6 7 Ab ¼ 4 matrix 5 Arox¼ to NP¼ 2 pb ðArox¼20 ! NP¼3 Þ 6 6 .. 6. 6 6 6 6 pb ðArox¼20 ! NP¼13 Þ ¼6 6 6 60 6 6. 6 .. 4 0

0 pb ðArox¼21 .. . pb ðArox¼21 .. . 

 .. ! NP¼4 Þ . .. . .. ! NP¼14 Þ . .. . 0

0 .. . 0 pb ðArox¼42 .. . pb ðArox¼42

3 7 7 7 7 7 7 7 7 7 7 ! NP¼25 Þ 7 7 7 7 5 ! NP¼35 Þ ð16Þ

If b = 0, Ab ¼ A0 ,

A0 ði; jÞ ¼ p0 ðAro¼ðiþ19Þ ! NP¼ðjþ2Þ Þ If b > 0,

For j > 3; Ab ði; jÞ ¼ ð1  bÞA0 ði; jÞ þ bðA0 ði; 2j þ 1Þ þ 2A0 ði; 2j þ 2Þ þ A0 ði; 2j þ 3ÞÞ

ð17Þ

Since the paraffins C3–C5 don’t crack, for j from 1 to 3 we have:

For j ¼ 2; 3; Ab ði; jÞ ¼ A0 ði; jÞ þ bðA0 ði; 2j þ 1Þ þ 2A0 ði; 2j þ 2Þ þ A0 ði; 2j þ 3ÞÞ For j ¼ 1; Ab ði; 1Þ ¼ A0 ði; 1Þ þ bð2A0 ði; 4Þ þ A0 ði; 5ÞÞ

    pðArox¼42 ! Arox¼7 Þ     .. . . . .

pðArox¼42 ! Arox¼8 Þ .. .

3 7 7 7 7 5

ð14Þ

pðArox¼20 ! Arox¼17 Þ     pðArox¼42 ! Arox¼17 Þ A normalised gamma distribution, Eq. (15), requires two parameters, u and h, to be estimated and gives a reasonable match. ðnC  6Þ or ðnC  10Þ is used as the independent variable because the lightest measured Aro1 and N1 is of size C7. ðnC  6Þ or ðnC  10Þ corresponds to the Aro or N carbon number minus 6 carbons for Aro1 and N1 or 10 carbons for Aro2 and N2, respectively. This distribution is applied globally, so each Aro1 or N1 produces the same product distribution i’.

ðnC  6Þu1 eðnC6Þ=h pðAro¼i ! Aro¼ðnC6Þ Þ ¼ P17 u1 ðz6Þ=h e z¼7 ðz  6Þ

ð15Þ

For the Aro2 and the N2 the probability matrices are the same as for the Aro1 and N1 respectively with the probabilities of forming components smaller than C10 set to zero and the column totals readjusted to one. Complementary matrices for Arox= and Nx= to NP= and IP= are constructed using the elements of the matrices described above to complete the mass balances. For example, p0 ðArox¼21 ! NP¼13 Þ ¼ pðArox¼21 ! Aro¼8 Þ. Most Aro and N are likely to have a multibranched structure, so there may be a simultaneous removal of short branches from these components. To test this effect simply, the model includes a

For cracking the NP= and the IP=, the probability matrices used in the model are derived from the experimental data and the literature. The measured product distributions for the paraffins, given in Fig. 1(d), show that they form a continuous distribution with some over-cracking, defined as products lighter than MD (i.e. lighter than C8 or C9). So, it is possible that all the components behave in a similar way, with little limitation on the product range. Weitkamp published a product distribution for ideal hydrocracking of n-hexadecane [8]. In molar terms, the curve is symmetrical and relatively flat except for removal of small branches, where no C1 or C2 was formed and only small quantities of C3. For the model, this curve has been translated into a probability distribution for the products of cracking C16 NP and extrapolated to all the other reacting NP and IP. For parents smaller than C16, the larger products which cannot be formed were removed from the distribution. For parents greater than C16 the probability of forming any product of size C9 or larger is said to be the same as for C8, so for very big components a mostly flat probability distribution is produced. All the distributions are normalised, so the sum of each column in the probability density matrix is one. The matrix is then converted into moles of carbon, so as to be compatible with Eq. (12). The GC  GC analysis allows consideration of the molecular structure of the components within each lump. For the smaller IP lumps (up to nC = 11), many isomers may be less reactive due to their conformation. For example an IP with 8 carbons can’t have branches of size C3 or C4. These molecules would be much less likely to undergo cracking and would reduce the observed reactivity of their lump. In the model, the reactivity of these lumps is calculated by dividing the number of isomers with removable

277

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284

CH3 CH3 H3C

CH3

(1)

(2)

CH3

Fig. 4. Possible products from cracking aromatic/naphthenic molecules (1) a single paraffin (2) two paraffins of similar size.

Table 3 List of kinetic parameters used in the model. Aro product C distribution N product C distribution Fraction of Aro and N producing 2 NP or 2 IP

uAro, hAro uN, hN b Reaction constants kdehyd khyd

Dehydrogenation Hydrogenation Cracking Aro1=

Aro1¼

kc

Aro2 kc ¼ Aro c N1¼ kc N2 kc ¼ N c NP to kc ¼ NP to kc ¼ IP to kc ¼ IP to kc ¼

Table 4 Parameter estimation results.

Reaction constants for hydrogenation reactions kdehyd 1.98  105 ± 1.30  107 m3 s1 khyd 1.45  105 ± 3.26  107 m6 mol1 s1

Cracking Aro2=

a

Fraction of Aro producing NP=

Reaction constants for cracking (m3 s1)

Cracking N1= Fraction of N producing NP= Cracking NP= to NP=

NP¼ IP¼

Cracking NP= to IP=

NP¼

Cracking IP= to NP=

IP¼

Cracking IP= to IP=

and

uN, hN.

Concerning the C distributions, for N and Aro, the analytical distributions are given and adapted to our case. Each distribution has two parameters. For the paraffins, the distribution parameters are not estimated, the distributions are extrapolated from those of Weitkamp. For the lower carbon numbers these distributions are multiplied by the fraction of isomers with removable branches. To avoid including too many parameters the reaction rate parameters for cracking, both Aro and N are assigned the same fraction of NP to IP products as follows: Aro1¼ to NP¼

Aro1¼ to IP¼

kc

Aro1¼

¼ kc

Aro1¼

¼ kc

3.74  105 ± 2.03  104 0.658 ± 0.0036 5.58  105 ± 3.60  107 1.99  106 ± 4.05  107

a

3.1.4. Parameter estimation The number of estimated parameters is 17 (Table 3) and the number of observables is 438 corresponding to components, >C7. The lumps representing the lighter components,
kc

7.94  105 ± 3.62  107

Aro2 kc ¼ Aro c N1 kc ¼ N2 kc ¼ N c NP to kc ¼ IP to kc ¼ IP to kc ¼ NP to kc ¼

a

branches by the total. This is an important factor in the model, without which the reactivity of these lumps is too high and it becomes impossible to fit the measured data.

– Relating to the feed: b. – Relating to the C distribution: uAro, hAro – Relating to the kinetic constants

Aro1¼

kc

Cracking N2=

a

3.84 ± 1.88  102 1.66 ± 5.30  103 2.86 ± 6.70  103 1.62 ± 4.60  103 5.19  101 ± 5.70  103

uAro hAro uN hN b

 aAro c

ð18Þ

 ð1  aAro c Þ

ð19Þ

The kinetic model described above was combined with a model of the reactor system which takes into account mass transfer and VLE. See Appendices A and B for a description of the reactor model and the mass transfer calculations.

NP¼

0 ± 0.0036 4.83  106 ± 3.69  108

NP¼

4.24  106 ± 3.61  108

IP¼

4.00  106 ± 3.59  108

IP¼

6.39106 ± 3.63108

3.2. Resolution of equations The kinetic model was written in FORTRAN 90. The simultaneous differential equations representing the reactor system mass balances are integrated using the solver DDASPG with a time interval of 1 s. Model inputs are the total quantities of material measured at the moment the reactor reaches its operating temperature of 400 °C. For this, the data used was obtained by analysis of cooled reactor contents for an experiment of zero duration. The initial vapour and liquid phase concentrations are set to equilibrium, calculated by flash calculation using ZBREN. The initial concentrations inside the catalyst particle are set to be the same as for the liquid phase. The measured quantities of hydrogen fed to the reactor are input to the simulation continuously via the mass balance. Parameter estimation is carried out in FORTRAN 90 using the DBCLSF algorithm. The objective function to be minimised is the sum of the squares of the error (SSE) between the calculated and measured mass fractions of each hydrocarbon lump measured by GC  GC. The precision of the estimated parameters was calculated using the following assumptions. First, errors associated with two successive measurements were independent and centred and follow a normal distribution. Moreover, the precision of each parameter was calculated as follows:

^i  diag p

 1  1 H H r t0:975 < pi < p^i þ diag r t0:975 2 2

ð20Þ

278

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284

Mass Fraction (%)

(a)

60

N data Aro data

50

N model Aro model

40

30

20

10

0

0

10

20

30

40

50

60

50

60

Conversion (%)

Mass Fraction (%)

(b) 60

NP data IP data

50

IP model NP model

40

30

20

10

0

10

20

30

40

Conversion (%) Fig. 5. Comparison between experimental and calculated mass fraction by hydrocarbon family for (a) Aro and N (b) NP and IP.

  ^i is the estimate of pi , diag H2 1 plays the role of the variwhere p ance, r is the standard deviation and t 0:975 is the student variable which corresponds to 97.5% probability in the confidence interval [48]. Data sets from three separate experiments are used, covering three time points, 28, 68 and 158 min. Results from a sample taken from the hot reactor after 38 min reaction time were used to validate the model.

4. Results and discussion The values of the estimated parameters and their confidence limits are given in Table 4. The results seem to be consistent with the experimental observations. With these parameters the SSE is at a minimum of 1.4098  103 and we calculate the model to be statistically significant at the 0.01 level. We find that the parameters for aNc and kc ¼ are not significant to 95% and could be removed in a future version of the model. As would be expected, there is quite high correlation between the parameters used to describe the Aro and N distributions and also between the reaction constants for hydrogenation and dehydrogenation. On the other hand, there is very little correlation between the reaction rate parameters for cracking. Aro2

Fig. 5 shows the results for each hydrocarbon family as a whole. The calculated values correspond well with the experimental data. The data points at the highest conversion are for the experiment of 158 min duration. For the olefins, the calculated total mole fraction remains below 0.5% throughout the simulation, in line with experimental observations. Fig. 6 compares the model results to the GC  GC analysis for the Aro and N hydrocarbon lumps of size C6 and greater after 28, 68 and 158 min. The C distributions can be seen to provide a reasonable fit for the product mass fractions across the lumps and time period considered. As the reaction progresses, there is a tendency for the measured N and Aro product mass fractions to move slightly to the left compared to the calculated values, which suggests a small amount of cracking does occur within these populations. Moreover, the evolution of the N and Aro parent distributions suggests cracking of aromatics and naphthenes occurs more and more slowly with time. Fig. 7 shows the breakdown between the mono- and di N, with the mono- being converted much more rapidly. This is logical because the N carbon numbers C21–C30 include some isomers with two and three naphthenic rings and only short branches. These are particularly difficult for hydrocracking as the rate is likely to be limited by both their chemistry and mass transfer rate. The results for the paraffins, given in Fig. 8 are also quite good with the shapes of the predicted mass fraction curves resembling the experimental data well. The model tends to under predict the accumulation of IP C8 and C9 and after 158 min it slightly over predicts the amount of C10–C14 NP. This could possibly be due to the same value of the parameter, b, the fraction of N and Aro producing 2 paraffins, being used for all lumps. If the di-N have more branches for potential removal and, as a result, produce relatively more light IP products than the mono-N, we might see the effect observed in the graphs. Fig. 9 compares the measured and calculated values for conversion and selectivity for MD against time. The simulated conversion and selectivity follow the experimental data well. As the experiments are carried out in a fed batch reactor, there is relatively little product at low conversions. This leads to a very wide margin of error when the selectivity is calculated from the measured data. As the conversion increases the selectivity calculation becomes much more accurate. For the hydrocarbon gases, a precise comparison between the model outputs and the experimental data is difficult. Some light materials are dissolved in the liquid phase and not measured by GC  GC, moreover the sample method is more accurate with increased conversion. Table 5 compares the total calculated quantities of C3 and C4 in the reactor gas phase against the measured values. The model estimates the amount of C3 well and the value for C4 seems reasonable given the expected losses. Fig. 10 compares the calculated mole fraction of hydrogen in the reactor vapour samples against the experimental data and also the equilibrium values. These are not exactly the same as the values in the reactor as there is a condenser on the vapour outlet line upstream of the sample point. The model results are almost identical with the equilibrium values, which are close to the experimental data. The dip in hydrogen mole fraction towards the end of the time period is due to increasing quantities of volatile hydrocracking products. The model results were validated using samples taken from the hot reactor after 38 min reaction time. Table 6 shows the measured and calculated conversion and selectivity at this time point and Fig. 11 compares the measured and calculated mass fractions for the lumps. The results are seen to be in good agreement and confirm the estimated parameter values. The results for the estimated parameters of Eq. (15) show that h is almost identical for Aro and N whereas / is greater for the Aro.

279

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284

Mass fraction (%)

4

5

naphthenes t0

3 2 1 0 2

10

15

aromatics

20 25 30 Carbon number

35

1

Mass Fraction (%) Mass Fraction (%)

5

4

t=28 min

3 5 2 4

t=28 min

1 3

0 40 2

10

15

10

15

35

40

35

40

1

0

10

15

20

25

30

35

40 0

Carbon number

20

25

30

Carbon number 5

t=68 min

10 15 20 25 30 aromacs Carbon number

35

Mass Fraction Mass (%) Fraction (%)

naphthenes

Mass fraction (%)

20 25 30 Carbon number

4

GC x GC data Model results

t=158 min

3 5 2 4 1 3

0 40 2

10

15

10

15

20 25 30 Carbon number

35

40

35

40

1

10

15

20

25

30

35

40

0

Carbon number

20

25

30

Carbon number

Fig. 6. Comparison between model outputs and measured data for Aro and N at different times.

4

5

naphthenes t=0

Mass Fraction (%)

Mass Fraction (%)

5

3 2 1 0

10

15

20

25

30

35

4

2 1 0

40

t=28 min

3

10

15

Carbon number

t=68 min

3 2 1 0

25

30

5

Mass Fraction (%)

Mass Fraction (%)

5 4

20

10

15

20

25

30

Carbon number

35

40

Carbon number

35

40

4

di-N+ mono-N

t=158 min

3 2 1 0

10

15

20

25

30

35

Carbon number

Fig. 7. Calculated mass fractions of mono-N and N with more than one ring at different times.

40

280

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284

5

4

iso-paraffins

t0 Mass fraction (%)

3 2 1 0 2

10

15

20 25 30 Carbon number

35

n-paraffins

Mass Fraction (%) Mass Fraction (%)

5

t=28 min

4 3 5 2 4 1 3 0 2

10

15

10

15

35

1

1 0

10

15

20 25 30 Carbon number

35

40

0

20

25

5

3 2 1

10 15 20 25 30 n-paraffins Carbon number

35

Mass Fraction (%) Mass Fraction (%)

iso-paraffins

0 2

35

40

GC x GC data Model results

t=158 min

t=68 min 4

30

Carbon number

5

Mass fraction (%)

20 25 30 Carbon number

1

4 3 5 2 4 1 3 0 2

10

15

10

15

20 25 30 Carbon number

35

1

0

10

15

20

25

30

35

40

0

Carbon number

20

25

30

35

40

Carbon number

Fig. 8. Comparison between model outputs and measured data for paraffins at different times.

Table 5 Comparison of measured and calculated quantities of C3 and C4 (mol) in the gas phase.

100

Conversion (%) Selectivity (%) 80

60

40

Measured

Calculated

28 min C3 C4

9.40  104 1.96  103

1.18  103 5.17  103

68 min C3 C4

2.87  103 5.63  103

2.71  103 1.11  102

158 min C3 C4

8.99  103 1.85  102

6.88  103 2.57  102

20

0

20

40

60

80

100

120

140

160

time (min) Fig. 9. Comparison between measured and calculated conversion and selectivity for MD.

Both distributions are skewed towards the lighter products and the effect of increasing / is to move the Aro distribution towards the heavier lumps. This result can be explained by the feedstock composition, since the initial fraction of Aro2, which have more than one ring, is much greater than that of the Aro1. For the N the

opposite is true, and so, we can expect the N products to be lighter than for the Aro. Regarding the coefficient b, approximately half of the cracking reactions remove a single branch from the Aro and N and the rest remove two branches simultaneously (see Fig. 4). b is required to allow the model to calculate the observed product distributions. The model’s limitation in Aro and N product size, combined with the fact that they undergo cracking only once, means that the number of carbon atoms which must be removed in a single step is high and so, with the value of b set to zero, only long chain paraffins are removed. We could imagine a scenario of consecutive hydrocracking, where a paraffin is formed from an Aro or N, then itself cracks to form two smaller components. However,

281

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284

N2, is always lower than for the mono-ringed species, N1. This is coherent with the experiments, as we observe a decrease in the naphthenes cracking rate with conversion. For cracking Aro, the paraffin products are distributed between NP and IP which suggests that Aro seems to crack with no preliminary isomerization, although the order of magnitude of the reaction parameters is the same as for the N1. The reaction rate parameters for the Aro2 are not significant to 95% and this is because the lump concentrations are too low to expect good precision. Finally, the parameters for cracking the paraffins are all of the same order of magnitude, estimated with good precision. There is very little difference between the values except for NP to IP, which is greater than for the others. This could be because the model uses NP to IP to account for an isomerisation step, which is not included separately. The values of the constants for the paraffins are about 10 times less than for removal of the Aro and N branches which implies that the rings have a substantial effect on the hydrocracking rate and mechanisms. Probably, the main pathway is formation of tertiary carbenium cation on the cycle with further beta scission and transfer of charge towards an atom of branched chain. In this respect it would be instructive to study the exact chain branching distribution within the lumps of naphthenes and aromatics.

1 0.9 0.8

Mole fraction, yH2

0.7 0.6 0.5 0.4 0.3 0.2 0.1 0

Experimental data Equilibrium 0

20

40

60

80

100

120

140

160

Time (min) Fig. 10. Comparison between experimental and model values for mole fraction of H2 in the vapour phase samples and equilibrium [55].

Table 6 Measured and calculated conversion and selectivity after 38 min.

Conversion (%) Selectivity (%)

Measured

Calculated

23.6 72.4

25.1 68.0

5. Conclusions A many lumped model for hydrocracking kinetics has been constructed based on experiments carried out in a semi-batch reactor. The model is able to fit the experimental data quite well. The use of GC  GC analysis results allowed lumping by hydrocarbon family and carbon number which gives the model a greater emphasis on molecular structure than in a continuous approach. In particular, to fit the experimental data, it was necessary to separate the N and Aro into mono-ringed and multi-ringed lumps. It was also important to include parameters for the number of branches shed by the Aro and N and the reactivity of the light IP. This method maintains the positive aspects of lumped models, such as relative simplicity and a direct relationship with measured data, whilst incorporating more kinetic detail. Furthermore, the estimated pseudo-reaction constant values respect the molecular reactivities observed in studies of hydrocracking on model molecules. This work shows the power of using a relatively simple reactor in conjunction with GC  GC analysis as a step towards understanding the intrinsic reaction kinetics of hydrocracking.

5

5

4 naphthenes

4 iso-paraffins

3 2 1 0 2

10

15

20

25

30

35

aromatics

Mass Fraction (%) Mass Fraction (%)

Mass fraction (%)

when this case was tested, the paraffin product distributions calculated did not fit the measured data at all. Also, we don’t exclude the idea that consecutive cracking, with removal of shorter branches from the Aro and N, one after the other, may take place within the catalyst. For this, the reaction rate would be very high, so that the intermediate Aro and N products are formed but not observed. The hydrogenation reaction constant is of the same order as that for dehydrogenation. This agrees with the lack of olefins measured in the system as the relatively high hydrogen concentration drives down the quantity of the olefin lumps. The reaction constants for cracking show that no NP are produced from the N and the rate constant for cracking N1 is relatively high. This agrees with the rapid isomerisation of naphthene branches reported by Souverijns [10] which would lead to much more IP being produced. It is also coherent with our assumption of no ring opening as this would result in NP formation from N. We observe that the cracking rate for the multi-ringed species,

1 0

GC x GC data model results

3 5 2 4 1 3 0 2

10

15

20

25

30

35

20

25

30

35

n-paraffins 1

10

15

20

25

30

Carbon number

35

40

0

10

15

Carbon number

Fig. 11. Comparison between model outputs and measured data after 38 min.

40

282

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284

The authors would like to thank IFP Energies nouvelles, France, for supporting this work, and especially Gerhard Pirngruber for the fruitful discussions.

those used by Henry et al. [53] and the materials were hydrogen in diesel oil at 25 °C and 340 °C. It has not been validated by testing with alternative solvents but is used in the model as no alternative suitable correlation exists. All resistance to mass transfer is assumed to be on the liquid side of the interface.

Appendix A. Reactor model

Shi ¼

Acknowledgements

2

dC i;z ¼ F iðs!zÞ dt

ðA:1Þ

X X dC j;z r j!g þ rh!j ¼ F jðs!zÞ þ c c dt g¼1 h¼1 Nps

Nimp dimp ql 2

For the dehydrogenated lumps in the zeolite phase,

Vz

ðB:1Þ

with

Component material balances are as follows: For each hydrogenated lump in the zeolite phase,

Vz

ðkL aÞi dimp ¼ 3:15:106 Re0:95 We1:55 Sc0:8 i DiB

Nps

ðA:2Þ

In the reaction rate terms in Eqs. (A.2) and (A.3) to (A.5), the stoichiometric coefficients are signed and correspond respectively to formation and consumption of lump j. The mass balances for component i, the corresponding dehydrogenated lump, j, and hydrogen in the metal phase are given by Eqs. (A.3)–(A.5):

Re ¼

ll

N2imp dimp ql ll ; Sci ¼ sl ql DiB 3

; We ¼

ðB:2Þ

Interface concentrations were estimated for each component from repartition coefficients found from flash calculations in Prosim Plus (version 3.3) with the Peng Robinson equation of state and Grayson Streed activity coefficient model [55]. Repartition coefficients were obtained for several sets of experimental data at different, known conversions and interpolation was used to find intermediate values. Diffusion coefficients for each component diffusing through the reaction mixture are estimated using the modified Wilke-Chang correlation given in Eq. (B.3) with the molar volume, v, in (cm3 mol1) [58].

7:4  1012 ð£M B Þ1=2 T

Vs

dC i;s ¼ F iðl!sÞ  F iðs!zÞ þ r i dt

ðA:3Þ

DiB ¼

Vs

dC j;s ¼ F jðl!sÞ  F jðs!zÞ þ r j dt

ðA:4Þ

The molar volume of each component at its normal boiling point is estimated by the Schroeder method [58] where 7⁄ is a reduction for each aromatic/naphthenic ring:

Vs

dC H2 ;s ¼ F H2 ðl!sÞ  F H2 ðs!zÞ þ rH2 dt

ðA:5Þ

The mass balance for the liquid phase,

dðV l C i;l Þ dC i;l dV l ¼ Vl þ C i;l ¼ F iðl!v Þ  F iðl!sÞ dt dt dt

ðA:6Þ

For the dehydrogenated components in the liquid phase the mass balance is identical to Eq. (A.6) with F jðl!v Þ = 0. Total liquid volume is calculated from the molar volumes and changing concentrations of the individual components:

X dV l v iCi 1 dt i¼1 Nps

!

¼ Vl

Nps X

vi

i¼1

dC i dt

ðA:7Þ

Total reactor volume is constant so,

ðB:3Þ

lB v 0:6 i

v i ¼ 7ðnC i þ nHi þ nDBi Þ  7i

ðB:4Þ

The viscosity is estimated by Glaso’s correlation as recommended by Korsten and Hoffmann [59]. The API used is that of the feedstock (calculated as API = 141.5/relative density – 131.5). Variation in the liquid viscosity due to changing liquid composition is small and not included in the model.

l ¼ 3:141  1010 ðT  460Þ3:444 ½log10 ðAPIÞ a

ðB:5Þ

a ¼ 10:313½log10 ððT  460ÞÞ  36:447

ðB:6Þ 00 ksL;i ,

A combined mass transfer coefficient for each lump, is used to represent the mass transfer resistance due to the catalyst particle – liquid phase boundary, ksL;i , and the resistance within the 0 pores of the catalyst metal phase, ksL;i .

Vr ¼ Vl þ Vv

ðA:8Þ

F iðl!sÞ ¼ ksL;i Ss ðC i;l  C i;s Þ

dV l dV v ¼ dt dt

ðA:9Þ

For each component lump, the combined mass transfer coefficients are found by summing the two mass transfer resistances:

The mass balances for the hydrogen and the other components in the vapour phase are given by Eqs. (A.10) and (A.11) respectively:

dðV v C H2 ;v Þ dC H2 ;v dV v ¼ Vv þ C H2 ;v ¼ F H2 ðl!v Þ þ F H2 ;in dt dt dt

ðA:10Þ

dðV v C i;v Þ dC i;v dV v ¼ Vv þ C i;v ¼ F iðl!v Þ dt dt dt

ðA:11Þ

Appendix B. Mass transfer coefficients for reactor model The gas/liquid mass transfer coefficient, kLa, is found for the H2 and each lump, i, from the correlation determined by Fongarland [57] and given in Eqs. (B.1) and (B.2) is used. The reactor, agitator and catalytic basket used to find this correlation were identical to

00

1 1 1 þ 0 00 ¼ ksL;i ksL;i ksL;i

ðB:7Þ

ðB:8Þ

The solid–liquid boundary mass transfer coefficient, ki,sL, is estimated using Magnico & Fongarland’s correlation [60], Eq. (B.9), for n-heptane at standard conditions in a 300 mL Parr reactor with Robinson Mahoney stationary catalyst basket. dks is the diameter of the sphere of equivalent surface to the catalyst and Reks depends on the superficial velocity within the catalyst basket. The value of Ss is adjusted to allow for zeolite crystallites at the catalyst surface. 0

Shi ¼

ksL;i dks 0:4 ¼ 1:1Re0:41 ks Sc i DiB

ðB:9Þ

To estimate the effects of the resistance to mass transfer in the 0 metal phase, ki;sL , a modified Gluekauf correlation is used. The correlation relates the mass transfer rate across the particle boundary

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284

to the difference between external concentration and average concentration in the catalyst pellet. Gluekauf develops the relation for a spherical particle [61] and Kim extends it to Eq. (B.10), for an infinite cylinder [62]. This is used in the model, as it is closer to the shape of the hydrocracking catalyst:

dC i;s 8De;i ¼  2 ðC i;l  C i;s Þ dcat dt

ðB:10Þ

2

0

Comparing with Eq. (B.7) it is deduced that the term, ksL;i in Eq. (B.11), can be applied at the liquid–solid boundary as a mass transfer coefficient accounting for resistance within the catalyst metal phase.

4De;i K 0sL;i ¼ d 

ðB:11Þ

cat

2

The effective diffusivity for hydrocarbons in the pores of the catalyst metal phase is estimated using Eq. (B.12). The measured value of the porosity, es, is adjusted to allow for the volume of zeolite crystallites within the particle and the tortuosity, s, is estimated as 3. The remaining term accounts for the hindering effects of small pore sizes [63], with dp the average pore diameter.

De;i ¼

Di es



ss

1

2rh;i dp

4 ðB:12Þ

The hydrodynamic radius, rh, is found from the Stokes–Einstein relation:

Di ¼

kB T 6plr h;i

ðB:13Þ

Mass transfer resistance within the zeolite is represented using the Gluekauf relation for a spherical particle:

F iðs!zÞ ¼ ksz;i Sz ðC i;s  C i;z Þ

ðB:14Þ

5De;i ksz;i ¼ d 

ðB:15Þ

cat

2

The zeolite structure is designed to facilitate mass transfer and contains both mesopores and micropores. Effective diffusivities in the mesopores were calculated using Eq. (B.15) and for the micropores a study was carried out using molecular dynamics. Self-diffusivity of different families of hydrocarbons from C4 to C40 in all-silica faujasite was determined by molecular dynamics using Cerius2 v. 4.7 software from Accelrys Inc. on a 32-node Linux workstation. The cff91 and Universal forcefields were applied. Hydrocarbon molecules were packed into an isosurface enclosed volume, formed by van der Waals spheres of zeolitre atoms. In order to avoid artifacts created by periodic boundary conditions for large hydrocarbon molecules, a faujasite supercell of 5  5  5 was taken for simulations. The simulations were carried out for times from 100 ps to 5 ns, with 1 fs step, using NPT ensemble at 673 K and 100 bar. Prior to diffusivity simulations the system was optimised and annealed in five cycles, 10 ps each, with temperature rising from 300 to 700 K. Flexibility of the zeolite framework was allowed to obtain a realistic picture. Mean square displacement of molecules was plotted against time to determine diffusivity. The resulting diffusion coefficients under the reaction conditions show values comparable to those of the corresponding liquids. Self-diffusion in the large channels of the faujasite is relatively rapid even for C30 and C40 molecules. The calculated values for NPs decrease from 1.15  108 m2 g1 for C4 to 0.035  108 m2 g1 for C40. The differences in diffusivity between various families of hydrocarbons are significant, but less than its variation as a function of nC. The diffusivity of IPs and Ns is similar

283

to that of NPs with the same nC, whereas it decreases for Aros and diaromatics. Thus the values for NP-C30, IP-C30, N-C30, Aro-C30 and diAro-C30 are respectively 0.12  108, 0.11  108, 0.15  108, 0.05  108 and 0.035  108 m2 g1. Here we give only a general estimate of the diffusivity. The results of simulations will be discussed in more detail elsewhere. The results discussed above showed that the values for the micropores were of a similar order to those in the mesopores. However, the distance for mass transfer within the micropores is much smaller. So, most of the resistance is in the mesopores and, in the model, mass transfer in the micropores is neglected.

References [1] I.E.A., Key World Energy Statistics, 2014. [2] M. Steijns, G.F. Froment, Hydroisomerization and hydrocracking. 2. Product distributions from n-decane and n-dodecane, Ind. Eng. Chem. Prod. Res. Dev. 20 (1981) 654–660. [3] G. Cui, J. Wang, H. Fan, X. Sun, Y. Jiang, S. Wang, D. Liu, J. Gui, Towards understanding the microstructures and hydrocracking performance of sulfide Ni–W catalysts: effect of metal loading, Fuel Process. Technol. 92 (2011) 2320– 2327. [4] M. Steijns, G.F. Froment, Hydroisomerization and Hydrocracking. 3. Kinetic analysis of rate data for n-decane and n-dodecane, Ind. Eng. Chem. Prod. Res. Dev. 20 (1981) 660–668. [5] M.A. Baltanas, H. Vansina, G.F. Froment, Hydroisomerization and hydrocracking. 5. Kinetic analysis of rate data for n-octane, Ind. Eng. Chem. Prod. Res. Dev. 22 (1983) 531–539. [6] G.F. Froment, Kinetics of the hydroisomerization and hydrocracking of paraffins on a platinum containing bifunctional Y-zeolite, Catal. Today 1 (1987) 455–473. [7] M. Sau, K. Basak, U. Manna, M. Santra, R.P. Verma, Effects of organic nitrogen compounds on hydrotreating and hydrocracking reactions, Catal. Today 109 (2005) 112–119. [8] J. Weitkamp, Catalytic hydrocracking – mechanisms and versatility of the process, ChemCatChem 4 (2012) 292–306. [9] V. Calemma, M. Ferrari, S. Rabl, J. Weitkamp, Selective ring opening of naphthenes: from mechanistic studies with a model feed to the upgrading of a hydrotreated light cycle oil, Fuel 111 (2013) 763–770. [10] W. Souverijns, R. Parton, J.A. Martens, G.F. Froment, P.A. Jacobs, Mechanism of the paring reaction of naphthenes, Catal. Lett. 37 (1996) 207–212. [11] J.W. Thybaut, C.S. Laxmi Narasimhan, J.F. Denayer, G.V. Baron, P.A. Jacobs, J.A. Martens, G.B. Marin, Acid-metal balance of a hydrocracking catalyst: ideal versus nonideal behaviour, Ind. Eng. Chem. Res. 44 (2005) 5159–5169. [12] K.P. de Jong, J. Zecevic, H. Friedrich, P.E. de Jongh, M. Bulut, S. van Donk, R. Kenmogne, A. Finiels, V. Hulea, F. Fajula, Zeolite Y crystals with trimodal porosity as ideal hydrocracking catalysts, Angew. Chem. Int. Ed. 49 (2010) 10074–10078. [13] G.G. Martens, G.B. Marin, Kinetics for hydrocracking based on structural classes: model development and application, AIChE J. 47 (7) (2001) 1607– 1622. [14] J. Ancheyta, S. Sanchez, M. Rodriguez, Kinetic modeling of hydrocracking of heavy oil fractions: a review, Catal. Today 109 (2005) 76–92. [15] V.W. Weekman, D.M. Nace, Kinetics of catalytic cracking selectivity in fixed, moving and fluid bed reactors, AIChE J. 16 (3) (1970) 397–404. [16] J.F. Mosby, R.D. Buttke, J.A. Cox, C. Nikolaides, Process characterisation of expanded-bed reactions in series, Chem. Eng. Sci. 41 (4) (1986) 989–995. [17] M.A. Callejas, M.T. Martinez, Hydrocracking of a Maya residue. Kinetics and product yield distributions, Ind. Eng. Chem. Res. 38 (1999) 3285–3289. [18] C. Botchwey, A.K. Dalai, J. Adjaye, Product selectivity during hydrotreating and mild hydrocracking of bitumen-derived gas oil, Energy Fuels 17 (2003) 1372– 1381. [19] S. Sanchez, M.A. Rodriguez, J. Ancheyta, Kinetic model for moderate hydrocracking of heavy oils, Ind. Eng. Chem. Res. 44 (2005) 9409–9413. [20] J. Singh, M.M. Kumar, A.K. Saxena, S. Kumar, Reaction pathways and product yields in thermal cracking of vacuum residues: a multi-lump kinetic model, Chem. Eng. J. 108 (2005) 239–248. [21] L. Pellegrini, S. Locatelli, S. Rasella, S. Bonomi, V. Calemma, Modeling of Fischer–Tropsch products hydrocracking, Chem. Eng. Sci. 59 (2004) 4781– 4787. [22] M.T. Martinez, A.M. Benito, M.A. Callejas, Thermal cracking of coal residues: kinetics of asphaltene decomposition, Fuel 76 (9) (1997) 871–877. [23] S. Sanchez, J. Ancheyta, Effect of pressure on the kinetics of moderate hydrocracking of Maya crude oil, Energy Fuels 21 (2007) 653–661. [24] J. Martinez, J. Ancheyta, Kinetic model for hydrocracking of heavy oil in a CSTR involving short term catalyst deactivation, Fuel 100 (2012) 193–199. [25] J.-M. Schweitzer, S. Kressmann, Ebulliated bed reactor modeling for residue conversion, Chem. Eng. Sci. 59 (2004) 5637–5645. [26] S. Sadighi, A. Ahmad, M. Rashidzadeh, 4-Lump kinetic model for vacuum gas oil hydrocracker involving hydrogen consumption, Korean J. Chem. Eng. 27 (4) (2010) 1099–1108.

284

B. Browning et al. / Chemical Engineering Journal 284 (2016) 270–284

[27] T.S. Nguyen, M. Tayakout-Fayolle, N. Ropars, C. Geantet, Hydroconversion of an atmospheric residue with a dispersed catalyst in a batch reactor: kinetic modelling including vapor–liquid equilibrium, Chem. Eng. Sci. 94 (2013) 214– 223. [28] S.M. Jacob, B. Gross, S.E. Voltz, V.W. Weekman Jr., A lumping and reaction scheme for catalytic cracking, AIChE J. 22 (4) (1976) 701–713. [29] I. Pitault, D. Nevicato, M. Forissier, J.-R. Bernard, Kinetic model based on a molecular description for catalytic cracking of vacuum gas oil, Chem. Eng. Sci. 49 (24) (1994) 4249–4262. [30] R. Krishna, A. Saxena, Use of an axial-dispersion model for kinetic description of hydrocracking, Chem. Eng. Sci. 44 (3) (1989) 703–712. [31] J. Sun, D. Li, R. Yao, Z. Sun, X. Li, W. Li, Modeling the hydrotreatment of full range medium temperature coal tar by using a lumping kinetic approach, React. Kinet. Mech. Catal. (2014), http://dx.doi.org/10.1007/s11144-014-0791-2. [32] B.E. Stangeland, A kinetic model for the prediction of hydrocracker yields, Ind. Eng. Chem. Process Des. Dev. 13 (1) (1974) 71–76. [33] C.S. Laxminarasimhan, R.P. Verma, Continuous lumping model for simulation of hydrocracking, AIChE J. 42 (9) (1996) 2645–2653. [34] A. Arefi, F. Khorasheh, F. Farhadi, Application of a continuous kinetic model for the hydrocracking of vacuum gas oil, Pet. Sci. Technol. 32 (2014) 2245–2252. [35] J.R. Hernandez-Barajas, R. Vazquez-Roman, Ma.g. Felix-Flores, A comprehensive estimation of kinetic parameters in lumped catalytic cracking reaction models, Fuel 88 (2009) 169–178. [36] I. Elizalde, M.A. Rodriguez, J. Ancheyta, Modeling the effect of pressure and temperature on the hydrocracking of heavy crude oil by the continuous kinetic lumping approach, Appl. Catal. A: General 382 (2) (2010) 205–212. [37] H.M.S. Lababidi, F.S. AlHumaidan, Modeling the hydrocracking kinetics of atmospheric residue in hydrotreating processes by the continuous lumping approach, Energy Fuels 25 (5) (2011) 1939–1949. [38] M. Adam, V. Calemma, F. Galimberti, C. Gambaro, J. Heiszwolf, R. Ocone, Continuum lumping kinetics of complex reactive systems, Chem. Eng. Sci. 76 (2012) 154–164. [39] G. Li, Y. Xia, W. Zeng, Kinetic mechanism research of an industrial hydrocracker based on strict calculation of stoichiometric coefficients, Fuel 103 (2013) 285–291. [40] T. Jansen, D. Guerry, E. Leclerc, M. Ropars, M. Lacroix, C. Geantet, M. TayakoutFayolle, Simulation of petroleum residue hydroconversion in a continuous pilot unit using batch reactor experiments and a cold mock-up, Ind. Eng. Chem. Res. 53 (41) (2014) 15852–15861. [41] K. Basak, M. Sau, U. Manna, R.P. Verma, Kinetics of bitumen-derived gas oil upgrading using a commercial NiMo/Al2O3 catalyst, Can. J. Chem. Eng. 82 (3) (2004) 478–487. [42] P.J. Becker, B. Celse, D. Guillaume, H. Dulot, V. Costa, Hydrotreatment modeling for a variety of VGO feedstocks: a continuous lumping approach, Fuel 139 (2015) 133–143. [43] G.F. Froment, Single event kinetic modeling of complex catalytic processes, Catal. Rev. 47 (2005) 83–124. [44] J.W. Thybaut, G.B. Marin, Single-event microkinetics: catalyst design for complex reaction networks, J. Catal. 308 (2013) 352–362. [45] T. Zhang, C. Leyva, G.F. Froment, Vacuum gas oil hydrocracking on NiMo/USY zeolite catalysts. Experimental study and kinetic modelling, Ind. Eng. Chem. Res. 54 (2015) 858–868.

[46] L.M. Chavez, F. Alonso, J. Ancheyta, Vapor–liquid equilibrium of hydrogen– hydrocarbon systems and its effects on hydroprocessing reactors, Fuel 138 (2014) 156–175. [47] G.D. Bellos, N.G. Papayannakos, The use of a three phase microreactor to investigate HDS kinetics, Catal. Today 79–80 (2003) 349–355. [48] G. Hoekstra, The effects of gas-to-oil rate in ultra low sulfur diesel hydrotreating, Catal. Today 127 (2007) 99–102. [49] I. Rossetti, C. Gambaro, V. Calemma, Hydrocracking of long chain linear paraffins, Chem. Eng. J. 154 (2009) 295–301. [50] J. Chen, V. Milgundmath, N. Wang, Accounting for vapor–liquid equilibrium in the modeling and simulation of a commercial hydrotreating reactor, Ind. Eng. Chem. Res. 50 (2011) 1571–1579. [51] L. Pellegrini, S. Gamba, V. Calemma, S. Bonomi, Modelling of hydrocracking with vapour–liquid equilibrium, Chem. Eng. Sci. 63 (2008) 4285–4291. [52] R. Henry, M. Tayakout-Fayolle, P. Afanasiev, F. Couenne, G. Lapisardi, L.J. Simon, V. Souchon, Methodology for the study of vacuum gas oil hydrocracking catalysts in a batch reactor – coupling of GC-2D data with vapour–liquid equilibrium, material sciences and technology, 1 & 2, book series, Adv. Mater. Res. 560–561 (2012) 207–213. [53] R. Henry, M. Tayakout-Fayolle, P. Afansiev, C. Lorentz, G. Lapisardi, G. Pirngruber, Vacuum gas oil hydrocracking performance of bifunctional Mo/Y zeolite catalysts in a semi-batch reactor, Catal. Today 220–222 (2013) 159– 167. [54] J. Stihle, D. Uzio, C. Lorentz, N. Charon, J. Ponthus, C. Geantet, Detailed characterization of coal derived liquids from direct coal liquefaction on supported catalysts, Fuel 95 (2012) 79–87. [55] B. Browning, R. Henry, P. Afanasiev, G. Lapisardi, G. Pirngruber, M. TayakoutFayolle, Vapor–liquid equilibrium of hydrogen, vacuum gas oil and middle distillate fractions, Ind. Eng. Chem. Res. 53 (2014) 8311–8320. [56] D. Leinekugel-le-Coq, M. Tayakout-Fayolle, Y. Le Gorrec, C. Jallut, A double linear driving force approximation for non-isothermal mass transfer modeling through bi-disperse adsorbents, Chem. Eng. Sci. 62 (2007) 4040–4053. [57] P. Fongarland, Etude Cinetique de l’Hydrodesulfuration de Composes Soufres de Gazoles en Reacteur Continu Parfaitement Agité, 2003, These UCBL1 No d’ordre 123-2003. [58] B.E. Poling, J.M. Prausnitz, J.P. O’Connell, The Properties of Gases and Liquids, fifth ed., McGraw-Hill, New York, 2004, http://dx.doi.org/10.1036/ 0070116822. [59] H. Korsten, U. Hoffmann, Three-phase reactor model for hydrotreating in pilot trickle-bed reactors, AIChE J. 42 (5) (1996) 1350–1360. [60] P. Magnico, P. Fongarland, CFD simulations of two stirred tank reactors with stationary catalytic basket, Chem. Eng. Sci. 61 (2006) 1217–1236. [61] E. Glueckauf, Theory of chromatography. Part 10.—Formulæ for diffusion into spheres and their application to chromatography, Trans. Faraday Soc. 51 (1955) 1540–1551. [62] D.H. Kim, Linear driving force formulas for unsteady-state diffusion and reaction in slab, cylinder and sphere catalyst, AIChE J. 55 (3) (2009) 834–839. [63] J.L. Anderson, J. Quinn, Restricted transport in small pores – a model for steric exclusion and hindered particle motion, Biophys. J. 14 (1974) 130–150.