- Email: [email protected]
Application of the Monte Carlo method to process calculations for a catalyst flow model
Application of the Monte Carlo Method to Process Calculations For a Catalyst Flow Model J. D. COLTHART,
D. J. COLLINS*,
AND
J. E. GWYN
Houston Research Laboratory, Shell Oil Company, P.O. Box 100, Deer Park, Texas 77536 (USA) (Received: 31 March, 1970; In final form 12 August, 1970)
ABSTRACT
which contained eleven probabilistic to evaluate computer requirements.
The Monte Carlo technique can provide a means to model chemical processes which are too complex for analytic solutions. The technique was jirst tested on simple reactor-regenerator catalyst J¶OWmodels for which analytic solutions exist. The relative errors in predicted average catalyst properties were about 11Jn where n is the number of particles involved in the controlling step. The technique was then applied to a model of a commercial catalytic cracking unit. One hundred particles were followed for ten weeks of process time. This required about ten minutes of Univac 1108 computer time during which about 10,000,OOOrandom numbers were generated. This is quite reasonable for technological calculations.
THE MONTE
METHOD
The Monte Carlo method simply employs random chance and the convergence to a true probability in the limit. For example, if we follow a catalyst particle through a reactor-regenerator system with a properly weighted random choice each and every time there is an option (residence time in a vessel, loss based on cyclone efficiency, addition and withdrawal, etc.), we obtain a single, non-representative history. If we repeat this n times, we obtain n random histories. As n becomes larger and larger, both the distribution of properties and average of particle properties approach the ‘true’ (or analytic) values. Thus, the principle of the Monte Carlo method is quite simple. Some of the detailed considerations in applying the method are treated where they occur in the illustrative problems. The first problem by which the technique itself was tested was that of the simple reactor-regenerator system. The reactor system and its analytic solutions are described first. The Monte Carlo solutions are then considered.
INTRODUCTION Although the Monte Carlo technique has found numerous applications in stochastic processes, especially queueing problems and thermodynamics, it has found little use in chemical process models described by continuous functions. The Monte Carlo technique is proposed here as a method which, by following the random history of enough individual catalyst particles, will give a fairly representative cross-section of what happens to the entire population in a reactor-regenerator catalyst flow system. The purposes of this project were (1) to determine the accuracy that could be achieved by this technique, and (2) to determine the computer requirements for reasonably complex problems. In the first instance we considered Petersen’s problem of temporary and permanent deactivation in a simple reactor-regenerator system. 1 Analytic solutions were obtained against which the Monte Carlo technique could be evaluated. In the second instance a complex process model,
THE REACTOR-REGENERATOR
PROCESS
Petersen’ presented a theoretical treatment for the catalyst activity in the simple system in Fig. 1. He gave a complete steady state solution for a system with no withdrawal or addition, no permanent deactivation, well-mixed vessels, and exponential (first order) deactivation and regeneration. He also proposed a steady state solution for a general case including permanent deactivation and additions and withdrawals. The authors of this paper later proposed some modifications to that general solution2.
* Georgia Institute of Technology. 63 Chemical Engineering (2) (1971)-Q
CARLO
events, was used
Elsevier Publishing Company Ltd. England-Printed
in Great Britain
J. D. COLTHART,
D. J. COLLINS AND J. E. GWYN
p = O.lh
Reactor
A=
A,=A,
kg-
k,
Regenerator hG
AR
_p-0.01 0.0
I I
0
I
I --. 2
I
I 3
I
( 4
! 5
kc/A
pi--b @ hR , A G =
At
Fig. 2.
Residence
time
p
=
Fractional
flow
,AR
=
Average
0.6
constants rate
catalyst
for
Activity level leaving the reactor.
.
p = O.lh
withdrawal
A =A,=A,
activities
k,=
5k,
(O
Fig. 1.
Reactor-regenerator
A, Table
B. C. D refer A-l
0.7
to
z \
I
flow system.
*
0.6 I\
The basic assumptions are that the temporary activity is expressed as a fraction of the permanent (uiz., maximum recoverable) activity, and the rate constants for temporary deactivation and reactivation are not affected by the level of permanent activity. The actual activity of a particle at some time, t, is thus the product of temporary and permanent activity : A(t) = R(t)K(t)
(1)
where A is actual activity; R is temporary activity and K is permanent (maximum recoverable) activity. (A summary of Petersen’s analytic solutions as modified by the authors is given in Table A-l of the Appendix.) Values for the analytic solutions of the mean activity leaving the reactor (equals mean activity in the reactor for CSTR) were calculated using the equations in Table A-l and are shown in Figs. 2 and 3. For convenience, residence time constants for the reactor and regenerator were assumed equal; ?L= I,
= I,
and the permanent deactivation rate constant, p, was taken as 0.1 times the residence time constant; p = O.llt
P=o.o1 ----
0.0
0
I
2
3
P
4
-
5
k,/A
Fig. 3.
Activity level leaving the reactor.
For no permanent deactivation p = 0. In Fig. 2 the temporary deactivation and reactivation rate constants, k, and k,, are equal, and in Fig. 3, k, = 5k,. The mean activities are plotted against the ratio of deactivation rate constant to residence time constant, E = k,/L The fractional withdrawal/addition, j3, (relative to circulation) is the parameter for the
APPLICATION
OF THE MONTE CARLO METHOD TO PROCESS CALCULATIONS
family of curves. These results served as the standard for the Monte Carlo calculations. Note that the problem would become considerably more complex for less ideal mixing. One might then be interested in the average activity in certain regions of the reactor or regenerator. Such a situation was not considered in this work.
MONTE CARLO SIMULATION: REACTORREGENERATOR SYSTEM WITH NO PERMANENT DEACTIVATION, ADDITIONS, OR WITHDRAWALS In this case the only random variables in the process are the residence times within vessels. For a wellmixed vessel the residence time distribution is F(t) = 1 - e-”
o
(2)
where F(‘(t) = the fraction of particles with residence times equal to or less than t; and 2 is the reciprocal of the mean residence time. F is a flat probability with values between 0 and 1. The random number generator used also has a flat probability between 0 and 1. Thus, when a particle is about to enter a vessel, its residence time in that vessel can be determined by obtaining a random number, RN, and substituting it for F(t) in eqn. (2). The residence time is then 1
1
t = n ln (1 - RN) and the activation or deactivation that occurs in that pass through the vessel is an algebraic function of the chosen time, t.
Note
: Each
point 1s the werage over 011 cycles of I2 monte carlo slmulatlons. wth different random number sequences. at
:
k,/h,=O5,
1.0,Z.O. 4.0,
Only one particle I” each c(1se
Fig. 4.
Mean absolute perczm; 7~
80
followed
over the range 0.5 5 kR/
65
The calculational approach for this system was to follow a single particle (and its activity) around the system for a given number of cycles, then compute the average activity leaving the reactor (or at any other point in the system) during those cycles. The average absolute percent errors resulting from this approach, basis the exact analytic solutions in Figs. 2 and 3, are shown in Fig. 4 as a function of the number of cycles followed. The relative error was approximately equal to the inverse of the square root of number of cycles. The error appeared to be independent of the value of k,ll, but varied slightly with the ratio of kg/k,.. One thousand cycles reduced the error to about 3 %.
MONTE CARLO SIMULATION: REACTORREGENERATOR SYSTEM WITH PERMANENT DEACTIVATION AND ADDITIONS AND WITHDRAWALS The approach used depends upon the information about the system that is desired. To obtain the steady-state conditions for continuous withdrawals and additions of catalyst, a number of particles must be circulated for a sufficient number of cycles that the particles as a whole have forgotten the initial conditions. For example, in this work the simulation was usually begun with the system assumed to contain only particles of unit activity. The choice was one of convenience only. The same steady state distribution resulted when using an initial state which described a distribution of activities similar to that which was eventually determined for the steady state. During these cycles particles are randomly withdrawn and replaced by fresh catalyst of unit activity in accordance with the replacement probabilities. After steady-state is obtained, the properties for each particle leaving a vessel in a given cycle are noted and appropriate average properties are calculated. In our case, we are interested in the properties of particles leaving a vessel during steady state. For this, the technique of averaging properties at equal numbers of cycles is correct in the ensemble sense. For the entire system, the particles could be examined wherever they occur at some specified time. However, it is not possible to examine particles leaving a vessel on a preset time basis. At some preset time a particle is always in a vessel and it is necessary to wait for it to leave before measuring its properties. Thus, looking at all particles at fixed points in ‘real’ or ‘process’ time biases the results: in fact, the weighting factor is proportional to their residence time. For example, consider the hypothetical situation illustrated in Fig. 5, which illustrates the history of one of a number of
J. D. COLTHART,
66
D. J. COLLINS AND J. E. GWYN
particles to be considered in a total recycle reactorregenerator system where there is a 50/50 probability of either a very long or very short residence time in either vessel. If the ‘bookkeeping’ is done at equal time intervals, it will appear that, when viewed, most of the particles leaving either vessel are experiencing
1
REAL(PROCESS)
Fig. 5.
TIME
TABLE 1 FUNCTIONS
No.
Section
kg=
Sk,
o
kg=
kr
I
I
I
200
300
400
Number
2 Regenerator
Residence Time
Coke burning/ reactivationfunction of time
3 Elutriation
Separation efficiency
Lost or returned to No. 1, regenerator inlet
4 Withdrawal
Route-based flow rates
on Routing only-withdraw or pass on to reactor inlet
5 Reactor Inlet
Predetermine route
exit Elutriate or pass on to stripper inlet
6 Reactor
Residence time
7 Elutriation
Separation Lost or returned to efficiency includes No. 5, reactor inlet clarifier
8 Stripper
Residence time and/or stripping efficiency
9 Riser
Contacting Temperature rise, efficiency (oxygen loss of surface area/ utilisation) permanent deactivation Attrition-age dependent only
of
/,!,, 600
particles,
II
N
10 Top of Riser
Fig. 6. Monte Carlo simulation for the model of Petersen with no withdrawals and no permanent deactivation.
20; 0
Process function
1 Regenerator Inlet Predetermine exit Exit routing omyroute; probability elutriated or direct to regenerator exit based on flow rates
0
100
Probabilistic function
2
a
11 Top of Riser 0
CRACKING
UNIT MODEL
Residence time history.
k,/A=
2.
IN A CATALYTIC
Particle sizeattrited fines
______Particle size of new particle
Coke laydown/ temporary deactivation-function of time
Dummy, not used at present
Accumulate fractional masses attrited-when unit mass obtained assign and store properties of last particle attrited, assign a fines particle size ___Make up lost particle-with attrited fines if available, use new particle if not
0
6-
21
100
I 200 Number
Fig. 7.
I
I
300
400
of particles,
I
I,,, 600
II
1
N
Monte Carlo simulation for the model of Petersen with fraction withdrawn = 0.10.
a long residence time (and high degree of deactivation or regeneration) even though we know that at any instant in time, one-half of the particles leaving either vessel will have experienced a short residence time. On the other hand, at the end of a preset number of cycles, the residence time of each particle is independent and distributed in accordance with the vessel residence time distribution and is thus unbiased. This cycle approach was the one used. The test calculations are presented in Figs. 6 through 8. When permanent deactivation and withdrawals were set to zero (Fig. 6), the average error
APPLICATION
OF THE MONTE CARLO METHOD TO PROCESS CALCULATIONS
decreased as the inverse of the square root of the number of particles. For withdrawals and permanent deactivation (Figs. 7 and 8) the errors also varied about as 1 Jn, and errors on the order of about 3 % were obtained with 1000 particles. Thus, in this case, accurate results can be obtained by considering a reasonable number of particles. 20 kr/A 0
Numbe of partlcles,
Fig. 8.
Add
= I
kg=kr
N
Monte Carlo simulation for the model of Petersen with fraction withdrawn = 0.30.
Elutriated and lost
APPLICATION
\
Clarified oil
Withdraw1
3
I
I
Regen.
Riser
Fig. 9.
COMPLEX
The Monte Carlo technique was applied to a realistic reactor-regenerator process model involving eleven random decisions and associated calculations for each cycle of a particle through a system. A flow diagram of the model is shown in Fig. 9, and the functions to be evaluated at each ‘decision point’ are described in Table 1. One hundred particles were followed for ten weeks of process time (- 10,000,000 random numbers) in ten minutes of Univac 1108 running time. The number of particles and the process time can be varied according to the accuracy desired and the time constants of the important process steps. As defined in this model, each ‘particle’ actually represented a unit mass of particles of the same size rather than an individual particle. This was strictly a practical consideration. Based on the particle size distribution by number, almost all particles are ‘fines’. The average residence time of fines in the system is relatively short. Thus, on a number basis, most of the computer time would be spent considering smaller particles which have very little effect on the properties of the inventory. Furthermore, cracking severity and coke-burning parameters are commonly expressed on a weight basis.
h-A 0
TO MORE PROCESSES
67
Catalytic cracking unit catalyst flow.
68
J. D. COLTHART,
Figure 10 shows a steady-state carbon-on-catalyst calculated after
o-2-
‘\\
Regenerated \
catalyst
spent
distribution simulating
D. J. COLLINS AND J. E. GWYN
of the
t = Mean value
integrated rate equations, and exponentiations in general. In summary, the Monte Carlo technique provided us with a simple-minded method for investigating a sophisticated process problem for which no analytic solution could be found. The computer time required was not found to be excessive.
APPENDIX TABLE A-l SUMMARY OF ANALYTIC SOLUTIONS FOR WELL-MIXED VESSELS OF EQUAL RESIDENCE TIMES Carbon on catalyst, %
Fig. 10.
w
Coke distributions on 100 catalyst particles by Monte Carlo simulation.
histories of 100 unit mass particles for ten weeks. Here, the simulation was begun at a point in time where it was assumed that all particles in the system had unit activity. The computer time required to define a steady-state distribution of properties could be somewhat reduced by assuming a more realistic initial state which more closely approximates the final distribution of properties. However, in all cases investigated, the model eventually came to the same steady state, regardless of the initial conditions chosen. This was expected since fresh catalyst is continually being added to compensate for losses from the regenerator cyclones, losses in ‘slurry oil’ from the reactor, and planned withdrawals.
Dejinitions and Assumptions: Assume both the reactor and regenerator mixed vessels, and let I = $
= (mean residence
time)-’
COMMENTS
Random number generators are numerous and specific to the computer used. The statistical properties of those used were not investigated. However, a random number generator which repeated its sequence about every 35,000 numbers gave systematic and biased results. The results often do not converge strongly and a large number of calculations and random numbers can be required. The generator used in the preceding example did not repeat in the 10,000,OOO numbers used. Speed in process calculations is also important for efficient application of the Monte Carlo method. Table lookup is much faster than calculations, especially in exponentiation and the extraction of logarithms. This approach was used in obtaining residence times from random numbers (RN was converted to integers from 1 to 1000 as an index in previously calculated residence times), in evaluating
(4)
m be set equal for both. Activity from temporary particle resides in the reactor
deactivation after for time, 8, is
R(8) = emkre
a
(5)
Activity from temporary reactivation after particle resides in the regenerator for time tl is G(or) = 1 - eek$L Activity (maximum deactivation
recoverable)
from permanent (7)
where t is the particle age. Define j3 = fractional makeup/withdrawal catalyst flowing through the reactor Summary Results : 1. Mean Maximum Regenerator,
Recoverable
a
(6)
K(t) = eCpr GENERAL
are well-
Activity
basis (8)
Leaving
(9) 2. Mean Maximum Regenerator (Leaving K
Recoverable Activity Entering Reactor, Withdrawn),
=
D p2
B4P + 4 -t 2pn + gn
(10)
3. Mean Maximum Recoverable Activity of Stream Leaving Reactor Which Came From Regenerator (Excludes First Pass of Makeup Catalyst),
PA”
K, =
(P + UP2
+ 2Pl +
fin21
(II)
APPLICATION OF THE MONTE CARLO METHOD TO PROCESS CALCULATIONS
4. Mean Maximum Recoverable Activity of Stream Leaving Reactor Which Came From Fresh Feed (Without Recirculation),
n
Krf =
K(t)
(12)
4
(recoverable) of time
activity
mean maximum recoverable vity leaving regenerator
*
(P +
permanent a function
KD
5. Mean Activity in Reactor Equals Mean Activity of Stream Leaving. Similarly for Regenerator, Since Both Are Assumed to be Well-Mixed. 6. Mean Activity of Catalyst Entering the Reactor,
mean maximum regenerator
69
of catalyst
(permanent)
recoverable
as
acti-
activity entering
mean maximum recoverable activity of portion of stream leaving reactor which came from regenerator
Ai A.
=
P + (1 - PVM’, - fg)
I
1 -
mean maximum recoverable of stream leaving reactor fresh catalyst
(13)
(I - /3)t,llt,l
where t’, =
deactivation of activity
rate constant
kg
reactivation
rate constant
n, N
number
P
deactivation of activity
R(t)
temporary permanent
RN
random number; distribution
t
time or particle
age
TWI
mean residence
time in a vessel
t ‘sl
Wg
+ P + A)
t ‘P t ‘R
Uk,
+ P + 2)
A
6%+
t’, = ___
P +
(14)
4
A
A
t’g= (kg + p of Catalyst
+ lj Leaving the Reactor,
A0 A, = A, 8. Mean generator
Activity
A
for temporary
loss
for permanent
loss
of particles
(15)
(P + 2)
7. Mean Activity
activity of portion which entered as
= Ait’,
(17)
of Catalyst
Leaving
the
=(Ai-B) G U-8)
Re-
(18)
ACKNOWLEDGEMENT The authors wish to thank Shell Oil Company permission to publish this article.
rate constant
activity activity
of catalyst
as fraction
0 I RN I
1
with
of flat
W(P + A>
for Greek Symbols
Note: This paper is based upon material presented at the Annual Meeting of the A.1.Ch.E. Chicago, December, 1970.
a
time in regenerator
B
ratio of catalyst withdrawal rate
&
ratio of deactivation time constant
A
residence
6
time in reactor
rate to circulation
rate constant
to residence
NOMENCLATURE actual activity time; A = RK
of catalyst
as a function
of
time constant,
time- I
mean activity of catalyst leaving the regenerator mean activity
of catalyst
entering
mean activity
of catalyst
leaving the reactor
mean activity
of catalyst
in reactor
residence
time distribution
the reactor
0 < F(t) 5 1
REFERENCES 1. PETERSEN, E. E., A.Z.Ch.E.
.Z. 1960 6
p. 488.
2. GWYN, J. E., AND COLTHART, J. D., A.Z.Ch.E. p. 932.
J. 1969 15
70
J. D. COLTHART,
D. J. COLLINS
RESUME La me’thode de Monte-Carlo peut permettre de calculer les processus chimiques trop complexes pour qu’on dispose de solutions analytiques. Elle a tout d’abord e’te’ testee sur des systemes simples reacteurregene’rateur avec circulation de catalyseur, pour lesquels ces solutions existent. Si n designe le nombre de particules de catalyseur intervenant dans P&ape controlant le processus, les erreurs relatives commises sur les valeurs estime’es des proprie’tes moyennes du catalyseur sont de I’ordre de l/Jn. La methode a et& ensuite applique’e au modele dune unite’ industrielle de craking catalytique. Cent particules ont ete suivies au tours de dix semaines de fonctionnement de I’unite’. Le temps de calcul sur Univac 1108 est environ de dix minutes et, pendant ce temps, environ 10 millions de nombre au hasard sont engendres, ce qui est tout d fait raisonnable pour des calculs techniques.
AND J. E. GWYN
Z USAMMENFASSUNG Die Monte-Carlo-Technik eignet sich, chemische Prozesse zu modellieren, die fiir eine analytische Liisung zu komplex sind. Diese Technik wurde zunachst auf einfache Reaktor-Regenerator-Katalysator-Strom-Modelle angewendet, fur welche analytische Liisungen bereits vorliegen. Die relativen Fehler in den vorausgesagten mittleren Katalysatoreigenschaften betrugen etwa I/l/n, wobei n die Zahl der Teilchen ist, die an dem bestimmenden Schritt beteiligt ist. Diese Technik wurde dann auf ein Model1 einer industriellen katalytischen Krackanlage angewendet. Einhundert Teilchen wurden fiir die Dauer von zehn Wochen Betriebszeit beobachtet. Dazu beniitigte man etwa 10 Minuten Rechenzeit an einem Univac110%Rechner, wiihrend der 10’ Zufallszahlen entwickelt wurden. Dieser Aufwand ist durchaus annehmbar fiir technologische Berechnungen.
D. J. COLLINS*,
AND
J. E. GWYN
Houston Research Laboratory, Shell Oil Company, P.O. Box 100, Deer Park, Texas 77536 (USA) (Received: 31 March, 1970; In final form 12 August, 1970)
ABSTRACT
which contained eleven probabilistic to evaluate computer requirements.
The Monte Carlo technique can provide a means to model chemical processes which are too complex for analytic solutions. The technique was jirst tested on simple reactor-regenerator catalyst J¶OWmodels for which analytic solutions exist. The relative errors in predicted average catalyst properties were about 11Jn where n is the number of particles involved in the controlling step. The technique was then applied to a model of a commercial catalytic cracking unit. One hundred particles were followed for ten weeks of process time. This required about ten minutes of Univac 1108 computer time during which about 10,000,OOOrandom numbers were generated. This is quite reasonable for technological calculations.
THE MONTE
METHOD
The Monte Carlo method simply employs random chance and the convergence to a true probability in the limit. For example, if we follow a catalyst particle through a reactor-regenerator system with a properly weighted random choice each and every time there is an option (residence time in a vessel, loss based on cyclone efficiency, addition and withdrawal, etc.), we obtain a single, non-representative history. If we repeat this n times, we obtain n random histories. As n becomes larger and larger, both the distribution of properties and average of particle properties approach the ‘true’ (or analytic) values. Thus, the principle of the Monte Carlo method is quite simple. Some of the detailed considerations in applying the method are treated where they occur in the illustrative problems. The first problem by which the technique itself was tested was that of the simple reactor-regenerator system. The reactor system and its analytic solutions are described first. The Monte Carlo solutions are then considered.
INTRODUCTION Although the Monte Carlo technique has found numerous applications in stochastic processes, especially queueing problems and thermodynamics, it has found little use in chemical process models described by continuous functions. The Monte Carlo technique is proposed here as a method which, by following the random history of enough individual catalyst particles, will give a fairly representative cross-section of what happens to the entire population in a reactor-regenerator catalyst flow system. The purposes of this project were (1) to determine the accuracy that could be achieved by this technique, and (2) to determine the computer requirements for reasonably complex problems. In the first instance we considered Petersen’s problem of temporary and permanent deactivation in a simple reactor-regenerator system. 1 Analytic solutions were obtained against which the Monte Carlo technique could be evaluated. In the second instance a complex process model,
THE REACTOR-REGENERATOR
PROCESS
Petersen’ presented a theoretical treatment for the catalyst activity in the simple system in Fig. 1. He gave a complete steady state solution for a system with no withdrawal or addition, no permanent deactivation, well-mixed vessels, and exponential (first order) deactivation and regeneration. He also proposed a steady state solution for a general case including permanent deactivation and additions and withdrawals. The authors of this paper later proposed some modifications to that general solution2.
* Georgia Institute of Technology. 63 Chemical Engineering (2) (1971)-Q
CARLO
events, was used
Elsevier Publishing Company Ltd. England-Printed
in Great Britain
J. D. COLTHART,
D. J. COLLINS AND J. E. GWYN
p = O.lh
Reactor
A=
A,=A,
kg-
k,
Regenerator hG
AR
_p-0.01 0.0
I I
0
I
I --. 2
I
I 3
I
( 4
! 5
kc/A
pi--b @ hR , A G =
At
Fig. 2.
Residence
time
p
=
Fractional
flow
,AR
=
Average
0.6
constants rate
catalyst
for
Activity level leaving the reactor.
.
p = O.lh
withdrawal
A =A,=A,
activities
k,=
5k,
(O
Fig. 1.
Reactor-regenerator
A, Table
B. C. D refer A-l
0.7
to
z \
I
flow system.
*
0.6 I\
The basic assumptions are that the temporary activity is expressed as a fraction of the permanent (uiz., maximum recoverable) activity, and the rate constants for temporary deactivation and reactivation are not affected by the level of permanent activity. The actual activity of a particle at some time, t, is thus the product of temporary and permanent activity : A(t) = R(t)K(t)
(1)
where A is actual activity; R is temporary activity and K is permanent (maximum recoverable) activity. (A summary of Petersen’s analytic solutions as modified by the authors is given in Table A-l of the Appendix.) Values for the analytic solutions of the mean activity leaving the reactor (equals mean activity in the reactor for CSTR) were calculated using the equations in Table A-l and are shown in Figs. 2 and 3. For convenience, residence time constants for the reactor and regenerator were assumed equal; ?L= I,
= I,
and the permanent deactivation rate constant, p, was taken as 0.1 times the residence time constant; p = O.llt
P=o.o1 ----
0.0
0
I
2
3
P
4
-
5
k,/A
Fig. 3.
Activity level leaving the reactor.
For no permanent deactivation p = 0. In Fig. 2 the temporary deactivation and reactivation rate constants, k, and k,, are equal, and in Fig. 3, k, = 5k,. The mean activities are plotted against the ratio of deactivation rate constant to residence time constant, E = k,/L The fractional withdrawal/addition, j3, (relative to circulation) is the parameter for the
APPLICATION
OF THE MONTE CARLO METHOD TO PROCESS CALCULATIONS
family of curves. These results served as the standard for the Monte Carlo calculations. Note that the problem would become considerably more complex for less ideal mixing. One might then be interested in the average activity in certain regions of the reactor or regenerator. Such a situation was not considered in this work.
MONTE CARLO SIMULATION: REACTORREGENERATOR SYSTEM WITH NO PERMANENT DEACTIVATION, ADDITIONS, OR WITHDRAWALS In this case the only random variables in the process are the residence times within vessels. For a wellmixed vessel the residence time distribution is F(t) = 1 - e-”
o
(2)
where F(‘(t) = the fraction of particles with residence times equal to or less than t; and 2 is the reciprocal of the mean residence time. F is a flat probability with values between 0 and 1. The random number generator used also has a flat probability between 0 and 1. Thus, when a particle is about to enter a vessel, its residence time in that vessel can be determined by obtaining a random number, RN, and substituting it for F(t) in eqn. (2). The residence time is then 1
1
t = n ln (1 - RN) and the activation or deactivation that occurs in that pass through the vessel is an algebraic function of the chosen time, t.
Note
: Each
point 1s the werage over 011 cycles of I2 monte carlo slmulatlons. wth different random number sequences. at
:
k,/h,=O5,
1.0,Z.O. 4.0,
Only one particle I” each c(1se
Fig. 4.
Mean absolute perczm; 7~
80
followed
over the range 0.5 5 kR/
65
The calculational approach for this system was to follow a single particle (and its activity) around the system for a given number of cycles, then compute the average activity leaving the reactor (or at any other point in the system) during those cycles. The average absolute percent errors resulting from this approach, basis the exact analytic solutions in Figs. 2 and 3, are shown in Fig. 4 as a function of the number of cycles followed. The relative error was approximately equal to the inverse of the square root of number of cycles. The error appeared to be independent of the value of k,ll, but varied slightly with the ratio of kg/k,.. One thousand cycles reduced the error to about 3 %.
MONTE CARLO SIMULATION: REACTORREGENERATOR SYSTEM WITH PERMANENT DEACTIVATION AND ADDITIONS AND WITHDRAWALS The approach used depends upon the information about the system that is desired. To obtain the steady-state conditions for continuous withdrawals and additions of catalyst, a number of particles must be circulated for a sufficient number of cycles that the particles as a whole have forgotten the initial conditions. For example, in this work the simulation was usually begun with the system assumed to contain only particles of unit activity. The choice was one of convenience only. The same steady state distribution resulted when using an initial state which described a distribution of activities similar to that which was eventually determined for the steady state. During these cycles particles are randomly withdrawn and replaced by fresh catalyst of unit activity in accordance with the replacement probabilities. After steady-state is obtained, the properties for each particle leaving a vessel in a given cycle are noted and appropriate average properties are calculated. In our case, we are interested in the properties of particles leaving a vessel during steady state. For this, the technique of averaging properties at equal numbers of cycles is correct in the ensemble sense. For the entire system, the particles could be examined wherever they occur at some specified time. However, it is not possible to examine particles leaving a vessel on a preset time basis. At some preset time a particle is always in a vessel and it is necessary to wait for it to leave before measuring its properties. Thus, looking at all particles at fixed points in ‘real’ or ‘process’ time biases the results: in fact, the weighting factor is proportional to their residence time. For example, consider the hypothetical situation illustrated in Fig. 5, which illustrates the history of one of a number of
J. D. COLTHART,
66
D. J. COLLINS AND J. E. GWYN
particles to be considered in a total recycle reactorregenerator system where there is a 50/50 probability of either a very long or very short residence time in either vessel. If the ‘bookkeeping’ is done at equal time intervals, it will appear that, when viewed, most of the particles leaving either vessel are experiencing
1
REAL(PROCESS)
Fig. 5.
TIME
TABLE 1 FUNCTIONS
No.
Section
kg=
Sk,
o
kg=
kr
I
I
I
200
300
400
Number
2 Regenerator
Residence Time
Coke burning/ reactivationfunction of time
3 Elutriation
Separation efficiency
Lost or returned to No. 1, regenerator inlet
4 Withdrawal
Route-based flow rates
on Routing only-withdraw or pass on to reactor inlet
5 Reactor Inlet
Predetermine route
exit Elutriate or pass on to stripper inlet
6 Reactor
Residence time
7 Elutriation
Separation Lost or returned to efficiency includes No. 5, reactor inlet clarifier
8 Stripper
Residence time and/or stripping efficiency
9 Riser
Contacting Temperature rise, efficiency (oxygen loss of surface area/ utilisation) permanent deactivation Attrition-age dependent only
of
/,!,, 600
particles,
II
N
10 Top of Riser
Fig. 6. Monte Carlo simulation for the model of Petersen with no withdrawals and no permanent deactivation.
20; 0
Process function
1 Regenerator Inlet Predetermine exit Exit routing omyroute; probability elutriated or direct to regenerator exit based on flow rates
0
100
Probabilistic function
2
a
11 Top of Riser 0
CRACKING
UNIT MODEL
Residence time history.
k,/A=
2.
IN A CATALYTIC
Particle sizeattrited fines
______Particle size of new particle
Coke laydown/ temporary deactivation-function of time
Dummy, not used at present
Accumulate fractional masses attrited-when unit mass obtained assign and store properties of last particle attrited, assign a fines particle size ___Make up lost particle-with attrited fines if available, use new particle if not
0
6-
21
100
I 200 Number
Fig. 7.
I
I
300
400
of particles,
I
I,,, 600
II
1
N
Monte Carlo simulation for the model of Petersen with fraction withdrawn = 0.10.
a long residence time (and high degree of deactivation or regeneration) even though we know that at any instant in time, one-half of the particles leaving either vessel will have experienced a short residence time. On the other hand, at the end of a preset number of cycles, the residence time of each particle is independent and distributed in accordance with the vessel residence time distribution and is thus unbiased. This cycle approach was the one used. The test calculations are presented in Figs. 6 through 8. When permanent deactivation and withdrawals were set to zero (Fig. 6), the average error
APPLICATION
OF THE MONTE CARLO METHOD TO PROCESS CALCULATIONS
decreased as the inverse of the square root of the number of particles. For withdrawals and permanent deactivation (Figs. 7 and 8) the errors also varied about as 1 Jn, and errors on the order of about 3 % were obtained with 1000 particles. Thus, in this case, accurate results can be obtained by considering a reasonable number of particles. 20 kr/A 0
Numbe of partlcles,
Fig. 8.
Add
= I
kg=kr
N
Monte Carlo simulation for the model of Petersen with fraction withdrawn = 0.30.
Elutriated and lost
APPLICATION
\
Clarified oil
Withdraw1
3
I
I
Regen.
Riser
Fig. 9.
COMPLEX
The Monte Carlo technique was applied to a realistic reactor-regenerator process model involving eleven random decisions and associated calculations for each cycle of a particle through a system. A flow diagram of the model is shown in Fig. 9, and the functions to be evaluated at each ‘decision point’ are described in Table 1. One hundred particles were followed for ten weeks of process time (- 10,000,000 random numbers) in ten minutes of Univac 1108 running time. The number of particles and the process time can be varied according to the accuracy desired and the time constants of the important process steps. As defined in this model, each ‘particle’ actually represented a unit mass of particles of the same size rather than an individual particle. This was strictly a practical consideration. Based on the particle size distribution by number, almost all particles are ‘fines’. The average residence time of fines in the system is relatively short. Thus, on a number basis, most of the computer time would be spent considering smaller particles which have very little effect on the properties of the inventory. Furthermore, cracking severity and coke-burning parameters are commonly expressed on a weight basis.
h-A 0
TO MORE PROCESSES
67
Catalytic cracking unit catalyst flow.
68
J. D. COLTHART,
Figure 10 shows a steady-state carbon-on-catalyst calculated after
o-2-
‘\\
Regenerated \
catalyst
spent
distribution simulating
D. J. COLLINS AND J. E. GWYN
of the
t = Mean value
integrated rate equations, and exponentiations in general. In summary, the Monte Carlo technique provided us with a simple-minded method for investigating a sophisticated process problem for which no analytic solution could be found. The computer time required was not found to be excessive.
APPENDIX TABLE A-l SUMMARY OF ANALYTIC SOLUTIONS FOR WELL-MIXED VESSELS OF EQUAL RESIDENCE TIMES Carbon on catalyst, %
Fig. 10.
w
Coke distributions on 100 catalyst particles by Monte Carlo simulation.
histories of 100 unit mass particles for ten weeks. Here, the simulation was begun at a point in time where it was assumed that all particles in the system had unit activity. The computer time required to define a steady-state distribution of properties could be somewhat reduced by assuming a more realistic initial state which more closely approximates the final distribution of properties. However, in all cases investigated, the model eventually came to the same steady state, regardless of the initial conditions chosen. This was expected since fresh catalyst is continually being added to compensate for losses from the regenerator cyclones, losses in ‘slurry oil’ from the reactor, and planned withdrawals.
Dejinitions and Assumptions: Assume both the reactor and regenerator mixed vessels, and let I = $
= (mean residence
time)-’
COMMENTS
Random number generators are numerous and specific to the computer used. The statistical properties of those used were not investigated. However, a random number generator which repeated its sequence about every 35,000 numbers gave systematic and biased results. The results often do not converge strongly and a large number of calculations and random numbers can be required. The generator used in the preceding example did not repeat in the 10,000,OOO numbers used. Speed in process calculations is also important for efficient application of the Monte Carlo method. Table lookup is much faster than calculations, especially in exponentiation and the extraction of logarithms. This approach was used in obtaining residence times from random numbers (RN was converted to integers from 1 to 1000 as an index in previously calculated residence times), in evaluating
(4)
m be set equal for both. Activity from temporary particle resides in the reactor
deactivation after for time, 8, is
R(8) = emkre
a
(5)
Activity from temporary reactivation after particle resides in the regenerator for time tl is G(or) = 1 - eek$L Activity (maximum deactivation
recoverable)
from permanent (7)
where t is the particle age. Define j3 = fractional makeup/withdrawal catalyst flowing through the reactor Summary Results : 1. Mean Maximum Regenerator,
Recoverable
a
(6)
K(t) = eCpr GENERAL
are well-
Activity
basis (8)
Leaving
(9) 2. Mean Maximum Regenerator (Leaving K
Recoverable Activity Entering Reactor, Withdrawn),
=
D p2
B4P + 4 -t 2pn + gn
(10)
3. Mean Maximum Recoverable Activity of Stream Leaving Reactor Which Came From Regenerator (Excludes First Pass of Makeup Catalyst),
PA”
K, =
(P + UP2
+ 2Pl +
fin21
(II)
APPLICATION OF THE MONTE CARLO METHOD TO PROCESS CALCULATIONS
4. Mean Maximum Recoverable Activity of Stream Leaving Reactor Which Came From Fresh Feed (Without Recirculation),
n
Krf =
K(t)
(12)
4
(recoverable) of time
activity
mean maximum recoverable vity leaving regenerator
*
(P +
permanent a function
KD
5. Mean Activity in Reactor Equals Mean Activity of Stream Leaving. Similarly for Regenerator, Since Both Are Assumed to be Well-Mixed. 6. Mean Activity of Catalyst Entering the Reactor,
mean maximum regenerator
69
of catalyst
(permanent)
recoverable
as
acti-
activity entering
mean maximum recoverable activity of portion of stream leaving reactor which came from regenerator
Ai A.
=
P + (1 - PVM’, - fg)
I
1 -
mean maximum recoverable of stream leaving reactor fresh catalyst
(13)
(I - /3)t,llt,l
where t’, =
deactivation of activity
rate constant
kg
reactivation
rate constant
n, N
number
P
deactivation of activity
R(t)
temporary permanent
RN
random number; distribution
t
time or particle
age
TWI
mean residence
time in a vessel
t ‘sl
Wg
+ P + A)
t ‘P t ‘R
Uk,
+ P + 2)
A
6%+
t’, = ___
P +
(14)
4
A
A
t’g= (kg + p of Catalyst
+ lj Leaving the Reactor,
A0 A, = A, 8. Mean generator
Activity
A
for temporary
loss
for permanent
loss
of particles
(15)
(P + 2)
7. Mean Activity
activity of portion which entered as
= Ait’,
(17)
of Catalyst
Leaving
the
=(Ai-B) G U-8)
Re-
(18)
ACKNOWLEDGEMENT The authors wish to thank Shell Oil Company permission to publish this article.
rate constant
activity activity
of catalyst
as fraction
0 I RN I
1
with
of flat
W(P + A>
for Greek Symbols
Note: This paper is based upon material presented at the Annual Meeting of the A.1.Ch.E. Chicago, December, 1970.
a
time in regenerator
B
ratio of catalyst withdrawal rate
&
ratio of deactivation time constant
A
residence
6
time in reactor
rate to circulation
rate constant
to residence
NOMENCLATURE actual activity time; A = RK
of catalyst
as a function
of
time constant,
time- I
mean activity of catalyst leaving the regenerator mean activity
of catalyst
entering
mean activity
of catalyst
leaving the reactor
mean activity
of catalyst
in reactor
residence
time distribution
the reactor
0 < F(t) 5 1
REFERENCES 1. PETERSEN, E. E., A.Z.Ch.E.
.Z. 1960 6
p. 488.
2. GWYN, J. E., AND COLTHART, J. D., A.Z.Ch.E. p. 932.
J. 1969 15
70
J. D. COLTHART,
D. J. COLLINS
RESUME La me’thode de Monte-Carlo peut permettre de calculer les processus chimiques trop complexes pour qu’on dispose de solutions analytiques. Elle a tout d’abord e’te’ testee sur des systemes simples reacteurregene’rateur avec circulation de catalyseur, pour lesquels ces solutions existent. Si n designe le nombre de particules de catalyseur intervenant dans P&ape controlant le processus, les erreurs relatives commises sur les valeurs estime’es des proprie’tes moyennes du catalyseur sont de I’ordre de l/Jn. La methode a et& ensuite applique’e au modele dune unite’ industrielle de craking catalytique. Cent particules ont ete suivies au tours de dix semaines de fonctionnement de I’unite’. Le temps de calcul sur Univac 1108 est environ de dix minutes et, pendant ce temps, environ 10 millions de nombre au hasard sont engendres, ce qui est tout d fait raisonnable pour des calculs techniques.
AND J. E. GWYN
Z USAMMENFASSUNG Die Monte-Carlo-Technik eignet sich, chemische Prozesse zu modellieren, die fiir eine analytische Liisung zu komplex sind. Diese Technik wurde zunachst auf einfache Reaktor-Regenerator-Katalysator-Strom-Modelle angewendet, fur welche analytische Liisungen bereits vorliegen. Die relativen Fehler in den vorausgesagten mittleren Katalysatoreigenschaften betrugen etwa I/l/n, wobei n die Zahl der Teilchen ist, die an dem bestimmenden Schritt beteiligt ist. Diese Technik wurde dann auf ein Model1 einer industriellen katalytischen Krackanlage angewendet. Einhundert Teilchen wurden fiir die Dauer von zehn Wochen Betriebszeit beobachtet. Dazu beniitigte man etwa 10 Minuten Rechenzeit an einem Univac110%Rechner, wiihrend der 10’ Zufallszahlen entwickelt wurden. Dieser Aufwand ist durchaus annehmbar fiir technologische Berechnungen.