FCCU transition-probability model

Mathematical and Computer Modelling 45 (2007) 241–251 www.elsevier.com/locate/mcm

FCCU transition-probability model Robbie J. Dixon, Maki Matsuka, Roger D. Braddock ∗ , Josh M. Whitcombe, Igor E. Agranovski Faculty of Environmental Sciences, Griffith University, Nathan, Qld 4111, Australia Received 9 September 2005; accepted 3 March 2006

Abstract The adequacy of the use of transition-probability matrices for modelling fluidised catalyst cracker unit emissions was investigated. A number of different-sized matrices that modelled the processes of attrition and agglomeration were used, and it was found that an 8 × 8 sized matrix provided the best results. The processes of attrition within the matrix were studied, indicating an oscillatory attrition curve, and may suggest a preferred attrition size range. Studies on the effects that the agglomeration parameters had on the model indicated that, as time and the size of the matrix increased, agglomeration became more important. The results of the modelling were compared with laboratory experiments, and indicated very good agreement between the model outputs and the observed emissions. c 2006 Elsevier Ltd. All rights reserved.

Keywords: Attrition; Agglomeration; Probability matrix

1. Introduction Fluidised Catalytic Cracking Units (FCCU) are used in refinery processes to convert long-chained hydrocarbons (crude oil) into more valuable short-chained molecules (petrol or gasoline) [1]. An FCCU typically consists of three parts: a rising main where the catalyst enhances the shortening (cracking) of the hydrocarbon chains, a reactor to separate the product from the catalyst, and a regenerator to recharge the spent catalyst. The regenerator contains a fluidised bed of catalyst which is coated in carbon. A controlled amount of air is supplied to the fluidised bed and the coke is burnt off the catalyst, although complete combustion is not achieved. The heat released from this reaction heats the catalyst, which in turn provides the necessary energy for the cracking process in the rising main. Flue gas exits the regenerator through a series of internal cyclones located at the top of the system [2]. Regenerated catalyst passes from the regenerator, back into the rising main, whilst smaller uncaptured particles are vented to the atmosphere [1,3]. The FCC catalyst is an essential part of the refining process and must be continually added to the system to overcome the loss of particulate material due to particle fracture and attrition. There is a growing awareness of pollution problems and emission levels from industrial processes. The increasing strictness of legislation around the world has led to a situation where oil refineries must continually reduce their ∗ Corresponding author.

E-mail address: [email protected] (R.D. Braddock). c 2006 Elsevier Ltd. All rights reserved. 0895-7177/$ - see front matter doi:10.1016/j.mcm.2006.03.019

242

R.J. Dixon et al. / Mathematical and Computer Modelling 45 (2007) 241–251

