Systematic Procedures for Interpretation and Modeling of Catalyst

Feb 23, 2011 - ABSTRACT: The method of decreasing conversion, constant space-time (DCCST) is sometimes used for empirical determination of catalyst ...
0 downloads 0 Views 2MB Size
ARTICLE pubs.acs.org/IECR

Systematic Procedures for Interpretation and Modeling of Catalyst Deactivation Using Integral Fixed-Bed Reactors John J. Birtill* Department of Chemistry, University of Glasgow, Glasgow G12 8QQ, Scotland, United Kingdom * ABSTRACT: The method of decreasing conversion, constant space-time (DCCST) is sometimes used for empirical determination of catalyst decay kinetics in fixed beds. The method of constant conversion, increasing space-time (CCIST) is used infrequently. Although it is well known that these methods are most readily applied to cases in which the deactivation is independent of the composition of the fluid phase, the ease with which they can be misapplied to concentration-dependent cases has not been studied systematically. In this work, the empirical determination of catalyst decay kinetics in fixed beds has been examined by numerical simulation. It has been shown that it is rather easy to generate an invalid concentration-independent decay model, either by fitting a dubious straight line to a limited range of data using the DCCST method or by fitting a good, or even perfect, straight line using the CCIST method. It has also been shown that both methods can, in principle, be used reliably for empirical determination of catalyst decay kinetics, whatever the dependence on composition, provided that multiple experiments are carried out over a range of conversion. The CCIST method has much greater potential for resolving reaction and decay kinetics than has been realized previously, and it can be applied to reversible as well as irreversible reactions. This method will generally give a linear plot, provided that the selected plot corresponds to the actual decay order. The plot is exactly linear for any first-order decay kinetics and almost linear for non-first-order decay kinetics. Exact analytical solutions have been derived for various first-order cases. These solutions can be used to determine the true decay kinetics from multiple parallel experiments. Deviation from ideal behavior in this method is a useful indication that the power-law decay kinetics has a limited range of validity, beyond which a more detailed, systematic investigation may be required.

1. INTRODUCTION The fundamental mechanistic description of the decay kinetics of a solid catalyst is usually rather complex, even if pore-diffusion effects can be ignored, and so extensive, detailed studies over the catalyst lifetime of both catalytic performance and material characterization are required for the determination of any fundamental decay model. In consequence, empirical activity power-law models, as described by Levenspiel,1,2 are widely used for practical purposes. Such power laws are, of course, simplifications of fundamental mechanistic expressions, which are valid under the normal range of operational conditions. This applies to deactivation by parallel poisoning3,4 and series poisoning,3 driven by reactants and products, respectively. In some cases of sintering, the fundamental mechanistic decay kinetics can be reduced to simple power laws with high deactivation order.5 The combination of empirical coking kinetics with empirical activity-coke expressions can also be expressed in terms of activity power-law decay kinetics.6,7 Hence, simple powerlaw expressions are generally the most useful and direct way of describing decay kinetics. Indeed, the fitting of complex mechanistic models with poorly determined coefficients or insignificant parameters may be regarded as misleading overelaboration. However, if the decay kinetics is too approximate or oversimplified, then this is also potentially misleading and an industrial catalyst operated outside the narrow range of research test conditions may suffer unexpectedly rapid deactivation. This is usually an expensive mistake in commercial operation. Hence, catalyst test procedures should always be designed so as to avoid both overelaborate and oversimplified decay kinetics. r 2011 American Chemical Society

Experimental techniques for the practical measurement of catalyst decay have been reviewed recently, including a discussion on the interrelation of the effects of pore-diffusion resistance and decay.7 The present paper is concerned with the influence of the composition of the fluid phase on the decay kinetics at the particle level, and it is assumed that the decay kinetics is not influenced by pore-diffusion resistance. For test purposes the catalyst will usually be in the form of small particles of uniform composition and with a suitable size range relative to the reactor dimensions so as to avoid deviation from plug flow. The scale-up of decay kinetics to describe large pellets with significant porediffusion resistance will require additional testing and/or modeling.7 It is also assumed that the deactivation and reaction kinetics are practically separable.1,2 This assumption might prove to be invalid, and so this point will be addressed again in the discussion. Decay kinetics, which is independent of the fluid composition, may be described in the form of a power law as in eq 1. -

da ¼ kd am dt

ð1Þ

The ‘activity’ a refers to the reaction rate at any point in the catalyst bed relative to fresh catalyst under identical conditions.1,2 Received: July 11, 2010 Accepted: December 23, 2010 Revised: December 10, 2010 Published: February 23, 2011 3145

dx.doi.org/10.1021/ie101481y | Ind. Eng. Chem. Res. 2011, 50, 3145–3158

Industrial & Engineering Chemistry Research

ARTICLE

The equivalent integral form of eq 1, the ‘time-on-stream’ model for catalyst activity, can also be used. If the decay kinetics is dependent on the concentration (Ci) of one or more components of the fluid phase, then eq 1 is obviously insufficient and it is necessary to incorporate a function to describe the concentration dependence as in eq 2. This can be an additional concentration power-law term1,2 or a more general function of composition.8 -

da ¼ kd ΨðCÞam dt

ð2Þ

Levenspiel described methods for investigating power-law decay kinetics.1,2 For concentration-dependent decay, the separate effects of fluid composition on the reaction kinetics and on the decay kinetics can be decoupled by following suitable experimental strategies in a gradientless reactor, so as to maintain constant concentration for selected components. However, due to the lack of suitable equipment and/or time, in practice many deactivation studies are carried out using fixed-bed plug-flow reactors, in some cases without sufficient checks for composition dependence. For true concentration-independent decay, the deactivation rate constant kd can be determined, using integral data from a plug-flow reactor, by plotting the appropriate linear function of activity against time. The integral conversion (X) at constant feed rate and constant temperature will decrease with time as the catalyst deactivates. For an irreversible reaction with no change in fluid density (ΔVR = 0), first-order reaction kinetics, and firstorder concentration-independent decay kinetics, it can be shown easily that the exit concentration of reactant A will decline according to eq 31,2 lnðlnðCA0 =CA ÞÞ ¼ lnðkr τÞ - kd t

ð3Þ

where kr = reaction rate constant for fresh catalyst, CA0 = inlet concentration, and (weight) space-time τ = WCA0/FA0. The left side of this equation can also be written as the equivalent function of conversion ln(ln(1/(1 - X))) or as ln(activity). A linear plot against time is consistent with concentration-independent firstorder deactivation and can be used to determine kd.1,2 This analysis can easily be extended to derive alternative plots for different combinations of irreversible reaction kinetics and decay order. For instance, eq 4 will apply for first-order reaction kinetics with constant fluid density and second-order decay kinetics. lnðCA0 =CA Þ ¼ kr τ=ð1 þ kd tÞ

ð4Þ

In comparative screening tests of catalyst life, the integral conversion is often held constant by gradually increasing the reaction temperature (deactivation compensation), especially when this is the normal practice in commercial operation. However, this method, which will be considered further in the discussion, introduces the additional complexity of the temperature responses of both the main reaction and the deactivation process, and so it is not suitable for determination of decay kinetics. Alternatively, for the purpose of research, both temperature and conversion can be held constant if the feed rate is reduced throughout the experiment (i.e., increasing space-time τ).1,2 Again assuming concentration-independent decay kinetics (no time-dependent axial activity gradient) but making no assumption regarding fluid density,

it can be shown that WCA0 -1 ¼ τ ¼ kr at FA0

Z

CA

CA0

dCA -I ¼ kr at ΛðCA Þ

ð5Þ

where at is the activity at time t, Λ(CA) describes the reaction kinetics, and I is the reaction integral for the corresponding composition function over the length of the catalyst bed. Equation 5 applies to reversible reactions as well as nonreversible reactions, although operation at less than equilibrium conversion would be necessary in practice to ensure an exactly constant integral. The expression for the variation of space-time over real time can then be generated by insertion of the appropriate definition of activity as a function of time on stream, the integrated form of eq 1. For first-order decay (m = 1) lnðτÞ ¼ kd t þ b

ð6Þ

