Residence Time Distribution
Residence Time Distribution
Residence Time Distribution
PERFORMANCE FOR
HOMOGENEOUS SYSTEMS INTRODUCTION
(CHE 512)
M.P. Dudukovic
Chemical Reaction Engineering Laboratory
(CREL),
Washington University, St. Louis, MO
ChE 512
Spring 2004
Updated 01/05
In dealing with homogeneous reaction systems we have seen that basically all problems encountered can be
classified into two broad categories:
I.
II.
Reactor design is required for a given reaction system, and "a priori" prediction of reactor
performance is therefore needed.
Assessment of performance of existing reactors is asked for, and suggestions for improvement are
solicited.
So far in dealing with these problems we have assumed two limits in mixing pattern (realizing that mixing
and flow patterns affect reactor performance). For systems with continuous flow at steady state they are:
For a single reaction at isothermal conditions, reactor performance is measured by conversion of the
limiting reactant and we know that for each of the ideal reactors it is a function of kinetics, feed
concentrations and thermodynamic limitations i.e
x A = f (Da,
C Bo
CAo
, K)
where
Da Damkohler Number =
(1)
characteristic process time
= P
R
characteristic reaction time
C Bo
reactant ratio in the feed
C Ao
K - equilibrium constant
For an n-th order irreversible reaction with constant density, a mass balance on the limiting reactant A
readily yields an explicit relationship for the Damkohler number as a function of conversion, which, for
each of the two ideal flow patterns, takes the following form:
xA
;
(CSTR)
(2a)
Dan =
(1 x A ) n
1
ChE 512
1
D an =
1 (1 x A) 1n
1n
1
Da1 = ln
1 xA
Spring 2004
; n1
(PFR)
(2b)
; n =1
where
Dan = kCAn1
o t
For example, for a first order irreversible reaction the ideal reactors design equations become:
xA
In a CSTR: Da1 = 1 x ,
A
In a PFR:
1
,
Da1 = ln
1 xA
Da1
1+ Da1
(3a)
x A = 1 e Da1 .
(3b)
xA =
Implicitly, in deriving the above expressions, we have assumed that the feed can be intimately premixed,
down to molecular level, without reaction (i.e the mixing time in generating a homogeneous feed is
extremely short compared to characteristic reaction time and no reaction occurs during that time).
Thus, we realize that besides the characteristic reaction and process time we have to worry about
characteristic mixing time (or mixing intensity which is proportional to its reciprocal) and about the scale of
mixing. Then
Reactor Performance = f (kinetics, flow pattern, mixing intensity and scale) (4)
We can intuitively recognize that the reaction system and its performance are characterized by four
parameters (in a simplistic approach):
p=
V
Q
R=
1
k C nao 1
m
lm
is small, ideal reactor concepts are useful but even then extremely high
conversion cannot be predicted accurately in real reactors unless we have more information on the flow
pattern and on the distribution of times that various elements of the fluid reside in the vessel. The concept
of Residence Time Distribution (RTD), to be introduced next, is useful here.
2
ChE 512
Spring 2004
When m >> R reactor performance completely depends on the scale of mixing l m and mixing intensity
(1 / m ). If we can describe mixing in the vessel we can describe reaction progress as well. There is
another way of looking at the concepts of scale and intensity of mixing by describing the degree of
segregation of the fluid and the nature of exchange, or mixing, of all fluid elements that entered together
with their environment comprised of all other fluid elements. The former approach is useful in very fast
reactions, the outcome of which is completely determined by the turbulent field. The second approach is
useful for those reactions the rate of which is comparable or slower to the rate of mixing. Then one can say
that
kinetics, flow pattern , state of aggregation
Reactor Performance =
of the fluid and earliness of mixing
(5)
We can now consider the description of the mixing process by two levels: macromixing and micromixing.
Macromixing is a global description of mixing yielding the information on how long various
elements of the incoming flow spend in the process vessel.
Micromixing is a more detailed description of mixing on a microscopic level providing the
information on the scale of aggregation of the fluid and earliness or lateness of mixing, i.e on the
environment that the individual fluid elements experience in the vessel.
CHARACTERIZATION OF
MACROMIXING
(CHE 512)
M.P. Dudukovic
Chemical Reaction Engineering Laboratory
(CREL),
Washington University, St. Louis, MO
Chapter 2
2.
CHARACTERIZATION OF MACROMIXING
2.1
In describing a
flow pattern in any flow system, reactors in particular, since the hydrodynamic
equations of flow are too complex to solve, it is useful to at least provide the information on what is the
distribution of residence times for the outflow. Here by residence time we mean the time that a fluid
element (particle) spends within the boundaries of the system (reactor). We will first define all the
functions that are customarily used to characterize the flow pattern.
2.1.1 Exit Age Density Function, E(t)
We define a probability density function or exit age density function by
(1)
The flowing fluid contains entities that are conserved. These entities may be molecules, atoms,
particles, etc. and from now on we will call these conserved quantities fluid elements.
2.
Every fluid element has some original entry point and final departure point from the system.
3.
The system consists of a volume in a three dimensional space and there is no ambiguity with
regard to its boundaries.
4.
Fluid elements have zero age as they enter and acquire age equal to the time spent in the system.
Aging stops if the fluid element leaves the system but resumes if the same element returns into the
system and stops completely when the element leaves never to return again. At that point its age
becomes the residence time of the element in the outflow.
(2)
E (t) dt = 1
(3)
Equation (2) simply reminds us that fractions of the outflow of any residence time must be non-negative
and that residence time can take only positive values. Equation (3) requires that the sum of all fractions
be unity.
(Updated 01/05)
Chapter 2
(4)
and is obtained by summing all the fractions of the outflow between residence time of 0 and t . In terms
of probability theory, F(t ) is the probability that the fluid element of the outflow has residence time
less than t .
2.1.3 Washout Function, W(t)
The so-called washout function, W (t ) , is the probability that the fluid element in the outflow has
residence time larger than t ; it is the fraction of the outflow of residence times larger than t.
W (t ) = 1 F (t ) = E( )d
(5)
The functions defined so far (i.e. E(t ), F(t ), W (t )) are based on the fluid elements of the outflow as
their sample space (population) and characterize the outflow. (How to determine these is a question that
we will address later).
2.1.4 Internal Age Density Function, I(t)
Let us now consider the fluid elements within the system at some actual time t a = 0 and consider how
they are distributed in their ages, the age of an element being the time that elapsed since its entry to the
system.
Let us define:
(6)
where I( ) is the internal age density function. To relate I( ) to the functions already defined we can
consider the system at time and seconds later. Make a mass balance on fluid elements around
age
(Fluid elements in the system of age about + )
(7)
Since elements of any age other than zero cannot be introduced to the system (since age can only be
acquired by residing in the system) the first term on the right hand side is zero. The other terms, using
the definitions introduced earlier, can be expressed as follows:
2
(Updated 01/05)
Chapter 2
VI ( + ) VI ( ) = Q E( )
3
144244
(7a)
E ( ) d
where + .
The limit process gives
Q
I ( + ) I( )
=
lim E( )
{
V O
O
dI
1
1 dF
1 dW
=
E( ) =
=
t
t d
t d
d
lim
{
(8)
The last two equalities in eq. (8) are obtained using the relationship between E, F,W defined by eqs. (4,
5). We also took the mean residence time to be t = V / Q . The boundary condition required to solve eq.
(8) is that at I = 0 and W = 0, F =1 so that
I ( ) =
1
1
[1 F ( a)] = W ( )
t
t
(9)
Now it only remains to be established that the mean residence time, which is indeed the mean or first
moment of the E function, is equal to V / Q .
(mean residence time) = (residence time t of a fluid element) x
t
V
V dI
t = t E (t ) dt = t
dt =
tI +
Q
Q o dt
o
o
V
I
dt
o = Q
(10)
t +
1
1
t I = t I(t) dt =
t dt E( ) d =
t o
t
o
t
1
E( ) d t dt = 2t
o
o
E( ) d
(11)
We will see later that we can put to good use the moments of the E curve defined by:
n = t n E(t )d t
(12)
= 1,
(Updated 01/05)
tI=
2
2t
Chapter 2
2
2 1
(13)
It is customary to use the variance, 2 , or the second central moment, which measures the spread of
the curve. The central moments are defined by:
nc
= (t 1 ) E( t )d t
(14)
2c
= 2 12 = 2 t 2 = 2
(15)
t 2 t
2
+
= [1 + D ]
2 2t 2
(15a)
so that
Then:
tI=
where
2
D
2
2
2
D = / t .
We can now ask the question as to what the above density and distribution functions look like for ideal
reactors.
For a PFR all elements of the outflow have the same residence time equal to the mean residence time.
E PFR (t ) = (t t )
(16a)
F (t ) = H (t t )
(16b)
W (t ) = 1 H (t t )
(16c)
I (t ) =
1
[1 H (t t )]
t
(16d)
t
since there is no spread of the E curve around the mean. Then
2
ECSTR (t ) = ICSTR (t )
(17a)
(Updated 01/05)
Chapter 2
(18a)
I (t ) dt = 1
(18b)
1 t t
e = ECSTR (t )
t
(19)
F (t ) = 1 e t t = 1 W (t )
(20)
tI = t
(21)
The mean age of the fluid elements in a perfectly mixed stirred tank is equal to the mean residence time
of the exiting fluid. This is to be expected since the assumption of perfect mixing requires that there is
no difference between the exit stream and the contents of the vessel.
2.1.7 Intensity Function, (t )
(t )dt = (fraction of fluid elements in the system of age between t and t + dt that will exit
during the next time interval dt ).
(t ) can be determined from the previously defined functions by the following procedure:
(t ) =
(22)
E (t )
E (t )
E (t )
=
=
t I (t ) W (t ) 1 F (t )
(23)
(Updated 01/05)
Chapter 2
It is customary to use dimensional time = t / t . Then the set of function defined over the domain
are E , F , I , . Since
E ( ) = E(t )
dt
= t E(t )
d
(24)
Equation (24) indicates that if in the functional form for E(t) we substitute t = t and
multiply the whole function with the mean residence time, t , we obtain the dimensionless E ( ).
Equation (24) is a well known relationship from probability and statistics where it is presented as the
general rule for an independent variable change in a p.d.f. We are in fact compressing the independent
variable by 1 / t and, hence, we have to expand the ordinate by t in order to preserve the area under the
E ( )d = 1 .
F ( ) = F (t )
(25a)
I ( ) = t I (t )
(25b)
W ( ) = W (t )
(25c)
( ) = t (t )
(25d)
SUMMARY
Review: E, I, W, F, and their inter-relationships
dI
1
1 dW
1 dF
= E =
=
dt
t
t dt
t dt
t
F = E(t )d t = 1 W
o
E (t )
=
t I( t )
E (t )
W (t )
(Updated 01/05)
V
Q
Chapter 2
t E (t) dt
0
n = t n E(t ) dt
= 2D =
2 12
12
when 0 = 1
t I (t) dt
PFR E = (t t )
2.2
t
[1 +
2
1
CSTR E(t) = e t / t
t
t
t
E(t)dt = E( ) d
Experimentally we can obtain directly the exit age density function, E( t ) , residence time distribution,
F( t ) , and internal age density function, I( t ) , by using tracers. By injecting, for example, a step input of
tracer at time t = 0 at the inlet of the system we can monitor the distribution of residence times of the
tracer elements in the outflow. This information can be inferred from some signal measured in the
outflow which is proportional to tracer concentration such as light absorption or transmission, reflection,
current, voltage, etc. This output signal can be interpreted in terms of the residence times of the tracer
only if:
a)
The system is closed, i.e the tracer enters and leaves the system by bulk flow only, i.e diffusion or
dispersion effects are negligible in the inlet and outlet plane.
b)
Tracer injection is proportional to flow, i.e. at the inlet boundary tracer injection rate is
proportional to the velocity component normal to the boundary at each point of the boundary.
c)
The total rate at which the tracer leaves the system is the integral of the product of the velocity
times concentration integrated in a vectorial sense over the whole exit boundary.
In addition, the residence time distribution of the tracer will yield the residence time distribution for the
carrier fluid, which is what we want, if and only if:
a)
the system is at steady state except with respect to (w.r.t.) the tracer concentration;
b)
the system is linear, i.e the response curve is proportional to the mass of tracer injected;
c)
the tracer is perfect, i.e . behaves almost identically to the carrier fluid;
d)
there is a single flowing phase and single homogeneous phase within the system;
e)
f)
Under the above set of conditions we can interpret the response to a step-up tracer injection to directly
obtain the F curve for the carrier fluid. Suppose we had no tracer in the inlet stream (white fluid), and
(Updated 01/05)
7
Chapter 2
then at time t = 0 we started injecting the tracer (red fluid) at such a rate that its concentration at the
inlet is CO . The quantity of tracer elements injected per unit time is Q CO . For each tracer element
there are K carrier fluid elements that , if tracer is perfect, behave identically to the red elements of the
tracer. Hence, KQ CO white fluid elements are entering the system per unit time. We monitor at the
outlet tracer concentration C . At each time t all the red elements that we see at the outlet, QC , have
residence times less than t because they only could have entered the system between time 0 (when
tracer injection started) and time t . For each red element of residence time t there must be K white
elements of carrier fluid that have the same residence time since they entered with these red elements
and have behaved in the same manner. By definition
F (t )
QC + KQC
(1 + K )QCO
The denominator above results from a simple mass balance. Per unit time we feed into the system Q Co
tracer elements and K Q Co carrier fluid elements. Total rate of input must be equal to the total rate of
output which, therefore, is (1 + K ) Q Co . Then the RTD of the carrier fluid is given by
F(t) =
C(t )
Co
(26)
The E, W, I curves can now be evaluated from the RTD (i.e, the F curve) by previously reported
relationships.
The residence time density function, E( t ) , can also be obtained directly from an impulse tracer
injection. During a short time interval dt at t = 0 we inject a pulse of m T of tracer elements (and of
course for each tracer element K carrier fluid elements entered). At the outlet we monitor tracer
response C . During time period dt at time t we collect QC dt tracer elements and KQC dt elements
of the carrier fluid of the same residence time. All these elements have residence time between t and
t + dt because they entered the system between time 0 and dt . By definition:
E(t) dt
(Updated 01/05)
Chapter 2
E(t) dt =
Q C dt + K Q Cdt
mT + Km T
Q C dt (1+ K)
m T (1 + K)
QC
mT
E(t) =
(27)
Other functions F, I, W can now be derived from the E curve using previously reported relationships.
The total mass balance on tracer in a pulse injection requires that all the tracer injected must eventually
mass balance and ensure that the experiment is executed properly. The formula can also be used (when
one is confident that tracer is indeed conserved) to determine the unknown flow rate:
Q =
(28)
C (t)dt
0
Whenever the mass balance for the tracer is not properly satisfied the tracer test does not represent a
proper way of determining the E curve. Various pitfalls were discussed by Curl, R. and McMillan, M.
L. (AIChE J., 12, 819, 1966).
The washout curve, W (t ) , can of course be obtained from the step-up tracer test by subtracting the F( t )
curve from unity, W = 1 F . This produces inaccurate results at large times because of subtraction of
numbers of the similar order of magnitude since
step-down tracer test.
lim F =1
t .
Imagine that at the end of the step-up test both the inlet and exit tracer
concentration are CO . Now at t = 0 we start the stop watch and reduce the inlet tracer concentration to
0. Then all the tracer elements appearing at the outflow at time t are older than t since they have
entered the system before time 0. Due to linearity and perfect behavior of the tracer, for each tracer
element there are K carrier fluid elements of the same residence time. By definition:
W (t ) =(fraction of the outflow that has residence times larger than t )=
(1+ K )QC C
=
(29)
(1+ K )QC o C o
The area under the washout curve gives the mean residence time:
(Updated 01/05)
V
=
Q
Chapter 2
W (t ) dt
(30)
Based on the previously derived relationships you should be able to prove the above.
The I( t ) curve can be evaluated readily using the step-down tracer test, its integral and eq (31):
I (t) =
1
W (t)
t
(31)
The I( t ) curve is often determined directly in biomedical applications by injecting a pulse of tracer mTo
at time 0, and by monitoring the response of the whole system (not the outflow) which is proportional to
the mass of tracer remaining at time t , m T . Then:
I (t ) =
m T (t)
m To
(32)
threaded ends with a self balancing bridge working at 1000 Hz to eliminate polarization).
Dyes - colorimetric detectors and spectrophometer may give nonlinear response.
2.3
We know what the age density functions for the two ideal reactors (PFR and CSTR) look like. Using
these we can derive the age density functions for a generalized compartmental model with time lags. By
this we mean that if we have evidence or reason to believe that a real flow system or reactor can be
represented by a set of CSTR's (well mixed compartments) and PFR's (time lags) in series or parallel we
can readily derive an E curve for any such combination.
Here, and in later applications, it is very useful to use Laplace transforms defined by
10
(Updated 01/05)
E (s ) =
Chapter 2
L {E (t )} =
s t
(33)
E (t ) d t
Any network of CSTR's and PFR's consists of: elements in parallel, elements in series, split points and
mixing points (where points are considered to have no volume). The rules for dealing with these are
explained below.
2.3.1 Systems in Parallel
V
(1 )Q
(1 - ) V
Consider two parallel branches represented above. The top branch contains the fraction of the total
volume of the system and the bottom branch has (1 ) of the total volume. The splitting point S and
mixing point M have no volume. Let the transfer function (Laplace transform of the unit impulse
response and hence the Laplace transform of the residence time density function, E ) of the top branch
be E 1 (s ) , if that branch had all the volume of the system and all the flow passed through it. Let
E 2 (s) be the Laplace transform (LT) of the impulse response of the bottom branch if all the volume
was in it and all the flow ran through it. Since only volume V is in the upper branch and Q flow
rate passes through it, the LT of the response is
E1
1
E2
s . The overall response is obtained by the weighted average of the two, where at point M
1
1
+ (1 ) E 2
s
1
(34)
j=1
jE
j s
j
(35)
where j is the fraction of the total flow rate going through branch j and j is the fraction of the total
volume of the system that is in branch j .
11
(Updated 01/05)
Chapter 2
We build block diagrams for flow pattern representation out of two ideal patterns (plug flow and
complete backmixing) or their combinations. We know the E ( ) curves for our building blocks with
= t / t where t is the total mean residence time. They are the impulse responses of a PFR and a
CSTR. If the building block is a subsystem then the mean residence time is t sub =
t ( = fraction of
total system's volume present in subsystem, - multiple or fraction of total system's throughput that
flows through the subsystem) so that
sub
= t /t
sub
t
= . Then the subsystem's dimensionless
t
More generally speaking we are stating that if the dimensionless response of a system with volume V
and flow rate Q that exhibits a certain flow pattern is E ( ) , then the impulse response of the same
system when it is a subsystem, i.e. a building block within a larger system) (containing fraction of the
volume of the whole system and with flow rate of Q going through it) is given by E .
The dimensional impulse response of the subsystem then is given by
E ( sub ) d sub = E
sub
(t ) dt
(36)
1
E
t
(t / t )
(37)
Thus, by replacing t by t / we get Esub (t) from E(t) . Now the previously stated relation for the
Laplace transforms follows.
L {E(t )}= e
st
L {Esub (t)}=
t
by replacing s with
1
E(t ) dt =
t
st
E (t / t )dt = E (s )
st
(38)
t
1 s u u
s (39)
E
dt = e E du = Esub
t
t O
t
s in E (s ) we get Esub .
Also:
L {E ( )} = e s E ( )d = E (s )
(40)
12
(Updated 01/05)
Chapter 2
L Esub ( ) = e
E d =
su
E (du) =
E
s
sub
(41)
_______________________________________________________________________
Example
1.
V1
S
Q
Q2
V2
Q1 + Q2 = Q
Q1
Q
V1 +V2 =V
V1
V
Now E (s) = E1
where t1 =
1
s + (1 )E2
s
1
(1 )
t , t2 =
t
(1 )
E1 (s ) = E2 (s) =
1
; i = 1 or 2
1+ ti s
E (s) =
1+
ts
+ (1 )
1
1
ts
1+
1
(1 ) t
2 t (1 ) 2 (1 )t
e
+
e
E (t ) = L 1 {E ( s)}=
t
(1 )t
2.3.2 Systems in Series
Let us now consider a system in series consisting of two elements as shown below.
13
(Updated 01/05)
Chapter 2
V1
V2
V1 = V
V2 = (1- ) V
Let the transfer function of the first one be E 1 (s ) and of the second one E 2 ( s ) . If the first one
contains a fraction
of the total volume of the system then the systems response in the Laplace
((1 )s)
(42)
E (t ) =
1
t
E1
E2
d
1
1 1
(43)
E (s ) = E j ( j s )
(44)
j =1
where E j (s ) is the Laplace transform of the impulse response of the j-th individual subsystem as if it
contained the whole volume of the system.
_____________________________________________________________________
Example:
1
1+ t s
where t = V / Q.
1
fraction of the total volume. Hence
N
1
s
E j ( s) = E j =
t
N
1+ s
N
14
(Updated 01/05)
Chapter 2
1
1
=
t
t N
j =1 1+
s
1+ s
N
N
E (s) =
The overall impulse response then is obtained by inversion of the Laplace transform
1
N 1
1
E(t ) = L1
N = L
t
N N
t
1+ s
s +
t
N
N N t N 1 Nt/t
=
e
t (N 1)!
The dimensionless E-curve is
E ( ) = t E (t ) =
NN
N 1e N
(N 1)!
We can now use the above to derive the responses of N-CSTR's of equal size in series, CSTR's and
PFR's in parallel and systems with recycle.
In the example on nonideal stirred tank that follows the section on systems with recycle we will show
how to determine the parameters of the model from the experimentally determined E-curve and how to
use this in assessing reactor performance.
15
(Updated 01/05)
Chapter 2
Once you have mastered the derivation of the transfer function for subsystems in series and subsystems
in parallel, you should be able to handle systems with recycle. The only new rule is the splitting rule
where by if you split a stream each of the outgoing streams possesses the same transfer function.
Let us consider a general recycle system depicted below:
BV
M (R+1)Q
G1
G2
V
(R+1)Q
S
Q
RQ
(1 )V
Flow rate Q flows through a recycle system (the system within the dashed box is the system with
recycle) of total volume V. Internally, at point M flow rate Q is joined by recycle flow rate, RQ, so that
the flow rate of (R+1)Q flows through the forward branch of the system that contains volume V . At
splitting point S, R, Q, is recycled through the recycle branch of volume (1 )V which flow rate Q
s
). The
leaves the system. The transfer function of the forward branch is G1 (we mean by it G1
R + 1
(1 )s
transfer function of the recycle flow branch is G2 i.e. G2
.
R
The transfer function of the total system is E (s ) .
Applying what we have learned so far, we note that for a normalized impulse injection of (t ) the
transfer function for the inlet is 1. The mass balance in the Laplace domain yields
(1 + R G E ) G = (R + 1) E
2
(45)
This response (R + 1)E is obtained by the product (i.e. convolution in the time domain) of the transfer
function for the forward path, G1 , and the transfer function for the inlet stream after mixing point M.
The later is the sum of the transfer function for the fresh inlet stream (i.e. 1) and the product (i.e.
convolution in time domain) of the transfer function for the recycle branch, G 2 , and the transfer
function for the exit stream, E , multiplied by R due to flow rate being RQ.
We can now solve eq (45) for the transfer function of the system E :
E=
G1
R + 1 R G1 G2
(46a)
16
(Updated 01/05)
Chapter 2
1
R +1
G1
R
1
G1 G2
R +1
(46b)
G1 n G2 n
n = 0 R + 1
(47)
1
E=
R
G1 n G2 n 1
n = 0 R + 1
(48)
To obtain the impulse response in the time domain one can now invert the series in equation (48) term
by term.
Example
Consider a recycle system with volume V in the forward branch and volume (1 )V in the recycle
branch as sketch below. Plug flow occurs in both branches so that
G1 = e t s ( R + 1)
G2 = e (1 ) t s R
(R + 1) Q
RQ
Plug Flow
Plug Flow
(1 ) V
t ns
R R +1
e
e
n = 1 R + 1
(1 ) t ( n 1)s
R
1
R
t n (1 ) t (n 1)
R
R +1
R
n = 1 R + 1
17
(Updated 01/05)
Chapter 2
(1 )t
R
t
2t
1
+
E (t ) =
t
t
2
R +1
R + 1 (R+ 1)
R +1
R
2(1 ) t
R2
3 t
+
t
3
R +1
R
(R + 1)
+ ...
E (t )
1
R +1
(R + 1)2
R2
(R + 1)3
t a1
t a2
t a3
t an
Where
t a1 =
t
R +1
t a2 =
(1 )t
2t
+
R +1
R
t a3 =
3 t
2(1 )t
+
R +1
R
t an =
(n 1)(1 )t
n t
+
R +1
R
Now, mentally perform a tracer experiment in which you injected normalized mass of tracer of unity
(i.e. mT Q = 1 ). The first time the tracer appears at the exit is at t a1 = t R + 1 , which is the time it
takes to transverse volume V at flow rate (R+1)Q. (Recall t = V Q always!). The tracer is carried
by (R+1)Q, only Q exits and RQ is recycled. So the amount of tracer carried out by the outflow at t a1 is
1/(R+1). Hence 1/(R+1) is the area under that delta function at t a1 .
For the tracer to appear the second time, it must transverse the volume in the recycle branch ( (1 )V at
flow rate RQ, which takes (1 ) V RQ and then transverse to forward branch again, which takes
18
(Updated 01/05)
Chapter 2
another t a1 . So t a2 = 2 t a1 +
(1 )V .
R
(remember 1 (R+ 1) left at the first passage time) and again only the fraction
1
is recycled. So the
R +1
area under the second delta function at t a2 is R (R + 1) . The reast follows by analogy.
2
We note that if we had only the response curve we could tell how much volume the recycle system has
associated with the forward and recycle branch. The difference
t an t an 1 =
Since t a1 =
t
R +1
(1 )t
R
= const
t
R +1
We know that if
t an t an 1 = t a1
All the volume is in the forward stream i.e. = 1 , Otherwise
t an t an 1 t a1 =
(1 )t
(*)
1
R
with the area under the second peak A2 =
R +1
(R + 1)2
we get
A1 R + 1
=
A2
R
(**)
We can use eqs (*) and (**) to get estimates for and R.
19
(Updated 01/05)
Chapter 2
A reactor of volume V = 25 m (25,000 lit ) with flow rate Q = 1000 lit / min = 1 m 3 / min which was
3
designed to operate as a CSTR and give very high conversion for a 2nd order irreversible reaction ( A
product) is operating poorly at x A = 0.75 . A pulse of mi = 250g of tracer is injected instantaneously
into the reactor. At the outlet the following exit concentration is measured for the tracer. Initially rapid
fluctuations within the first five seconds of very high tracer concentration are observed. Afterwards the
following data is obtained:
t (min)
10
20
30
40
2.15 1.10
50
60
0.70
70
80
a)
b)
c)
What volume of a perfect CSTR do we need to get conversion that currently is produced by the
well mixed region?
The rapid initial rise of tracer concentration in the outflow seems to suggest bypassing. The slope of
ln C vs t is
slope = -
ln (1.06 / 0.01)
= 0.0549
85
We cannot a priori discard the possibility of stagnancy either. Use Cholette and Cloutier model
schematically presented in the figure below.
(1 ) Q
Q
S
Q
Vm
Vd
1+ t s
(Updated 01/05)
Chapter 2
2 t
e
E(t ) = (1 ) (t ) +
t
where t =
V
.
Q
mi E (t )
, then:
Q
t
m
2 t
C(t ) = i (1 ) (t ) +
e
Q
t
mi ( )2
= 10.7 (mg / L ) = 10.7x10 3 (g / L)
Q t
mi ( )
Q
I
S
Q I
mi S
Q Cexp ( 0)
1000 10.7x10 3
=
= 0.778
mi
250
S
0.055
0.78
1 = 0. 22 of the flow bypasses the vessel
25
= 25 min
1
Hence, we can find now the fraction of the total volume, , that is actively mixed.
= 0.055
t
0.055 t
0.78
= 0.567 0.57
0.055 x 25
0.57
V active = 0.57 x 25 = 14.25 m
V dead = 10.75 m
Now we need to set the CSTR design equation in order to find the unknown rate parameters:
21
(Updated 01/05)
Chapter 2
V
CAo x Ar
V active
=
=
2
2
Qactive
( )Q
kC Ao (1 x Ar )
Now x Ar is the actual conversion produced by the reactor found in the stream leaving the active section
of the reactor before mixing with the bypass stream
kC Ao =
x Ar
2
(1 x Ar ) t
However, we must first relate this conversion x Ar to the conversion produced by the reactor as a whole.
This requires a balance around the mixing point M. The stream Q arriving from the reactor has
conversion x Ar , the stream bypassing the reactor has conversion of zero.
represented by:
1 x Ar = 1 x A so that x Ar =
xA
0.75
= 0.96
0.78
0.96
0.78
=
= 32.8 (min 1 )
2
(1 0.96) 0.57 x 25
Recall that V =
b)
c)
kC Ao
QxA
, Then
(1 x A )2
V new (x A = 0.75) =
1,000 x 0.75
3
2 = 366 (lit) = 0.37 m
32.8 (1 0.75)
V new (x A = 0.96) =
1000 x 0.96
= 18,293 lit = 18.3 m3
2
32.8 (1 0.96)
________________________________________________________________________
22
(Updated 01/05)
Chapter 2
Let us consider now a single CSTR with bypassing. If the portion of the flow that bypasses is (1 )
the schematic of the system is as follows:
(1 )Q
Q
S
and the impulse response and its transform are given below by eqs. (48) and (49). We expand the
transform (see eq (49)) to evaluate its moments, and we tabulate below the dimensionless variance of
the system as a function of the fraction of flow that bypasses which is (1 ) .
Clearly, the
dimensionless variance is larger than one for a CSTR with bypassing indicating pathological behavior.
E(t ) = (1 ) (t ) +
2
t
e t / t
(49)
t
t2
t2 2
= 1 + 1 s + 2 s 2 = 1 t s +
s + O(s3 ) (50)
t
1 +
s
E ( s) = 1 +
By a simple technique that is illustrated later, we have evaluated the moments from the above Laplace
transform expansion as:
0 = 1
1 = t
2 =
2
E
2t 2
2t 2
2
E
t 2 (2 )
0.1
0.2
2
E
2E
t
0.3
0.4
2.33
0.5 0.6
3
0.7
5.67
0.8 0.9
9
19
Now let us consider a single CSTR with 'dead' volume. Such dead volume cannot be reached by tracer
or by reactant molecules.
Q
Q
Vm
Vd
23
(Updated 01/05)
Chapter 2
Let V d = V , so is the fraction of the total reactor volume that is for all practical purposes
inaccessible. Now the impulse response of the system is:
t
1
e (1 )t
E(t ) =
(1 )t
(51)
()
1
2 2 2
3
= 1 (1 )t s + (1 ) t s + O s
1 + (1 )t s
(52)
Recall that
n
(
1)
E (s ) =
n
n=0
n!
sn
(53)
Hence,
o = 1
(54a)
1 = (1 )t
(54b)
2 = 2(1 ) t 2
2
E2 =
(54c)
2 1 2
=1
12
(54d)
Now we detected the presence of dead volume not from the value of the variance but from the fact that
1 < t , i.e the central volume principle is violated. We may note, however, that the central volume
principle is never violated and, while a fraction of the volume of the system may be difficult to
reach, i.e is relatively "stagnant", as long as that volume is a physical part of the flow system under
consideration it is never "dead", i.e it will be reached by at least a few elements of the fluid and,
hence, elements of the tracer if not by flow then at least by diffusion. The fascinating feature of the
central volume principle, which is little known, is that the zeroth moment of the tracer impulse
concentration at any point of a closed system is constant and equal to mT/Q where mT is the mass of
the instantaneous tracer injection and Q is the volumetric flow rate through the system. This means
that the area under the concentration response to an impulse injection is constant anywhere in the
system (i.e at all points of the system) and it implies that at points where the tracer concentration
response rises rapidly, high values will be of short duration, while where the response is barely
detectable, it lasts seemingly forever! This makes sense, as it says that if a point is readily accessible
and easy to get to, it is also easy to get out of and the converse is also true - if it is hard to get to a
point it is equally hard to get out. With this in mind let us say that our system has a stagnancy
24
(Updated 01/05)
Vd = V
Chapter 2
which exchanges its contents very slowly with the main well mixed region Vm at a rate
Q with 1 .
By setting up differential mass balances for the tracer in region Vm and Vd, normalizing tracer
concentrations in Vm and Vd by multiplying with Q/mT to get the impulse response, and by applying
the LaPlace transform and solving the equations in the LaPlace domain you should be able to show
that
ts
E (s) =
(1 ) 2 2
1+ 1 + t s +
t s
1+
(55)
By expanding the denominator via binomial series, multiplying the result with the numerator and
grouping items with equal power of s you should be able to get:
2 2 2
E (s) = 1 t s + 1+
t s + 0(s3 )
(56)
From the above series we readily identify the moments of the E-curve
0 = 1
(57a)
1 = t
(57b)
+ 1 t 2
2 = 2
(57c)
2 12
2
=
=1 + 2
>1
12
2
E
(58)
Clearly no matter what volume fraction is relatively stagnant, the dimensionless variance is always
greater than one. The quantity E2 1 is proportional to the square of the stagnant volume fraction
and inversely proportional to the ratio of the exchange flow rate Q between the ideally mixed
region, Vm, and stagnant region, Vd, and the flow rate Q through the system, i.e inversely proportional
to . However, while the dimensionless variance is larger than one, indicating pathological behavior,
the central volume principle is satisfied and 1 = t . Hence, regarding stagnancy we can adopt a
simple empirical rule. We will get an indication of stagnancy (dead volume) either from the fact that
1 < t , in which case the variance can be anything, or if our data are accurate enough and 1 = t , one
can detect stagnancy from E2 > 1.
25
(Updated 01/05)
Chapter 2
You should note that if , that is if the exchange flow rate Q between Vm and Vd is infinitely
faster than flow rate through the system itself, Q, we recover E2 = 1 and as evident from the
expression for E ( s) we recover the perfect mixer response. After all infinitely fast (instantaneous)
mixing between all regions of the system is the definition of the perfect CSTR, so we should not be
surprised by this result.
It should be clear from the above discussion, however, that just from the fact that the dimensionless
variance is larger than one we cannot tell whether we deal with stagnancy or bypassing. To be able to
distinguish between the two we should examine the shape of the E-curve and especially of the
intensity fraction .
It is tempting to talk about bypassing, if pathology is detected at small times, (e.g. a peak at small
times) and of stagnancy, if pathology is present at large times (e.g. a peak at very large times).
However, small and large time are ill defined. It may be tempting to say that if the problem is caused by
fluid residing less than the characteristic reaction time, R , in the vessel we have a bypassing problem,
and if the problem is caused by fluid of residence times order of magnitude larger than r then
stagnancy is the culprit.
This argument can be summarized as follows:
k CnAo1 t >> 1 stagnancy model is justified
k CnAo1 t << 1 bypassing model is justified
One can argue that stagnancy is justified if there is a significant fraction of the outflow with the
residence times t much larger than characteristic reaction times i.e
t =
1
kCo n
(59a)
1
> 0.1
1 F
kCo n 1
(59b)
Bypassing is then considered as a model if there is a significant fraction of fluid that emerges at times
much less than the characteristic reaction time
26
(Updated 01/05)
t=
t=0
Chapter 2
1
kCO n 1 E (t ) dt >
(60a)
0.1
(60b)
However, if one adopts these definitions then even PFR can exhibit bypassing if t < t r , and CSTR can
exhibit both bypassing or stagnancy as illustrated below:
ECSTR (t ) =
CSTR
1 t / t
e
t
(61)
(t )dt = e t / t t = e t
r
/t
> 0.1
(62)
tr
indicated
whenever
e t/t
tr
0
= 1 e t r /t > 0.1
0. 9 > e
t r /t
tr
1
> ln
.
t
0.9
To preclude considering stagnancy and bypassing for a perfect mixer (CSTR), and to reserve stagnancy
and bypassing for pathological systems only,
2
> 1, or
d
> 0 exhibit either stagnancy or bypassing depending on the criteria above.
dt
Example
Let us consider again the ill fated mixed reactor of our Example problem that was represented by the
Cholette-Cloutier model.
Characteristic reaction time is:
r =
1
1
=
= 0.03 min = 1.8 s
kC Ao
32.8
V
25
=
= 25 min
Q
1
27
(Updated 01/05)
Chapter 2
tp =
An actual ideal CSTR has in the outflow the following fraction of fluid of residence times between
0 and R
R
0
ECSTR (t )dt =
1
t
t/ t
dt = 1 e
R / t
= 1 e 0 .03/ 25 = 1 e 0. 0012
= 1 0.9988 = 0.001
Our flow pattern in the above example problem yields an outflow with the fraction of residence times
between 0 and R of:
E (t )d t =
(1 ) (t ) d t + 0
= 1 + 1 e R / t
= 1 - 0.78 + 0.78 1 - e
2 t/ t
e
dt
t
2 t / t
dt
t e
= 0 + e
R
t
bypassing or stagnancy) is by
(63)
28
(Updated 01/05)
Chapter 2
For a CSTR
1
t
Additional ways of dealing with stagnancy via the introduction of the holdup function have been
discussed by Tester and Robinson in the AIChE Journal in the 1980s.
________________________________________________________________________
Example of Determination of System Parameters by Matching Experimentally Determined E
Curve with Model Predicted One via the Method of Moments.
We illustrate here the Method of Moments which we have used above in our examples of CSTR with
stagnancy and bypassing.
Recall E (s ) =
e
O
st
E(t ) dt
dn E
sn
(
s
=
0
)
n
n!
n = O ds
29
(Updated 01/05)
doE
o
ds
Chapter 2
= E (s = o) =
s=o
E (t )dt
= o
dE
ds
= te st E (t ) dt
=
O
s = o
s=o
t E (t ) dt
= - 1
Hence,
dn E
n
ds
n = (1)
(1)
E (s) =
n= O
E (s) =
s= o
s = 0 1 s +
n
n!
dnE
sn
(
0
)
n!
dsn
n= O
(A)
(B)
By comparison of (A) and (B) moments can be obtained. For our stirred tank with bypassing and
stagnancy example,
E = 1 +
1+ t s
t
2
E = 1 + 1
s + 2 t 2 s 2
E = 1 + t s +
2 2 2
t s
2 2 2
t s
1 =t
E =1t s +
0 =1
2=
2=
2
D
2 2 t
2 2 t
2t
2t 2
[2 ]
2 2 t 2 [2 ] 2 2 0. 78
= 2 =
=
=
= 1.56
0. 78
1
2t2
exp
=t
2t
exp
2
exp
= 1.56
30
(Updated 01/05)
S=
The best estimate of parameters is now the choice of parameter vector k that minimizes the
error criterion, S, in same sense.
In RTD theory and practice most often y(t) = E measured (t), h(t, x) = E predicted (t,k)
Parameter Estimation Methods
1.
S=
[E
(1)
2.
[E
min k S = min k
(t) E p (t,k)]2 dt
(2)
3.
min k S =
4.
1
2
{[Re
E m ( ) Re E p (, k )] + [Im E m ( ) Im E p (, k )]2 d
2
(3)
min k S =
E m ( s) E p ( s, k ) ds
(4)
5.
min k S =
[E
( s) E p ( s, k )] ds
2
(5)
6.
min k S =
[E
(6)
The method of moments places too much weight and over-emphasizes too much the tail of
the curve which is usually inaccurately measured. It is popular due to simplicity.
Least squares method in time domain is the equivalent of the least squares in frequency
domain due to Parsevals theorem.
F 2 (t )dt =
2
1
F ( )
2 x
F ( ) = F ( F ( t )) =
F (t )e
(7a)
tot
dt = Re F ( ) + i / mF ( )
(7b)
Least squares in transform domain weight heavily the errors at small times. Integral
absolute errors in transform domain corresponds to least squares in time domain with
Q(t) =
es,t es1 t
t
n =
t n E(t)dt = (1)n
o
dn E
/s= o
dsn
(8)
=
c
n
(t )
1
E(t)dt
(9)
E (s) =
st
E w (s) =
(t )s
= s1c +
s2 c s3 c
2 + 3
2
3
nc = (1) n
d n s
[e E (s)]
dsn
Task 2
Task 3
Task 4
Task 5
Task 1)
Every tracer measurement also includes the response of auxiliary lines in front
and after the equipment of interest and of the measuring device itself. Deconvolute the
system response from the auxiliary system response as per the sketch below:
(t )
GI
Gs
Em
The normalized impulse tracer response at the exit Em(t) = Q Cm(t)/mT is obtained from the
measured impulse tracer response at the exit Cm(t) divided by the ratio of the tracer mass, mT,
that was instantaneously injected at the inlet and the steady flow rate, Q. This measured
impulse response represents the convolution of the injection system, GI(t), our system of
interest, E(t), and sampling system GS(t). In the Laplace domain this is:
E m (s) = GI (s)E (s)Gs (s) = E (s) Ga (s)
(1)
where Ga (s) = GI Gs (s) is the transfer function for the auxiliary system.
The above eq (1) is only true if we have taken precautions to satisfy the conditions of the
proper tracer test, (i.e. with respect to tracer injection and choice of tracer, if we have a
monitoring, closed system, etc.) and if we have demonstrated invariance of Em(t) to the
mass of tracer injected, i.e. Cm doubles as mT doubles!
Assuming that the above conditions are satisfied, we perform deconvolution of eq (1) to
extract E(t) which is of interest to us
E (s) =
E m (s)
Ga (s)
(2)
(3)
In order to perform the operations required by eq (1) we need to know Ga (s). This
knowledge can be obtained in several ways:
a) -
obtaining a unit impulse tracer response (i.e. transfer function) for the auxiliary system
itself
b) -
Gs
GI
Ga
independent experiments.
The most popular models are:
i)
ii)
iii)
iv)
v)
1
1+ t a s
et a1 s
1+ t a 2 s
The less spreading of the response occurs in the auxiliary system, the more justified the use
of simple time lag et 2 s where t 2 =
Va
+ t instrument
Q
This is OK for many plant data. In laboratory or pilot work often an exact Ga is needed.
Note: From E(t) data that is experimentally obtained we can get E (s) =
st
E(t)dt by using
a number of routine numerical programs. No difficulty involved in this step. In the same
2
st
G(t)dt . However,
performing eq (2) is non-trivial due to the inherent noise of the routine deconvolution
package programs. (Routinely s = i is used and Fourier transform is utilized to perform the
inversion). This step can create a completely distorted E (s) . Thus, filtering algorithms must
be used (see Mills and Dudukovic, AIChE J 34(10), 1752 (1988)).
Inversion back to time domain by eq (3) can also be tricky
+i
E(t) =
1
e st E (s)dt
2i i
(4)
If all of the above is done correctly then we have the experimentally determined E(t) = Ee(t)
that we can try to match to model predicted one Ep(t).
An alternative, which allows us to avoid deconvolution, is to convolute the measured Ga(t)
with model predicted Ep(t) to generate the predicted Emp(t) which can be matched to
experimental observations directly. This approach is to be preferred.
Task 3)
In this step we suggest a model that is simple but seems to be able to capture
the essence of experimental observations. This model then has a predicted impulse response
In general then we will have a difference at each sampling time ti between measured value
Ep (ti) and predicted values Ep(ti,). We want to minimize some function g of that difference,
g[E e (t i ) E p (t i ,P)] weighted with appropriate weights Qi over all times by selecting P in
or
= min P
g[E (t) E
e
(t,P ]Q(t)dt
Our result for the parameters that provide best match depends on the choice of g and Q.
The most frequently used methods are:
A) The Method of Moments
Q(t) = t n ; g = E e (t) E p (t,P)
[E
= n (P)
(t)t n dt
o 4243
1
e
(values from exp data) (moments of the system as function of system parameter)
This is a favored method due to simplicity but requires accurate data and tends to be
inaccurate for moments larger than the second moment.
b) Weighted moments
Q(t) = t n [e at e bt ] g = E e (t) E p (t,P)
[E
Tail effects are de-emphasized. Guidance for selection of a and b is provided. Analytical
solution possible.
c) Least squares in time domain
Q = 1 g = (E e E p ) 2
min P
[E (t) E
e
(t,P ] dt
2
[E
E p
dt = 0
pi
(t) E(t,P)]
for i = 1,3...N
Q = 1 g = E e ( ) E p (,P
1
min P
2
[Re
Really the same as least squares in time domain due to Parsevals theorem
F 2 (t)dt =
1
1F ( ) 2 d
F(t)e
F ( ) = F [F(t)] =
it
dt = Re F( ) + i ImF( )
MinP S
s2
E e (s) E p (s,P ds
s1
E = L(E(t))
es,t es2 t
Q(t) =
t
F) Integral square error in transform domain
min P S S =
[E (s) E
e
(s,P)] ds
2
will be bypassed. By selecting simple models often model testing can be done in a simpler
way. For example
-
semilog plots will reveal CSTR behavior or two CSTRs with distinctly different
t i' s
from the peak, width or height of the curve model parameters can be evaluated for
N CSTRs
N 1
N
2
=
max
N 1
Qmax =
at E imf
etc
When using the method of moments it is the easiest to evaluate the moments of the model
from its Laplace transform
n =
E(t)dt
E=
st
E(t)dt
dE
= est tE(t)dt
ds
o
dn E
= (1) n
dsn
dn E
dsn
t E(t)dt
n
s= 0 = (1)
Thus n = (1) n
st n
dn E
dsn
E(t)dt = (1) n n
s= 0
Often it is much easier to expand E (s) by Taylor theorem around s = 0 and identify the
moments from such an expansion
sn
= (1) n m s n
s= o
n! n= o
n!
dn E
E ( s) = n
n= o ds
= o 1s +
2
2!
s2
s
3!
s3
1
N
t
1+ s
N
Bimodal theorem
(1+ x )
= N
= 1+ x +
x=
( 1)
2!
x2 +
( 1)( 2)
3!
t
s
N
x3
2
N(N 1) t
t
E N (s) = 1 N s +
s
N
2!
N
E N (s) = 1 t s +
(n + 1) 2
t + 0(s 2 )
N2!
E (s) = 0 1s +
By comparison
2
2!
s2 + +0(s2 )
0 = 1
1 = t
2 =
N +1 2
t
N
N + 1 2 2 N + 1 N 2 1 2
t t =
t = t
N
N
N
2
2
1
t
2 = 2 = 2 =
1 Nt
N
2 = 2 12 =
Task 5) Plot model prediction vs data, evaluate mean absolute and relative error, judge by
Mills, P.L. and M.P. Dudukovic, Deconvolution of Noisy Tracer Response Data
by a Linear Filtering Method, AIChE J., 34(10), 1752-1756 (1988).
2.
Mills, P.L. and M.P. Dudukovic, Convolution and deconvolution of Non ideal
Tracer Response Data with Application to Three-Phase Packed Beds, Computers
Chem. Engng., 13(8), 881-898 (1989).
1a.
Q
Q
xA1
PFR
V/2
1b.
CSTR
V/2
PFR
xA1'
CSTR
V/2
V/2
Q
xA2
Q
xA2'
2( t t / 2 )
t
0
where H(t a) =
1
t
H t ;
2
E ( s) =
t
s
2
1 +
t
s
2
t < a
is the Heaviside's unit step function.
t > a
2
_
-t
exponential decay
E(t)
_t
2
EXAMPLE 2:
2a.
2b.
Q
Q
(1 - ) Q
V
Q
(1 - ) V
1 t / t
e
t
E ( s) =
1
1 + ts
1_
tE(t)
exponential
If the RTD of the reactor determines its performance uniquely, then systems 1a and 1b
should yield the same conversion for all reaction orders. Systems 2a and 2b should also
behave alike. We suspect that this will not be the case because we recall that reaction orders
higher than one will give higher conversions if mixing is delayed, and hence we suspect that
system 1a will perform better than 1b and system 2b better than 2a. For a first order
process conversion is entirely determined by the residence time of the individual reactant
elements in the reactor, not by their surroundings in the reactor, and hence we suspect that
systems 1a and 1b will perform the same, and systems 2a and 2b will perform alike. For
reaction orders less than one we suspect that the behavior will be the opposite to that for
orders larger than one.
Let us consider first a 1st order irreversible reaction at constant temperature and of constant
density in systems 1a and 1b.
SYSTEM 1a:
dCA
1
C
= ln A0
kCA
k
CA1
t
=
2
C CA 2
t
= A1
kC A2
2
Hence, 1 x A 2 =
CA 0
C A1
CA 2
e kt /2
=
kt
CA 0
1 +
2
SYSTEM 1b:
First reactor CSTR:
C C'
t
= A 0 ' A1
kCA 1
2
t
=
2
1 x A'2 =
CA'2
CAo
C'A1
C'A2
dC A
1
C'
= ln A' 1
kCA
k
CA2
kt / 2
x A2 = x A2' = 1
e
1 + kt / 2
e Da / 2
1 + Da / 2
where Da = Da1 = kt
Just as we expected the performance of the system is entirely determined by its RTD for a
first order reaction. This can readily be generalized to an arbitrary network of first order
processes.
Now consider a 2nd order irreversible reaction taking place in system 1.
SYSTEM 1a:
First reactor PFR:
t
=
2
dCA
1 1
1
=
2
kC A
k CA1
CAo
CAo
C A1
C C
t
= A1 2 A 2
kC A2
2
1 +
1 xA2 =
CA 2
=
C Ao
2kt CAo
1
t
1 + k CAo
2
kt CAo
SYSTEM 1b:
First reactor CSTR:
CAo CA1'
t
=
2
kCA'1
2
1 1
1
t
=
'
k CA 2
CA1'
2
1 x
'
A2
1 + 2kt CAO 1
CA '2
=
= kt
CAo
2 + 1 + 2kt CAo 1
2
Note that both expressions for conversion are a function of the Damkohler number for the
system, where Da = Da2 = kt CAo .
x A2 = 1
1
2Da
1+
1
Da
1+ Da / 2
x A2' = 1
2 1 + 2 Da 1
Da 1 + 1 + 2 Da
The table below compares the results for various values of the Damkohler number.
x'A 2
Da x A 2
1 0.472 0.464
2 0.634 0.618
4 0.771 0.750
You should plot now x A 2 and x A'2 vs Da. You will find that x A 2 > x 'A2 . When is this
difference larger, at low or high values of Da ?
xA2
Da
3Da2
Da
=
+
2
16
2
Da
Da2
1
+
2
8
Da2
5Da2
16
16
x A'2 = Da 1 +
for system 1a
for system 1b
x A'2
Da x A 2
0.5 0.424 0.426
1 0.708 0.718
2 0.957 0.986
You should check the above expressions and draw your own conclusions by plotting x A 2
and x A'2 vs Da = Da1/2 = kt / CAo1/ 2 .
Da
1 + Da
SYSTEM 2b:
CA1 =
CA2 =
CAo
is the exit reactant concentration from the 1st CSTR.
1 + kt
CA1
1 + kt
CAo
is the exit reactant concentration from the
(1 + kt ) (1 + kt )
2nd CSTR.
The balance around the mixing point where the two streams join yields:
CAb = C A1 + (1 ) CA2 = CAo
1
+
=
1 + kt
(1 + kt )(1 + kt )
CAo ( + kt + 1 )
CAo
CAo
=
=
1
+
kt
1
+
kt
1
+
kt
1
+ Da
(
)(
)
Hence,
x Ab = 1
CAb
Da
=
CAo
1 + Da
x Aa = x Ab
For an exercise consider System 2 and a 2nd order reaction and then a zeroth order
reaction. Are the conversions now different? Why?
3.1
The above examples demonstrate that the RTD of the reactor determines its performance
uniquely in case of first order processes. In a first order process conversion is determined
by the time spent by individual fluid reactant elements in the reactor, not by their
environment in the reactor. Let us generalize this conceptually in the following way.
Assume that all the elements of the inflow can be packaged into little parcels separated by
invisible membranes of no volume. Mixing among various parcels is not permitted, i.e no
elements can cross the membranes (parcel walls) except at the very reactor exit where the
membranes vanish and the fluid is mixed instantaneously on molecular level. Consider the
life expectation density function of the inflow Li and define it as
Li (t ) dt = (fraction of the elements of the inflow with life expectancy between
t & t + dt )
(54)
Clearly, due to steady state, if we look at the inflow at time t = 0 and if we could substitute
among white fluid a parcel of red fluid for the elements of the inflow of life expectancy, ti ,
then ti seconds later the red parcel would appear at the outlet and would represent the
elements of the outflow of residence time ti . Hence, at steady state the life expectation
density function for the inflow equals the exit age density function
Li (t ) = E (t )
(55)
Hopefully, it is clear to everyone that since the flow rate is constant, and there can be no
accumulation, the fraction of the fluid entering with life expectancy of ti must equal the
fraction of the fluid exiting with residence time ti . This means that we can consider that
each of our hypothetical parcels formed at the inlet engulfs the fluid of the same life
expectancy. When parcels reach the exit, each parcel will contain fluid of the same
residence time specific to that parcel. Since, there is no exchange between parcels during
their stay in the reactor, each can be considered a batch reactor. The reactant concentration
at the outflow is then obtained by mixing all the parcels of all residence times in the right
proportion dictated by the exit age density function and can be expressed as:
CAout =
reactant concentration
CAout = lim
C
Abatch
(t i ) E (ti ) ti
i =1
ti 0
N
CAout =
C
Abatch
( ) E ( ) d
(56)
where CAbatch is obtained from the reactant mass balance on a batch reactor, which for an nth order reaction takes the following form:
dCAbatch
d
n
= kCAbatch
; at = 0; CA batch ( 0) = CAO
(57)
(58)
so that
k
CAout = CAo e E ( ) d
(56a)
The above expression is nothing else but the Laplace transform of the E function evaluated
at s = k , i.e, Laplace transform variable takes the numerical value of the rate constant.
Conversion can then be calculated as:
xA = 1
CAout
CAo
= 1
e
k
E ( ) d = 1 E (s ) s =k
(56b)
xA = 1
e kt / 2
e Da / 2
= 1
1 + kt / 2
1 + Da / 2
(59a)
For system 2:
xA = 1
1
kt
Da
=
=
1 + kt
1 + Da
1 + kt
(59b)
xA = 1
e
Da
E ( ) d = 1 E ( s) s = Da
(56c)
The model that we have just developed is called the segregated flow model because it
assumes that the fluid elements that enter together always stay together and are surrounded
at all times by the fluid elements of the same age, except at the outlet where they finally mix
intimately with elements of all ages in proportion dictated by the residence time density
function E( t ) . A fluid for which this model is applicable behaves as a macro fluid and has
the tendency of travelling in clumps. As we have seen the segregated flow model gives the
exact prediction of performance for 1st order reactions, but formula (56) is general and
predicts exit concentration for any reaction order for a macro fluid. However, for nonlinear
rate forms, n 1 , as we have seen from our examples, reactor performance does not only
depend on the RTD but also on the details of the mixing pattern. Hence, formula (56)
cannot predict for n 1 the exact performance but perhaps predicts a bound on the
performance. In dimensionless form, the prediction of the segregated flow model is:
xA =
x
Abatch
( ) E ( ) d
(60)
= Dan (1 x A ) ; = 0, x A = 0
(61)
For an n-th order reaction eq. (60) gives either the upper or lower bound on conversion
depending on the concavity (convexity) of the x A vs curve. We always expect to
obtain higher conversion in a system where the reaction rate on the average is higher. Eq.
(60) requires that mixing between fluid elements of various ages occurs only at the exit. Let
us consider only two fluid elements of equal volume but different age, and hence different
reactant concentration, and let us examine how mixing or lack of it affects the reaction rate
for the system comprised of these two elements. We consider the rate obtainable if the
elements are first intimately mixed, rm , and the rate if they remain unmixed, rav . These two
rates are:
CA + CA2
rm = r 1
2
rav =
( )
r CA1
( )
+ r CA 2
2
10
(62)
r
r2
rm
rav
rm
rav
rav
rm
r1
CA1 CAm CA2
FIGURE 5:
CA
CAm CA2 CA
CA
CA
CAm CA2
CA
( )
( )
= r1 , r CA2
react each at its own rate, the average rate lies on the chord connecting the points
r + r2
CA1 ,r1 and CA2 ,r2 and is rav = 1
. If the two fluid elements are mixed first, the
2
1
concentration in the mixed element of double volume is CAm =
CA1 + CA2 , the rate at
2
this new concentration is rm = r CAm and lies at the rate vs concentration curve evaluated
( )
at the abcissa of CAm . Clearly then, whenever the cord lies above the curve (the curve is
concave up) rav > rm and late mixing, or fluid segregation by age, leads to increased rate
and larger conversion. When the cord is always below the curve rav < rm , then late mixing
or fluid segregation leads to reduced rates and reduced conversion. For first order reactions
micromixing, i.e. earliness or lateness of mixing plays no role, it is only the RTD of the
system that counts.
The above discussion can be generalized (see E.B. Nauman and B.A. Buffham, Mixing in
Continuous Flow Systems, Wiley, 1983) and a proof can be given that for all monotonic
d 2r
rate forms for which
the segregated flow model gives the upper bound on
2 > 0
dCA
d 2r
conversion, while for the monotonic rate forms for which
the segregated flow
2 < 0
dCA
11
model gives the lower bound on conversion. For first order processes
d 2r
2 0 and the
dCA
x
xA =
Abatch
( ) 2e
1
2
2
dx A batch
d
= Da 1 x A batch
1
H d
2
(63)
(64)
Da
1 + Da
(65)
e 2
1 1 + Da d
2
(66)
2
(1 + Da )
Da
(67)
12
we get
2
2 1 + Da
e
xA = 1
Da
2
1 +
Da
e x
dx
x
(68)
Again it is instructive to plot x A for segregated flow vs Da and compare it to the two
models based on ideal reactor concepts for system 1a and 1b. You can find the values of
the above integral tabulated as values of one of the family of exponential integrals.
The segregated flow model represents a useful limiting behavior of the system and gives a
bound on performance for monotonic rate forms. It's requirements are: i) that the fluid be
segregated by age, and ii) that mixing between elements of various ages occurs at the latest
possible time, i.e. only at the reactor exit. Thus, every point in the system has its own age
i.e. the point age density function in the system is a delta function
I p ( ) = ( p ) where p is the mean age for the point under consideration.
Consider the sketch below as a representation of the segregated flow model:
+
Q
VI ( )
VI ( + )
Q E ( )
Q
V
FIGURE 6a:
13
We can write:
where + . Now with the help of age density functions the above can be
written as
+
VI ( + ) = VI ( )
+
Q E ( ) d
d
(69)
+
lim
0
I ( + ) I( )
=
lim E( )
(69)
O
dI
= E ( )
d
= I = 0
t
(70)
(70a)
We have derived, based on the above, the already well known relationship between the I and
E function.
Now, based on the above model framework, we focus on the reactant A (see Figure 6B). We
make a balance on the fluid elements of reactant between age and + .
Q
C Ao
VI ( ) CA ( )
VI ( + )C A ( + )
VI ( ) rA
QC A ( ) E( )
Q
C Aout
FIGURE 6b:
14
(71)
Here r A (CA ( )) = ( RA ) is defined as positive for the rate of disappearance of reactant
A . It is multiplied by since this is the time during which the reaction takes place.
2
Dividing by and taking the limit as 0 we get:
t lim
O
lim
I ( + ) CA ( + ) I ( ) CA ( )
=
O
{E( ) CA( )}
lim
O
{I ( )r (C
A
d
( ICA ) = EC A t IrA
d
dI
dCA
t
C A + tI
= ECA t IrA
d
d
dCA
ECA + t I
= EC A t IrA
d
dC
t I A + rA = 0
d
( ))}
(72)
(73)
(74)
= 0,
C A = CAo
which is the reactant mass balance in a batch system!
15
(75)
The outflow is now represented by summing and averaging the composition of all side
streams that form the exit stream (see Figure 6).
CAout
C ( )
A
E( ) d
(76)
16
3.2
Besides classifying the fluid elements according to their age, one can also group them
according to their life expectation, where life expectancy is the time that will elapse between
the time of observation of the system and the moment when the fluid element leaves the
system. In this manner each fluid element could be characterized by its residence time t so
that t = + . At the entrance, i.e. in the inflow each fluid element has a life
expectancy equal to its residence time. The life expectation density function for the
elements of the inflow, as seen before, equals the exit age density function Li (t ) = E (t ).
Each fluid element of the outflow has zero life expectation and an age equal to its residence
time, and the age density function is the density function of residence times. The fluid
elements in the system are distributed in their age, I( ) , and life expectancy, Ls ( ) . The
number of elements in the system of life expectancy about , VLs ( ) d , must equal the
number of elements of age = , VI ( ) d , if steady state at constant volume is to be
maintained. Thus, the life expectation density function of the elements in the system equals
the internal age density function I( ) = Ls ( ).
We can now view micromixing, i.e. mixing process on small scale down to molecular
level, to be a transition of the fluid elements from segregation by age to segregation by life
expectation. Namely, at the reactor entrance all fluid elements are segregated by age and
have the same age of = 0 and are distributed in life expectation with the density function
E( ) . At the reactor exit all fluid elements are segregated by life expectation and have life
expectation = 0 while they are distributed by age with the density function E( ) . The
details of this transition from segregation by age to segregation by life expectation could
only be described with the complete knowledge of the flow pattern and turbulent fields.
However, just like it was possible to describe one limit on micromixing (the latest possible
mixing) by the segregated flow model by maintaining segregation by age until zero life
expectancy was reached for each element of the system, it is also possible to describe the
other limit on micromixing. This other limit of the earliest possible mixing, or the
maximum mixedness model, is obtained by requiring that all elements of various ages are
assembled (segregated) by their life expectancy. Zwietering suggested and proved (Chem.
Eng. Sci. ll, 1 (1959)) that the following conditions should be satisfied for the maximum
mixedness model.
a)
All fluid elements within any point of the system have the same life expectancy, i.e.
the point life expectancy density function is a delta function Lp ( ) = ( p ) ,
and all elements of the point exit together.
17
All points with equal life expectation p are either mixed or at least have identical
b)
v
+
VI ( + )C A ( + )
()
VI rA
VI ( )C A ( )
Q
C A ( = 0)
()
QC Ao E
Q
CAo
FIGURE 7:
The life expectation scale runs opposite to the volume scale of the system which is counted
from the inlet. Thus dv = VI ( ) d represents the element of reactor volume occupied
by the fluid of life expectancy between and + d because at = , v = 0 and at
= 0, v = V.
The model schematically shows that we have distributed the inlet flow according to its life
expectation and are mixing it with the fluid of other ages but with the same life expectation.
We can now write a balance on all fluid elements between life expectation and + :
(Elements of the system with life expectancy about ) =
(Elements of the system with life expectancy about + ) +
(Elements of the inflow of life expectancy about added to the system during time
).
18
where + .
()
(77)
O
t
I ( + ) I ( )
=
()
lim E
O
+
dI
= E( )
d
(78)
I=0
(78a)
Again we have derived the known relationship between the I and E function.
Let us make now a balance on all fluid reactant elements between life expectancy
and + (Refer to Figure 7).
(Reactant elements in the system of life expectancy about ) =
(Reactant elements in the system of life expectancy about + ) +
(Reactant elements of life expectancy about added to the system during time
)(Reactant elements of life expectancy about that have reacted during time )
VI ( ) CA ( ) = VI ( + ) CA ( + )
+ Q E C Ao - VI rA CA
()
()
O
()
O
( ) ( ( ))
d ( ICA )
t
= CAo E t IrA
3
142d4
19
( ( ))
(79)
(79)
dI
dC
t CA t I A = CAo E t IrA
d3
d
12
E
ECA t I
dCA
= CAoE t Ir A
d
E
dCA
+
(CAo CA ) rA = 0
d
tI
(80)
or
E ( )
dCA
= rA
(CAo CA )
d
W( )
(80a)
This is the governing differential equation for the maximum mixedness model. For
elements of extremely large life expectation the change in reactant concentration is small.
Therefore, the proper boundary condition is
dCA
= 0
d
(80b)
W ( ) CAo CA ( = )
=
E ( )
rA (CA ( = ))
(81)
x = e ;
= ln
1
x
(82a)
or
(b)
x =
1
;
1 +
=
1x
x
(82b)
20
(83)
2
where = / (cm / s) is the kinematic viscosity of the fluid and is the energy
dissipated per unit mass. In water k = O (10 m) . The characteristic diffusion time then
is
D
10 3 )
k2
(
=
=
D
10 5
= 10 1 s
(84)
Water and lower chain hydrocarbons are likely to behave like a microfluid, and a CSTR at
sufficient energy dissipation will operate at maximum mixedness condition except for very
fast reactions and unpremixed feeds, i.e. for reactions with characteristic reaction time
< 102 s .
However in highly viscous polymer systems k = O (100 m)
or larger and
2 2
(10 )
10 8
21
is likely to behave as in a segregated flow mode. Remember, however that if the feed is not
premixed segregated flow leads to no conversion.
For a CSTR
t =
W ( )
= t and eq (5a) leads to the condition:
E ( )
CAo CA
rA (CA )
(85)
1
is independent of , eq (80a) is satisfied by eq (85)
t
for all . ThereforeCA ( = 0 ) = CA and is given by solution of eq (81) which
However, since E( ) / W ( ) =
represents the classical ideal reactor design equation for a CSTR mixed on a molecular
level. Naturally, the reactant concentration in a perfectly mixed CSTR in all fluid elements
is the same and independent of their life expectation.
Thus, for a CSTR, the segregated flow model yields:
1
CA =
t
e
t /t
CAbatch (t )dt =
e
CAbatch ( )d
(86)
C Ao CA
.
r A (C A )
(87)
On the other hand for a PFR, the segregated flow model gives
CA =
(t
t ) CAbatch (t ) dt = CAbatch (t )
(88)
Maximum mixedness model for a PFR, upon reflection, shows that condition (81) is not
valid, and must be replaced, due to the presence of a (t t ) (CAo CA ) in eq (80a), by
= t , CA = CAo
(89)
(90)
22
(91)
t = 0
(92)
CA = CAo
CA =
C
Ao
e kt E (t )dt
(93)
(94)
dCA
= 0
d
(95)
One should note that the function appearing in the above formulation of the maximum
mixedness model, i.e. E / W has special significance, is called the intensity function and is
very useful in evaluation of stagnancy or bypassing within the system as described in the
previous section.
23
For a PFR
For a CSTR
1
t
3.3
Attempts have been made to characterize the intermediate levels of micromixing between the
extremes of segregated flow and maximum mixedness by appropriate models. Early on
Danckwerts (1951) and later Zwietering (1952) introduced the degree of segregation as a
measure of micromixing in the system and proposed ways to evaluate this measure
experimentally or to calculate it from the models.
Key Classical References:
1.
2.
3.
24