emissions [4]. In addition, the catalyst is expensive, and any reduction in emissions will improve the economics of the refinery. In circulating through an FCCU, the particles of catalyst may undergo fracture (or shattering), whereby the particle splits into two or more parts, each of a significant size [5]. The catalyst particle may also undergo attrition, whereby small pieces, or “fines”, are removed or abraded from the particle [5]. Mechanical forces in various forms of collision are usually assumed to be the cause of the change in the particle-size distribution over time. Recent research suggests that thermal effects can lead to fracture of the catalyst particles [6,7]. Most of the research on attrition and fracture has been carried out in specific laboratory tests, and the results or models cannot predict the nature of attrition (or fracture) inside an industrial fluidised bed or regenerator, where the stress mechanisms are exceedingly complex [8]. There is still the scaling problem of trying to apply the results of laboratory-based experiments to full-scale refinery operations. Process models have been developed for various forms of attrition, i.e. bubble, jet and cyclone attrition, and Werther and Reppenhagen [8] have derived a differential-equation model to describe the total loss rate of particles from an FCCU. Whitcombe [9] attempted to apply the Werther and Reppenhagen model to the observations of emissions from a scale experimental model of an FCCU. The fit of the differential-equation model to the laboratory-scale FCCU was disappointing. Industrial observation of particle dynamics is difficult due to the FCCU’s operating temperatures of 600–700 ◦ C. There are other less prescriptive modelling methods available, such as the use of transition-probability matrices which have found extensive use in Ecology and the dynamics of populations [10,11]. The aim of this paper is to apply a transition-probability model to the refinery problem. This model will be calibrated using the observations from the laboratory plant of Whitcombe [9]. The resulting transition-probability values will then be interpreted and related to the dominant fracture and attrition processes occurring in the laboratory plant. 2. Transition-probability model Consider a collection of particles such as the catalyst in the laboratory plant. This collection is categorised by size into eight size classes, labelled from size class 1, the largest particles, to size class 8, which is the smallest. The number of particles in each size class is denoted by P (t) = Pi (t) for i = 1, 2, . . . , 8 at time t. Time is also discrete ˜ and t = 0, 1, 2, . . ., while P (0) is the initial distribution of particles. The fracture and attrition processes within the ˜ laboratory plant are partly discrete (i.e. collision) and partly continuous (i.e. abrasion). It is convenient to use a discrete representation of time in terms of minutes so that the properties of the transition-probability model can be utilised. The vector P (t) then describes the state of the system at discrete time t. Note that the particle-size classes are discrete, ˜ P i (t) is the number of particles in a continuous range of physical “diameters” represented by the class and that each ˜ the ranges of “diameter” associated with each size class will be discussed later. i. The choice of Now consider the transitions of the particles between the size classes in one discrete time step from t to t + 1. (See Fig. 1.) A particle in size class i at time t may remain in size class i at time t + 1, with a probability denoted by ai,i . The particle in size class i at time t may also undergo fracture or attrition, and generate a smaller particle of size class j > i at time t + 1. The probability of this transition is denoted by a j,i . These events are not mutually exclusive, and attrition or abrasion may produce a small particle and a large particle at the same time. Particle fracture may produce two approximately equally sized particles. Agglomeration is also handled by the model, in that a1,i is the probability of a large particle of size class 1 being formed from other size classes. The smallest particle sizes require special consideration, as these states are needed to model the emission of particles by the cyclones (see Fig. 2). Size class 8 is the smallest, and does not suffer any further attrition. However, the catalyst particles do pass through the cyclones, may be removed from the FCCU, and emitted to the atmosphere. The return efficiency of the cyclone, given by the parameter a8,8 (expressed as a ratio) is particle-size dependent, can be measured and also estimated theoretically [12]. The situation for size class 7 depends on the typical particle diameter for that class. These particles may be large enough to undergo total capture by the cyclone which thus has a capture efficiency of 100%. In this case, this size class is not emitted, and this size class can be handled as in Fig. 1. Where the particles in size class 7 are small enough, these particles may escape capture in the cyclones, and may be emitted. The fraction of particles of size class 7 emitted from the cyclone is assumed to be 1 − a7,7 .

R.J. Dixon et al. / Mathematical and Computer Modelling 45 (2007) 241–251

243

Fig. 1. Transition probabilities for particles in size class i at time t to other size classes at time t + 1.

Fig. 2. State diagram for size classes 7 and 8, without agglomeration.

Then, for each size class i, the number of particles in this size class at time t + 1 is given by Pi (t + 1) =

8 X

ai, j P j (t),

(i = 1, 2, . . . , 8),

(1)

j=1

where ai, j P j (t) is the expected number of particles transferring from size class j to size class i in the time step. This may be written in vector form as P (t + 1) = A P (t), (2) ˜ ˜ where A = {ai, j } is the probability transition matrix. Note that the diagonal elements of A, i.e. ai,i , i = 1, 2, . . . , 8, represent particles which stay in the same size class from time t to t + 1. The elements ai, j , with j > i, represent agglomeration effects, and the elements ai, j , with j < i, represent fracture and attrition effects. The total number of particles emitted during each transition, E(t), is given by E(t) = E 7 (t) + E 8 (t),

