- Email: [email protected]
Thermal runaway analysis of a three-phase reactor for LCO hydrotreatment
Chemical Engineering Science 65 (2010) 313 -- 321
Contents lists available at ScienceDirect
Chemical Engineering Science journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / c e s
Thermal runaway analysis of a three-phase reactor for LCO hydrotreatment Jean-Marc Schweitzer ∗ , Clementina López-García, Daniel Ferré Institut Français du Pétrole, Process Development and Engineering Division, PO Box 3, 69360 Solaize, France
A R T I C L E
I N F O
Article history: Received 30 June 2008 Received in revised form 1 July 2009 Accepted 14 July 2009 Available online 21 July 2009 Keywords: Thermal stability Gasoils hydrotreatment Packed bed Dynamic simulation Light cycle oil
A B S T R A C T
Safety is a high-priority topic for the chemical industry to minimize the frequency and severity of accidents while keeping the productivity and quality of the production. The processes that may undergo thermal runaways due to exothermic reactions are at the heart of the risks of accidents. The study of such highly reactive systems is essential to achieve a safe and productive operation of existing processes and to ensure inherently safe new designs. It is well known that the stationary analysis with the van Heerden criterion must be satisfied, however, this is not sufficient to ensure reactor stability. Only dynamic analysis can provide an accurate answer concerning the safe operation of the reactor. Most of the reported stability studies are carried out for relatively simple systems (pseudo-homogeneous models, simple reaction schemes). This work presents the dynamic thermal stability study of a refining process (hydrotreatment of light cycle oils) carried out with real gasoils at industrial operating conditions. A 1D dynamic model that accurately represents the reactive system (gas–liquid–solid) was developed and validated with experimental pilot plant data. This mathematical model was used to perform the thermal stability analysis of the dynamic system using a perturbation method. The effect of the variation of the heat transfer coefficient on the thermal stability is presented. The spectral analysis of the eigenvalues indicating the stable/unstable behavior of the reactive system was compared with dynamic simulations. An excellent agreement was found between the simulations and the stability analysis. The case of oscillating behavior is also described. The frequencies of oscillation determined by the stability analysis are compared to the frequencies calculated by Fourier transform applied to the simulated signal. The reactor behavior and the oscillations features are accurately predicted with this stability analysis method. © 2009 Elsevier Ltd. All rights reserved.
1. Introduction One of the security priorities involving highly reactive systems is the risk of thermal runaway. The temperature increase for the reactions that follow an Arrhenius law, induces the rise of the heat generation that further increases the reaction temperature and so on; this situation may result in a thermal runaway. The consequences of a runaway can be: the early deactivation of the catalyst (coking), the loss of selectivity, conversion or the operability of the unit, the onset of secondary reactions and in some cases even the reactor explosion. The design of industrial reactors must rely on an accurate thermal stability study. Indeed, the design of the whole reactive system (reactor and the devices used for heat input/removal) should ensure a priori the thermally stable operation of the reactor. It is also important to determine the regions of operating conditions where the reactor has an unreliable behavior (runaway regions). In practice,
∗ Corresponding author. Tel.: +33 4 78 02 25 98; fax: +33 4 78 02 20 08. E-mail address: [email protected] (J.-M. Schweitzer). 0009-2509/$ - see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2009.07.012
most of the industrial designs of reactive systems are based only on the van Heerden (1953) stationary criterion. This criterion imposes that the slope of heat generated by the reactions (dQgen /dT) must be lower than the slope of the heat extracted from the system (dQexc /dT). Since this comparison is carried out only under stationary conditions the dynamic behavior of the system is not taken into account and the thermal stability cannot be guaranteed. Indeed, only the dynamic analysis will provide an accurate answer concerning the safe operation of the reactor. An extensive number of literature works have tackled the reactors thermal stability topic. The objective is to dispose of theoretical methods mainly to determine the boundaries of runaway/nonrunaway. There are two foremost approaches for the thermal stability assessment: parametric sensitivity and application of dynamic systems theory. Parametric sensitivity consists of the determination of the operation regions where a small change of the inlet parameters produces high changes of the variables of the system such as temperature. This concept was introduced in the field of chemical reactors by Bilous and Amundson (1956). Since then, many authors have extensively applied the parametric sensitivity theory to study different reactive systems. Among the studies
314
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
concerning packed bed catalytic reactors we can cite the work of Morbidelli and Varma (1986, 1987), Balakotaiah and Luss (1991), Balakotaiah et al. (1995). These authors present the study of the role of inter and intraparticle resistances in the criteria to determine the runaway boundaries as well as the role of the interphase gradients. Wu et al. (1999) used the parametric sensitivity to determine the critical ignition conditions to achieve the steady-state auto-thermal operation of reactors with reverse flow. Concerning the dynamic stability analysis, the development of an intrinsic criterion based on the sum of Lyapunov exponents is presented in several literature works (Alós et al., 1996; Zaldívar et al., 2003) where the instabilities are characterized with chaotic attractors. These studies have been performed in batch and semi-batch reactors. The perturbation stability analysis ant the chaos theory has more recently been extended to deal with the determination (1D or 2D) hot spots in packed-bed reactors (Viswanathan and Luss, 2006; Luss and Sheintuch, 2005; Agrawal et al., 2007). Most of the models study the instabilities in homogeneous or pseudo-homogeneous systems which are not the case of this work. The objective of this study is to achieve the dynamic thermal stability analysis of a complex refining process. The hydrotreatment (HDT) of an industrial light cycle oil (LCO) gasoil in a fixed bed reactor is used as a case study. Experiments were carried out at industrial operating conditions in a representative pilot plant. A complete dynamic model of the three-phase system was developed, fitted and validated on experimental data. The dynamic simulation model was used to assess the thermal stability analysis. This work presents this reliable stability method. The effect of the variation of the heat transfer coefficient is presented. In the case of unstable behavior the significance of the eigenvalues is also illustrated (oscillations). In the case of oscillations occurrence, their frequency is determined. 2. LCO hydrotreatment process 2.1. Introduction Commercial diesels are constituted of a blend of gasoils coming from different refining processes or operations. These different gasoils can be among others: straight run (SR) issued from the atmospheric distillation, coker gasoil produced in coking units, light cycle oil (LCO) issued from fluid catalytic cracking (FCC), hydrocracking (HDC). These gasoils must follow a hydrotreatment (HDT) process to meet commercial specifications. The main hydrotreatment reactions are: hydrogenation of aromatics and olefins, hydrodesulfurisation of organic sulfur compounds and hydrodenitrogenation. The chemical composition of the feedstocks depends directly on the provenance of the crude oil as well as the process from which it is produced and its severity (SR, LCO, coker . . . ). Table 1 presents an example of some characteristics of a SR gasoil and a LCO. The current European onroad specification of gasoils (NF EN 590) is also indicated. As illustrated, the aromatics content is completely different between both gasoils. The LCO are generally characterized by very high aromatics contents that result in a very low cetane number and high density. The aromatics and olefins content in LCO gasoils is reduced by the Table 1 Comparison of some physico-chemical characteristics of a SR and a LCO gasoil with the European on-sroad gasoil specification (NF EN 590). Description 3
Density 15/4 (kg/m ) Olefins (wt%) Sulfur (wt ppm) Total aromatics (wt%) Diaromatics+ (wt%) Cetane index
SR gasoil
LCO A gasoil
NF EN 590 specification
857 – 11 400 32 19 56
951 11 17 400 79 57 18
820–845 – 50 – 11 51
TRI DI MONO SULF OLEF
+ + + + +
2H2 2H2 3H2 5H2 H2
DI MONO SAT MONO + H2S SAT
(reaction 1) (reaction 2) (reaction 3) (reaction 4) (reaction 5)
Fig. 1. LCO hydrotreatment scheme.
HDT hydrogenation reactions which are highly exothermic. Therefore, if the temperature of the reactor is not properly controlled a runaway event may occur. At industrial scale, LCO gasoils are commonly diluted in other gasoils such as SR before hydrotreatment. One of the reasons is the high exothermicity of LCO HDT. Only a fraction of the total production of LCO is therefore hydrotreated, the rest is commonly integrated to the fuel oil which is not a very valuable product. Since there is an increasing demand of on-road diesel, the incorporation of a more important fraction of LCO in the commercial diesel pool should be interesting, particularly in the European Union where diesel is imported. However, it is important to ensure the safe operation of the unit if higher fractions of LCO are integrated to the feeds. For this reason, the hydrotreatment of LCO feeds was chosen as case study to perform a thermal stability study. The next paragraphs describe the development of the dynamic model of the reactor used for this stability study. 2.2. Kinetics of the reactive system An apparent kinetic model of LCO hydrotreatment was developed to represent the source term of the material and thermal balances. The purpose of this kinetic model is to accurately describe the heat generated by the reactions. All experiments were carried out in pilot plants (fixed bed reactors) at industrial operating conditions with different real feeds. These experiments were used to fit the parameters of the kinetic model. The analytical characterizations of all the feeds and their hydrotreated effluents are as follows. The chemical composition of the gasoils was analyzed using standardized methods. The reaction scheme is illustrated in Fig. 1. It includes eight chemical families: six hydrocarbon families and two pure compounds (hydrogen and hydrogen sulfide produced from HDS reactions). The hypotheses considered are as follows. Triaromatics family (TRI) is hydrogenated via the external cycle to form diaromatics (DI). Diaromatics mainly constitute condensed doubled ring molecules. Reversibility of hydrogenation reactions is neglected for the moment since our operating conditions (high H2 partial pressures and T < 380 ◦ C) do not promote the reversibility but this should be taken into account in a further work. For HDS reactions, dibenzothiophenes were considered to be the only representative sulfur compounds in the cut since they are very refractory to HDS. Dibenzothiophenes are supposed to be desulfurized via the hydrogenation pathway to form monoaromatics (Whitehurst et al., 1998). This implies a higher level of hydrogen consumption. The hydrodenitrogenation reactions are not taken into account since the nitrogen content in feeds is very low ( < 0.05 wt%) and the heat generated by these reactions can be therefore neglected. The hydrocarbon families illustrated in Fig. 1 are present all along a boiling point interval of the gasoil feeds (150–380 ◦ C). In order to respect the material balance and to describe the different volatilities of the compounds along the reactor, three boiling point (BP) cuts were distinguished as illustrated in Table 2. Triaromatics are not present in the IBP-200 ◦ C cut since the boiling point of phenantherene/anthracene is about 340 ◦ C, therefore these compounds do not exist in lighter cuts. The same reasoning can be made for diaro-
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
a rj = kj0 · Cl j m
b · Cl j H2
Ej exp − R
1 1 − T Tref
(1)
where j is the reaction number, kj0 is the rate constant at a reference temperature Tref , Cim is the liquid concentration of the m hydrocarbon family, ClH is the liquid concentration of hydrogen and Ej is the 2 activation energy. The reactivity of each hydrocarbon lump was supposed to be independent of the BP cut. Therefore, according to the reaction scheme, the kinetic model constitutes five rate expressions that include five rate constants (kj0 ), 10 reaction orders (ai and bi ) and five activation energies Ej . The results of the optimized parameters are illustrated in Fig. 2. The total aromatics content measured in the effluents is compared to the predicted total aromatics content in Fig. 2a. The dotted lines indicate ± 5 wt%. The points labeled as “all feeds” correspond to the prediction of all the effluents obtained from the six different industrial feeds used for the parameters estimation. As illustrated, a very good agreement is obtained for the prediction of this property: 41 from 49 experimental points are predicted in Table 2 Description of the LCO pseudo-components. BP cut
IBP-200 ◦ C
200–300 ◦ C
300 ◦ C+
Chemical family
OLEFINS SATurates MONOaromatics – – SULFur compounds
OLEFINS SATurates MONOaromatics DIaromatics – SULFur compounds
OLEFINS SATurates MONOaromatics DIaromatics TRIaromatics SULFur compounds
Other compounds: H2 and H2 S.
Kinetic model prediction (wt%)
90
All feeds LCO A 5%
80 70 60 50 40 30 20 10 0
the ± 5 wt% interval; the other points are close to this precision. In the same figure, the points labeled “LCO A” are the effluents obtained from the LCO A (Table 1) at different HDT operating conditions. These points were distinguished since the thermal stability study was performed based on HDT experiments carried out with this LCO. The prediction of the aromatic sub-families illustrated in Fig. 2b, shows that a good prediction is also obtained for each sub-family. These results can be considered as validation of the kinetic model. 2.3. Experimental data of dynamic temperature profiles To carry out the thermal stability study of the gasoils hydrotreatment, an experimental test was carried out in a pilot plant. The objective of these experiments was to collect dynamic data for the development and fitting of the dynamic reactor model. The test was carried out in a reactor working in up-flow mode. The reactor scheme is illustrated in Fig. 3a. The internal temperatures in the reactor are measured and controlled using an internal multipoint thermowell (four thermowells TW1 to TW4 ). The first thermowell (TW1 ) is located in the inlet bed of glass beads and the other three thermowells (TW2 to TW4 ) in the catalytic bed. The internal reactor temperature is controlled by the temperature of the air layer located between the reactor and the heater (isolated device). The controller regulates the heat power of each electric heater to warm up the mentioned air layer. It should be pointed out that in case of excessive temperature, the control only diminishes or shuts-down the heater power, this means that the system is not designed to cool the reactor with another fluid. The external temperatures of the reactor are also measured at four different points (TC1 to TC4 ). These external temperatures will be later referred to as the cooling temperatures of the reactor. The dynamic hydrotreatment test was performed with two different gasoils. First, a SR gasoil (Table 1) was hydrotreated. The lower content of aromatics in the SR gasoil (32 wt%) compared to the LCO (79 wt%) implies a lower exothermicity due to less hydrogenation reactions. In this case, the isothermal temperature profile was easily controlled with the heater devices. The LCO hydrotreatment was carried out in the same test, just after the SR hydrotreatment. Before the shift of feed from SR to LCO, the internal temperatures (TW ) in the reactor were diminished to 300 ◦ C to prevent an excessive increment of temperature. The LCO was injected when the TW setpoints were reached. The experimental profiles (TW exp ) are illustrated in Fig. 3b, where time = 0 corresponds to the injection of the LCO feed.
80 Kinetic model prediction (wt%)
matics in the IBP-200 ◦ C cut. The same reaction scheme (Fig. 1) was applied to each BP cut to obtain 17 lumps (15 hydrocarbon pseudocomponents, H2 and H2 S). The kinetic parameters were estimated on the basis of 49 experimental pilot plant points with six different industrial feeds. The range of operating conditions is as follows: LSHV: 0.5–4 h−1 , temperature: 320–370 ◦ C, total pressure: 30–120 bar, H2 /HC ratio: 350–1300 NL/L. Temperature profile was maintained isothermal for the experiments intended to the kinetic fit. The gas–liquid phases were considered to be in equilibrium. The Grayson–Streed method was used to perform the flash calculations. The pilot plant was operated in up-flow mode. The catalyst was therefore considered to be completely wet and the reactions were supposed to be carried out only in the liquid phase. The general rate expression for all the reactions is
MONOaromatics DIaromatics TRIaromatics 5%
70 60 50 40 30 20 10 0
0
10 20 30 40 50 60 70 80 90 Measured aromatics (wt%)
315
0
10 20 30 40 50 60 70 80 Measured aromatics (wt %)
Fig. 2. Kinetic model results: (a) parity plot of total aromatics content and (b) parity plot of aromatics families.
316
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
outlet Legend:
heater 1
TC 4
TW 3
Experimental
350
Model
340 TC 3
TW 2 TC 2
TW 1 TC 1
Dint Dext
H2
Temperature (°C)
L2 L1
inlet
TW 3
heater 2
L3
TW : internal thermowell TC : external thermowell L : length Dint: internal reactor diameter Dext : external reactor diameter
heater 4
L4
catalyst bed
TW 4
heater 3
L5
360
glass beads
TW 4
330 320
TW 2
310 300 290
TW 1
280 0:00
0:30
Gasoil
1:00 1:30 time (hh:mm)
2:00
2:30
Fig. 3. (a) Experimental setup: pilot plant reactor scheme and (b) HDT reactor model. Comparison of experimental and calculated dynamic temperature trends.
The increasing temperature profile in the reactor was mainly caused by the heat released by the hydrogenation reactions. This transient behavior was recorded and used as experimental data to validate the dynamic reactor model.
Since gas and liquid phases are considered to be in equilibrium, the concentrations in both phases for each lump are linked as follows:
2.4. Dynamic model of the HDT reactor
The combination of Eq. (3) with the material balance (Eq. (2)), results in the following expression:
The model mainly includes the material and thermal balances for all the pseudo-components in all phases. The heating/cooling devices are also included. The following assumptions were taken into account to develop the model. As for the kinetic model, all the compounds and pseudo-components in the gas and liquid were considered to be in equilibrium. This is a conservative hypothesis for thermal runaway analysis since in the presence of mass transfer limitations, it is well known that the reactor is more stable. The liquid flow is described by an axial dispersion model while the gas flow is considered to be in plug flow. The gas holdup is assumed to be constant all along the reactor due to the low consumption of hydrogen compared to the inlet flowrate (hydrogen is added in large excess). The pressure drop is negligible compared to the total pressure, therefore, the reactor was considered to be isobaric. Liquid–solid external mass transfer limitations are also neglected. Internal diffusion limitations are taken into account with the apparent kinetic model. The upflow operation ensures a complete catalyst wetting. The superficial liquid velocity is considered to be constant along the reactor. Since heat transfers between all phases are very high, gas, liquid and solid were considered to be at the same temperature at a given reactor position. The transient material balance for each lump is given by
*(g · Cig + (l + s · por )Cil ) *C l *(ug · Cig ) = − ul · l i − g *t *z *z * Cil + ij · rj · s · s *z2 j
*Cil
g
*t
(3)
Cil
Hi + l + s · por = − · g R·T R
(2)
*
Hi T
*t
2
* Cil + Dlax · l *z2 Cil · Hi * ug R·T
*C l − ul · l i − g *z + ij · rj · s · s
*z (4)
j
The gas velocity changes all along the reactor due to the temperature profile, H2 consumption and H2 S production and the partial vaporization of the lumps. To be able to calculate this variable, the transient material balances (Eq. (2)) are summed for all species. At each time step and reactor axial coordinate, the constant pressure balance must be satisfied. Then, the following differential equation is obtained for the gas velocity.
*
ug T
*z
⎛ 2 * i Cil l * i Cil * i Cil R ⎝ −(l +s ·por ) = +Dax ·l −ul · l g ·Pt *t *z 2 *z ⎞ g · Pt 1 *T (5) + + ij · rj · s · s ⎠ R T 2 *t i
2
+ Dlax · l
H Pi g = Cil ⇒ Ci = Cil i Hi R·T
j
Concerning the thermal balance, the gas, the liquid and the solid phases were considered to be at the same temperature. According to this assumption, the following transient thermal balance is obtained,
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
where the heat generated by all the reactions is also taken into account. *T S g · Cpg · g + l · Cpl · l + s · Cps · s + w · Cpsteel · steel Sr *t * ¯ *T *T ax = l − (l · Cpl · l · ul + g · Cpg · g · ug ) *z *z *z + rj · s · s (−Hj ) − Uwall · Awall (T − TC ) (6) j
air · Cpair
In a vectorial formalism this equation becomes
*y = f (y ) *t
*TC *TC Sr = − air · Cpair · uair + Uwall · Awall (T − TC ) Svoid *t *z − Uloss · Awall
Sr Qk (TC − Tamb ) + Svoid Svoid · hheater,k
(7)
Finally, the control of the internal the fours heaters are implemented through Eq. (8), each controller having its own temperature setpoint in relation to each corresponding internal thermowell (TW1 –TW4 ).
*TTW,k * Qk = G(Tset,k − TTW,k ) − cont *t *t
(11)
Substituting Eq. (11) into Eq. (10), the reactor model becomes
*(y + x) = f (y + x) *t
(12)
If the perturbations are very small, the first order Taylor expansion can be applied to linearize the reactor model around the operating point. *(yi + xi ) *fi = fi (y1 , y2 , . . . , y3 , . . . , yn ) + x1 *t *y1 op *fi *fi *fi + x2 + · · · + xi + · · · + xn (13) *y2 op *yi op *yn op Subtracting Eq. (9) from Eq. (13), the perturbation model can be expressed as follows: *xi *fi *fi *fi *fi = x1 + x2 +· · ·+ xi + · · · + xn *t *y1 op *y2 op *yi op *yn op (14)
(8)
The results obtained with the dynamic model of the HDT reactor were compared to the experimental trends observed in the pilot plant test carried out with the shift of SR feed to LCO. The comparison is illustrated in Fig. 3b where the dotted lines represent the experimental temperatures and the uninterrupted lines the model predictions. As described before, after the SR hydrotreatment, the internal temperatures of the reactor (TW1 to TW4 ) were set to 300 ◦ C. When this setpoint was reached, the feed was shifted to the LCO. This moment corresponds to the beginning of the curves (time = 0) in Fig. 3b. After the LCO injection (about 15 min), an important increase of the temperature is observed. The simulation results accurately predict these dynamic temperature trends. Since the simulated temperature trends are in very good agreement with the experimental profiles, these results were used to validate the transient reactor model. 3. Dynamic thermal stability analysis The stability analysis is performed according to the perturbation of the operating variables (Perlmutter, 1972) and consists in the following main steps: (a) integration of the dynamic equations along the z axis using a numerical scheme, (b) perturbation of the dynamic reactor model, (c) linearization of the perturbation model and (d) resolution of the perturbation model and analysis of the solution. The condition of stability implies that all the perturbations must tend to zero when time tends to infinity. The next paragraphs describe these steps. The integration along the z axis is carried out using an upwind scheme for the convective terms and a centered scheme for the diffusion–dispersion terms. The reactor model can be written into a generalized formalism as follows:
*yi = fi (y1 , y2 , . . . , yi , . . . , yn ) *t
(10)
where y is the variables vector and f is the vector of the second member of partial differential equations. The perturbation vector x corresponds to the perturbation of the variable vector y around an operating point and is defined as follows: x = y p − y ⇒ y p = y + x
The accumulation in the gas, the liquid and the solid catalyst are taken into account as well as the inertia of the metallic mass of the reactor that plays an important role in the dynamic behavior. A heat transfer term is also considered to describe the heat loss through the reactor wall. The four electric heaters (Fig. 3a) were also integrated in the reactor model. As indicated before, the internal reactor temperature is controlled by the temperature of the thin air layer located between the reactor and the heaters. This air can be warmed up by an electrical resistance. The transient thermal balance of the air layer is described by Eq. (7). The heat losses through the isolated device are also taken into account.
317
(9)
Using a matrix formalism, the perturbation model becomes
*x = J · x *t
(15)
where J is the Jacobian of the application f . The resolution of this linear first order differential equation system allows to know if the perturbation of at least one of the whole set of variables diverge with time and propagates it. This means that the reactor has reached an unstable state that will lead to a runaway. The solution of the perturbation around the stationary point is given by Eq. (16) where Vi,j is the jth rank of the eigenvector i and i is the corresponding eigenvalue. xi = Vi,j · aj · ej ·t (16) j
Eigenvalues can be real (R) or complex (C). It clearly appears that xi → ∞ with time only when at least one of the n eigenvalues has a positive real part. Hence, the condition for stability is indicated in Eq. (17). ∀i,
lim xi = 0
t→∞
if i ∈ R ⇒ i < 0 and if i ∈ C ⇒ Re(i ) < 0
(17) (18)
To satisfy Eq. (17), the conditions indicated in Eq. (18) must be verified. This means that after the solution of the perturbation model, if the eigenvalues are only real (R), when all i are negative, perturbations are exponentially attenuated. If the eigenvalues are only real (R), but at least one is positive, runaway conditions are achieved and perturbations are exponentially amplified. However, for complex eigenvalues, oscillatory behavior of the perturbation is observed. This analysis was applied on the complete reactor model described before. The equations are discretized into 50 spatial steps
318
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
1.5
0.010 0.008
1.0
0.006 0.004 Im (λi)
Im (λi)
0.5 0.0 -0.5
0.002 0.000 -0.002 -0.004 -0.006
-1.0
-0.008 -1.5 -8.0
-6.0
-4.0 Re (λi)
-2.0
0.0
-0.010 -0.020 -0.015 -0.010 -0.005 0.000 Re(λi)
Fig. 4. Eigenvalue spectrum for U = 500 W/m2 /K: (a) all the eigenvalues and (b) zoom of values close to zero.
(nz). The dimension of the y , * y / *t and f vectors is 900 (each).
460
Temperature
In the y vector, the first 850 lines correspond to the liquid concentrations of the 17 lumps at each axial coordinate (nz = 50). The last 50 lines correspond to the temperature at each spatial step. An ana-
440
log reasoning can be made for the transient vector * y / *t. The first
420
Controller OFF
Temperature (°C)
850 lines of the f vector correspond to the second member of the material balance (MB) for each lump at each axial coordinate in the reactor. The last 50 lines contain the second member of the thermal balance (TB) at each spatial step. Therefore, the Jacobian (J) of the application f is a matrix of 900×900 dimension. For the thermal stability analysis, the effect of the variation of the heat transfer coefficient was studied. The desired operating reactor temperature was set to 380 ◦ C (isothermal profile in the catalytic zone). Three cases were studied: 500, 150 and 200 W/m2 /K. Hence, the eigenvalue spectrum of the Jacobian was determined for each case. Fig. 4 shows the eigenvalue spectrum obtained for the 500 W/m2 /K case. The abscissa corresponds to the real part of the eigenvalues while the ordinate is the imaginary part of the eigenvalues. All the eigenvalues are illustrated in Fig. 4a; a zoom near zero (Fig. 4b) shows that there are no eigenvalues with a positive real part. The stability condition of Eq. (18) is satisfied, therefore, the reactor must be stable. A simulation was carried out with the dynamic model in order to confirm the stable behavior of the reactor at 380 ◦ C and a heat transfer coefficient of 500 W/m2 /K. The simulation (Fig. 5) was started with the temperature controller activated; after 0.5 h of operation, stationary state conditions were reached and the controller was shut down. As illustrated, even without a control system, the reactor temperature remains at the same value (380 ◦ C). The stability analysis and the simulation lead to the same conclusions: for this heat transfer coefficient, the reactor is stable. The same analysis was applied for the case of a heat transfer coefficient of 150 W/m2 /K. The other parameters are the same compared to the first case. Fig. 6 shows, as for the previous case, the corresponding eigenvalue spectrum. The zoom of the eigenvalues close to zero (Fig. 6b) show that there are two complex eigenvalues of the spectrum which have a positive real part: 674 =8.51×10−4 +2.53×10−3 ·i and 675 = 8.51 × 10−4 − 2.53 × 10−3 · i. This means that the reactor should be unstable with an oscillatory behavior. The simulation of the case of a heat transfer coefficient of 150 W/m2 /K was done to confirm these results (Fig. 7a). As for the first case, the temperature control was activated to attain a temperature of 380 ◦ C and after 0.5 h the control was turned-off. As indicated by the simulation
Controller ON
Stop controller 400 380 360 340 320 300 0:00
1:00
2:00
3:00 time (h)
4:00
5:00
6:00
Fig. 5. Dynamic simulation of LCO HDT with a heat transfer coefficient of 500 W/m2 /K.
results, the reactor starts to oscillate about 3.5 h after the controller shut-down. The minimum and maximum temperature oscillations are close to 340 and 440 ◦ C, respectively. A Fourier transform was applied on the temperature signal in order to find the period of the oscillation. The frequency spectrum of the signal obtained by simulation is shown in Fig. 7b. One of the main characteristic frequencies is F = 3.54 × 10−4 Hz that corresponds to a period of T = 47 min. The solution of the perturbation (Eq. (16)) shows that when an eigenvalue is complex the exponential term can be expressed as indicated ei ·t = eRe(i )·t [cos(Im(i )t) + i sin(Im(i )t)]
(19)
The eRe(i )·t term clearly reveals that if the real part is positive, there is an exponential amplification of the perturbation. The imaginary part indicates that the perturbation leads to an oscillating behavior. Therefore, the period of the oscillation can also be determined
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
1.5
0.005 0.004
1.0
λ674
0.003 0.002 Im(λi)
0.5 Im (λi)
319
0.0
0.001 0.000 -0.001
-0.5
λ675
-0.002 -0.003
-1.0
-0.004 -1.5 -7.5
-5.5
-3.5 Re (λi)
-1.5
0.5
-0.005 -0.005 -0.004 -0.003 -0.002 -0.001 0.000 0.001 Re (λi)
Fig. 6. Eigenvalue spectrum for U = 150 W/m2 /K: (a) all the eigenvalues and (b) zoom of values close to zero.
Temperature (°C)
440
3.0E+06 Temperature Controller ON Controller OFF
420 400
Stop controller
380 360 340
Fourier coefficients
460
F = 3.54 × 10−4 Hz
2.0E+06 1.0E+06 0.0E+00 -1.0E+06
320 300 0:00 2:00 4:00 6:00 8:00 10:00 time (h)
-2.0E+06 0
0.0005 0.001 0.0015 0.002 Frequency (Hz)
Fig. 7. (a) Dynamic simulation of LCO HDT with a heat transfer coefficient of 150 W/m2 /K and (b) Fourier transform frequency spectrum for a heat transfer coefficient of 150 W/m2 /K.
using the imaginary part of this eigenvalue. T=
2· 2· = = 42 min Im(i ) 2.53 × 10−3
This period should correspond to one characteristic frequency obtained from the Fourier transform of the temperature signal. This value (42 min) is very close to the value obtained from the frequency spectrum (47 min). The small difference is due to the accuracy of the linearization when oscillations are far from the stationary point. To confirm this assertion, another simulation was performed with a heat transfer coefficient of 200 W/m2 /K. This value corresponds to the stability limit. In this case two eigenvalues of the whole spectrum were complex with a positive real part: 9.62 × 10−5 + 2.64 × 10−3 · i and 9.62 × 10−5 − 2.64 × 10−3 · i. The period of oscillation calculated with the eigenvalues is 39.68 min. The Fourier transform calculations result in a frequency value of F = 4.2 × 10−4 Hz that corresponds to a period of 39.68 min. The period of oscillation calculated with both methods is in excellent agreement. For this case, the real part of the eigenvalues is one order of magnitude lower than the one obtained from the previous case with a heat transfer coefficient of 150 W/m2 /K. In other words, for the U = 200 W/m2 /K case, the temperature oscillations are composed in values of 10−3 ◦ C while for the
U = 150 W/m2 /K case, the temperature oscillations were of about 110 ◦ C. Then, the time constant of the amplification is = 1/Re(i ) = 173.3 min compared to 19.6 min for the U = 150 W/m2 /K case. As shown in this work, the analysis of the eigenvalue spectrum allows to determine if the system is stable or unstable. Moreover, the interpretation of the eigenvalues provides information about the propagation of the perturbation which may be of different types: 1—exponential attenuation, 2—oscillations with exponential attenuation, 3—oscillations with exponential amplification, and 4—exponential amplification. If the system is oscillating the frequency of the oscillation can be calculated with the imaginary part of the eigenvalues. 4. Conclusions This work presents a thermal stability study of a refining process (hydrotreatment of gasoils). This process consists on the hydrogenation of unsaturated compounds which can be present in high amounts in some gasoil feeds like LCO ( ≈ 90 wt%). The hydrogenation reactions are highly exothermic. The hydrotreatment of this great content of aromatics and olefins produces a high heat release that may lead to a runaway of the reactor. Therefore, a
320
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
thermal stability study of this process was carried out. The dynamic stability study is performed on the basis of a mathematical model that describes the reactive system. For this reason, a dynamic model of the hydrotreatment of gasoils was developed. The kinetics was fitted using experimental points obtained with industrial catalyst and industrial feeds in the classical range of operating conditions of this process. A dynamic model of HDT reactor was developed to simulate an existing pilot plant. The model takes into account a detailed description of the reactive system and the configuration of the pilot plant: three phases (gas–liquid–solid), the kinetics, the flow characterization in the gas and liquid phases, the heater devices and the control system. An experimental test was carried out in the pilot plant to collect dynamic data for the validation of the model. The dynamic temperature profiles of the reactor are in good agreement with the simulation results and the model was validated with this comparison. The thermal stability analysis of the LCO hydrotreatment was performed with a perturbation method. The criterion of stability imposes that after the solution of the perturbation system, all the resultant eigenvalues must have a negative real part, otherwise the reactor will be unstable. The influence of the heat transfer coefficient on the reactor stability was evaluated. Different cases were studied in this work, a stable case and two unstable cases. The conclusions obtained via the thermal stability study were compared to the conclusions obtained by model simulations. For all cases, the conclusions were always in good agreement. For the unstable cases, an oscillatory behavior was observed by simulation and the stability analysis successfully predicts the occurrence of the oscillations. Furthermore, the frequency and the period of these oscillations were compared to those obtained by Fourier transform signal treatment of the simulation results. In the case of small oscillations where the 1st order Taylor expansion is applicable, the results obtained by both methods were the same. The stability analysis presented in this work is therefore a reliable method to establish if the reactor is stable or not. In further work, the impact of the variation of other parameters such as the feed composition and flowrates of reactants will be achieved. The localization of the runaway ignition will be carried out by analyzing the eigenvectors. In addition, the impact of Jacobian accuracy will also be studied, the mesh and the numerical scheme have to be also evaluated. Finally, this stability analysis will be applied in transient conditions in order to determine the dynamic unstable regions. Notation a A b C Cp Dax E f
f
F G h H
H J k0 nz P
reaction order of hydrocarbon pseudo-component in the liquid phase heat exchange surface, m2 /m3 reaction order of hydrogen in the liquid phase Concentration, mol/m3 specific heat capacity, J/kg/K axial dispersion coefficient, m2 /s activation energy, J/mol second member of partial differential equations vector of the second member of partial differential equations frequency, Hz controller gain heater height, m Henry coefficient heat of reaction, J/mol Jacobian matrix reaction rate constant at the reference temperature number of steps for the reactor discretization pressure, bar
t T T Tamb Tref u U
total pressure, bar electrical power, W reaction rate, mol/s/kgcatalyst ideal gas constant, J/mol K reactor section surface, m2 section surface between the heaters and the reactor wall, m2 time, s temperature, K period, min ambient temperature, K reference temperature, K superficial liquid velocity, m/s heat transfer coefficient, W/m2 /K
V x
eigenvector perturbation
x
vector of perturbations
y z
variables vector axial coordinate, m
Pt Q r R Sr Svoid
Greek letters
¯ ax cont
holdup eigenvalue effective thermal conductivity, W/m/K stoichiometric coefficient density, kg/m3 time constant of the amplification, min controller time derivative constant
Subscripts and superscripts air C g heater i j k l loss n op p por s set steel wall W
air between the reactor wall and the heater, K cooling gas heater pseudo-component index or eigenvalue index reaction index number of heater index liquid relative to heat loss total number of variables operating point perturbed porosity solid setpoint reactor steel reactor external wall internal reactor thermowell
Abbreviations C DI IBP Im MB MONO R Re
complex diaromatics initial boiling point imaginary part material balance monoaromatics real real part
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
SAT SULF TB TRI
saturates sulfur compounds thermal balance triaromatics
Acknowledgments Authors are grateful to Mr. D. Hudebine for supplying a part of the theoretical data of LCO hydrotreatment. This work was partially supported by the “Axelera” French Research Project. References Agrawal, R., West, D.H., Balakotaiah, V., 2007. Chemical Engineering Science 62, 4926–4943.
321
Alós, M.A., Strozzi, F., Zaldívar, J.M., 1996. Chemical Engineering Science 51, 3089–3094. Balakotaiah, V., Luss, D., 1991. A.I.Ch.E. Journal 37 (12), 1780–1787. Balakotaiah, V., Kodra, D., Nguyen, D., 1995. Chemical Engineering Science 50 (7), 1149–1171. Bilous, O., Amundson, N.R., 1956. A.I.Ch.E. Journal 2 (1), 117–126. Luss, D., Sheintuch, M., 2005. Catalysis Today 105, 254–274. Morbidelli, M., Varma, A., 1986. A.I.Ch.E. Journal 32 (2), 297–306. Morbidelli, M., Varma, A., 1987. A.I.Ch.E. Journal 33 (12), 1949–1958. Perlmutter, D.D., 1972. Stability of Chemical Reactors. Prentice-Hall, New Jersey. van Heerden, C., 1953. Industrial and Engineering Chemistry 45, 1242–1247. Viswanathan, G.A., Luss, D., 2006. Industrial and Engineering Chemistry Research 45, 7057–7066. Whitehurst, D.D., Isoda, T., Mochida, I., 1998. Advances in Catalysis 42, 345–471. Wu, H., Rota, R., Morbidelli, M., Varma, A., 1999. Chemical Engineering Science 54 (20), 4579–4588. Zaldívar, J.M., Cano, J., Alós, M.A., 2003. Journal of Loss Prevention in the Process Industries 16, 187–200.
Contents lists available at ScienceDirect
Chemical Engineering Science journal homepage: w w w . e l s e v i e r . c o m / l o c a t e / c e s
Thermal runaway analysis of a three-phase reactor for LCO hydrotreatment Jean-Marc Schweitzer ∗ , Clementina López-García, Daniel Ferré Institut Français du Pétrole, Process Development and Engineering Division, PO Box 3, 69360 Solaize, France
A R T I C L E
I N F O
Article history: Received 30 June 2008 Received in revised form 1 July 2009 Accepted 14 July 2009 Available online 21 July 2009 Keywords: Thermal stability Gasoils hydrotreatment Packed bed Dynamic simulation Light cycle oil
A B S T R A C T
Safety is a high-priority topic for the chemical industry to minimize the frequency and severity of accidents while keeping the productivity and quality of the production. The processes that may undergo thermal runaways due to exothermic reactions are at the heart of the risks of accidents. The study of such highly reactive systems is essential to achieve a safe and productive operation of existing processes and to ensure inherently safe new designs. It is well known that the stationary analysis with the van Heerden criterion must be satisfied, however, this is not sufficient to ensure reactor stability. Only dynamic analysis can provide an accurate answer concerning the safe operation of the reactor. Most of the reported stability studies are carried out for relatively simple systems (pseudo-homogeneous models, simple reaction schemes). This work presents the dynamic thermal stability study of a refining process (hydrotreatment of light cycle oils) carried out with real gasoils at industrial operating conditions. A 1D dynamic model that accurately represents the reactive system (gas–liquid–solid) was developed and validated with experimental pilot plant data. This mathematical model was used to perform the thermal stability analysis of the dynamic system using a perturbation method. The effect of the variation of the heat transfer coefficient on the thermal stability is presented. The spectral analysis of the eigenvalues indicating the stable/unstable behavior of the reactive system was compared with dynamic simulations. An excellent agreement was found between the simulations and the stability analysis. The case of oscillating behavior is also described. The frequencies of oscillation determined by the stability analysis are compared to the frequencies calculated by Fourier transform applied to the simulated signal. The reactor behavior and the oscillations features are accurately predicted with this stability analysis method. © 2009 Elsevier Ltd. All rights reserved.
1. Introduction One of the security priorities involving highly reactive systems is the risk of thermal runaway. The temperature increase for the reactions that follow an Arrhenius law, induces the rise of the heat generation that further increases the reaction temperature and so on; this situation may result in a thermal runaway. The consequences of a runaway can be: the early deactivation of the catalyst (coking), the loss of selectivity, conversion or the operability of the unit, the onset of secondary reactions and in some cases even the reactor explosion. The design of industrial reactors must rely on an accurate thermal stability study. Indeed, the design of the whole reactive system (reactor and the devices used for heat input/removal) should ensure a priori the thermally stable operation of the reactor. It is also important to determine the regions of operating conditions where the reactor has an unreliable behavior (runaway regions). In practice,
∗ Corresponding author. Tel.: +33 4 78 02 25 98; fax: +33 4 78 02 20 08. E-mail address: [email protected] (J.-M. Schweitzer). 0009-2509/$ - see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.ces.2009.07.012
most of the industrial designs of reactive systems are based only on the van Heerden (1953) stationary criterion. This criterion imposes that the slope of heat generated by the reactions (dQgen /dT) must be lower than the slope of the heat extracted from the system (dQexc /dT). Since this comparison is carried out only under stationary conditions the dynamic behavior of the system is not taken into account and the thermal stability cannot be guaranteed. Indeed, only the dynamic analysis will provide an accurate answer concerning the safe operation of the reactor. An extensive number of literature works have tackled the reactors thermal stability topic. The objective is to dispose of theoretical methods mainly to determine the boundaries of runaway/nonrunaway. There are two foremost approaches for the thermal stability assessment: parametric sensitivity and application of dynamic systems theory. Parametric sensitivity consists of the determination of the operation regions where a small change of the inlet parameters produces high changes of the variables of the system such as temperature. This concept was introduced in the field of chemical reactors by Bilous and Amundson (1956). Since then, many authors have extensively applied the parametric sensitivity theory to study different reactive systems. Among the studies
314
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
concerning packed bed catalytic reactors we can cite the work of Morbidelli and Varma (1986, 1987), Balakotaiah and Luss (1991), Balakotaiah et al. (1995). These authors present the study of the role of inter and intraparticle resistances in the criteria to determine the runaway boundaries as well as the role of the interphase gradients. Wu et al. (1999) used the parametric sensitivity to determine the critical ignition conditions to achieve the steady-state auto-thermal operation of reactors with reverse flow. Concerning the dynamic stability analysis, the development of an intrinsic criterion based on the sum of Lyapunov exponents is presented in several literature works (Alós et al., 1996; Zaldívar et al., 2003) where the instabilities are characterized with chaotic attractors. These studies have been performed in batch and semi-batch reactors. The perturbation stability analysis ant the chaos theory has more recently been extended to deal with the determination (1D or 2D) hot spots in packed-bed reactors (Viswanathan and Luss, 2006; Luss and Sheintuch, 2005; Agrawal et al., 2007). Most of the models study the instabilities in homogeneous or pseudo-homogeneous systems which are not the case of this work. The objective of this study is to achieve the dynamic thermal stability analysis of a complex refining process. The hydrotreatment (HDT) of an industrial light cycle oil (LCO) gasoil in a fixed bed reactor is used as a case study. Experiments were carried out at industrial operating conditions in a representative pilot plant. A complete dynamic model of the three-phase system was developed, fitted and validated on experimental data. The dynamic simulation model was used to assess the thermal stability analysis. This work presents this reliable stability method. The effect of the variation of the heat transfer coefficient is presented. In the case of unstable behavior the significance of the eigenvalues is also illustrated (oscillations). In the case of oscillations occurrence, their frequency is determined. 2. LCO hydrotreatment process 2.1. Introduction Commercial diesels are constituted of a blend of gasoils coming from different refining processes or operations. These different gasoils can be among others: straight run (SR) issued from the atmospheric distillation, coker gasoil produced in coking units, light cycle oil (LCO) issued from fluid catalytic cracking (FCC), hydrocracking (HDC). These gasoils must follow a hydrotreatment (HDT) process to meet commercial specifications. The main hydrotreatment reactions are: hydrogenation of aromatics and olefins, hydrodesulfurisation of organic sulfur compounds and hydrodenitrogenation. The chemical composition of the feedstocks depends directly on the provenance of the crude oil as well as the process from which it is produced and its severity (SR, LCO, coker . . . ). Table 1 presents an example of some characteristics of a SR gasoil and a LCO. The current European onroad specification of gasoils (NF EN 590) is also indicated. As illustrated, the aromatics content is completely different between both gasoils. The LCO are generally characterized by very high aromatics contents that result in a very low cetane number and high density. The aromatics and olefins content in LCO gasoils is reduced by the Table 1 Comparison of some physico-chemical characteristics of a SR and a LCO gasoil with the European on-sroad gasoil specification (NF EN 590). Description 3
Density 15/4 (kg/m ) Olefins (wt%) Sulfur (wt ppm) Total aromatics (wt%) Diaromatics+ (wt%) Cetane index
SR gasoil
LCO A gasoil
NF EN 590 specification
857 – 11 400 32 19 56
951 11 17 400 79 57 18
820–845 – 50 – 11 51
TRI DI MONO SULF OLEF
+ + + + +
2H2 2H2 3H2 5H2 H2
DI MONO SAT MONO + H2S SAT
(reaction 1) (reaction 2) (reaction 3) (reaction 4) (reaction 5)
Fig. 1. LCO hydrotreatment scheme.
HDT hydrogenation reactions which are highly exothermic. Therefore, if the temperature of the reactor is not properly controlled a runaway event may occur. At industrial scale, LCO gasoils are commonly diluted in other gasoils such as SR before hydrotreatment. One of the reasons is the high exothermicity of LCO HDT. Only a fraction of the total production of LCO is therefore hydrotreated, the rest is commonly integrated to the fuel oil which is not a very valuable product. Since there is an increasing demand of on-road diesel, the incorporation of a more important fraction of LCO in the commercial diesel pool should be interesting, particularly in the European Union where diesel is imported. However, it is important to ensure the safe operation of the unit if higher fractions of LCO are integrated to the feeds. For this reason, the hydrotreatment of LCO feeds was chosen as case study to perform a thermal stability study. The next paragraphs describe the development of the dynamic model of the reactor used for this stability study. 2.2. Kinetics of the reactive system An apparent kinetic model of LCO hydrotreatment was developed to represent the source term of the material and thermal balances. The purpose of this kinetic model is to accurately describe the heat generated by the reactions. All experiments were carried out in pilot plants (fixed bed reactors) at industrial operating conditions with different real feeds. These experiments were used to fit the parameters of the kinetic model. The analytical characterizations of all the feeds and their hydrotreated effluents are as follows. The chemical composition of the gasoils was analyzed using standardized methods. The reaction scheme is illustrated in Fig. 1. It includes eight chemical families: six hydrocarbon families and two pure compounds (hydrogen and hydrogen sulfide produced from HDS reactions). The hypotheses considered are as follows. Triaromatics family (TRI) is hydrogenated via the external cycle to form diaromatics (DI). Diaromatics mainly constitute condensed doubled ring molecules. Reversibility of hydrogenation reactions is neglected for the moment since our operating conditions (high H2 partial pressures and T < 380 ◦ C) do not promote the reversibility but this should be taken into account in a further work. For HDS reactions, dibenzothiophenes were considered to be the only representative sulfur compounds in the cut since they are very refractory to HDS. Dibenzothiophenes are supposed to be desulfurized via the hydrogenation pathway to form monoaromatics (Whitehurst et al., 1998). This implies a higher level of hydrogen consumption. The hydrodenitrogenation reactions are not taken into account since the nitrogen content in feeds is very low ( < 0.05 wt%) and the heat generated by these reactions can be therefore neglected. The hydrocarbon families illustrated in Fig. 1 are present all along a boiling point interval of the gasoil feeds (150–380 ◦ C). In order to respect the material balance and to describe the different volatilities of the compounds along the reactor, three boiling point (BP) cuts were distinguished as illustrated in Table 2. Triaromatics are not present in the IBP-200 ◦ C cut since the boiling point of phenantherene/anthracene is about 340 ◦ C, therefore these compounds do not exist in lighter cuts. The same reasoning can be made for diaro-
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
a rj = kj0 · Cl j m
b · Cl j H2
Ej exp − R
1 1 − T Tref
(1)
where j is the reaction number, kj0 is the rate constant at a reference temperature Tref , Cim is the liquid concentration of the m hydrocarbon family, ClH is the liquid concentration of hydrogen and Ej is the 2 activation energy. The reactivity of each hydrocarbon lump was supposed to be independent of the BP cut. Therefore, according to the reaction scheme, the kinetic model constitutes five rate expressions that include five rate constants (kj0 ), 10 reaction orders (ai and bi ) and five activation energies Ej . The results of the optimized parameters are illustrated in Fig. 2. The total aromatics content measured in the effluents is compared to the predicted total aromatics content in Fig. 2a. The dotted lines indicate ± 5 wt%. The points labeled as “all feeds” correspond to the prediction of all the effluents obtained from the six different industrial feeds used for the parameters estimation. As illustrated, a very good agreement is obtained for the prediction of this property: 41 from 49 experimental points are predicted in Table 2 Description of the LCO pseudo-components. BP cut
IBP-200 ◦ C
200–300 ◦ C
300 ◦ C+
Chemical family
OLEFINS SATurates MONOaromatics – – SULFur compounds
OLEFINS SATurates MONOaromatics DIaromatics – SULFur compounds
OLEFINS SATurates MONOaromatics DIaromatics TRIaromatics SULFur compounds
Other compounds: H2 and H2 S.
Kinetic model prediction (wt%)
90
All feeds LCO A 5%
80 70 60 50 40 30 20 10 0
the ± 5 wt% interval; the other points are close to this precision. In the same figure, the points labeled “LCO A” are the effluents obtained from the LCO A (Table 1) at different HDT operating conditions. These points were distinguished since the thermal stability study was performed based on HDT experiments carried out with this LCO. The prediction of the aromatic sub-families illustrated in Fig. 2b, shows that a good prediction is also obtained for each sub-family. These results can be considered as validation of the kinetic model. 2.3. Experimental data of dynamic temperature profiles To carry out the thermal stability study of the gasoils hydrotreatment, an experimental test was carried out in a pilot plant. The objective of these experiments was to collect dynamic data for the development and fitting of the dynamic reactor model. The test was carried out in a reactor working in up-flow mode. The reactor scheme is illustrated in Fig. 3a. The internal temperatures in the reactor are measured and controlled using an internal multipoint thermowell (four thermowells TW1 to TW4 ). The first thermowell (TW1 ) is located in the inlet bed of glass beads and the other three thermowells (TW2 to TW4 ) in the catalytic bed. The internal reactor temperature is controlled by the temperature of the air layer located between the reactor and the heater (isolated device). The controller regulates the heat power of each electric heater to warm up the mentioned air layer. It should be pointed out that in case of excessive temperature, the control only diminishes or shuts-down the heater power, this means that the system is not designed to cool the reactor with another fluid. The external temperatures of the reactor are also measured at four different points (TC1 to TC4 ). These external temperatures will be later referred to as the cooling temperatures of the reactor. The dynamic hydrotreatment test was performed with two different gasoils. First, a SR gasoil (Table 1) was hydrotreated. The lower content of aromatics in the SR gasoil (32 wt%) compared to the LCO (79 wt%) implies a lower exothermicity due to less hydrogenation reactions. In this case, the isothermal temperature profile was easily controlled with the heater devices. The LCO hydrotreatment was carried out in the same test, just after the SR hydrotreatment. Before the shift of feed from SR to LCO, the internal temperatures (TW ) in the reactor were diminished to 300 ◦ C to prevent an excessive increment of temperature. The LCO was injected when the TW setpoints were reached. The experimental profiles (TW exp ) are illustrated in Fig. 3b, where time = 0 corresponds to the injection of the LCO feed.
80 Kinetic model prediction (wt%)
matics in the IBP-200 ◦ C cut. The same reaction scheme (Fig. 1) was applied to each BP cut to obtain 17 lumps (15 hydrocarbon pseudocomponents, H2 and H2 S). The kinetic parameters were estimated on the basis of 49 experimental pilot plant points with six different industrial feeds. The range of operating conditions is as follows: LSHV: 0.5–4 h−1 , temperature: 320–370 ◦ C, total pressure: 30–120 bar, H2 /HC ratio: 350–1300 NL/L. Temperature profile was maintained isothermal for the experiments intended to the kinetic fit. The gas–liquid phases were considered to be in equilibrium. The Grayson–Streed method was used to perform the flash calculations. The pilot plant was operated in up-flow mode. The catalyst was therefore considered to be completely wet and the reactions were supposed to be carried out only in the liquid phase. The general rate expression for all the reactions is
MONOaromatics DIaromatics TRIaromatics 5%
70 60 50 40 30 20 10 0
0
10 20 30 40 50 60 70 80 90 Measured aromatics (wt%)
315
0
10 20 30 40 50 60 70 80 Measured aromatics (wt %)
Fig. 2. Kinetic model results: (a) parity plot of total aromatics content and (b) parity plot of aromatics families.
316
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
outlet Legend:
heater 1
TC 4
TW 3
Experimental
350
Model
340 TC 3
TW 2 TC 2
TW 1 TC 1
Dint Dext
H2
Temperature (°C)
L2 L1
inlet
TW 3
heater 2
L3
TW : internal thermowell TC : external thermowell L : length Dint: internal reactor diameter Dext : external reactor diameter
heater 4
L4
catalyst bed
TW 4
heater 3
L5
360
glass beads
TW 4
330 320
TW 2
310 300 290
TW 1
280 0:00
0:30
Gasoil
1:00 1:30 time (hh:mm)
2:00
2:30
Fig. 3. (a) Experimental setup: pilot plant reactor scheme and (b) HDT reactor model. Comparison of experimental and calculated dynamic temperature trends.
The increasing temperature profile in the reactor was mainly caused by the heat released by the hydrogenation reactions. This transient behavior was recorded and used as experimental data to validate the dynamic reactor model.
Since gas and liquid phases are considered to be in equilibrium, the concentrations in both phases for each lump are linked as follows:
2.4. Dynamic model of the HDT reactor
The combination of Eq. (3) with the material balance (Eq. (2)), results in the following expression:
The model mainly includes the material and thermal balances for all the pseudo-components in all phases. The heating/cooling devices are also included. The following assumptions were taken into account to develop the model. As for the kinetic model, all the compounds and pseudo-components in the gas and liquid were considered to be in equilibrium. This is a conservative hypothesis for thermal runaway analysis since in the presence of mass transfer limitations, it is well known that the reactor is more stable. The liquid flow is described by an axial dispersion model while the gas flow is considered to be in plug flow. The gas holdup is assumed to be constant all along the reactor due to the low consumption of hydrogen compared to the inlet flowrate (hydrogen is added in large excess). The pressure drop is negligible compared to the total pressure, therefore, the reactor was considered to be isobaric. Liquid–solid external mass transfer limitations are also neglected. Internal diffusion limitations are taken into account with the apparent kinetic model. The upflow operation ensures a complete catalyst wetting. The superficial liquid velocity is considered to be constant along the reactor. Since heat transfers between all phases are very high, gas, liquid and solid were considered to be at the same temperature at a given reactor position. The transient material balance for each lump is given by
*(g · Cig + (l + s · por )Cil ) *C l *(ug · Cig ) = − ul · l i − g *t *z *z * Cil + ij · rj · s · s *z2 j
*Cil
g
*t
(3)
Cil
Hi + l + s · por = − · g R·T R
(2)
*
Hi T
*t
2
* Cil + Dlax · l *z2 Cil · Hi * ug R·T
*C l − ul · l i − g *z + ij · rj · s · s
*z (4)
j
The gas velocity changes all along the reactor due to the temperature profile, H2 consumption and H2 S production and the partial vaporization of the lumps. To be able to calculate this variable, the transient material balances (Eq. (2)) are summed for all species. At each time step and reactor axial coordinate, the constant pressure balance must be satisfied. Then, the following differential equation is obtained for the gas velocity.
*
ug T
*z
⎛ 2 * i Cil l * i Cil * i Cil R ⎝ −(l +s ·por ) = +Dax ·l −ul · l g ·Pt *t *z 2 *z ⎞ g · Pt 1 *T (5) + + ij · rj · s · s ⎠ R T 2 *t i
2
+ Dlax · l
H Pi g = Cil ⇒ Ci = Cil i Hi R·T
j
Concerning the thermal balance, the gas, the liquid and the solid phases were considered to be at the same temperature. According to this assumption, the following transient thermal balance is obtained,
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
where the heat generated by all the reactions is also taken into account. *T S g · Cpg · g + l · Cpl · l + s · Cps · s + w · Cpsteel · steel Sr *t * ¯ *T *T ax = l − (l · Cpl · l · ul + g · Cpg · g · ug ) *z *z *z + rj · s · s (−Hj ) − Uwall · Awall (T − TC ) (6) j
air · Cpair
In a vectorial formalism this equation becomes
*y = f (y ) *t
*TC *TC Sr = − air · Cpair · uair + Uwall · Awall (T − TC ) Svoid *t *z − Uloss · Awall
Sr Qk (TC − Tamb ) + Svoid Svoid · hheater,k
(7)
Finally, the control of the internal the fours heaters are implemented through Eq. (8), each controller having its own temperature setpoint in relation to each corresponding internal thermowell (TW1 –TW4 ).
*TTW,k * Qk = G(Tset,k − TTW,k ) − cont *t *t
(11)
Substituting Eq. (11) into Eq. (10), the reactor model becomes
*(y + x) = f (y + x) *t
(12)
If the perturbations are very small, the first order Taylor expansion can be applied to linearize the reactor model around the operating point. *(yi + xi ) *fi = fi (y1 , y2 , . . . , y3 , . . . , yn ) + x1 *t *y1 op *fi *fi *fi + x2 + · · · + xi + · · · + xn (13) *y2 op *yi op *yn op Subtracting Eq. (9) from Eq. (13), the perturbation model can be expressed as follows: *xi *fi *fi *fi *fi = x1 + x2 +· · ·+ xi + · · · + xn *t *y1 op *y2 op *yi op *yn op (14)
(8)
The results obtained with the dynamic model of the HDT reactor were compared to the experimental trends observed in the pilot plant test carried out with the shift of SR feed to LCO. The comparison is illustrated in Fig. 3b where the dotted lines represent the experimental temperatures and the uninterrupted lines the model predictions. As described before, after the SR hydrotreatment, the internal temperatures of the reactor (TW1 to TW4 ) were set to 300 ◦ C. When this setpoint was reached, the feed was shifted to the LCO. This moment corresponds to the beginning of the curves (time = 0) in Fig. 3b. After the LCO injection (about 15 min), an important increase of the temperature is observed. The simulation results accurately predict these dynamic temperature trends. Since the simulated temperature trends are in very good agreement with the experimental profiles, these results were used to validate the transient reactor model. 3. Dynamic thermal stability analysis The stability analysis is performed according to the perturbation of the operating variables (Perlmutter, 1972) and consists in the following main steps: (a) integration of the dynamic equations along the z axis using a numerical scheme, (b) perturbation of the dynamic reactor model, (c) linearization of the perturbation model and (d) resolution of the perturbation model and analysis of the solution. The condition of stability implies that all the perturbations must tend to zero when time tends to infinity. The next paragraphs describe these steps. The integration along the z axis is carried out using an upwind scheme for the convective terms and a centered scheme for the diffusion–dispersion terms. The reactor model can be written into a generalized formalism as follows:
*yi = fi (y1 , y2 , . . . , yi , . . . , yn ) *t
(10)
where y is the variables vector and f is the vector of the second member of partial differential equations. The perturbation vector x corresponds to the perturbation of the variable vector y around an operating point and is defined as follows: x = y p − y ⇒ y p = y + x
The accumulation in the gas, the liquid and the solid catalyst are taken into account as well as the inertia of the metallic mass of the reactor that plays an important role in the dynamic behavior. A heat transfer term is also considered to describe the heat loss through the reactor wall. The four electric heaters (Fig. 3a) were also integrated in the reactor model. As indicated before, the internal reactor temperature is controlled by the temperature of the thin air layer located between the reactor and the heaters. This air can be warmed up by an electrical resistance. The transient thermal balance of the air layer is described by Eq. (7). The heat losses through the isolated device are also taken into account.
317
(9)
Using a matrix formalism, the perturbation model becomes
*x = J · x *t
(15)
where J is the Jacobian of the application f . The resolution of this linear first order differential equation system allows to know if the perturbation of at least one of the whole set of variables diverge with time and propagates it. This means that the reactor has reached an unstable state that will lead to a runaway. The solution of the perturbation around the stationary point is given by Eq. (16) where Vi,j is the jth rank of the eigenvector i and i is the corresponding eigenvalue. xi = Vi,j · aj · ej ·t (16) j
Eigenvalues can be real (R) or complex (C). It clearly appears that xi → ∞ with time only when at least one of the n eigenvalues has a positive real part. Hence, the condition for stability is indicated in Eq. (17). ∀i,
lim xi = 0
t→∞
if i ∈ R ⇒ i < 0 and if i ∈ C ⇒ Re(i ) < 0
(17) (18)
To satisfy Eq. (17), the conditions indicated in Eq. (18) must be verified. This means that after the solution of the perturbation model, if the eigenvalues are only real (R), when all i are negative, perturbations are exponentially attenuated. If the eigenvalues are only real (R), but at least one is positive, runaway conditions are achieved and perturbations are exponentially amplified. However, for complex eigenvalues, oscillatory behavior of the perturbation is observed. This analysis was applied on the complete reactor model described before. The equations are discretized into 50 spatial steps
318
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
1.5
0.010 0.008
1.0
0.006 0.004 Im (λi)
Im (λi)
0.5 0.0 -0.5
0.002 0.000 -0.002 -0.004 -0.006
-1.0
-0.008 -1.5 -8.0
-6.0
-4.0 Re (λi)
-2.0
0.0
-0.010 -0.020 -0.015 -0.010 -0.005 0.000 Re(λi)
Fig. 4. Eigenvalue spectrum for U = 500 W/m2 /K: (a) all the eigenvalues and (b) zoom of values close to zero.
(nz). The dimension of the y , * y / *t and f vectors is 900 (each).
460
Temperature
In the y vector, the first 850 lines correspond to the liquid concentrations of the 17 lumps at each axial coordinate (nz = 50). The last 50 lines correspond to the temperature at each spatial step. An ana-
440
log reasoning can be made for the transient vector * y / *t. The first
420
Controller OFF
Temperature (°C)
850 lines of the f vector correspond to the second member of the material balance (MB) for each lump at each axial coordinate in the reactor. The last 50 lines contain the second member of the thermal balance (TB) at each spatial step. Therefore, the Jacobian (J) of the application f is a matrix of 900×900 dimension. For the thermal stability analysis, the effect of the variation of the heat transfer coefficient was studied. The desired operating reactor temperature was set to 380 ◦ C (isothermal profile in the catalytic zone). Three cases were studied: 500, 150 and 200 W/m2 /K. Hence, the eigenvalue spectrum of the Jacobian was determined for each case. Fig. 4 shows the eigenvalue spectrum obtained for the 500 W/m2 /K case. The abscissa corresponds to the real part of the eigenvalues while the ordinate is the imaginary part of the eigenvalues. All the eigenvalues are illustrated in Fig. 4a; a zoom near zero (Fig. 4b) shows that there are no eigenvalues with a positive real part. The stability condition of Eq. (18) is satisfied, therefore, the reactor must be stable. A simulation was carried out with the dynamic model in order to confirm the stable behavior of the reactor at 380 ◦ C and a heat transfer coefficient of 500 W/m2 /K. The simulation (Fig. 5) was started with the temperature controller activated; after 0.5 h of operation, stationary state conditions were reached and the controller was shut down. As illustrated, even without a control system, the reactor temperature remains at the same value (380 ◦ C). The stability analysis and the simulation lead to the same conclusions: for this heat transfer coefficient, the reactor is stable. The same analysis was applied for the case of a heat transfer coefficient of 150 W/m2 /K. The other parameters are the same compared to the first case. Fig. 6 shows, as for the previous case, the corresponding eigenvalue spectrum. The zoom of the eigenvalues close to zero (Fig. 6b) show that there are two complex eigenvalues of the spectrum which have a positive real part: 674 =8.51×10−4 +2.53×10−3 ·i and 675 = 8.51 × 10−4 − 2.53 × 10−3 · i. This means that the reactor should be unstable with an oscillatory behavior. The simulation of the case of a heat transfer coefficient of 150 W/m2 /K was done to confirm these results (Fig. 7a). As for the first case, the temperature control was activated to attain a temperature of 380 ◦ C and after 0.5 h the control was turned-off. As indicated by the simulation
Controller ON
Stop controller 400 380 360 340 320 300 0:00
1:00
2:00
3:00 time (h)
4:00
5:00
6:00
Fig. 5. Dynamic simulation of LCO HDT with a heat transfer coefficient of 500 W/m2 /K.
results, the reactor starts to oscillate about 3.5 h after the controller shut-down. The minimum and maximum temperature oscillations are close to 340 and 440 ◦ C, respectively. A Fourier transform was applied on the temperature signal in order to find the period of the oscillation. The frequency spectrum of the signal obtained by simulation is shown in Fig. 7b. One of the main characteristic frequencies is F = 3.54 × 10−4 Hz that corresponds to a period of T = 47 min. The solution of the perturbation (Eq. (16)) shows that when an eigenvalue is complex the exponential term can be expressed as indicated ei ·t = eRe(i )·t [cos(Im(i )t) + i sin(Im(i )t)]
(19)
The eRe(i )·t term clearly reveals that if the real part is positive, there is an exponential amplification of the perturbation. The imaginary part indicates that the perturbation leads to an oscillating behavior. Therefore, the period of the oscillation can also be determined
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
1.5
0.005 0.004
1.0
λ674
0.003 0.002 Im(λi)
0.5 Im (λi)
319
0.0
0.001 0.000 -0.001
-0.5
λ675
-0.002 -0.003
-1.0
-0.004 -1.5 -7.5
-5.5
-3.5 Re (λi)
-1.5
0.5
-0.005 -0.005 -0.004 -0.003 -0.002 -0.001 0.000 0.001 Re (λi)
Fig. 6. Eigenvalue spectrum for U = 150 W/m2 /K: (a) all the eigenvalues and (b) zoom of values close to zero.
Temperature (°C)
440
3.0E+06 Temperature Controller ON Controller OFF
420 400
Stop controller
380 360 340
Fourier coefficients
460
F = 3.54 × 10−4 Hz
2.0E+06 1.0E+06 0.0E+00 -1.0E+06
320 300 0:00 2:00 4:00 6:00 8:00 10:00 time (h)
-2.0E+06 0
0.0005 0.001 0.0015 0.002 Frequency (Hz)
Fig. 7. (a) Dynamic simulation of LCO HDT with a heat transfer coefficient of 150 W/m2 /K and (b) Fourier transform frequency spectrum for a heat transfer coefficient of 150 W/m2 /K.
using the imaginary part of this eigenvalue. T=
2· 2· = = 42 min Im(i ) 2.53 × 10−3
This period should correspond to one characteristic frequency obtained from the Fourier transform of the temperature signal. This value (42 min) is very close to the value obtained from the frequency spectrum (47 min). The small difference is due to the accuracy of the linearization when oscillations are far from the stationary point. To confirm this assertion, another simulation was performed with a heat transfer coefficient of 200 W/m2 /K. This value corresponds to the stability limit. In this case two eigenvalues of the whole spectrum were complex with a positive real part: 9.62 × 10−5 + 2.64 × 10−3 · i and 9.62 × 10−5 − 2.64 × 10−3 · i. The period of oscillation calculated with the eigenvalues is 39.68 min. The Fourier transform calculations result in a frequency value of F = 4.2 × 10−4 Hz that corresponds to a period of 39.68 min. The period of oscillation calculated with both methods is in excellent agreement. For this case, the real part of the eigenvalues is one order of magnitude lower than the one obtained from the previous case with a heat transfer coefficient of 150 W/m2 /K. In other words, for the U = 200 W/m2 /K case, the temperature oscillations are composed in values of 10−3 ◦ C while for the
U = 150 W/m2 /K case, the temperature oscillations were of about 110 ◦ C. Then, the time constant of the amplification is = 1/Re(i ) = 173.3 min compared to 19.6 min for the U = 150 W/m2 /K case. As shown in this work, the analysis of the eigenvalue spectrum allows to determine if the system is stable or unstable. Moreover, the interpretation of the eigenvalues provides information about the propagation of the perturbation which may be of different types: 1—exponential attenuation, 2—oscillations with exponential attenuation, 3—oscillations with exponential amplification, and 4—exponential amplification. If the system is oscillating the frequency of the oscillation can be calculated with the imaginary part of the eigenvalues. 4. Conclusions This work presents a thermal stability study of a refining process (hydrotreatment of gasoils). This process consists on the hydrogenation of unsaturated compounds which can be present in high amounts in some gasoil feeds like LCO ( ≈ 90 wt%). The hydrogenation reactions are highly exothermic. The hydrotreatment of this great content of aromatics and olefins produces a high heat release that may lead to a runaway of the reactor. Therefore, a
320
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
thermal stability study of this process was carried out. The dynamic stability study is performed on the basis of a mathematical model that describes the reactive system. For this reason, a dynamic model of the hydrotreatment of gasoils was developed. The kinetics was fitted using experimental points obtained with industrial catalyst and industrial feeds in the classical range of operating conditions of this process. A dynamic model of HDT reactor was developed to simulate an existing pilot plant. The model takes into account a detailed description of the reactive system and the configuration of the pilot plant: three phases (gas–liquid–solid), the kinetics, the flow characterization in the gas and liquid phases, the heater devices and the control system. An experimental test was carried out in the pilot plant to collect dynamic data for the validation of the model. The dynamic temperature profiles of the reactor are in good agreement with the simulation results and the model was validated with this comparison. The thermal stability analysis of the LCO hydrotreatment was performed with a perturbation method. The criterion of stability imposes that after the solution of the perturbation system, all the resultant eigenvalues must have a negative real part, otherwise the reactor will be unstable. The influence of the heat transfer coefficient on the reactor stability was evaluated. Different cases were studied in this work, a stable case and two unstable cases. The conclusions obtained via the thermal stability study were compared to the conclusions obtained by model simulations. For all cases, the conclusions were always in good agreement. For the unstable cases, an oscillatory behavior was observed by simulation and the stability analysis successfully predicts the occurrence of the oscillations. Furthermore, the frequency and the period of these oscillations were compared to those obtained by Fourier transform signal treatment of the simulation results. In the case of small oscillations where the 1st order Taylor expansion is applicable, the results obtained by both methods were the same. The stability analysis presented in this work is therefore a reliable method to establish if the reactor is stable or not. In further work, the impact of the variation of other parameters such as the feed composition and flowrates of reactants will be achieved. The localization of the runaway ignition will be carried out by analyzing the eigenvectors. In addition, the impact of Jacobian accuracy will also be studied, the mesh and the numerical scheme have to be also evaluated. Finally, this stability analysis will be applied in transient conditions in order to determine the dynamic unstable regions. Notation a A b C Cp Dax E f
f
F G h H
H J k0 nz P
reaction order of hydrocarbon pseudo-component in the liquid phase heat exchange surface, m2 /m3 reaction order of hydrogen in the liquid phase Concentration, mol/m3 specific heat capacity, J/kg/K axial dispersion coefficient, m2 /s activation energy, J/mol second member of partial differential equations vector of the second member of partial differential equations frequency, Hz controller gain heater height, m Henry coefficient heat of reaction, J/mol Jacobian matrix reaction rate constant at the reference temperature number of steps for the reactor discretization pressure, bar
t T T Tamb Tref u U
total pressure, bar electrical power, W reaction rate, mol/s/kgcatalyst ideal gas constant, J/mol K reactor section surface, m2 section surface between the heaters and the reactor wall, m2 time, s temperature, K period, min ambient temperature, K reference temperature, K superficial liquid velocity, m/s heat transfer coefficient, W/m2 /K
V x
eigenvector perturbation
x
vector of perturbations
y z
variables vector axial coordinate, m
Pt Q r R Sr Svoid
Greek letters
¯ ax cont
holdup eigenvalue effective thermal conductivity, W/m/K stoichiometric coefficient density, kg/m3 time constant of the amplification, min controller time derivative constant
Subscripts and superscripts air C g heater i j k l loss n op p por s set steel wall W
air between the reactor wall and the heater, K cooling gas heater pseudo-component index or eigenvalue index reaction index number of heater index liquid relative to heat loss total number of variables operating point perturbed porosity solid setpoint reactor steel reactor external wall internal reactor thermowell
Abbreviations C DI IBP Im MB MONO R Re
complex diaromatics initial boiling point imaginary part material balance monoaromatics real real part
J.-M. Schweitzer et al. / Chemical Engineering Science 65 (2010) 313 -- 321
SAT SULF TB TRI
saturates sulfur compounds thermal balance triaromatics
Acknowledgments Authors are grateful to Mr. D. Hudebine for supplying a part of the theoretical data of LCO hydrotreatment. This work was partially supported by the “Axelera” French Research Project. References Agrawal, R., West, D.H., Balakotaiah, V., 2007. Chemical Engineering Science 62, 4926–4943.
321
Alós, M.A., Strozzi, F., Zaldívar, J.M., 1996. Chemical Engineering Science 51, 3089–3094. Balakotaiah, V., Luss, D., 1991. A.I.Ch.E. Journal 37 (12), 1780–1787. Balakotaiah, V., Kodra, D., Nguyen, D., 1995. Chemical Engineering Science 50 (7), 1149–1171. Bilous, O., Amundson, N.R., 1956. A.I.Ch.E. Journal 2 (1), 117–126. Luss, D., Sheintuch, M., 2005. Catalysis Today 105, 254–274. Morbidelli, M., Varma, A., 1986. A.I.Ch.E. Journal 32 (2), 297–306. Morbidelli, M., Varma, A., 1987. A.I.Ch.E. Journal 33 (12), 1949–1958. Perlmutter, D.D., 1972. Stability of Chemical Reactors. Prentice-Hall, New Jersey. van Heerden, C., 1953. Industrial and Engineering Chemistry 45, 1242–1247. Viswanathan, G.A., Luss, D., 2006. Industrial and Engineering Chemistry Research 45, 7057–7066. Whitehurst, D.D., Isoda, T., Mochida, I., 1998. Advances in Catalysis 42, 345–471. Wu, H., Rota, R., Morbidelli, M., Varma, A., 1999. Chemical Engineering Science 54 (20), 4579–4588. Zaldívar, J.M., Cano, J., Alós, M.A., 2003. Journal of Loss Prevention in the Process Industries 16, 187–200.