For second-order decay (m = 2) lnðτÞ ¼ lnð1 þ kd tÞ þ b

ð7aÞ

τ ¼ ckd t þ c

ð7bÞ

or

where b = ln(-I/kr) = ln(τi) and c = (-I/kr). In general, for m 6¼ 1 (based on the general solution of Wojciechowski for activity against time on stream 9) lnðτÞ ¼ ð1=ðm - 1ÞÞlnð1 þ ðm - 1Þkd tÞ þ b

ð8aÞ

τ ¼ cð1þðm - 1Þkd tÞ1=ðm - 1Þ

ð8bÞ

Hence, a linear plot of ln(τ) against time on stream is consistent with concentration-independent first-order decay1,2 and a linear plot of τ against time is consistent with concentration-independent second-order decay. At t = 0, either plot should pass through the intercept, which is defined by the reaction kinetics and reaction conditions. For first-order reaction kinetics and constant fluid density, Λ(CA) = CA b ¼ lnðlnðCA0 =CA Þ=kr Þ ¼ lnðlnð1=1 - X A Þ=kr Þ

ð9Þ

If the fluid density is not constant (ΔVR 6¼ 0), then linear plots will still apply.10 However, the reaction integral I in eq 5 above will need to incorporate the change in fluid density. The problem with these methods is that a linear plot derived from integral conversion data over a fixed catalyst bed may be consistent with concentration-independent decay, but it is not proof of this.7 If an axial activity gradient is actually developing over time, the effect may be obscured by the resultant change in the concentration gradient, so that the estimated ‘activity’ is effectively a mean value over the whole catalyst bed. It has previously been reported that eq 6 will give perfectly linear plots for some cases of concentration-dependent decay.7 The problem of the hidden activity gradient has previously been recognized or inferred in the literature,11-13 but little systematic examination has been reported. Froment and Bischoff simulated cases of series and parallel coking kinetics, combined with empirical activity-coke relationships, and examined the effect of an axial coke profile on the observed average bed activity.14 They warned that “the distance-averaged rate coefficient 3146

dx.doi.org/10.1021/ie101481y |Ind. Eng. Chem. Res. 2011, 50, 3145–3158

Industrial & Engineering Chemistry Research at a given carbon level is not the same as the rate coefficient evaluated at the distance-averaged carbon level”. It has been noted already that the activity-power-law approach and the activity-coke approach are usually interchangeable, and so they share the same weakness when misapplied to plug flow reactors with timedependent activity gradients.7 Santos et al. simulated two cases of parallel, concentrationdependent deactivation, based on power-law decay kinetics, for a fixed-bed reactor.13 They concluded that reasonable values for kd could be determined by the decreasing conversion method, apparently from differential examination of the simulated axial conversion profiles over time, although the methodology was not entirely clear. They noted that incorrect values were determined when the methodology was applied to the deactivation-compensation method. The study did not include practical methodology for determination of the axial conversion profiles over time. The composition dependence of deactivation can be identified by carrying out tests with controlled variation of space-time and/or fluid composition.7,15,16 However, the possibility of hidden deactivation gradients due to composition dependence is sometimes overlooked altogether. The first aim of this work was to force fit equations for concentration-independent decay to simulated cases of concentration-dependent decay and thereby examine how easy it would be to determine invalid concentration-independent decay kinetics (or the equivalent integrated time on stream models) from concentration-dependent deactivation data. The second aim was to design and evaluate easy, rapid, and robust procedures for determining approximate decay kinetics. Some initial results have been presented previously.7,17,18 A comprehensive description is presented in this paper.

2. MODELING METHOD The simulations were carried out in multiple cycles of small, sequential steps in space-time along an isothermal plug flow reactor with each complete reactor cycle separated by a short step in real time. For tests with constant space-time and decreasing conversion, the reaction was modeled over 100 consecutive reactor segments, assuming first- or second-order irreversible reaction kinetics. The initial activity was set at 1.0, and either the reaction rate constant kr or the flow rate was adjusted to set the initial conversion to the target value. In most cases the deactivation rate was modeled by power laws of the form described by eq 2 but with power-law terms for reactant CA and/or product CB. During each time interval, the deactivation in each reactor segment was calculated using the mean concentration values in that segment. For the next time interval, the activity values were adjusted accordingly and the concentration in each segment was recalculated, followed again by the deactivation. The calculation cycle was repeated over 100 successive time intervals for each complete simulation. For irreversible poisoning by an impurity, the deactivation was modeled with a constant, very low inlet poison concentration CP0 and with progressive poison consumption and hence progressive deactivation of individual catalyst sites along the catalyst bed. This may be compared to a slow titration process for active sites, with a sharp deactivation front for a fast reacting poison (converted immediately on contact with an active site), or a front that is smeared out over the bed for a moderately reactive poison (kr ≈ kp) or a slowly reacting poison (kr > kp). These simple models are

ARTICLE

valid for any ratio of sites deactivated per poison molecule provided that this ratio is constant but may not be appropriate for some more complex cases of poisoning. In this work, the main interest is in the diagnosis of catalyst deactivation for the main reaction. The breakthrough behavior of a moderately reactive poison can also be used for direct characterization of poisoning kinetics.19 For tests with constant conversion and increasing spacetime, similar procedures were used, except that the flow rate was decreased after each time interval so as to maintain constant overall conversion at the end of the catalyst bed. An iterative procedure was used to determine the correct adjustments to flow rate after each time interval. The time intervals were adjusted so as to achieve the target degree of deactivation in 100 cycles. Case 15 was simulated with smaller bed segments at the front and rear in order to reduce the step error in the region where the activity gradient was most severe. The data were plotted in accordance with the equation selected for examination. The departure from linearity was at first tested by visual inspection of each plot. The data were then fitted by linear regression to a straight line, constrained to pass through the initial data point. The data were simulated without any experimental uncertainty, and so a rigorous analysis for curvature was inappropriate. The value of R2, the ‘goodness of fit’, was used as an indicative index. By definition, R2 = 1 Σ(Yobs - Yline)2/Σ(Yobs - Ymean)2 over all the data points. Basically, the deviation of the line is compared to the deviation of a horizontal line drawn through the mean Y value (the ‘null hypothesis’). When R2 = 1.0, all points lie exactly on the straight line with no scatter. The higher the value of R2 for these simulated data, the more difficult it would be to detect curvature in real data, which are subject to experimental error and usually more limited in quantity. For these noise-free simulated data, very slight curvature could be detected by visual inspection even for R2 as high as 0.995. As an approximate guide, the curvature was generally obvious for R2 ≈ 0.98 and very obvious for R2 < 0.97 and so should be evident with real data with moderate scatter, although this also depends on the exact shape of the curve.

3. DECREASING CONVERSION, CONSTANT SPACETIME (DCCST) METHOD Regardless of the decay kinetics which had actually been used to generate the integral conversion-time data, most data sets were tested against eq 3, which is valid only for a first-order reaction with concentration-independent decay kinetics and decay order m = 1. A single set of second-order reaction data with decay order m = 1 was tested against the corresponding plot of ln(1/CA - 1/CA0) against time. The results of the simulations are shown in Table 1. Provided that the correct equation has been used for the decay order and reaction kinetics, true concentration-independent decay should produce a straight line, with a negative gradient in the case of eq 3. The data were fitted by linear regression with the intercept constrained to the theoretical initial value, ln(krτ). The degree of curvature and its direction were recorded for a standard drop in conversion from 98% to 50%. The examples in Figure 1a illustrate a range of behavior from slight to a high degree of curvature. The corresponding activity profiles at the end of the simulation (conversion = 50%) are shown in Figure 1b and 1c. It is noteworthy and cautionary that different types of decay kinetics can produce near-linear plots in Figure 1a even 3147

dx.doi.org/10.1021/ie101481y |Ind. Eng. Chem. Res. 2011, 50, 3145–3158

Industrial & Engineering Chemistry Research

ARTICLE

Table 1. Appearance of Trial Plots for the DCCST Methoda decay order case

decay model

am, m =

Cn, n =

deviation from linearity for

final % conversion

conversion drop from 98% to 50%

for linear plot

severity

first-order reaction

direction

R2