(3)

where E 7 (t) = (1 − a77 )P7 (t), E 8 (t) = (1 − a88 )P8 (t), where, in this case, size classes 7 and 8 only are emitted. Note that the size-class ranges selected to represent the particle-size distribution will be reflected in the emissions. Cyclones used in oil refineries are generally very efficient at retaining particles above about 25 µm. The number of size classes below 25 µm will affect the structure of E(t). Note that this model will not conserve mass. The size classes i are based on the particle diameter, and each class incorporates a range of diameters of the particles. Thus a range of diameters from 10 µm to 20 µm encompasses a corresponding eightfold range of particles by volume (or mass). Obviously, a particle of 20 µm diameter can fracture

244

R.J. Dixon et al. / Mathematical and Computer Modelling 45 (2007) 241–251

Fig. 3. Graph showing 0.3 µm raw and ‘smoothed’ data.

into two particles which remain in this size class. The accounting by the number of particles does not conserve mass because of this effect. Also note that the number of size classes, say eight, can be selected to suit the data available. 3. Emissions data To test the effectiveness of the model, operational data obtained from an experimental fluidised bed were used. The emission results were taken from a scale model FCCU used in experiments performed by Whitcombe [9]. The laboratory scale equipment was filled to the same level as the industrial FCCU, so that the cyclones and diplegs operated in the same way. Sampling points were connected to a real-time particle counter, providing total and sizespecific results in terms of particle concentrations from each sampling point. An Autocounter 300A (Malvern, UK) was used with sampling nozzles being designed to allow USEPA Method 5 isokinetic sampling conditions to be maintained. Catalyst used in the laboratory was sourced from the industrial FCCU and filled to the standard operating level used in the industrial FCCU, ensuring standard results. In total, 250 data points were measured at 0.03 h intervals, representing about 7.5 h of data. The catalyst particles emitted from the FCCU cyclones were grouped into eight discrete particle diameters (i.e. 0.3, 0.8, 2.0, 4.0, 8.0, 15.0, 20.0 and 25 µm), and the raw data output for 0.3 µm diameter particles is shown in Fig. 3. Note that a number of particles within each diameter size were emitted from the FCCU, and therefore must be represented in the smallest size class of the model, as discussed in Section 2. The data was highly variable, as shown in Fig. 3, with a lot of higher frequency noise. The data was mathematically filtered three times using the following formula F=

1 (Y−1 + 2Y0 + Y1 ), 4

where Y0 is the data point at time t Y1 is the data point at time t + 1 Y2 is the data point at time t − 1 F is the ‘smoothed’ data point. This filtering yielded the ‘smoothed’ graph which is also shown in Fig. 3. The raw data for the other sizes of particles were also filtered in the same manner. The particle counts in these ranges are approximately 10 times less than for the 0.3 µm case, and the number of particles >20 µm emitted after one hour is almost negligible. A full description of the experimental techniques used and the results obtained is given in [9]. 4. Size classes The particle-size distribution (PSD) of the catalyst in the laboratory equipment was divided up into different size classes, to represent the catalyst circulating in the equipment. For the example used in Section 2, an 8 size-class

R.J. Dixon et al. / Mathematical and Computer Modelling 45 (2007) 241–251

245

Fig. 4. Equilibrium catalyst (e-cat) particle-size distribution (PSD) used in the FCCU.

Table 1 List of size class systems (ranges in µm) Size class P1 P2 P3 P4 P5 P6 P7 P8 P9 P10

Size class system 4×4

6×6

8×8

10 × 10

100–300 30–100 5–30 0–5 – – – – – –

160–300 100–160 60–100 20–60 0.5–20 0–0.5 – – – –

160–300 130–160 100–130 70–100 40–70 10–40 1–10 0–1 – –

