ChemicdEngineering
Science,
1972, Vol. 27, pp. 1497-U
13.
Pergamon
Press.
Printed in Great Britain
Monte Carlo methods of simulating micromixing in chemical reactors C. R. TRELEAVENt
and A. H. TOBGYS
Department of Chemical Engineering, University of Exeter, England. (Received 16 July 1971; Accepted 25 January 1972)
Abstract-The random coalescence model, simulated by Monte Carlo techniques, has been extended to the case of reactors for which the two separate reactant feed streams have differing residence time distributions. The parameter characterising the micromixing is the average number of collisions each fluid element experiences in its sojourn in the reactor vessel. This parameter has been varied for a simple second order, isothermal, homogeneous reaction for two arbitrarily chosen sets of residence time distributions, and it has been shown that an infhtite collision rate corresponds to one of the maximum mizedness models presented by the authors in a previous paper. Methods by which the micromizing parameter can be measured, and the value subsequently used for predicting the outcome of any given reaction are discussed. Comments are also offered on the reliability of Monte Carlo techniques, in terms of the parameters affecting the averaged output conversions, and on possible ways of reducing computation times.
1. INTRODUCTION
the pioneering work of Danckwerts [ 11and Zwietering [2], who described the extreme conditions of complete segregation and maximum mixedness for chemical reactors having arbitrary residence time distributions, a considerable amount of attention has been devoted to the characterisation of intermediate states of mixing and their relation to chemical conversion. Using the population balance approach as a basis, a number of authors have introduced a variety of mixing history parameters, each of which allow the degree of micromixing (i.e. mixing between fluid elements of different ages and equal residual life-time) to be varied for a given residence time distribution. Particular attention has been focussed on the ideal CSTR residence time distribution[3-161. One of the inherent dithculties of these micromixing studies is relating any particular model parameter to the flow patterns existing in the reactor vessel of interest. Single parameter models (e.g. [5, 11, 13, 151) have the advantage of computational simplicity, in that the values of the parameter can be readily determined from FOLLOWING
input-output measurements of conversion, (or in certain circumstances, tracer concentrations), but can be limited in their predictive powers because of their ill-defined relationship to the actual complexities of the fluid motion. Multiparameter models (e.g.[6, 17,181) may purport to describe the vessel fluid mechanics more realistically, but considerable ambiguity can arise when attempting to evaluate the parameters from residence time or conversion measurements. It should be pointed out, however, that in some instances, mixing history models have been proposed not so much as to describe an actual reactor, but more to establish whether micromixing is of importance in determining the outcome of various reaction types. At the present time, the most successful approach towards relating a mixing history parameter to the flow conditions within a reactor appears to be that involving the use of a random coalescence model, simulated by Monte Carlo [7,12,19-221, analytic or numerical procedures [3,4,8-10,13,14]. In this model, the reactor contents are depicted as consisting of a large constant number of equally sized elements (or
tPresent address: Armfield Engineering Ltd, Ringwood, Hampshire, England. SOn leave from Department of Chemical Engineering, Cairo University, Giza, Cairo, U.A.R.
1497
C. R. TRELEAVEN
droplets for the dispersed phase of liquid-liquid systems). Each element consists of a large number of molecules, and it is a basic postulate of the model that the contents of each of these elements are well mixed. The mechanism of mixing consists of coalescing two elements at a time, whereupon their intensive properties are assumed to be averaged instantaneously. In order to maintain a constant inventory of elements in the vessel, re-dispersion into two identical elements takes place after every coalescence. Whatever method is used for selecting elements for this coalescence-redispersion process, and its rate of occurrence, one restriction must always be imposed. Each coalescence implies the transfer of molecular material between two elements, and since each element is considered well-mixed, each re-dispersed element contains molecules drawn from both coalescing elements. Whatever the distribution of ages, chemical species, temperatures etc. within each element, at the vessel outlet every molecule must have a residual life-time of zero. (Otherwise a spontan: eous de-mixing of molecules of non-zero residual life-times would have to occur; a process which is contrary to the Second Law of Thermodynamics). In order to ensure this condition, elements from the feedstreams associate only with those elements already in the reactor having the same residual life-time and, furthermore, each pair-wise coalescence and re-dispersion process is restricted to elements with the same residual life-time. The flow of an element into or out of the reactor and the coalescence process are assumed to have no effect on the probability of entrance or exit of other elements-this probability being determined by the known residence time functions of the system. Provided these conditions are fulfilled, the way and the rate at which elements are selected for coalescence is a matter of choice. The basis of the random coalescence model, however, is that elements are selected at random from populations of elements having the same residual lifetime. The residual life-time distribution is of course found from the residence time distribu-
and A. H. TOBGY
tion, according to the following relationship [21:
W) = + [ 1 -F(h)].
well known
(1)
Although it is possible, in principle, to vary the number of pairs, h, coalescing in a given time interval, T, with such parameters as residual life-time, age, etc. it has proved more convenient to use a single value for the rate of coalescence for the whole reactor population. This implies that all elements in the vessel have an equal chance of coalescing in the given time interval T, i.e. the coalescence is a Poisson process. In this time interval between successive coalescences, elements behave as independent batch reactors which are assumed to have no effect on the flow characteristics. Flow through the vessel is simulated by entering and exiting a number of elements, I, over the time interval T; the number being proportional to the flow rate. Thus, since the average residence time Q-is the ratio of the hold-up of elements in the vessel, n, to the throughput, it follows that
7=& The destination of these entering elements, in terms of their assignment to the residual lifetime distribution, is governed by the residence time functions for the system. The parameter of the random coalescence model that characterises micromixing is I, the average number of coalescences undergone by a single element during its stay in the vessel. As each element has an equal chance of coalescence, the average rate of a single element is Z/r. Within the vessel, therefore, the number of pairs of elements coalescing per unit time is nZ/2r. The number of pairs of elements, h, coalescing in the time interval T is given by
The independent parameter of this micromixing model is I, and it is the task of the modelling procedure to assign a value for Z to a reactor of
1498
Monte Carlo methods of simulating micromixing
known residence time distribution that is of use in correlating experimentally determined conversion results and/or to evaluate the effect of various degrees of micromixedness on the reactor performance. Evangelista ef al.11 31 have shown that the parameter expressing the rate of coalescence I can be related to the properties of a homogeneous, isotropic turbulence field, in a way that correlates experimentally determined data for a stirred vessel. Similarly, Kattan and Adler [201, Harris and Srivastava[211 and Rao and Dunn[221, have shown that their coalescence parameters account for the experimental mixing and reaction data obtained by Vassitalos and Toor[23] for a plug flow reactor having unpremixed feed streams. There is a prima fucia case, therefore, for investigating whether the random coalescence model of micromixing can provide a useful approach to reactors having arbitrary residence time distribution, particularly those in which turbulent flow conditions exist. Kattan and Adler [ 121 have applied this model to homogeneous reactors that are characterised by a single residence time distribution, although no correlation with experimental data has yet been attempted. In an earlier paper[24], the present authors considered the extreme mixing states of complete segregation and maximum mixedness as applied to reactors having two un-premixed reactant feed streams with different residence time distributions. This paper is concerned with applying the random coalescence model to this type of reactor arrangement, using a Monte Carlo simulation procedure, in order to calculate conversions for intermediate states of mixing. 2. THE
RANDOM
COALESCENCE
MODEL
response of the reactor vessel is required and on the computer facilities available. The basic features of both simulation methods are as follows:(a) Stochastic method
simulation
In essence, this method involves adapting the work of Kattan and Adler [ 121 to reactors having un-premixed feed streams whose residence time distributions differ. The model and associated computer programme are constructed in the following order: (i) The two feed streams are divided into two large numbers of discrete elements m, and m2, the ratio of these two populations being in proportion to their flow rates Q 1 and Q 2. Thus: m,:m,:m
= QI:Qz:Q.
(4)
(ii) The theoretical or experimentally determined residence time functions fi(t) and f&t), are re-formulated on a discrete basis, by segmenting the functions over a large number, K, of small time intervals At. Thus, for the kth residence time interval for thejth feedstream, =f,(k . At). At = fjle. At = F,[(k+ l)At]
A(t
-Fj(k.
At);
j=
1,2
k=1,2,...K
(5)
where k . At corresponds to residence time t,. Elements of the two feed populations are then assigned to each residence time interval in numerical proportion to the frequency for each interval. For example, the number of elements assigned to the kth residence time interval, is given by: mjl,
Although the essential feature of the model is the stochastic simulation of the micromixing process, it is possible to simulate the flows in and out of the reactor vessel either on a stochastic or deterministic basis. Both methods are presented here, where it will be seen that each offers certain advantages and disadvantages: which is preferable depends on whether the dynamic
jlow and micromixing
,
mjl, = l?‘ljfJk. At;
j= k=
1,2 1,2,.
. . K.
(6)
In view of the necessity of maintaining an integer number Of elements in each interval, Iltl]r and mk are always taken as the nearest integer number given by these equations.
1499
C. R. TRELEAVEN
(iii) The discrete versions of the residual life-time frequency functions of the contents of the vessel, $,(A) and &(A) are formed from the residence time functions by virtue of the relationship given in Eq. (1) and the considerations of more than one feed stream as outlined in Ref. [24]. Thus:
=i [l -Fj(A)l; j=
1,2
and discretising A in the same way as t:
Jlj@tc) =
$[l-if;x+At]=i[l-Fj,l; i=l
3
j= 1,2 k= 1,2,...K
(7)
where Akis any specified residual life-time, and
The reactor contents are divided into two populations of a large number of equally sized elements, n, and n2, such that the ratio n, : n2 is equal to the ratio of the hold-ups from each stream, V, : V,. In an earlier paper [24], it is shown that: I’j = Q~T~z QjJr
(1sFj(t))dt;
j=
1,2.
(8)
Each of the n elements is given an identification number and assigned to a residual life time interval, k, the number in each being proportional to the residual life-time frequency of that interval. Thus njk = nj$j(A,) AA;
j= 1,2. k=l,2,...K,
(9)
and 2 n(A,) = jzl n,, ;
k= 1,2,. . . K.
(10)
(iv) Flow into the vessel is simulated by drawing elements at random from the two feed
and A. H. TOBGY
populations, ml and m,, in the ratio Q I : Q2, and according to their residence time identification numbers, entering them into the intervals containing those elements in the reactor having the same residual life-times. (The residual lifetime of an element at the vessel inlet is equal to its residence time.) The average numbers of elements entering per unit time from each stream are H~/T~and n2/r2 respectively. Any particular flow simulation step consists of entering 1 elements from the feed population m. If the ratio of the flow rates Q 1: Q 2 involves a non-integer number, 1 (which must be an integer) has to be varied over a series of flow steps in order to maintain the correct average flow ratio over the whole series. The variation in 1 should be such that in the ratio of the number of elements entering from each stream does not deviate from the correct ratio by more than 10 per cent. Within the reactor, the same number of elements as enter each interval are selected at random from those intervals and their residual life-times decreased by one interval. This process of random selection and movement to a lower life-time interval continues until the appropriate number of elements leave the first residual life-time interval and exit from the vessel. In order to maintain a constant inventory of both populations representing V, and V,, elements only associate with those that entered from the same feed stream when moving from one life-time interval to the next. This does not imply that the elements from each stream do not associate with each other from the mixing point of view, as the next section makes clear. The flow simulation is depicted in Fig. 1. (v) Micromixing is simulated by selecting pairs of elements at random from the same (randomly selected) residual life-time interval, whatever their origin. Coalescence of each pair of elements results in an instantaneous averaging of the concentration of each species (and temperature, if non-isothermal cases are being considered), so that re-dispersion produces two identical elements. The number of pairs of elements, h, chosen for coalescence per micromixing step is obtained from Eq. (3) as applied to a two
1500
Monte Carlo methods of simulating micromixing
Feed I
Feed 2 m2 *
Z
I 1
-f,(t)
----I
-Jr,(X,-)
-
---Jr,(X)-=
v
fz(i)---)
Exit m Fig. 1. Diagrammatic representation of stochastic flow simulation. (Horizontal arrows across centre line represent micromixing only.)
feed stream situation, i.e.
(11) where I is the average number of coalescences experienced by a single element during its sojourn in the vessel. (If h is not an integer number, then a compensated, nearest integer number of micromixing simulations are performed in each micromixing step, the compensation being such as to obtain the correct average h over a series of micromixing steps.) Z is therefore the independent variable characterising micromixing in the vessel, and it follows that when I equals infinity, the condition of maximum mixedness between material of different ages and species, Model I of the authors’ previous paper 1241, should apply. (vi) The flow and micromixing simulation steps are repeated every T units of time, where, from Eq. (2):
T=k n,+n,’
(12)
In this time interval T, each element behaves as an independent batch reactor, and integration
of the appropriate kinetics is carried out to give the concentrations in each element before the next flow and micromixing simulation steps. However, if Z&z1+ 4) is very small, considerable computation time can be saved by performing the integration of the batch kinetics for each element after K flOW and miCrOkXing simulation steps (i.e. K can be greater than 1). The value of KT is selected by trial and error to give a negligible effect on the final averaged conversions of elements leaving the reactor. This simulation model allows the dynamic response of a reactor vessel to be calculated, by arranging for the initial concentrations of the elements in the two reactant feed streams to be a function of time. (b) Deterministic flow, stochastic micromixing simulation method Experimentally determined or theoretical residence time distributions are usually presented on a deterministic basis (although fluctuations about the mean distributions may well have significance in relation to micromixing, as will be discussed elsewhere). It is thus possible to simulate the flow process in a non-random manner while still retaining the features of the
1501 Cl33 VoL 27 No. 8-B
C. R. TRELEAVEN Feed a
I
Feed
and A. H. TOBGY
where N, is the reactor population (i.e. hold-up) originating from the jth feed stream. For the two streams, where j = 1 and j = 2 respectively,
2 _
Ml
(1% The function gj(al, hK)A~Ah represents the fraction of the hold-up population originating from the jth feed stream which has the earliest age (0-Acr) and longest residual life-time (XK_-I-AK). From this, it follows that: r
gj(al,A,-l)AaAA = [~~(~~-~)-~~(~~)lA~; j=
1,2.
(16)
In practice, it is most convenient to arrange for the time intervals over which the various frequency functions are discretised to equal each other, i.e. At=Aa=AA. Fig. 2. Diagrammatic representation of deterministic flow, stochastic micromixing simulation method.
random coalescence model described above. The following steps are involved (see Fig. 2). (i) As before, the normalised residence time frequency functions are arranged in discretised forms. While the total areas $ fikAt and i fileAt
(iii) The total number of elements with the longest residual life-time, [NII(II(AK)+ N2&(AK)] X AA, each of which has an identification number, forms the population for the tirst micromixing and reaction simulation step. In the Kth time interval A?, the number of pairs of elements due for coalescence, H(A,), is given by:-
k=l
k=l
the same, where K is the index of the longest residence time of either of the two feed streams, the number of elements within each total area, M, and M,, are given by are
H(h) =
I. N(A,)At 2T
M represents the number of elements entering the vessel in a time interval At. (ii) All the elements with the longest residence time, M,f,,At and M&At are entered into the last residual life-time interval K. Thus:
MjK = Mj e&KAt = Njgj(al, A,)Acx . AX; j=
1,2.
(14)
+ N,$&)IAA. (19)
(13)
Q,’
(18)
where N(A,) = N~K+ N,, = [N,JI,(U
Ml Q, -=Mz
(17)
As in the method (a) outlined previously, this micromixing is simulated by random selection, and is considered to take place instantaneously and to produce identical pairs of elements. Each of the elements is then treated as a batch reactor, and chemical conversions calculated for a time interval At . (iv) At the end of the Kth residual life-time interval, all of the elements in the Kth residual life-time interval are of age Aa, and therefore
1502
Monte Carlo methods of simulating micromixing
have residual life-times of (K- 1)AA. These elements are therefore moved into the (K - I)th residual life-time interval, where they are joined by all the feed elements entering from the (K - 1)th residence time interval, [M,f,((K - 1) x At) + M,f,((K - 1)At)lAt. The total population in the (K- I)th residual life time interval N(X,_,), is then
After the final micromixing simulation and conversion calculations, all the M elements then exit from the reactor, and the average composition of the exit stream is calculated from suitably selected samples of these leaving elements. The averaging procedures used for calculating the exit conversions are discussed below. 3. SIMULATIONS
[N,JI,@,-d
+ N&&+-1)IAA
In order to test the validity of the simulation procedures and associated computer programmes, a number of initial trials were undertaken for a residence time distribution of one CSTR and a simple second order, isothermal reaction following the stoichiometry:
= N . $FI(b--l)A~
where #(A) = &
*+,(A) +$&
*&(A).
A+B+
(20)
P.
For these trials: The micromixing simulation, and the subsequent batch kinetic calculations, are then carried out for the N(A,_,) elements in a similar way to that used for the N(A,) elements. (v) The flow, micromixing and reaction simulation procedure is repeated for each decreasing residual life-time index, K - 2, K - 3, etc. In general, for the kth sequence of simulation steps, the following relationships apply: steps, the following relationships apply: Mjk = MjJkAt;
j= 1,2. k=K,K-I
Njk = Nj$j(Ak)AA; N(A,J=Nlk+Nzk; Ho,
k
j= 1,2. k=K,K-l,...,
27
(21)
1,
(22)
9K-1
(23) ,*-*,
1.
= : e+.
Both premixed and un-premixed reactant feed stream situations were studied: in the former case, each stream contained both A and B, while in the latter, A alone was present in stream 1 and B alone in stream 2. The results showed good agreement with those of other authors [7, 13,191, with the exception of two apparent misprints in certain of these reported resultst. Computations by both methods of simulating the random coalescence model were then performed, for the same reaction and two pairs of residence time distributions as were used in an earlier paper[241; i.e. A+B+
(24)
P
with
(The computer programme Eqs. (2 l)-(24) incorporated errors caused by rounding to integers.) When the first reached, the total population + N&&WA
1,
k=K,K--l,...,l,
)=z’N(Ak)*At;k=K
[N,MQ
,...r
. h(t) =fi(t)
used for calculating compensation for fractional numbers life time interval is is
= W, + M-3 = M.
(25)
rA =
k,.CACp
tin the work of Spielman and Levenspiel[7], the values of their dimensionless reaction parameter kc,,7 presented in their Fig. 3 should be multiplied by 100, whilst in Fig. 2 of the paper of Evangelista et al. [13], the curve labelled X07= 1 would appear in fact to represent the results for /rCoT= 0.5.
1503
C. R. TRELEAVEN
For Case 1: Stream 1:
and A. H. TOBGY
tisation times and 8 were chosen, are discussed in a subsequent section of this paper.
t fi(t) = 16 eetj4,
4. RESULTS
Stream 2:
_MO=~e.
p
For Case 2, the distributions Stream 1:
-t/2
.
were 0 G t d O-75 set
F,(t) = 10-9t; = l-exp-
[
-t-o.75 7.25
1 *9
0.75 < t < 59.25 set, stream 2:
F&) = 10-9; = l-exp-
0 G t s 4.5 set [
t-4.5 -. 7_5
I
,
4-5 c t < 59.25 sec. . Reactant A was fed in stream 1 at a concenttation of cAo while B was fed in stream 2 at cBo, and in all cases, the volumetric feed rates Q, and Q2 were equal. The progress of the reaction for any given element was calculated from the integration of the batch kinetic expression, namely, for CA(t) # cg(t): cAt)[c.4(t) - cew1 “(‘+ @= CA(t)- cc(t) exp [(k,fI) . (cc(t) -CA(t))]
Although Monte Carlo simulation techniques are attractive in solving problems not amenable to analytical solutions, they suffer from the disadvantage that the results so produced are only known in a probability sense. For the reactor simulations carried out in this work, the factors affecting what might be termed stochastic convergence will be discussed below, but it is first worth discussing the averaged results. Using the stochastic feed and micromixing method, it was possible to observe the dynamic response (in effect, the start-up) of the reactor, and steady state conditions were reached after some 3-4N elements had passed through the vessel. Results for the Case 1 residence time distributions, (2 and 6 CSTRs respectively for stream 1 and stream 2), are shown in Figs. 3-5 for values of the stoichiometric feed ratio, p, equal to O-25, 1.0 and 4.0, and for I = 4 and 50. For the results reported here, k,.iio~ was maintained as unity, where co is the initial mean concentration of the limitng reactant. Thus:-
co = Eaoif ?A0 < cBo,and 4 = cBoif cao 3 cBo, where
(26) and for CA(t)= cg(t):
(29) and (30)
k-@c‘4(01
(27)
/3 is defined as ~~o/&o, and the conversion X, as &,/~o. For a sample size of L elements,
CB(t+ e) - CB(t)= CA(t+ e) - CA(t).
(28)
(31)
For the stochastic flow and micromixing simulation method, 8 = KT, while for the deterministic flow but stochastic micromixing method, 8 = At = AX. The results were calculated on the basis of m=n= 600 element (first method), and M = 200 (second method). The reason why these particular values and also those for the discre-
Computation time for each of the quadruple dynamic simulation tests took approximately 14 minutes at the Atlas Computer Laboratory (Berkshire, England). Each point on the fluctuating traces shown in Figs. 3-5 represents an average conversion, calculated from four replications of the whole conversion-time run and lumping the con-
c‘4(t+ 0) = cA(t)l[ 1+ and
1504
Monte Carlo methods of simulating micromixing
versions of every 10 consecutive leaving elements together; viz. the average conversion of 40 leaving elements. The time-averaged smoothed curves were obtained as follows: (i) lumping together the results of 5 points from the
fluctuating traces, such that each point on the smoothed curve represents the average conversion of 200 leaving elements, distance r/ 12 apart, (ii) fitting a smooth curve through these equidistance results using a least squares, third
Simulation of start-up behaviour of Case 1 unpremixed feed reactor, as a function of I and /3. [k&r = 1; Ar=0.35sec;m=n=600;1=K=2.1
0 640 X .E E 2 s 0
0.460
0320
0.160
Dimensbnless
time,
t/r
Fig. 3. Simulation for p = O-25.
Dimensonless
time,
t/r
Fig. 4. Simulation for p = 1.O.
C. R. TRELEAVEN
and A. H. TOBGY
oaoo -
0640X .5
0.480 -
f u"
0.320-
O-160-
Oimensionless time,
t/T
Fig. 5. Simulation for p = 4.0.
degree,
seven point, moveable
strip technique
M. It may be noticed that at t > 37, where it might be expected that a stationary random conversion to be approached asymptotically, a certain amount of periodicity exists. This is probably due to a defect in the random number generator used. The second method of simulation, i.e. deterministic flow and stochastic micromixing, while losing the advantage of giving the dynamic response of the reactor vessel, allowed considerable saving in computer time, in that each simulation run of six replicates, took only approximately 5 min on an I.C.L. 4/50 computer. Averaged results for the two sets of residence time distributions pairs for various values of k&,7 and the stoichiometric feed ratio p are presented in Figs. 6-l 1. Each point represents the average of between 4 and 12 replications, the number of replications being dependent on Z (q.u.). Although only a few of the results of the first simulation method are presented here (Figs. 3-Q the averaged final steady state conversions were within the same confidence limits as those of the second method. In effect, therefore, Figs. 6- 11 represent
the averaged outlet conversions calculated by both methods. As expected, all results for increasing values of the micromixing parameter I, asymptote to the maximum species and age mixedness model (Model I) postulated in the authors’ earlier paper[24]. The greatest effect of variations in the micromixing on conversion is seen to be in the region of Z = l-20, whatever the value of the stoichiometric ratio ~3 or the parameter k&r. 5. STATISTICAL
CONSIDERATIONS
In view of the inherently stochastic nature of Monte Carlo techniques, it is necessary to consider the factors that alfect the variations in the reactor output conversions, and the degree of confidence that can be placed on averaged results. These factors fall into two classes; those characterising the operation of the reactor, such as the residence time distributions, reaction velocity constant, micromixing parameter etc., and those relating to the simulation techniques themselves, namely, the number of elements chosen to represent the feed and reactor populations, time intervals selected for discretising the
1506
Monte Carlo methods of simulating micromixing Steady state conversions for unpremixed feed reactors, as a function of I, p and k&7. (Values at I = m, shown as horizontal dotted lines, are taken from a previous paper[24] for Model I maximum mixedness.) 10 k, F,, T -50 08-
O-6 -
x
x
I
5
o-4 Q2 02
-I 01
Fig. 6. Case 1 residence time functions, p = 0.25.
I
Fig. 7. Case 1 residence time functions, p = 1.0.
2I
I
I11111l 5
20 1
IO
50 I1111111
IO0
5
10
20
100
I
,
m
Fig. 10. Case 2 residence time functions, /3 = 1-O.
I
Ol1
2
Fig. 9. Case 2 residence time functions, B = 0.25.
4 I
--
0.1
k&r-JO/
ik
I
Fig. 11.Case 2 residence time functions, /3 = 4-O.
Fig. 8. Case 1residence time functions, fl= 4-O. 1507
m
C. R. TRELEAVEN
and A. H. TOBGY
various continuous residence time functions, and the integration time of the batch kinetics. (i) Reactor-parameters For all the simulation runs reported in this work, calculations of the conversion variances, var X, of the elements leaving the reactor were carried out. For a total of M elements, var X is defined in the usual manner:
2 (Xi-X)" varx=‘=lM_l
(32)
X
$
>
where
X=
5 Xi/M.
i=l
(33)
Although all the results of the statistical calculations are too lengthy to tabulate heret, a number of conclusions were drawn from these, as selected examples given below will indicate. The most important, from the point of view of simulation procedure, was that the variances calculated from the steady state results of the stochastic flow and micromixing technique (e.g. those for t > 37 in Figs. 3-5) showed an identical pattern to those calculated from replications of the deterministic flow and stochastic micromixing method, for any values of the parameters fi and fi, p, k,.?,,rand I. Figure 12 shows the conversion variances for the case of 2 and 6 CSTR residence time distributions, for various values of /I and 1. The confidence limits for these variances were obtained by plotting the variances of samples of 200 elements leaving the reactor on log-normal probability co-ordinates (e.g. Fig. 13) which gave approximate straight lines of almost identical gradients for given values of Z and the functions fi and&. These plots suggested that although the variance about the mean conversion was a function of p and k&r, the coefficient of variation (defined as 100 m/X) should be virtually independent of these parameters. This was in fact found to be the case as the selected examples tAvailable on request.
Mean
confiina3
Fig. 12. Conversion variances for Case 1 unpremixed feed reactor, for /C&T= 1-O.
presented in Fig. 14 indicate. Apart from the strong dependence of the coefficient of variation on I, the residence time distribution functions fi and fi were also found to affect the coefficient. An example of this interaction is shown in Fig. 15, for a single, un-premixed feed CSTR and also for the Case 2 distributions referred to in Section 3 above, namely, PFR followed by a CSTR. These studies indicate that, for a given reactor configuration, the region where the lowest degree of confidence can be placed on average conversion results, based on a given sample size L of the leaving elements, is that in which a low value of the micromixing parameter Z prevails. In view of the sensitivity of the mean conversion, and its variance to the parameters characterising any chemical reactor, it is necessary to consider how the Monte Carlo simulation parameters are chosen to give results within the
1508
Monte Carlo methods of simutatmg mlcromlxmg
b
k&T
I.0
I
I.0 I.0 0.25
2 IO I
I
I 002
'oa5
I
I
var
A V q
IIIII 006
0.04
0.03
0
006
O-I
x
Fig. 13. Log-normal plot of conversion variance for Case 1 unpremixed feed reactor, for I = 4.
desired confidence limits, while at the same time not allowing the computation time to rise to uneconomic proportions.
200 1
6
E .Q
95% k,CObr
p
Mean value
0.25 I.0 I.0
0 0 v
2
2
t
0”
I
Fig. 15. Coefficient of variation of steady state conversions for two different unpremixed feed reactors. (k&g = 1; /3= 1.)
I I IO
2
5
(ii) Simulation parameters From the preliminary trials on the system of a single CSTR, the stochastic simulation of flow and micromixing showed that the selection of m, n, At, I and KT may cause bias in the conversion results. Too large a value of At results in loss of resolution in the discretisation of the residence time functions, while small values of m and/or 12may cause truncation of these functions. Similarly, a large number of flow and micromixing simulations before integration of the batch kinetics results in incorrect values of the conversion. Although m and n are functionally related by
_I\
confidence litits ----
IO
I 20
I1111111 50
loo
n=r.m
(34)
I Fig. 14. Coefficient of variation of steady state conversions for Case 1 unpremixed feed reactor.
practice, it may be useful (for computer timesaving) to ignore this relation, and to select
in
1509
C. R. TRELEAVEN
and m independently. This arises because if m is selected large enough to reduce the discretisation error of the residence time frequency functions and r % 1, the resulting value of n would be prohibitively large. On the other hand, if T Q 1, n would be too small to adequately represent the reactor frequency functions. In order to avoid these difficulties, the following points were found to be important in selecting the size of m/n. The flow of n particles takes place in r time units, and if m is chosen to be much larger than n, then only a few elements of the m feed population can be selected for entry to the vessel during each time interval over which the dynamic response is being studied. Thus, selection from the residence time frequency functions may be biased, leading to large errors in the dynamic simulation results. The lower limit of m is, of course, dictated by the adequacy of the discretisation of the frequency functions. In the present work, the range O-5 B n/m 3 2-O was found to give reasonably accurate results. One method of reducing m without truncating the residence time functions to an unacceptable level that proved useful, was to determine the number of particles in each feed subpopulatiOn, mjk, from the relationship; n
mjk = mjFj(hk)
k-l z mji;
-
j=
i=l
k=
1,2. 1,2,. . .K
(35)
rather than from Eq. (6). Similarly, the reactor contents subpopulations are determined from: k-l %k
=
%*dhk)
-
2
j=
?ii ;
i=l
1,2.
k= 1,2,. . .K
(36)
rather than from Eq. (9). In the Monte Carlo method of simulating the flow deterministically and the micromixing stochastically, the same sources of error as mentioned above mav arise. However, for this method, N and M must necessarily obey the relationship: N=x
MT
(37)
and A. H. TOBGY
The selection of At is therefore more crucial in determining the bias in the results as well as the computation time. For systems strongly affected by the interaction of mixing and reaction, At may have to be so small that the required length of the computation time may suggest that the other Monte Carlo method as economic. It should be noted that equations analogous to (35) and (36) may advantageously be used to determine MJk and Njk, rather than equations (21) and (22). For the simulations carried out in this work, i.e. for second order reactions with O-25 < p G 4 and 0.1 s k&r c 1000, 100 particles proved to be near the lowest tolerable values for ml, m2, nl, n2, M, and M2, while At had to be in the range (r/40) G At < (r/10), in order to prevent any considerable bias in the conversion results. So far in this discussion, attention has been focussed on the factors that may cause systematic biases in the results. Of equal importance is the selection of sample sizes from the reactor output population, in a manner that minimises the random errors inherent in any averaging procedure. Unfortunately, the convergence of a sample mean to that of an infinite population is a slow process, in view of the fact that the standard deviation of the mean value of a sample of size L is inversely proportional to fi. Thus, to decrease the random error by a factor of 10, L must be increased 100 fold. However, the Central Limit Theorem allows access to the confidence limits that can be placed on the mean conversions obtained by the two Monte Carlo methods used in this work, in terms of the sample size L[27]. Thus, the 95 per cent confidence limits of the mean conversion x of samples of size L drawn from a population of mean conversion and variance var X, are
IX-(/
s
where to.oz5,~-1 is
$_
to.o!25,r,-1
*.
(38)
the parameter used in “Student’s” t-test, at the O-025 level with L- 1 degrees of freedom. Different methods may be employed to obtain
1510
Monte Carlo methods of simulating micromixing
the appropriate sample size L that yields mean conversion results within a required degree of confidence. For steady state conversions, the stochastic flow and micromixing Monte Carlo method can be employed in such a way that the length of the simulation run is sufficiently large to give an adequate value of L. The deterministic flow method will yield sufhciently high values of L if a large value of M is selected and/or a number of replicate runs are made. For dynamic tests, where only the first Monte Carlo simulation method is appropriate, either n can be increased or a number of replicate runs can be made. 6. PRJZDICTION OF REACTOR PERFORMANCE
The present work has extended the random coalescence model to allow for the effects of micromixing on chemical reactions to be studied for vessels having two feed streams of differing residence time distributions. From the practical point of view, the important task is to devise some means, either theoretical or by simple experiments, to evaluate Z for a given reactor configuration, and to discover whether this micromixing parameter does in any sense represent the fluid mixing characteristics of the reactor. Implicit in the present work is the assumption that Z is a hydrodynamic parameter, and that providing the reaction itself does not a.iIect the flow patterns within a reactor, a value for Z in given hydrodynamic circumstances should be independent of parameters characterising the reaction, such as /3 or k,. If this assumption is obeyed in practice, then a value for Z can be deduced from one reaction experiment, and thereafter used to simulate the reactor behaviour for any other reaction of interest to be carried out under the same hydrodynamic conditions. This is essentially the same approach as that of Kattan and Adler 1201, although the present work indicates that it is not necessary to select an infinitely fast reaction to determine I. There is some evidence that Z can be determined from experimental measurements of the
fluctuations in concentration of an injected, nonreacting tracer. For example, the work of Evangelista et aZ.[131 suggests that the micromixing parameter that characterises random coalescence model is related to the properties of a homogeneous, isotropic turbulence field, in such a way that scale-up of stirred reactor vessels is possible. Their parameter (denoted here as y to avoid confusion with the stoichiometric ratio p) is related to Z used in the present work by: (39) Similarly, Toor[26] and Kattan and Adler[20] have proposed methods by which the micromixing characteristic can be determined from the statistical properties of tracer concentration measurements in a turbulent field. In a subsequent paper, the single chemical reaction method of determining I, and its efficacy in predicting conversions for different reactions for an experimental, two feed stream, concentric jet reactor arrangement will be reported, while comments will be offered elsewhere on the relation between measurements of concentration fluctuations at the reactor outlet and the micromixing parameter. The basis for these latter comments is the observation that the nature of the output fluctuations found in the simulation work (e.g. those shown in Figs. 3-5 in this paper), appear to bear the same relationship to Z as those to be reported in the experimental work.
NOTATION
chemical species, reactants concentration, moles vol.-’ continuous stirred tank reactor residence time frequency function, time-’ F cumulative residence time distribution &z, A) joint age and residual lifetime frequency function, time-’ h number of pairs of elements chosen for coalescence per micromixing step A, B
1511
c CSTR f
C. R. TRELEAVEN
and A. H. TOBGY
H
number of pairs of elements chosen for coalescence during At as defined in I micromixing parameter, text k arbitrary identification number k second order reaction velocity constant, vol mole-’ time-’ K identification number of the interval of longest residence time 1 average number of elements entering reactor in one flow simulation step L sample size m number of elements representing feed population, time-’ M number of elements representing feed population number of elements representing reacn, N tor population P chemical species, product PFR plug flow reactor Q flow rate, vol time-’ rate of consumption by reaction of rA reactant A, mole vol-’ time-’ t residence time, time T time interval between successive flow and micromixing steps, time variance Var
V X Greek
Reactor volume, vol estimated conversion
symbols
age, time stoichiometric ratio of feed streams Y average rate of coalescence of each element (definea as /3 in Ref.[ 131) number of flow and micromixing steps K performed before integration of batch kinetics residual life-time, time ; residual life-time frequency function, time-’ residual life-time distri0 cumulative bution average residence time, time e’ time, time 5 population conversion SuJixes 0 initial value 192 stream 1,2. (does not apply to (Y,h or t) species’A, B, P -4,&P i ith interval or element jth feed stream kth or Kth interval k,;
REFERENCES DANCKWERTS P. V., Chem. Engng Sci. 1958 8 93. ZWIETERING Th. N., Chem. Engng Sci. 1959 111. 131HARDA M., ARIMA K., EGUCHI W. and NAGATA S., Mem. Fat. Engng. Kyoto Univ. 1962 24 43 1 [41 CURL R. L., A.Z.Ch.E. Jl1963 9 175. PI NG D. Y. C. and RIPPIN D. W. T., Third European Symposium on Chemical Reaction Engineering, Amsterdam, Sept. 1964 p. 161. Pergamon Press, Oxford 1965.
[email protected].,Paper 10.5a,A.Z.Ch.E.-Z.Chem.E. JointMtg. London 1965. ;; SPIELMAN L. A. and LEVENSPIEL O., Chem. Engng Sci. 1965 20 247. CHEN MING-HENG, CHEN LIANG-HENG and YUAN WEI-KANG, Scientia Sinica 1966 15 123. ;z; SHAIN S.,A.Z.Ch.E. Jl1966 12 806. [lOI CURL R. i., Chem. Engng Sci: 1967 22 353. [ill WEINSTEIN, H. and ADLER R. J. Chem. Engng Sci. 1967 22 65. r121KATTAN A. and ADLER R. J., Paper presented at 62nd Nat.A.Z.Ch.E. Mtg, Salt Lake City, Utah. [I31 EVANGELISTA J. J.. KATZ S. and SHINNAR R.. A.Z.Ch.E. JI 1969 15 843. [I41 HULBERT H. M. and’AILYAMA T., Znd. Engng &em. Fundls 1969 8 3 19. VILLERMAUX J. and ZOULALIAN A., Chem. Engng Sci. 1969 24 1513. t::; METHOT J. C. and ROY P.-H., Chem. Engng Sci. 197126 569. [171 ADLER R. J., LONG W. M., ROOZE J. and WEINSTEIN H., Proc. 6th World Pet. Congress, Sect. VII, Paper 19, Frankfurt-Maine 1963. 1181NISHIMURA Y. and MATSUBARA M., Chem. Engng Sci. 1970 25 1785. KATTAN A., Ph.D. Thesis, Case Instituteof Technology, 1967. $3 KATTAN A. and ADLER R. J., A.Z.Ch.E. Jl1967 13 580. Pll HARRIS I. J. and SRIVASTAVA R. D., Can. J. Chem. Engng 1968 46 66. WI RAO D. P. and DUNN I. J., Chem. Engng Sci. 1970 25 1275.
‘,:;
1512
[23] [24] [25] [26] [27]
VASSILATOS G. andTOORH. L.,A.Z.Ch.E.JI 1965 11666. TRELEAVEN C. R. and TOBGY A. H., Chem. Engng Sci. 197126 1259. HERSHEY H. C., ZAKIN J. L. and SIMHA R., Ind. Engng Chem. Fundls 1967 6 413. TOOR H. L.,A.Z.Ch.E. Jll962 8 70. HALD A., Srarisrical Theory with Engineering Applications. Wiley, New York 1967. R&me-Les auteurs dtendent le modele de coalescence irregulibre, simule par les techniques de Monte Carlo, aux reacteurs dans lesquels les deux courants distincts d’alimentation de corps reagissant ont des distributions differentes du temps de residence. Le parametre qui caracterise le micromelange est le nombre moyen de collisions auxquelles chaque element fluide est soumis au tours de son sejour dans le reacteur. Les auteurs ont fait varier le parametre pour le cas d’une reaction simple, isothermique, homogene et de deuxitme ordre, avec deux ensembles arbitraires des distributions des temps de residence; ils montrent qu’un taux de collision infini correspond a I’un des modtles du melange maximal present& anttrieurement par les auteurs. Une discussion suit sur les methodes d’apres lesquelles on peut mesurer le parametre du micromelange et la valeur utilisie ulterieurement pour pr6dire le r&hat de n’importe quelle reaction donnee. Les auteurs commentent tgalement la fiabilitd des techniques de Monte Carlo par rapport aux parametres qui atfectent les rendements moyens et les possibilitts de reduire les temps de calculs. Zusanunenfassung-Das regelloser Model1 Koaleszenz, nachgebildet mit Hilfe von Monte Carlo Methoden, wurde auf den Fall des Reaktors erweitert, bei welchem zwei separate Zufuhrstriime der Reaktionsteilnehmer verschiedene Verweilzeitverteilungen aufweisen. Der die Mikromischung kennzeichnende Parameter stellt die durchschnittliche Anzahl von Zusammenstiissen dar, die jedes Fliissigkeitselement im Laufe seines Aufenthaltes im Reaktionskessel erfiihrt. Dieser Parameter wurde fib eine enfache, isothermie, homogene Reaktion zweiter Ordmmg fur zwei willkiirlich gewtilte Gruppen von Verweilzeitverteilungen variiert, und es konnte gezeigt werden, dass eine unendliche Stossgeschwin einem der durch die Autoren in einem frilheren Artikel dargelegten Modelldigkeit maximaler Mischung entspricht. Es werden Methoden erortert mittels welcher die Mikromischparameter gemessen und der Wert hemach zur Voraussage des Ergebnisses einer bestimmten Reaktion verwendet werden kann. Im iibrigen wet-den Betrachtungen ilber die Verliisslichkeit der Monte Carlo Methoden hinsichtlich der Parameter, die den sich ergebenden mittlem Unisatz beeinflussen, angestellt, sowie ilber Miiglichkeiten zur Verringenmg der Rechenzeiten.
1513