R2 = 0.990

R2 = 0.980

trial plot: ln ln (CA0/CA) vs t

4

parallel

1

1

high

down

0.959

87.7

76.8

9

parallel

2

1

slight

up

0.993

46.1

38.4

5

parallel

1

2

high

down

0.934

89.0

83.0

11

series

1

1

high

up

0.955

74.6

64.4

13

series

2

1

very high

up

0.485

94.2

91.9

12

series

1

2

high

up

0.803

86.6

81.1

14.1 14.2

mixed ParþSer (h = 1.0)b mixed ParþSer (h = 0.5)b

1 1

1 1

none slight

none down

0.996

5.3

14.3

mixed ParþSer (h = 0.2)b

1

1

high

down

0.983

69.1

42.2

15

mixed Par0.5Ser0.5 b

1

1

slight

downc

0.996

12.4

9.2

irreversible poisoning titrate active sitesd

16

fast poison

high

down

0.924

90.1

83.4

17

moderate poison

high

down

0.957

88.3

78.0

18

slow poison

none

none

1.000

0.87

95.9

94.3

0.75

96.6

95.7

trial plot: ln(1/CA - 1/CA0) vs te

second-order reaction 19

parallel

1

1

very high

down

trial plot: ln ln(CA0/CA)) vs t

f

All cases: initial conversion =98.17%, 100 time intervals and 100 bed intervals. Intercept constrained at ln krτ. b Cases 14.1-14.3 describe mixed parallel/series decay with additive dependence with the same order on reactant A and product B but different relative values for kd, with h = kdB/kdA. Case 14.1 is concentration pseudoindependent, and so gives a linear plot. Case 15 describes the multiplicative half-order concentration dependence on both reactant and product with overall composition decay order n = 0.5  2 = 1. c Curve changes direction at a higher degree of decay for case 15. d Poisoning cases 16-18 were modeled as progressive poison consumption along catalyst bed with deactivation of individual catalyst sites, c.f., a titration process: sharp for a fast poison (converted immediately on contact with an active site), smeared out over the bed for a moderate (kr ≈ kp) or a slow (kr . kp) poison. e Correct reaction kinetics. f Incorrect reaction kinetics. a

though the resulting axial decay profiles in Figure 1b and 1c are very different. The falls in conversion corresponding to two selected target values for R2 (0.990 and 0.980) were also recorded. In order to be confident in the validity of a linear fit, the decreasing conversion must be followed over a sufficiently wide range so that any divergence from nonlinearity can be detected, as explained in Table 1 and shown by the comparative plots in Figure 1a. In many cases, such as fast poisoning, series decay, and parallel firstorder decay, the severity of the curvature (Table 1, column 5) is high, and so the deviation from linearity should become apparent even with noisy real data after a moderate drop in conversion over time (see values for R2 = 0.980 in column 9). In other cases, the severity of the curvature is slight, and so the deviation from linearity may not be obvious even after a large drop in conversion. Hence, it is not possible to discriminate with certainty between concentration-independent decay and all types of concentrationdependent decay with this type of plot. For parallel decay with order m = 2 (case 9), the fall in conversion would need to be followed down to below 40% in order to see the deviation. Note that in this case the plot is invalid for two reasons: the decay is not independent of concentration and the decay order m 6¼ 1. The plot of eq 3 also appeared linear over wide ranges for cases of mixed parallel/series decay, driven by both reactant and product. This is not surprising for additive parallel/series decay because the axial decay gradient is small, but it is surprising for the

multiplicative case 15 which develops a shallow activity minimum in the midbed region, as shown in Figure 1b. The simulation of fast, irreversible poisoning (case 16) produced a curved plot after quite a small drop in conversion. For a drop from 98% to 50% the plot was similar to that of parallel decay (case 4) in Figure 1a, although with sharper curvature, but the deactivation profile in Figure 1c was obviously much more abrupt. Fast poisoning at a constant rate reduces the effective amount (W) of catalyst in the bed at a constant rate, and so the alternative simulated plot of ln(ln(1/(1 - Xi))) against axial catalyst charge (W) at t = 0 produced a matching result. The simulation of moderate poisoning produced a plot which was similar in shape to parallel decay, and as expected, the match to parallel decay was exact when the rate constants for conversion of reactant and poison were set equal. The simulation of slow, irreversible poisoning gave a near linear plot because the axial decay gradient was shallow, and so concentration-independent decay is a reasonable approximation in this case, provided that CP0 remains constant. The nature of the curvature in this type of plot provides a clue to the nature of the decay kinetics. Strong downward (concave) curvature indicates first-order parallel decay or feed poisoning, whereas strong upward (convex) curvature indicates series decay. Alternatively, the change in gradient in the latter case could also indicate a change in the decay mechanism after an initial period of rapid deactivation.11 Equation 3 is based on the assumption that kd is constant in all parts of the bed. The downward deviation of 3148

dx.doi.org/10.1021/ie101481y |Ind. Eng. Chem. Res. 2011, 50, 3145–3158

Industrial & Engineering Chemistry Research

ARTICLE

Figure 1. Decreasing conversion, constant space-time. Drop in conversion from 98% to 50%. (a) Deviation from the linearity of eq 3. Plot is equivalent to ln(Æaækrτ) against time. Deviation from the initial straight line, degree, and direction (Table 1): case 4 (par, m = n = 1) high, down; case 15 (mixed par0.5  ser0.5) very slight, down; case 9 (par, m = 2, n = 1) slight, up; case 11 (ser, m = n = 1) high, up; case 13 (ser, m = 2, n = 1) very high, up. (b) Final activity profile after conversion drop to 50%. (c) Final activity profile for parallel decay, fast and moderate poisoning decay. Note: Moderate poisoning case 17c has been adjusted to avoid exact overlap with parallel decay case 4 by setting kp = 0.8kr.

eq 3 for parallel decay is expected because the ‘effective’ value of kd in eq 3 also depends on the local reactant concentration CA. The quantity kdCA increases in each segment as the catalyst decays, and so the overall bed function, ln ln(1/(1 - XA))  ln(Æaækrτ), decreases faster than expected, leading to downward deviation from linearity. For similar reasons, for series decay the function decreases more slowly than expected. For fast poisoning, the deactivation rate does not decrease with the decreasing overall bed activity. The rate of loss of active sites remains constant, although the deactivation front moves progressively down the bed. Hence, the function decreases faster than expected.

The overall conclusion is that a single plot of eq 3, which is equivalent to a plot of ln(Æaækrτ) against time, is unreliable for determination of concentration-independent decay kinetics or the corresponding ‘time-on-stream’ expressions for activity.

4. DECLINING CONVERSION METHOD WITH VARIATION OF INITIAL CONVERSION The DCCST method can be used for diagnostic and comparative purposes by running multiple tests (at least 3) in parallel 3149

dx.doi.org/10.1021/ie101481y |Ind. Eng. Chem. Res. 2011, 50, 3145–3158

Industrial & Engineering Chemistry Research

ARTICLE

Figure 2. Decreasing conversion, constant space-time with variation of initial conversion. (a) Variation of initial ‘effective’ kd with initial conversion. (b) Variation of initial ‘effective’ kd with flow rate or 1/τ. For both figures, the actual kd values used in simulations have been selected so that all plots lie on the same scale. The plots show approximate values for the simulated effective kd from the early near-linear region of eq 3. At t = 0, the exact value of the effective kd = kdi for power law decay kinetics with any m. Key: Mix PþS denotes additive parallel and series decay with h = kdB/kdA = 0.2 (case 14.3). For fast poison, s denotes simulated values and c denotes exactly calculated values.

on each catalyst, covering a wide range of initial conversion values by setting different values of space-time W/F (by variation of catalyst charge or feed rate). Any composition dependence of the decay kinetics will in due course become evident. Each plot of eq 3 will have a different gradient (or ‘effective’ kd), and the values will vary systematically with the initial conversion. The initial values of the ‘effective’ kd can be measured approximately in the early near-linear region. Simulated values are plotted against initial conversion in Figure 2a and against flow rate or reciprocal space-time in Figure 2b. These plots show that it should be possible to discriminate different types of decay kinetics to some extent. The main types, concentration-independent, parallel, series, and fast poisoning, show different characteristic profiles. However, the decay order m is not indicated by these plots, as shown by the similar curves for parallel decay with m = 1 or 2. The plot profiles do not provide a definitive indication of the composition dependence of the decay kinetics. Additive dependence on reactant and product concentrations (mixed parallel þ series) can resemble parallel, concentration-independent, or series decay with decreasing value