160–300 140–160 120–140 100–120 60–100 40–60 20–40 5–20 1–5 0–1

system was considered. The PSD of the fresh (e-cat) catalyst [9], as determined using a Mastersizer S, laser particlesize analyser with a small-volume sample-handling unit attached (Malvern Instruments, UK), is shown in Fig. 4. The PSD can easily be divided into any number of size classes. For example, a 4 size-class system (i.e. i = 4) could consist of 0–5, 5–30, 30–100, 100–300 µm (diameter) size ranges. These ‘divisions’ may be altered, but it is important to include all the particles. Further, a class of small-sized particles is required to allow the effects of the cyclone emissions to be modelled. The model was implemented for 4, 6, 8 and 10 size-class systems, and each system was broken down into various size-class ranges as shown in Table 1. The size classes emitted by the cyclones were determined on the basis of the midpoint diameters of the size-class ranges, and a cyclone efficiency cut-off at 25 µm. These size-class ranges are shown in bold in Table 1. The particles were assumed to be spherical in shape and have a uniform density of 1400 kg/m3 . The number of particles within each size range in each size-class system was found from the initial amount of catalyst present in the FCCU at t = 0. The cyclones in the FCCU are not efficient at removing the particles in the final two (highlighted bold in Table 1) size classes in each size-class system. Therefore, the model will have two emission outputs for the 4, 6 and 8 size-class systems. In the 10 × 10 system, a larger number of ranges is available, which allows a ‘finer’ division of the emitted size classes, with a commensurate change to E(t) in Eq. (3). The raw, or observed, emissions were collected on a different set of particle sizes (i.e. 0.3, 0.8, 2.0, 4.0, 8.0, 15.0, 20.0 and 25 µm). All the observed emissions need to be matched with the relevant size classes in each size system in Table 1. For example, in the 4 × 4 system, there are four size classes, P1 , P2 , P3 and P4 . The particles in size classes P1 and P2 are assumed to be 100% collected by the cyclones, and therefore are not represented in the emissions. However, P3 and P4 are not 100% collected, and therefore must incorporate all observed emissions from the test rig.

246

R.J. Dixon et al. / Mathematical and Computer Modelling 45 (2007) 241–251

Fig. 5. Model schematic.

This is done by summing the emissions data within each size-class range P3 and P4 . Then, O3 (t) = observed emissions corresponding to size class P3 (see Table 1), at time t = sum of number of observed emitted particles in sizes 8 µm, 15 µm, 20 µm and 25 µm at time t. and O4 (t) = observed emissions corresponding to size class P4 (see Table 1) at time t = sum of number of observed emitted particles in sizes 0.3 µm, 0.8 µm, 2 µm and 4 µm at time t. This groups the observed emissions to fit to the size classes in the 4 × 4 transition matrix model. Similar definitions are readily obtained for the observed emissions for the 6 × 6, 8 × 8 and 10 × 10 models. 5. Precise fitting of parameter values and numerical experiments The fitting of the transition matrix model to the emissions data obtained from the test rig is achieved by leastsquares analysis. Let F=

n2 m X X

(E i (t j ) − Oi (t j ))2 ,

j=1 i=n 1

where m is the number of time steps, and i = n 1 , . . . , n 2 , and n 1 and n 2 delimit the size classes emitted by the model. In Eq. (3), n 1 = 7 and n 2 = 8, and similar definitions hold for the other 4 × 4, 6 × 6 and 10 × 10 models. The error function F = F(A) is to be minimised with respect to the coefficients, or parameters, in the transition-probability matrix A. This minimisation was performed using a genetic algorithm. The fitting of the model to observed emissions consists of three main parts: a genetic algorithm, the function to be optimised, and the transition-probability black-box model. The way in which the different sections are linked is shown in Fig. 5. The genetic algorithm generates ‘test’ probability parameter values for each element of A. These parameters are then fed into the transition-probability model, and the model is run over 250 time steps, which gives the modelled emissions. The modelled emissions and observed emissions are then fed into the objective function and the result stored in the GA. The GA then generates a new set of parameters based on the input from the objective function, and the process starts again. The number of generations controls the number of passes or iterations by the GA. 5.1. Parameter bounds Bounds are required for the parameters, or elements of A, in the model, for use in the genetic algorithm. These bounds are used to restrict the search region by the GA, and thus reduce the computation time. Some pilot experiments were run to obtain crude bounds on the parameter values. These bounds are to be used in more exact fitting of the model, and the following example relates to the 8 × 8 model.

R.J. Dixon et al. / Mathematical and Computer Modelling 45 (2007) 241–251

247

5.1.1. Probability of particles staying in the same size class For example, in Fig. 1, the probability that the particles in size class i will not undergo attrition is quite high, i.e. a greater than 95% chance. In the bounds, aii can range from 0.95 ∼ 1 for i = 1 to 8 (i.e. a 95% to 100% chance of not undergoing attrition). 5.1.2. Attrition probability bounds After some experimentation, it was found that the bounds for the attrition terms should lie in the range 0–0.1. These parameters correspond to the elements ai j , i > j. 5.1.3. Agglomeration probability bounds The probabilities of agglomeration occurring were found to be much lower than for attrition. From running the experimental model, the probabilities were shown to be in the range from 0 to 1 × 10−5 . These parameters correspond to the elements ai j , j > i. 5.1.4. Emission parameters The parameters a77 and a88 relate to the efficiency of the cyclone. Based on testing and the collection efficiency curves for the worst-case scenario in Fig. 2, the range varied from 0.1 to 0.5 and this range was kept constant for each size system. 5.2. Numerical experimentation The numerical experimentation consisted of modelling the emissions data using a 4, 6, 8 and 10 size-class system. In the first set of experiments, each system included both attrition and agglomeration terms. Once the transition matrices were developed, a plot of attrition through each size class was made. In the second set of experiments, the agglomeration terms were left out, and the model was run in order to determine the effect the agglomeration terms had on the model. 6. Results and discussion 6.1. Results for the 8 × 8 system The parameter-transition matrix obtained for the 8 × 8 system, with both attrition and agglomeration, is shown in Eq. (4). The agglomeration terms are in italics and the attrition terms are in bold.

(4)

The diagonal elements aii , i = 1, . . . , 6 are large and nearly unity, indicating the probability that these particles will remain in their size classes. Note that the final two highlighted diagonal numbers relate to the efficiency of the cyclone. Two were required because the particles in those size classes are not efficiently removed by the cyclone and are therefore involved in the emissions (i.e. size classes 7 and 8). The values of a7,7 and a8,8 show the decreasing ability of the cyclones to capture the smaller particles. The attrition terms in A (in bold) are up to 104 times greater than the agglomeration terms. The attrition terms are also at least an order of magnitude less than the diagonal terms. This tends to support literature in that the attrition process is more pronounced than agglomeration [9]. The modelled and observed emissions for size classes 7 and 8, are shown in Fig. 6. The model is predicting the emissions of these two size classes very well. The only feature missing from the modelled results is the higher

248

R.J. Dixon et al. / Mathematical and Computer Modelling 45 (2007) 241–251

Fig. 6. Size class 8 × 8, modelled and observed emissions.

Fig. 7. Size class 6 × 6, modelled and observed emissions.

frequency oscillations in the observed data. The source or reason for this feature of the observed data is not clear, but it may arise from fluctuations in the air supply, or from flow resonance in the laboratory rig, or the operation of the particle counter. 6.2. Results for the 4 × 4, 6 × 6 and 10 × 10 systems The results for the 6 × 6 system were similar to the results for the 8 × 8 system. The observed and modelled emissions for size classes 5 and 6 (see Table 1) are shown in Fig. 7. The emissions in size class 6, i.e. 0–0.5 µm, are well modelled, with good agreement between observed and modelled values. Note that the emissions from size class 5 are not modelled as well as for size class 7 in the 8 × 8 model. The transition-probability matrix, A6 , is given by