of parameter h, the ratio of the kd values for product and reactant. If a poison is moderately reactive (kp ≈ kr) then its profile will resemble parallel decay, whereas a slowly reacting poison will resemble concentration-independent decay. In some cases, the definitive composition dependence might require systematic variation of individual concentration terms in a well-mixed reactor.1,2 The initial values for the ‘effective’ kd in eq 3 can be calculated for fast, irreversible poisoning as this is essentially equivalent to a reduction in the quantity W of active catalyst caused by a poison front, moving progressively down the bed. The slight discrepancy between simulated and calculated values of ‘effective’ kd for fast poisoning in Figure 2a and 2b is due to slight inaccuracy in the simulations. In Figure 2b, fast poisoning gives a straight line with gradient QP for the plot against 1/τ, as shown by eq 10 in Appendix 1.a. The concept of describing catalyst deactivation in terms of a decreasing quantity of catalyst has also been used recently for a model of parallel decay of zeolite catalyst in the conversion of methanol to hydrocarbons, based on the breakthrough curves of methanol over time.20 3150

dx.doi.org/10.1021/ie101481y |Ind. Eng. Chem. Res. 2011, 50, 3145–3158

Industrial & Engineering Chemistry Research

ARTICLE

Table 2. Appearance of Trial Plots for the CCIST Methoda plot vs t

decay order

case

model

am, m

Cn, n

conv %

m = 1.5, source

m = 2, source eq 7b

ln τ, R2 c

eq 8b √ 2c τ, R

m = 1, source eq 6

Ff/Fib

m = 0.5, source eq 8b √ 1/ τ, R2 c

τ, R2c

notesd

main reaction: first-order irreversible 1.1

independent

0.5

0

98.17

0.902

1.0000

0.9999

0.9996

1.2

independent

0.5

0

98.17

0.25

1.0000

0.985

0.942

0.9991 0.88

2.2

independent

1

0

98.17

0.25

0.981

1.0000

0.984

0.941

3.1 6.1

independent parallel

2 0.5

0 1

98.17 98.17

0.25 0.50

0.76 0.998

0.906 0.989

0.979 0.972

1.0000 0.95

6.2

parallel

0.5

1

98.17

0.25

0.999

0.987

0.949

0.89

6.4

parallel

0.5

1

22.12

0.25

1.0000

0.985

0.942

0.88

7.1

parallel

0.5

2

98.17

0.50

0.998

0.991

0.977

0.957

7.2

parallel

0.5

2

98.17

0.25

0.994

0.996

0.969

0.92

i

7.3

parallel

0.5

2

98.17

0.071

0.938

0.999

0.946

0.83

i

7.4

parallel

0.5

2

22.12

0.25

1.0000

0.984

0.941

0.88

22.1 22.1

series series

0.5 0.5

2 2

98.17 22.12

0.25 0.25

0.9997 0.9991

0.980 0.980

0.935 0.937

0.87 0.88

any n

any conv

any concentration dependence

8.1

parallel

1

0.75

0.9992

1.0000

0.9992

0.997

ii

series or

1

0.5

0.995

1.0000

0.996

0.984

ii

mixed P/S

1

0.25

0.981

1.0000

0.984

0.941

ii

parallel

1.5

0.25

0.90

0.975

0.9998

0.985

1

98.17

8.2

parallel

1.5

1

22.12

0.25

0.91

0.979

1.0000

0.984

9.1 9.2

parallel parallel

2 2

1 1

98.17 22.12

0.25 0.25

0.73 0.76

0.891 0.905

0.974 0.979

0.9997 1.0000

10.1

parallel

2

2

98.17

0.455

0.87

0.932

0.972

0.994

10.3

parallel

2

2

98.17

0.25

0.71

0.882

0.971

0.999

13.2

series

2

1

98.17

0.25

0.75

0.901

0.977

0.99997

13.3

series

2

1

63.21

0.25

0.75

0.899

0.976

0.99992

20

parallel

7

1

98.17

0.61

0.17

0.33

0.47

0.59

iii

21

series

7

1

98.17

0.60

0.19

0.35

0.49

0.61

iii

16.1 16.2

fast poisoning fast poisoning

98.17 98.17

0.378 0.25

0.991 0.981

1.0000 1.0000

0.992 0.984

0.969 0.940

iv iv

17

mod poisoning

98.17

0.25

0.981

1.0000

0.984

0.940

iv

18

slow poisoning

98.17

0.23

0.981

1.0000

0.983

0.940

iv

main reaction: first-order reversible 23

parallel

1

1

92.96

0.25

0.981

1.0000

0.984

0.940

v

24

series

1

1

92.96

0.25

0.981

1.0000

0.984

0.940

v

25

parallel

2

1

92.96

0.25

0.750

0.900

0.978

0.99992

v

98.21

0.25

0.981

1.0000

0.984

0.941

main reaction: second-order irreversible 19

parallel

1

1

All cases: ΔVR = 0, P = parallel, S = series (CB= CA0 - CA). b Ff/Fi: Final flow/initial flow. c Intercept constrained to the initial value for all plots. R2 values in bold correspond to correct m. d Other notes: (i) m = 0.5, n = 2. Wrong order m = 1 gives a better fit than the correct order at high conversion and high decay. Doubling the number of bed steps in the simulation made no difference. (ii) Applies to all cases in Table 1 with m = 1. (iii) m = 7. R2 > 0.999 for a plot of τ6 vs t. (iv) Modeled as progressive poison consumption along catalyst bed with deactivation of individual catalyst sites (c.f. titration). Fast poison: kp . kr. Moderate poison: kp = kr. Slow poison kp , kr (v) Reversible reaction held at 99.0% of equilibrium conversion. a

5. CONSTANT CONVERSION, INCREASING SPACETIME (CCIST) METHOD The progressive deactivation data were simulated, with decay order set in the range 0.5 e m e 2.0 and in two cases with m = 7. Many cases were examined. Some selected results are shown in

Table 2. The main reaction kinetics was, as indicated, irreversible first or second order or reversible first order. In most cases the deactivation was followed over a wide range, as shown by the ratio of the final to initial flow rate (Ff /Fi). The deactivation data were tested against eqs 6-8b, which are valid 3151

dx.doi.org/10.1021/ie101481y |Ind. Eng. Chem. Res. 2011, 50, 3145–3158

Industrial & Engineering Chemistry Research

ARTICLE

Figure 3. CCIST method. Case 4: Parallel decay, m = 1. Trial plots for m = 0.5, 1.0, 1.5, 2.0. Notes: 98.17% conversion, Ff/Fi = 0.25, least-squares fits constrained to pass through starting values. For 22.1% conversion, a similar degree of model discrimination is achieved in the shorter time of 27.5 units.

Figure 4. Simulated final activity profiles for constant conversion, increasing space-time. Notes: 98.17% conversion, Ff/Fi = 0.25.

for concentration-independent decay kinetics for particular values of the decay order m. The values for ‘goodness of fit’ R2 are shown for the plots against time of four different functions of space-time τ which apply to m = 0.5, 1.0. 1.5, and 2.0, the R2 values in bold corresponding to the ‘correct’equations for decay order m. In each case the best fit was constrained to pass through the intercept at t = 0 for the respective function of space-time. Once again, in the absence of experimental error, even moderate deviation from linearity was evident by visual inspection. As expected, the plots for concentration-independent decay were perfectly linear. However, the plots for cases of concentration-dependent decay also appeared to be linear over a wide range of deactivation, provided that the ‘correct’ equation was used for the simulated decay order m. This can be seen in each case by comparing the R2 value in bold for the ‘correct’ order equation with the R2 values for the other equations. It was also found that for m = 1 the deviation from linearity for wrong-order plots was similar in all cases when followed over the same degree of deactivation. The results for parallel decay with