(5)

R.J. Dixon et al. / Mathematical and Computer Modelling 45 (2007) 241–251

249

Table 2 Results of the sum of least-square error values Size-class system

Error value

A4 A6 A8 10 × 10–3 output 10 × 10–2 output A∗8

1.33 × 1010 1.37 × 109 1.26 × 109 1.57 × 109 1.35 × 109 1.29 × 109

and displays a structure similar to that of A8 . The values of a5,5 and a6,6 again show that the cyclones are not very efficient at removing small particles from the exhaust streams. The particle-size classes for the 6 × 6 and 8 × 8 models are different (see Table 1), but the cyclone efficiencies are comparable. The results for the 10 × 10 model showed that the fits between the observed and modelled emissions were similar to the 8 × 8 model. In Table 1 there are three size classes, i.e. 8, 9 and 10, which are emitted through the cyclones. The emissions for size class 8 were not as well fitted by the model. However, size classes 9 and 10 provided a good fit between the modelled and observed values. The full results are not given here in the interests of brevity. The transition-probability matrix shows a similar structure to A8 , and the pattern in Fig. 6 was repeated. The 10 × 10 model was reapplied, using just two size classes for the emissions. The results were improved, as indicated in Table 2, where the sum of the least-square error is given, summed over time. The two-output model shows better performance than the three-output model. A 4 × 4 system model was also applied to the emissions data, and the genetic algorithm did produce A4 . However, the fit to the data was much worse, as is indicated by the least-squares error in Table 2. Note that the size classes used in the modelling are different in number and range, and a discrete comparison is not possible. However, the 4 × 4 system is an order of magnitude worse than all the other models, in terms of the least-squares error. This indicates that 16 parameters are not sufficient to represent the emissions data. 6.3. Sensitivity and errors The results in Figs. 6 and 7 were obtained from Eq. (2), using the best-fit matrix, A8 , given in Eq. (4), and this incorporated both attrition and agglomeration. The sensitivity of the emissions can be tested by (a) replacing A8 with A¯ 8 , where A¯ 8 is constructed from A8 by setting a¯ i j = 0 for i < j and retaining all of the remaining values; (b) by refitting the data to A∗8 , where ai∗j = 0, i < j, and recalculating the least-squares fit to the data to obtain new elements. The first method removes the agglomeration terms from the model, and the results are given in Fig. 8 for the number of emitted particles in size class 8. The results showed little difference for size class 7, and these are not shown. For size class 8, the difference becomes discernible at about four hours, and increases to about 5% difference at seven hours (the limit of experimental data). The mathematical solution of Eq. (2) can be expressed in terms of the eigenvalues of the matrix, raised to a power of time. This difference between the emissions will continue to increase in time. Refitting the data, using A∗8 in the genetic algorithm, gives

(6)

250

R.J. Dixon et al. / Mathematical and Computer Modelling 45 (2007) 241–251

Fig. 8. 8 × 8 Comparison of agglomeration + attrition versus attrition.

Fig. 9. Emission curves for A∗8 transition matrix.

and the corresponding emission curves are also shown in Fig. 9. The values of a7,7 and a8,8 for attrition only are less than for the values where attrition and agglomeration are considered (Eq. (4)). When agglomeration is permitted, small particles can combine to form larger particles, and have the opportunity of capture by the cyclones at least temporarily. The least-squares error summed over time (see Table 2) is almost the same as using the model with the matrix A8 . The diagonal elements of A8 and A∗8 generally differ by small amounts, less than 0.015, except for the two emission elements a7,7 and a8,8 , where the difference is as high as 0.076. The elements in columns 3, 4, 5, 6 and 7, beneath the diagonal, show changes of the same order and pattern, indicating a similarity in the attrition process. The subdiagonal elements in columns 1 and 2 of A8 and A∗8 show similar changes, but the pattern is less consistent. The break-up pattern of the particles is less well-defined by the model. Similar features were observed in the transition-probability matrices for the 6 × 6 and the 10 × 10 model. The break-up pattern of the particles can be displayed by plotting the subdiagonal columns of the transitionprobability matrix as a function of the mean diameter of the size class. Fig. 10 shows the break-up pattern of the largest sized particles for the 6 × 6, 8 × 8 and 10 × 10 models, and is essentially a plot of the first subdiagonal column of A6 , A8 and A10 . Generally, the large particles produce more particles by attrition, i.e. produce more small particles than by fracture.

R.J. Dixon et al. / Mathematical and Computer Modelling 45 (2007) 241–251

251

Fig. 10. Break-up probabilities for A6 , A8 , and A10 transition matrices.

7. Conclusion The transition-probability matrix model provides a novel solution to the difficult task of predicting and categorising emission trends from fluidised catalytic cracker units. The 6 × 6, 8 × 8 and 10 × 10 models have captured the initial emissions resulting from the quick break-up of relatively weak particles in the laboratory test rig. The 4 × 4 model contained insufficient parameters for the model to adequately fit the data. The processes of attrition are more important than fracture in the laboratory experiments, and agglomeration is relatively insignificant. In hightemperature environments, as in industrial FCCU’s, agglomeration may be more important. There is insufficient data currently available from an industrial FCCU to be able to apply the model with any confidence. This is an avenue for further experimentation, and the data requirements will imply extensive monitoring over start-up and recharge periods. The model will then provide information for use in predicting the probability size distribution of the catalyst and provide a tool for managing the catalyst. References [1] D. Geldart, Characterisation of fluidised powder, in: D. Geldart (Ed.), Gas Fluidisation Technology, John Wiley and Sons, UK, 1986, pp. 33–52. [2] M. Stittig, Petroleum Refining Industry: Energy Saving and Environmental Control, Noyes Data Corporation, USA, 1978. [3] D. Kunii, O. Levenspiel, Fluidisation Engineering, 2nd edition, Butterworth-Heinemann, USA, 1991. [4] S.T. Eagleson, E.H. Weaver, Controlling FCCU emissions, Int. J. Hydrocarbon Eng. 4 (1999) 76–80. [5] S.Y. Wu, J. Baeyens, C.Y. Chu, Effect of grid-velocity on attrition in fluidised beds, Can. J. Chem. Eng. 77 (1999) 738–744. [6] O.V. Barsukov, R.K. Nasirov, N.B. Artemova, L.A. Toplygina, Mechanical strength, activity and integrity of cracking catalysts under service conditions, Chem. Technol. Fuels Oils 34 (1998) 235–240. [7] J.M. Whitcombe, I.E. Agranovski, R.D. Braddock, F. Gandola, A. Hammond, Catalyst attrition and fracture due to thermal shock in fluidised catalytic cracker units, Chem. Eng. Commun. 191 (2004) 1259–1274. [8] J. Werther, J. Reppenhagen, Catalyst attrition in fluidised bed systems, AIChE J. 45 (1999) 2001–2010. [9] J.M. Whitcombe, Study of catalyst emissions from a fluidised catalytic cracker unit, Ph.D. Thesis, Griffith University, Australia, 2003. [10] M. Kot, Elements of Mathematical Ecology, C.U.P., Cambridge, UK, 2001. [11] F. Brauer, C. Castello-Ch´avez, Mathematical Models in Population Biology and Epidemiology, Springer, New York, 2001. [12] C.D. Cooper, F.C. Alley, Air Pollution Control: A Design Approach, Waveland Press, Inc., Prospect Heights, IL, 1994, pp. 127–150.