m = 1 are plotted in Figure 3 with degree of deactivation Ff /Fi = 0.25 (i.e., apparent activity loss ≈ 75%). Only the correct plot for m = 1 is clearly linear. The plots were similar over this range for other cases with m = 1, i.e., series, parallel, or mixed composition dependence, regardless of conversion and reaction order. This is because all the plots for m = 1 are perfectly linear, and so the effect of fitting the wrong equation over the same range of values (Ff/Fi) is the same in all cases. The discrimination of decay order was less obvious when the catalyst had deactivated to a smaller extent, but for concentration-dependent decay, the deactivation rate can be increased in order to improve the discrimination of decay order, in the case of parallel decay, by operation at low conversion (see notes to Figure 3). The activity profiles within the catalyst beds with the same apparent activity loss are shown for selected cases in Figure 4. These all show severe axial deactivation profiles despite the linearity of the plots corresponding to concentration-independent decay. On closer inspection, √ √the apparently linear plots against time for the functions 1/ τ, τ, or τ, appropriate to decay order m = 0.5, 1.5, and 2, respectively, were found to be not exactly linear. This is indicated for cases with m 6¼ 1 in the respective columns of Table 2 by slight deviations in R2 from exactly 1.000 for the ‘correct’ plots. These slight deviations from linearity could not be corrected by simulation with smaller increments of catalyst mass or time. However, in most cases the deviation was so small that it would in practice be obscured by experimental error. Hence, the CCIST method of deactivation-compensation effectively produces almost linear plots even though there is an activity gradient within the catalyst bed. The mathematical reason for this outcome is investigated later in this section. The CCIST method can therefore, in principle, be used to determine the decay order m. In a wide range of cases with m = 2, followed over a similar range of deactivation (Ff/Fi = 0.25), the deviation from linearity for any plot for the wrong decay order was similar in each case and the discrimination of (0.5 in m was 3152

dx.doi.org/10.1021/ie101481y |Ind. Eng. Chem. Res. 2011, 50, 3145–3158

Industrial & Engineering Chemistry Research

ARTICLE

Table 3. CCIST Method: Some Analytical Solutions for First-Order Concentration-Dependent Decaya m

case

a

n

analytical solution for gradient = ‘effective’ kd

reaction order

4

parallel

1

1

1

-kd(CA0 - CA)/ln(CA/CA0)

5

parallel

1

2

1

-kd(CA02 - CA2)/(2 ln(CA/CA0))

11

series

1

1

1

kd(1 þ XA/ln(1 - XA))

12

series

1

2

1

kd[CA02 ln(CA/CA0) - 2CACA0 þ 3CA02/2 þ CA2/2]/ln(CA/CA0)

14

Par þ Ser CA þ hCB

1

1

1

-kdACA0[(1 - h)XA/ln(1 - XA) - h]

15

Par  Ser CA0.5  CB0.5

1

2  0.5 =1

1

-kdCA0[Arcsin(XA0.5) - (XA - XA2)0.5]/ln(1 - XA)

19

parallel

1

1

2

-kd ln(CA/CA0)/(1/CA-1/CA0)

23 cf 16

parallel fast poison

1

1

1 reversible 1

(kdCA0/(kr þ kr0 ))*[(kr þ kr0 )(CA - CA0)/ln(((kr þ kr0 )CA - kr0 CA0)/krCA0) þ kr0 CA0] CPNAV/SWi

Notes: Solutions to eq 6. All cases: ΔVR = 0. Irreversible reactions except case 23; ‘effective’ kd = kdi; intercept = ln τi. Case 14: h = kdB/kdA.

similar to that shown in Figure 3 for m = 1. However, for parallel decay with m = 0.5 and n = 2 (case 7), the discrimination between m = 0.5 and 1 was unreliable at high conversion. In this single case the fit with m = 0.5 became less good with increasing deactivation. The deviation became visually evident at high deactivation (Ff/Fi < 0.25), and the incorrect equation for m = 1 thereafter gave a better fit. This problem did not occur in case 7.4 at low conversion or in case 22 with series decay. The severe concentration dependence (n = 2) of decay in case 7 should still allow model discrimination by more sophisticated analysis over a range of conversion. In other cases examined, the discrimination of the correct decay order m became more obvious with increasing degree of deactivation. The analysis was also extended to higher decay order, m = 7, for parallel and series decay. Once again, the correct plot for the decay order, τ6 against t, was almost perfectly linear in both cases, although the discrimination of the exact decay order was more like (1 at Ff/Fi = 0.25. In general, deactivation compensation by the CCIST method has the effect of suppressing the influence of concentration on the linearity of the corresponding space-time plot for any decay order. The reason why plots with decay order m = 1 are perfectly linear can be understood from the following argument, based on site numbers rather than site density. The degree of conversion in any segment δW of a fixed bed, centered at point z along the bed, is given by FA0δXA = -rA0 δW = -krΛ0 (XA)azδW. FA0 is progressively decreased in order to compensate for deactivation and thereby maintain constant XA. If the activity az in any segment of the bed has decayed more than average, then the quantity of catalyst δW required to achieve the conversion element δXA must increase to compensate. Although the proportion of the mass of the catalyst bed required for conversion element δXA increases, the proportion of the residual activity that is required remains constant. In other words, the plot of conversion against fractional residual activity is constant, even though the axial conversion profile changes with degree of deactivation. Hence, the axial fractional residual activity profile remains under the influence of a constant composition profile during catalyst deactivation. Analytical solutions for the gradients were found on this basis for cases with m = 1. The same solutions can also be found by evaluation of the mean integral activity Æaæ in terms of the initial mean bed value of the decay kinetic function Ψ(C) followed by combination with the appropriate reaction kinetics. A general mathematical solution, based on active site density and valid for any value of m, is presented in Appendix 2. The general solution is always perfectly linear for m = 1 but it is not linear for m 6¼ 1. The solutions for any combination of reaction kinetics Λ(CA)

and decay kinetics Ψ(CA) with m = 1 can be derived from eq A8. The solutions for some selected cases are presented in Table 3. The agreement between calculated and simulated gradients was generally good (see Figure 5 in section 6), although in some cases, due to the stiffness of the equations, it was necessary to reduce the time intervals in the simulation by a factor of 10 in order to improve the accuracy. A linear plot is also expected for slow, irreversible poisoning because the axial deactivation profile is effectively uniform and therefore equivalent to concentration-independent decay, provided of course that the feed poison concentration CP0 is constant. The linear plot can easily be derived for fast, irreversible poisoning, which is once again equivalent to a reduction in the quantity W of active catalyst caused by a poison front moving progressively down the bed (section 4). See Appendix 1.b for details. It follows from the above analysis that the simple CCIST method can be used reliably to determine the decay order m, but it is still potentially misleading because plots, which appear to be linear, can be used to fit invalid concentration-independent decay kinetics and thereby derive invalid values for deactivation rate constants. In order to overcome this problem it is necessary to carry out tests in which composition is varied in a systematic manner, as described in the next section.

6. CCIST METHOD WITH VARIATION OF INITIAL CONVERSION The CCIST method can be used for diagnostic and comparative purposes by running multiple tests in parallel on each catalyst, covering a wide range of constant conversion values by selecting different initial values of space-time. Each test will produce a linear plot for m = 1 or a near linear plot for m 6¼ 1, provided that the correct equation for the actual decay order has been selected from eqs 6-8b. For first-order decay, the general nature of the concentration dependence will then be indicated by the change in the gradient (the effective kd) in the linear plots of eq 6 with change in either conversion or the initial value of space-time. This is illustrated in Figure 5a and 5b for a range of different composition-dependent decay kinetics. This analytical approach should provide a more straightforward means of model selection and discrimination than the purely statistical (trial and error) approach. It should certainly be possible to choose between the main types of decay kinetics. Provided that the decay order m has been correctly determined, the general nature of the composition dependence can then be inferred from the plot of ‘effective’ kd against conversion X. For 3153

dx.doi.org/10.1021/ie101481y |Ind. Eng. Chem. Res. 2011, 50, 3145–3158

Industrial & Engineering Chemistry Research

ARTICLE

Figure 5. Variation in effective kd from a plot of eq 6 for different test conditions. (a) Plot against conversion. (b) Plot against initial flow rate or 1/τi. ‘Effective’ kd from plots of eq 6. Actual kd values selected so that all plots lie on the same scale. c = calculated from analytical solution. s = estimated from numerical simulation. Shape of case 14 in Figure 5a for additive ParþSer decay will lie somewhere between concentration-independent and either parallel or series, depending on parameter h = kdB/kdA, and the effective kd will not approach zero at either X = 0 or 1. kd for case 9 was ‘wrongly’ estimated using eq 6 in the early near linear region, Ff/Fi = 0.75. As t f 0, this plot will approach that of case 4, i.e., ‘effective’ kd = kdi.

first-order deactivation (m = 1), if the ‘effective’ kd values are plotted against trial functions of X (listed in Table 3), then the correct function will give a linear plot. In practice, selection of the most nearly linear plot might be the best outcome, for example, if the concentration dependence is not a simple function or if the simplified power-law decay kinetics is not valid over the full range of composition. If the decay kinetics is wrongly determined to be first order, by following deactivation for too short a time or by drawing straight lines though curved or scattered plots, then the kd vs X plot should still indicate the nature of the composition dependence but the overall decay kinetics will be incorrect. This is shown by the example of case 9 (second-order parallel decay), incorrectly plotted in Figure 5a. The values of ‘effective’ kd were determined

at a low degree of decay (Ff/Fi = 0.75) such that a linear fit appeared plausible. The shape of the kd vs X plot indicates parallel decay, and a plot (not shown) of ‘effective’ kd against the function (CA0 - CA)/ln(CA/CA0), corresponding to first-order parallel decay from Table 3, also appears linear. Hence, this procedure might still be useful for determining the composition dependence even for second-order deactivation, but for correct overall decay kinetics, it is essential that the decay order m is correctly determined.

5. DISCUSSION It is evident from this study, and not really surprising, that concentration-independent decay kinetics cannot reliably be 3154

dx.doi.org/10.1021/ie101481y |Ind. Eng. Chem. Res. 2011, 50, 3145–3158

Industrial & Engineering Chemistry Research inferred from integral activity-time data based on any single experiment in a fixed catalyst bed, either with decreasing conversion (DCCST) or with constant conversion (CCIST). However, the DCCST and CCIST methods can both be used to investigate concentration-dependent decay kinetics by means of a systematic experimental strategy, varying space-time so as to cover a range of initial conversion values. Accurate flow control over a wide range is nowadays possible, and so, the CCIST method should be quite practicable. The CCIST method has a number of advantages, summarized in section 6. The decay order m and concentration dependence of deactivation can be resolved by separate plots, whereas they are likely to be confounded with each other and also with the reaction kinetics if different models are fitted statistically by trial and error. Statistical refinement can still be carried out as the final stage of the analytical approach if this appears to be justified. Ψ(C) might be an indirect function of conversion, depending, for instance, on the removal of some protective component such as O2 or H2, but even so it should be possible to rearrange the composition dependence to be a function of conversion. Any observed deviation from ideal behavior is potentially a useful indication of more complex, possibly even nonseparable, decay kinetics. However, even an approximate decay model is a useful starting point, which can provide direction for a more detailed systematic investigation, such as the effect of the concentration of individual feed or product components. The systematic procedure for diagnosis of the nature of decay kinetics using the DCCST method, as described in section 4, is based on the separate analysis of several conversion-time plots to generate a range of initial ‘effective’ kd values. Other methods of analyzing data from multiple tests over a range of values of space-time have been described previously.7,15,16 Corella et al. showed how point activity-time values in a fixed bed can be derived by differential analysis of multiple tests carried out over a range of different values of space-time.15 In the ‘parallel difference’ method, the changing performance over time in each integral sector of a full catalyst bed is determined by comparison with simultaneous tests in partially filled beds.16 The preferred decay model can be determined by statistical fitting of trial models to the data generated by any decreasing conversion method, but there is no visual indication of the decay order, and its value might be confounded to some extent with other model parameters. Provided that catalyst decay is followed over a wide range, the curvature of the conversion-time plot from a single test can provide some indication of the decay kinetics, as shown in Figure 1a. However, it is difficult to discriminate unambiguously between concentration-independent decay and various concentration-dependent decay models using a single test, even if the deactivation is followed over a wide range. However, analysis of conversion-time curves may be suitable in some cases for rapidly decaying catalysts. Analysis of the methanol breakthrough curves from single tests has been used to model the deactivation of zeolite catalysts for methanol conversion to hydrocarbons, assuming first-order parallel decay kinetics.20 An equivalent solution for the conversion-time dependence for first-order parallel decay has also been described previously, using arguments based on mean integral activity.21 This solution matches the simulated conversion-time plot for case 4 in Table 1. For true concentration-independent decay the observed kd, determined by either the DCCST or the CCIST method, should be invariant with conversion or indeed with any change in feed

ARTICLE

composition. Decay kinetics may still appear falsely to be independent of composition in certain cases. If the deactivation is caused by the separate, additive effects of reactants and products, then additional tests can be carried out with inert diluent to determine if the ‘effective’ kd varies with the total concentration of reactants and products. Reliable diagnosis of poisoning obviously requires feedstock of constant composition (preferably a single batch supply) and avoidance of any contamination from equipment. Deactivation by a slow-reacting poison resembles independent decay due to the low poison conversion across the bed. Some variation of the ‘effective’ kd might be detected by periods of operation at near total reactant conversion in order to force some variation in the axial poison concentration profile. If the poison is converted at a similar rate to the reactants, then the decay kinetics will resemble parallel decay. The influence of a poison can also be inferred by other means, such as a change in ‘effective’ kd when the feed is passed through absorbent guard beds and/or postmortem analysis of catalyst samples. It was noted in the Introduction that the combination of empirical coking kinetics with empirical activity-coke expressions can usually be re-expressed in terms of activity power-law decay kinetics. Hence, the DCCST and CCIST methods can be used to study decay kinetics caused by coking without measuring the coke content. In principle, if required, the activity-coke dependence can additionally be estimated by comparing the observed values for mean coke content for catalyst discharged after CCIST tests, carried out over a range of conversion values, with values estimated by numerical simulation, combining trial activity-coke relationships with the observed empirical powerlaw decay kinetics. One complication that has not been considered so far is the shift from one dominant decay mechanism to another under different reaction conditions or at different stages of deactivation, such as the “fast” and “slow” coke” suggested by Shum et al. for deactivation of Pt-Re/Al2O3 reforming catalysts using the decreasing conversion method.11 However, provided that deactivation is studied over a sufficient range of conditions and a sufficient degree of deactivation, then such complications should become evident. Deactivated samples should in general be subjected to suitable physical and chemical characterization in order to confirm the nature and consistency of the decay mechanism. For practical purposes, it is sometimes preferred to disregard an initial period of rapid deactivation (‘lining out’) if the main interest is in the subsequent long, slow deactivation over the active life of the catalyst. However, this approach will be invalid for the DCCST and CCIST methods, as described above, if the rapid deactivation process is composition dependent, because an activity gradient will already have been generated in the catalyst bed by the start of the slow deactivation process. Any such ‘lining out’ would need to be carried out as a batch pretreatment under gradientless conditions. Alternatively, a short period of rapid deactivation followed by a long period of slow deactivation might correspond to a single decay process with high decay order m, such as metal particle sintering, in which case the DCCST and CCIST methods are valid over the entire range of deactivation. It is beyond the scope of this work to consider in detail the effect of temperature on catalyst deactivation. Obviously, separate experiments can be conducted at different temperatures and then, if possible, the kd values can be fitted to an Arrhenius expression in the deactivation model, defining Ed as the ‘activation energy for decay’. Hence, the least number of experiments 3155

dx.doi.org/10.1021/ie101481y |Ind. Eng. Chem. Res. 2011, 50, 3145–3158

Industrial & Engineering Chemistry Research

ARTICLE

required for a simple deactivation model that incorporates both composition and temperature dependence is 5, i.e., 3 concentration values, one of them at 3 temperatures. Obviously, more experiments will increase confidence in the conclusions. If an additional parallel experiment is carried out with constant conversion and space-time, increasing temperature, then this will confirm the general validity of the deactivation model, regardless of the decay path. A further challenge is to characterize catalyst samples, which have deactivated to the same extent by different decay paths, and show that they have reached the same average physicochemical state. It is worth noting at this point that an activity-time relationship for a test at constant conversion, increasing temperature has been derived by Ho, based on the assumption that composition can be treated as constant.22 However, this relationship is valid for composition-dependent decay only if the decay order m = 1. This requirement has been confirmed by some further simulation. For parallel decay, defining ÆCAæi as the mean bed value at t = 0 and q = Ed/Er, the equation can be adapted to the form aðtÞ ¼ ½1 - kdi hCA ii qt - 1=q

ð12Þ

For m 6¼ 1, the assumption that composition can be treated as constant appears to be flawed in that the local decay rate is dependent on the product of CAz and azm, and so the mean bed concentration ÆCAæ does not apply. One potential difficulty of the DCCST and CCIST methods is the need to study deactivation to such a high degree, at least in some tests. This may take a long time. However, some test conditions will naturally lead to more rapid deactivation than others. For instance, with the CCIST method, the decay order for parallel decay will become evident most rapidly in a test at low conversion and for series decay in a test in which some reaction product has been added to the feed. Likewise, more rapid deactivation will be expected in experiments conducted at increased temperature. Hence, the decay order can be determined from tests with rapid deactivation, and other tests with slower deactivation need only be followed for a short time. These ‘accelerated decay’ conditions16 would thereby become part of an overall test strategy in which equipment is used with maximum test productivity. Intrinsic selectivity decay is not reliably observed in the decreasing conversion method because selectivity may change with conversion due to consecutive reactions of the products. This is not a problem when operating at constant conversion (and constant temperature), and so the direct observation of intrinsic selectivity decay is clearly an advantage for the CCIST method. Catalyst deactivation models usually do not need to be perfect, just fit for their purpose, e.g., assessing the impact of catalyst composition during catalyst development. The models derived using these methods for tests on small particles can be combined, if necessary, with equations for intra- and interparticle diffusion for a deactivation model at pellet scale and then used for optimization of reactor design and overall process performance. The ideas presented in this paper for application of the decreasing conversion and constant conversion methods are mainly theoretical, based on numerical simulation and analytical solutions for certain well-defined cases of decay kinetics. No suitable range of experimental data has been found in the literature to test these ideas in practice, but hopefully this deficiency will in due course be corrected by some dedicated research programs.

6. CONCLUSIONS It is easy to generate incorrect decay kinetics from limited tests in integral reactors. However, reliable decay kinetics can be determined by either the DCCST or the CCIST method, provided that tests are carried out over a suitable range of conversion. There are few examples of the CCIST method in the literature to date, but it has several potential advantages over the DCCST method, namely, (i) clear indication of the decay order, (ii) applicability to reversible reactions, (iii) exact analytical solutions for concentration-dependent first-order decay, (iv) indication of the range of validity of the decay kinetic model, (v) acceleration of the overall test procedure, and (vi) reliable indication of intrinsic selectivity decay. ’ APPENDIX 1. FAST POISONING 1.a. DCCST Method. Fast, irreversible poisoning is essentially equivalent to a reduction in the quantity of active catalyst caused by a poison front, moving progressively down the bed, with activity a = 0 in the deactivated zone and a = 1 elsewhere. The initial values for the ‘effective’ kd in eq 3 can be calculated as follows   kr W NAv - kr CP0 f ðXÞ ¼ ln ln ð1=ð1 - XÞÞ ¼ ln t F SWi df ðXÞ NAv 1 ¼ - CP0 : dt SWi ðW=F - CP0 ðNAv =SWi ÞtÞ At t ¼ 0,

df ðXÞ F NAv ¼ - QP where QP ¼ CP0 dt W SWi ð10Þ

1.b. CCIST Method. In order to maintain constant conversion, FA0 must be adjusted correspondingly so that the quantity Wactive/FA0 is constant and equal to W/FA0i.

Wactive þ dWactive Wactive - FQP dt W ¼ ¼ where QP FA0 þ dFA0 FA0i FA0 þ dFA0 Z CP0 NAv WCA0 dFA0 ¼ and FA0 ¼ FCA0 Hence SWi FA0i FA0 Z FA0 ¼ - QP dt, and so τi ln ¼ - QP t FA0i This can be rearranged to give ln τ ¼ ðQ P =τi Þt þ ln τi

ð11Þ

This equation reproduces the simulated plot for fast, irreversible poisoning.

’ APPENDIX 2. CCIST METHOD: GENERAL SOLUTION Consider a fixed catalyst bed, mass W, and site density SW per unit catalyst mass. Hence, -((da)/(dt)) = kdam can be reexpressed as -((d(SW/SWi))/(dt)) = kd(SW/SWi)m or (dSW)/ (dt) = -kdSSWm, where kdS = kd/SWim-1 and kdS = kd for m = 1. Let Sz = the variation in the number of sites with the fractional bed length z, where Szi = πr2FSWi. Note that Sz is directly proportional to the true site density. Hence, the total number of R sites in the bed at any time t = ST = 10Sz dz. The total fractional conversion is XA at the bed exit (z = 1) for all values of t and X(z,t) for any other point in the bed. XA is 3156

dx.doi.org/10.1021/ie101481y |Ind. Eng. Chem. Res. 2011, 50, 3145–3158

Industrial & Engineering Chemistry Research

ARTICLE

maintained constant by reducing FA0 to compensate for deactivation. If the decay rate varies with composition in any way that can be written as a function of CA and kdz is the deactivation rate constant per unit bed fraction (kdz = kdW(πr2LF)1-m), then at any point in the bed and at any time dSz ¼ - kdz ΨðCA ÞSz m ðA1Þ dt The reaction rate at any time is given by FA0 dX ¼ ks ΛðCA ÞSz dz Rearranging eq A2 at any fixed time t Z Z XA dX ks 1 ks ¼ Sz dz ¼ ST ΛðC Þ F F A A0 0 A0 0

ðA2Þ

ðA3Þ

Since CA = CA0(1 - X), the integral on the left side of eq A3 is solely a function of X. Since XA is constant for all t, then the integral from X = 0 to XA is constant for all t and is hereby designated I. Note also that this holds true even if the fluid density changes during the reaction (ΔVR 6¼ 0), although the integral needs to incorporate this effect. Hence ðA4Þ FA0 ¼ k1 ST where k1 ¼ ks =I Note that FA0 depends on the total number of sites and not on their distribution. If necessary, the active site can be defined as a multicenter cluster of surface atoms. Also, at any time t Z Z 1 dFA0 dST ∂ 1 ∂Sz ¼ k1 ¼ k1 dz ðA5Þ Sz dz ¼ k1 ∂t 0 dt dt 0 ∂t This follows from application of the Leibniz integral rule to the specific case when the integral limits are constant values. Substituting eq A1 into eq A5 Z 1 dFA0 ¼ k1 - kdz ΨðCA ÞSmz dz ðA6Þ dt 0 Rearranging eq A2 dX ks FA0 dX ¼ ΛðCA ÞSz or Sz ¼ dz FA0 ks ΛðCA Þ dz Substituting for Sz from eq A7 into eq A6  m Z 1 m dFA0 FA0 ∂X ¼ k1 - kdz ΨðCA Þ m dz m ∂z dt k ΛðC Þ A 0 s   Z 1 ΨðCA Þ ∂X m m dz ¼ k2 FA0 m ∂z 0 ΛðCA Þ

ðA7Þ

ðA8Þ

m-1 where k2 = (-k1kdz)/km I), substituting for k1 s = -kdz/(ks from eq A4. Note also from eqs A5 and A4 that R 1 ∂Sz dz m Z 1 dFA0 FA0 ∂Sz m 0 ∂t ¼ k1 m dz ¼ k1 FA0 ðk1 ST Þm dt FA0 0 ∂t R 1 ∂Sz dz dFA0 1 - m m 0 ∂t ¼ k1 and so FA0 ðA9Þ m dt ST

From eq A9, (dFA0)/(dt) can have the form (dFA0)/(dt) = m k3Fm A0 only if (((∂Sz)/(∂t))/ST ) is not a function R of t. This is identical to the condition from eq A8 that 10((Ψ(CA))/

(Λ(CA)m))((∂X)/(∂z))m dz is not a function of t. RFor first-order decay (m = 1) the above integral becomes I2 = 10((Ψ(CA))/ R1 (Λ(CA))((∂X)/(∂z))dz = 0g(X)((∂X)/(∂z))dz. LetRG(X) be the antiderivative of g(X),R i.e., ((dG)/(dX)) = g(X). I2 = 10((dG)/ (dX))((∂X)/(∂z))dz = X0 A((dG)/(dX))dX = constant. Hence, the condition is met for all cases for which m = 1, and for which Ψ, Λ can be written as functions of conversion of A. However, for m 6¼ 1, it is not generally true that the integral I2 is a constant. Hence, the deviation from linearity is explained for the numerical simulation of cases for which m 6¼ 1. Setting m = 1 in eq A8 and integrating, it follows that ln(FA0/FA0i) = -k3t or ln τ = k3t þ ln τi. This has the same form as eq 6 for concentration-independent decay in section 1. The proof is general for any reaction kinetics, even with a change in the 6 0)).ThevalueoftheconstantkdR=k3 dependson fluiddensity((ΔVR ¼ the functions Ψ and Λ in the integral k2 X0 A((Ψ(CA))/ (Λ(CA)))dX. This can easily be evaluated, and solutions for some cases with constant fluid density are shown in Table 3.

’ AUTHOR INFORMATION Corresponding Author

*Address: 15 Portman Rise, Guisborough, TS14 7LW, U.K. Tel: 0044 1287 638265. E-mail: [email protected].

’ ACKNOWLEDGMENT The University of Glasgow granted an Honorary Senior Research Fellowship to the author. The author is grateful to Paul A. Kirwan for deriving a general mathematical solution which has been adapted for Appendix 2. ’ KEY TO SYMBOLS a = activity (reaction rate at any point in bed)/(rate under equivalent conditions at t = 0) Æaæ = mean activity value across catalyst bed b, b0 , c = algebraic constants CA = concentration of reactant A (mol L-1) ÆCAæ = mean value across catalyst bed CP = concentration of poison P (mol L-1) h = kdB/kdA = ratio of deactivation rate constants for product B and reactant A kr = forward reaction rate constant per unit mass catalyst (units depend on rate law) kr0 = reverse reaction rate constant per unit mass catalyst (units depend on rate law) ks = reaction rate constant per active site (units depend on rate law) kd = catalyst deactivation rate constant (units depend on rate law) kd0 = site-loss rate constant (units depend on rate law) kdw = catalyst deactivation rate constant based on catalyst mass Er, Ed = activation energy for reaction and decay, respectively F = total molar flow rate (mol h-1) FA = molar flow rate of component A L = bed length NAv = Avogadro number m = decay order on activity n = decay order on concentration p = site order in reaction rate equation QP = CP0NAv/NS r = radius of catalyst bed S = number of active sites 3157

dx.doi.org/10.1021/ie101481y |Ind. Eng. Chem. Res. 2011, 50, 3145–3158

Industrial & Engineering Chemistry Research ST = total number of active sites in catalyst bed SW = number of active sites per unit mass of catalyst W = mass of catalyst (g) X = fractional conversion ΔVR = volume change per mole of reaction F = packing density of catalyst (mass per unit volume) Ψ(CA), Λ(CA) = undefined functions of CA τ = (weight) space-time = WCA0/FA0 = W/F

ARTICLE

(20) Janssens, T. V. W. A New Approach to the Modeling of Deactivation in the Conversion of Methanol on Zeolite Catalysts. J. Catal. 2009, 264, 130. (21) Ostrovskii, N. M. Problems in the Study of Catalyst Decay Kinetics. Kinet. Catal. (Engl.) 2005, 46, 693. (22) Ho, T. C. Some Aspects of the Constant-Conversion Policy Dealing with Catalyst Deactivation. J. Catal. 1984, 86, 48.

Subscripts

A, B = reactant A and product B 0 = inlet value to reactor z = position along the reactor bed i, f = initial and final values over a period of time

’ REFERENCES (1) Levenspiel, O. Experimental Search for a Simple Rate Equation to Describe Deactivating Porous Catalyst Particles. J. Catal. 1972, 25, 265. (2) Levenspiel, O. Chemical Reaction Engineering, 2nd ed.; Wiley: New York, 1972; Chapter 15. (3) Wolf, E. E.; Petersen, E. E. On the Kinetics of Self-Poisoning Catalytic Reactions. J. Catal. 1977, 47, 28. (4) Ostrovskii, N. M. General Equation for Linear Mechanisms of Catalyst Deactivation. Chem. Eng. J. 2006, 120, 73. (5) Sehested, J. Sintering of Nickel Steam-Reforming Catalysts. J. Catal. 2003, 217, 417. (6) Reiff, E. K., Jr.; Kittrell, J. R. Use of Active Site Balances for Catalyst Deactivation Models. Ind. Eng. Chem. Fund. 1980, 19, 126. (7) Birtill, J. J. Measurement and Modeling of the Kinetics of Catalyst Decay in Fixed Beds: the Eurokin Survey. Ind. Eng. Chem. Res. 2007, 46, 2392. (8) Corella, J.; Asua, J. M.; Bilbao, J. Kinetics of Catalyst Deactivation. Can. J. Chem. Eng. 1981, 59, 647. (9) Wojciechowski, B. W. The Kinetic Foundations and the Practical Application of the Time on Stream Theory of Catalyst Decay. Catal. Rev. 1974, 9, 79. (10) Sapre, A. V. Catalyst Decay kinetics from Variable SpaceVelocity Experiments. Chem. Eng. Sci. 1997, 52, 4615. (11) Shum, V. K.; Sachtler, W. M. H.; Butt, J. B. A Deactivation Correlation for Platinum-Rhenium/ Alumina in n-Hexane Conversion. Ind. Eng. Chem. Res. 1987, 26, 1280. (12) Haynes, H. W., Jr. A Parametric Study of Isothermal Fixed Bed Catalytic Reactors Undergoing Concentration-Dependent Deactivation. Chem. Eng. J. 1990, 43, 1. (13) Santos, A.; Garcia-Ochoa, F.; Romero, A. Determination of Deactivation Kinetic Parameters. II. Data from Integral Reactors. React. Kinet. Catal. Lett. 1989, 40, 163. (14) Froment, G. F.; Bischoff, K. B. Kinetic Data and Product Distributions from Fixed Bed Catalytic Reactors Subject to Catalyst Fouling. Chem. Eng. Sci. 1962, 17, 105. (15) Corella, J.; Asua, J. M.; Bilbao, J. Kinetics of the Deactivation of a 10% Copper-0.5% Chromia/Asbestos Catalyst for Benzyl Alcohol Dehydrogenation. Chem. Eng. Sci. 1980, 35, 1447. (16) Birtill, J. J. But will it last until the shut-down? Deciphering Catalyst Decay. Catal. Today. 2003, 81, 531. (17) Birtill, J. J. Some Pitfalls in Modeling Catalyst Decay using Integral Reaction Data. Poster 33, extended abstracts, 10th International Symposium on Catalyst Deactivation, Berlin, Dechema, 2006. (18) Birtill, J. J. Unpublished lectures (a) Some Pitfalls in Modeling Catalyst Decay using Integral Reaction Data. SURCAT 2007, Glasgow, (b) Progress in Catalyst Deactivation. Eurokin 10th anniversary symposium, IFP Lyon, 2008. (19) Wentrcek, P. W.; McCarty, J. G.; Ablow, C. M.; Wise, H. Deactivation of Alumina-Supported Nickel and Ruthenium Catalysts by Sulfur Compounds. J. Catal. 1980, 61, 232. 3158

dx.doi.org/10.1021/ie101481y |Ind. Eng. Chem. Res. 2011, 50, 3145–3158