Residence Time Distribution

Download as pdf or txt
Download as pdf or txt
You are on page 1of 126

EFFECT OF MIXING ON REACTOR

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

Chemical Reaction Engineering Notes


M. P. Dudukovic
1.

EFFECT OF MIXING ON REACTOR PERFORMANCE FOR


HOMOGENEOUS SYSTEMS - INTRODUCTION

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:

CSTR - perfectly mixed (on molecular level) tank


PFR - plug flow reactor (no mixing in direction of flow, piston like flow)

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

- characteristic process time

R=

1
k C nao 1

- characteristic reaction time

m
lm

- characteristic mixing time


- characteristic mixing scale

When  R >>  m , and l m

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

ChE 512 Course Notes

Chapter 2

2.

CHARACTERIZATION OF MACROMIXING

2.1

Residence Time Distributions, What are they?

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

fraction of the outflow that has resided


E(t ) dt =
in the system between time t and t + dt

(1)

Clearly then E(t ) has units of (time -1 )


The key concepts associated with the above are:
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.

The rules for a probability density function (p.d.f.) require that:


E (t) 0 on t (0,)

(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)

ChE 512 Course Notes

Chapter 2

2.1.2 Residence Time Distribution, F(t)

The RTD, or residence time distribution, can now be defined by


t

F (t ) = ( fraction of the outflow of residence time less than t) = E (t )d t

(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:

fraction of the fluid elements in the reactor


I( )d =
that has an age between and + d

(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 + )

(Fluid elements in the system of age about )


(Fluid elements fed to the system of age to + during time )

(7)

(Fluid elements leaving the system of age to + during time )

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)

ChE 512 Course Notes

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)

2.1.5 Mean Residence Time and Mean Age

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

(fraction of the elements of residence time t)

V
V dI
t = t E (t ) dt = t
dt =
tI +
Q
Q o dt
o
o

V
I
dt
o = Q

(10)

Using eq (8) a proof given above can readily be established. Indeed t = V / Q.


From the definition of the RTD it follows that lim
I (t ) = 0, lim F (t ) = 1
1
424
3
1
424
3
t

t +

The mean age of the fluid in the vessel by definition is:

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)

We have already seen that

= 1,

1 = t and, therefore, in terms of the moments of the E-curve

the mean age is

(Updated 01/05)

ChE 512 Course Notes

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

is the dimensionless variance,

2
2
2
D = / t .

Often we will use 2 for 2D .

2.1.6 Ideal Reactors

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

The mean age in the system is t I =

(16d)
t
since there is no spread of the E curve around the mean. Then
2

the mean age is equal to half of the mean residence time.


For a CSTR, the age density function is the same as the residence time (i.e exit age) density function of
the outflow since there is perfect mixing, and the probability of exiting does not depend on the age of
the fluid element.

ECSTR (t ) = ICSTR (t )

(17a)

(Updated 01/05)

ChE 512 Course Notes

Chapter 2

Using eq (8) this implies


dI 1
= I
dt t

(18a)

We also know that I (t ) must e a p.d.f. so that

I (t ) dt = 1

(18b)

Equations (18a) ad (18b) yield readily:


I CSTR (t ) =

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 )

The conditional probability density, or the intensity function, (t ) , is defined as follows:

(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:

(Elements in the system that are of aget to t + dt ) x


(Fraction of elements of age between t and t + dt that will exit in next interval dt ) =
(Elements in the outflow of residence time between t and t + dt collected during dt )
V I (t )d t ( t )d t = (Q d t )E (t ) dt

(t ) =

(22)

E (t )
E (t )
E (t )
=
=
t I (t ) W (t ) 1 F (t )

(23)

(Updated 01/05)

ChE 512 Course Notes

Chapter 2

2.1.8 Dimensionless Representation

It is customary to use dimensional time = t / t . Then the set of function defined over the domain
are E , F , I , . Since

(Fraction of the fluid of residence time t to t + dt ) =

(Fraction of the fluid of residence time to + d )


E(t )dt = E ( )d
This implies

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

curve to be unity i.e

E ( )d = 1 .

Similarly, using well known relations from the theory of

distributions we can show that

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)

ChE 512 Course Notes

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

How to Obtain RTDs or Age Density Functions Experimentally

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)

the system has one inlet and outlet;

f)

tracer injection does not perturb the system.

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

ChE 512 Course Notes

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 )

= (fraction of the outflow of residence times less than t )

(elements of the outflow of residence times less than t )


(total elements of the outflow )
(tracer elements) + (carrier fluid elements of the same residence time)
=
(total elements of the outflow )
=

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

= (fraction of the outflow of residence time between t and t + dt)

elements of the outflow of residence time


between t & t + dt collected during interval dt
=
(total elements collected during dt )

(Updated 01/05)

ChE 512 Course Notes

Chapter 2

(tracer elements in the outflow collected during time dt at t) +


carrier fluid elements of the same
residence time as tracer elements

(total elements collected during dt )

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

emerge, which is the same as requiring that

E (t) dt = 1 . This formula is used to check the tracer

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 .

Therefore, W (t ) can be obtained directly by a

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)

ChE 512 Course Notes

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)

Some Other Items of Interest:


Perfect tracer -detectable, yet same behavior as carrier fluid and at infinite dilution.
Radioactive tracers - half life >> t (only exception positron emitters).
Electrolytes - conductivity meters (cast epoxy tubular body carbon ring electrodes and female pipe

threaded ends with a self balancing bridge working at 1000 Hz to eliminate polarization).
Dyes - colorimetric detectors and spectrophometer may give nonlinear response.

- thermal conductivity detectors.


- flame ionization detectors (organics)
- R ( CO2 , SO2 , NH3 in mixture of diatomic gases) two channel design preferred.

2.3

How To Derive Age Density Functions

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)

ChE 512 Course Notes

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

and similarly for the bottom branch

1
E2
s . The overall response is obtained by the weighted average of the two, where at point M
1

the weighting is accomplished proportionally to flow.


E (s) = E

1
+ (1 ) E 2
s
1

(34)

This can readily be generalized to M - branches in parallel


E (s) =

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)

ChE 512 Course Notes

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

impulse response is: E ( ) .

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)

When = = 1 we get the whole system's response.


E (t) =

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)

ChE 512 Course Notes

Chapter 2

L Esub ( ) = e

E d =

However the mean of E sub is not 1 but

su

E (du) =


E
s

sub

(41)

_______________________________________________________________________
Example

1.

2-CSTR's of different volume in two parallel branches:


Q1

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)

ChE 512 Course Notes

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

domain (transfer function) is given by


E ( s) = E 1 ( s) x E

((1 )s)

(42)

In time domain this is:


t

E (t ) =

1
t

E1
E2

d
1
1 1

(43)

This is the well known convolution theorem for linear systems.


The above rule can be readily generalized to N subsystems in series, each containing a fraction j of
the total volume of the system. The overall transfer function is given by
S

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:

Find the E-curve of N-equal sized CSTRs in series.


We know that for a single stirred tank of volume V and with flow rate Q the impulse response is
given by
1 t/ t
E1 (t ) = e
t
so that the transfer function is
E1 (s) =

1
1+ t s

where t = V / Q.

For N equal size stirred tanks in series each tank contains j =

1
fraction of the total volume. Hence
N

1
s
E j ( s) = E j =
t
N
1+ s
N

14

(Updated 01/05)

ChE 512 Course Notes

Chapter 2

The overall transfer function is obtained by eq (44):


N

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)

ChE 512 Course Notes

Chapter 2

2.3.3 Systems with Recycle

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)

ChE 512 Course Notes

Chapter 2

Which can also be written as:


E=

1
R +1

G1
R
1
G1 G2
R +1

(46b)

We recall that G1 , G2 are functions of ( s R + 1) and (1 ) s R , respectively.


To get the impulse response in the time domain, inversion of equation (46a) or (46b) can be attempted
using the usual rules for Laplace transform inversion.
Often it is necessary to expand equation 46(b) by binomial theorem to get:
1
E=
G1
R +1

G1 n G2 n

n = 0 R + 1

(47)

Which can be rearranged to the following form:

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

According to equation (48)


1
E=
R

t ns

R R +1
e
e

n = 1 R + 1

(1 ) t ( n 1)s
R

Inversion, term by term, yields


E (t ) =

1
R

t n (1 ) t (n 1)
R

R +1
R
n = 1 R + 1

17

(Updated 01/05)

ChE 512 Course Notes

Chapter 2

Let us write out explicitly the first three terms

(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)

+ ...

This response looks as shown below:

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)

ChE 512 Course Notes

Chapter 2

another t a1 . So t a2 = 2 t a1 +

(1 )V .
R

The amount of tracer that now arrives at the exit is (R R + 1)

(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

(*)

By dividing the areas under the first peak A1 =

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)

ChE 512 Course Notes

Chapter 2

Example of a Nonideal Stirred Tank Reactor

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

c (mg / lit) 6.21 3.52

30

40

2.15 1.10

50

60

0.70

70

80

0.40 0.23 0.13

a)

How can we model the old reactor?

b)

If we had a perfect CSTR what volume do we need for x A = 0.75 ?

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

V m = volume of a perfectly mixed region (CSTR)


V d = 'dead' volume, not accessible to flow
Let V m = V where V = V m + V d .
Fraction of the flow that bypasses is (1 ) . Then the transfer function for the above system is given
by
E = 1 +

1+ t s

and the impulse response is:


20

(Updated 01/05)

ChE 512 Course Notes

Chapter 2

2 t
e
E(t ) = (1 ) (t ) +
t
where t =

V
.
Q

If we recognize that C(t ) =

mi E (t )
, then:
Q
t

m
2 t
C(t ) = i (1 ) (t ) +
e

Q
t

If we plot on semilog paper


m ( )2
t , then slope = 0.055 = slope = S
ln C = ln i
t
Q t
t
Extrapolation of the exponential to t = 0 gives
I = Cexp (0 ) =
Cexp (0 )
S

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

From the information given we have the mean residence time:


t=

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)

ChE 512 Course Notes

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.

This balance can be

represented by:

(1 )FAo + FAo (1 x Ar ) = FAo (1 x A )


1 + (1 x Ar ) = 1 x A

1 x Ar = 1 x A so that x Ar =

xA

We are told that x A = 0.75 . Then, for x A = 0.75 and = 0.78


x Ar =
kC Ao

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)

ChE 512 Course Notes

Chapter 2

2.3.4 Bypassing and Stagnancy

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

1 1.22 1.5 1.86

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)

ChE 512 Course Notes

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)

and its Laplace transform is


E(s ) =

()

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)

ChE 512 Course Notes

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)

Hence the dimensionless variance of the E-curve is

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)

ChE 512 Course Notes

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

E( t )dt > 0.1

(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)

ChE 512 Course Notes

t=

t=0

Chapter 2

1
kCO n 1 E (t ) dt >

(60a)

0.1

F(tr ) > 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

which certainly is true whenever


tr
< ln 10 which would then indicate "stagnancy". Similarly, bypassing would be
t

indicated

whenever
e t/t

tr
0

= 1 e t r /t > 0.1

0. 9 > e

t r /t

; i.e. for all

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

we say that only systems that are pathological

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

Characteristic design process time


t =

V
25
=
= 25 min
Q
1

27

(Updated 01/05)

ChE 512 Course Notes

Chapter 2

Actual characteristic process time


Va
0.57 x 25
=
= 18.3 min
0.78 x 1
Qa

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

0.78 x 0.03/ 0.57 x 25

= 0.22 + 0.78 (1 - e 0.001642 )


= 0.22 + 0.78 (1 - 0.9984 ) = 0.22 + 0.001 = 0.221
Since this is much larger than 0.001 , and we know that the system is pathological, i.e 2 > 1 , it is clear
that bypassing contributed significantly to our problem of low conversion.
Now let us examine the effect of stagnancy, which we suspect since 1 < t

E(t )dt = (1 ) (t )dt +

2 t / t
dt
t e

= 0 + e

R
t

= 0.78e 0.001642 = 0.779

Indeed stagnancy also contributes to the problem.


Ultimately, the best way to check for pathological behavior (e.g.

bypassing or stagnancy) is by

examination of the intensity function, (t )


Lack of undesirable behavior is indicated by
d2
0
dt 2

(63)
28

(Updated 01/05)

ChE 512 Course Notes

Chapter 2

Fluid passes through the vessel in a regular fashion. The


longer the fluid element has been in the vessel (the larger
its age) the more probable that it will exit in the next time
interval.

Bypassing. Some fluid has increased probability of


leaving at a very young age. Rest of the fluid behaves
normally.

Stagnancy. After main flow leaves, the probability of the


fluid caught in the dead zones to exit decreases until very
large ages are reached at which point all fluid must
ultimately leave.

For a CSTR

Note that the above definition of stagnancy via the


function will not detect an entirely 'dead' volume as used
initially as our example for a single CSTR with dead
volume on page 29. This we can only detect from the fact
that the central volume principle is not satisfied
(1 < t ) .
For a PFR

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

Taylor Series Expansion around s = 0 yields:


E (s) =

dn E
sn
(
s
=
0
)
n
n!
n = O ds

29

(Updated 01/05)

ChE 512 Course Notes

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

From the tracer experiment we get experimentally


t

exp

=t

2t

exp

2
exp

= 1.56

30

(Updated 01/05)

ChE 512 Addendum 1 to Chapter 2


Models for Nonideal Reactors
Approach to Model Development and Parameter Estimation
Two tasks arise:
1. Formulation of a proper configurational model based on: the observed RTD,,
physical description of the reactor, and, partial knowledge of the presumed flow
pattern.
2. Evaluation of model parameters by matching the model predicted RTD to the
experimental RTD.
Then the model can be used to predict reactor performance at given operating conditions.
Use of the model at different conditions requires the knowledge of how the model parameters
change with the change in operating parameters.
Matching experimentally determined RTDs to model predicted ones is just one of the many
problems that arise in the area of parameter estimation. Parameter estimation is then in
general based on minimization of some error criterion, where the error is defined as the
difference between the actual and the model predicted value.
The most general error criterion for a system whose response, or certain property, was
continuously measured as y(t) and the model predicted response, or property,
h(t, x (t,k)), where x(t,k) are the n-state variables of the system while k is the vector of punknown parameters, can be stated as:

S=

G[y(t) h(t, x(t,k ))]Q(t)dt


o

or if data is available at only R discrete points in time:


R

S = G[y(t i ) h(t i x(t i ,k)]Qi


i=1

Q(t) is the appropriate weighting function.

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

Most often used methods for parameter estimation are:


Matching of moments

1.

S=

[E

(t) E p (t,k )]t n dt, n = 0,1,2...

(1)

2.

Unweighted least squares in time domain

[E

min k S = min k

(t) E p (t,k)]2 dt

(2)

3.

Frequency response least squares (matching of transfer function and transformed


data for s = i )

min k S =

4.

1
2

{[Re

E m ( ) Re E p (, k )] + [Im E m ( ) Im E p (, k )]2 d
2

(3)

Integral absolute error in the transform domain

min k S =

E m ( s) E p ( s, k ) ds

(4)

5.

Integral square error in the transform domain

min k S =

[E

( s) E p ( s, k )] ds
2

(5)

6.

Weighted moments method

min k S =

[E

( t ) E p ( t, k )t r ][eat eby ]dt

(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

Other Important Relations


Noncentral moments are defined as

n =

t n E(t)dt = (1)n
o

dn E
/s= o
dsn

(8)

Central moments are defined as

=
c
n

(t )
1

E(t)dt

(9)

E (s) =

st

E(t)dt - moment generating function (transfer function)

E w (s) =

(t )s

E(t)dt = L{E(t ) = e s E (s)

= s1c +

s2 c s3 c
2 + 3
2
3

nc = (1) n

d n s
[e E (s)]
dsn

ChE 512: Addendum 2 to Chapter 2


Evaluation of Flow Model Parameters
Task 1

Perform a well controlled tracer study for the system of interest

Task 2

Eliminate end effects

Task 3

Suggest a flow model for the system

Task 4

Evaluate model parameters from tracer data

Task 5

Evaluate model suitability and either change the model or use it

Task 1)

Proper flow tagging and mixing cup concentration monitoring is essential.

(More about this later when discussing laminar flow)


Task 2)

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)

E(t) = L1{E (s)}

(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) -

by removing the system of interest out of the experimental configuration and

obtaining a unit impulse tracer response (i.e. transfer function) for the auxiliary system
itself

b) -

Gs

GI

Ga

by using a model for Ga (s) and attempting to quantify the model by

independent experiments.
The most popular models are:
i)

simple time lag Ga = et a s

ii)

first order process Ga =

iii)

combination of first order and time lag

iv)

second order response plus time lag, etc.

v)

actual physical model

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

way from the experimentally determined G(t) we can get G (s) =

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

E p (t,P) where P is the vector of the parameters of the model.


Task 4)

Our task now is to match E p (t,P) with E e (t) or E mp (t,P) to E m (t) .

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

the best possible way.


The problem can be stated as
N

min p S = min p g[E e ((t i ) E e (t i P ]Qi


i= 1

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

(t) E p (t,P)]t n dt = O for 0,1,2,...N

where N is the dimensionality of the vector P

= n (P)

(t)t n dt
o 4243
1
e

n ( valluefrom exp erimentaldatai. e . numbers

(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

(t) E p (t,P)]t n [e at e bt ]dt = 0

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

Run a nonlinear parameter estimation package

[E

E p
dt = 0
pi

(t) E(t,P)]

for i = 1,3...N

In general N nonlinear eqs to solve for the N pi s!


In the above Q 1 can be introduced.
D) Least squares in frequency domain

Q = 1 g = E e ( ) E p (,P
1
min P
2

[Re

E e ( ) Re E p (,P)] + [Im E c ( ) Im E p (,P ]dt


2

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( )

E) Integral absolute error in transform domain


S=

MinP S

s2

E e (s) E p (s,P ds

s1

This corresponds to least squares in time domain with

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

emphasizes errors at small times.


Note: One should keep in mind that in most practical situations the above general techniques

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

Example: Find 0, 1, 2 , 2 for N CSTRs in series :


E N (s) =

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

requirements on the model whether the fit is suitable.


For an Introduction to parameter estimation of reactor models from tracer data see:
1.

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).

3. RESIDENCE TIME DISTRIBUTION AND REACTOR PERFORMANCE


Let us assume for the moment that we can calculate from the knowledge of the flow pattern
the RTD of the system or that we can readily obtain it experimentally on a reactor prototype.
Is the reactor performance uniquely determined by its RTD? This seems to be a pertinent
question to ask.
We already know that every flow pattern has a unique RTD. Unfortunately a given RTD
does not determine the flow pattern uniquely. Two examples to illustrate this are shown
below.
EXAMPLE 1:

1a.
Q

Q
xA1

PFR
V/2

1b.

CSTR
V/2

PFR

xA1'

CSTR
V/2

V/2

Q
xA2

Q
xA2'

FIGURE 1: Schematic of two systems with the same RTD.

The E function for both system 1a and 1b is:


2
E( t ) = e
t

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

FIGURE 2: Impulse response of both systems 1a and 1b.

EXAMPLE 2:
2a.

2b.

Q

Q
(1 - ) Q

V
Q
(1 - ) V

FIGURE 3: Schematic of two flow systems with the exponential E curve.

The E-curve for both systems 2a and 2b is:


E(t ) =

1 t / t
e
t

E ( s) =

1
1 + ts

1_
tE(t)

exponential

FIGURE 4: Impulse response for systems 2a) and 2b).

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

First reactor PFR:

t
=
2

Second reactor CSTR:

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

Second reactor PFR:

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

Second reactor CSTR:

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

Second reactor PFR:

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 ?

For a reaction of half order (n = 1 2) we get

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 .

Let us now consider system 2 and a lst order irreversible reaction.


SYSTEM 2a:
x Aa =

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

SEGREGATED FLOW MODEL

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

  after batch reaction time


all ti

  fraction of the outflow 


x
ti   of residence time around t i 

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)

Now for a 1st order irreversible reaction:


 k
CAbatch ( ) = CAoe

(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)

Using the above formula we get:


For system 1:

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)

In terms of dimensionless quantities, t E (t ) = E ( ) , and for a first order reaction we can


write:

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)

Where for n-th order irreversible kinetics


dx A batch
d

= 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

Illustration of the Effect of Mixing or No mixing on the rate for reaction


of orders n = 1, n > 1 and n < 1 .

Rate vs reactant concentration is plotted for reaction orders of n = 1, n > 1, n < 1 in


Figure 5. The two fluid elements have concentrations CA1 , CA2 , respectively, and the

( )

corresponding rates are r CA1

( )

= r1 , r CA2

= r2 . If the fluid elements do not mix, but

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

segregated flow model gives the exact prediction.


In terms of our examples, the segregated flow model prediction for a 2nd order reaction for
Scheme 1 is:


x

xA =

Abatch

( ) 2e


1
2  

2


dx A batch
d

= Da 1  x A batch

1

H    d
2

(63)
(64)

So that upon integration we get


x Abatch =

Da
1 + Da

(65)

Substituting this into the above equation (63) yields:



 e 2
x A = 2e Da 
d = 2e  e 2  d 
1
+
Da

1
 1
2
2



e 2
1 1 + Da d


2


(66)

Upon substitution in the second integral of


x =

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:

Schematic of the Segregated Flow Model for All the Fluid.

An element of volume of the system occupied by fluid elements of age between


 and d is dV = VI ( ) d . Thus, the volume scale, V, where volume is measured
from the reactor entrance, and the age scale are proportional and point in the same direction.
We want now to make a balance on all fluid elements of age between  and  +  , i.e.
in the differential volume, dV = VI ( ) d
as sketched in the figure.

13

We can write:

(Elements in the reactor of age about  +  )


(Elements in the reactor of age about  ) 

(Elements of the outflow of age about

 collected during time interval  )

where      +  . Now with the help of age density functions the above can be
written as
 + 

VI ( +  )  = VI ( )   




 + 

Q  E ( ) d 
d 
 


(69)

Using the mean value theorem we get:


VI ( +  )   VI ( )   =  Q E(  ) 

     + 

Taking the limit as    0 yields


V
Q

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:

Schematic of the Segregated Flow Model for Reactant A.

14

(Reactant elements in the reactor of age about  +   ) =


(Reactant elements in the reactor of age about  ) (Reactant elements in the outflow of age about  collected during   ) (Reactant elements in the reactor of age about  that react into product during time
interval   ).
VI ( +  ) CA ( +  )  = VI ( ) CA ( ) 

(71)

QE( ) CA ( )  VI ( )rA  ( CA ( )) 

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)

This is satisfied for all  only if


dCA
+ rA  = 0
d

(74)

Reorganizing this we get:


dCA
= rA
d




 = 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)

which is the segregated flow model.


It was Zwietering (Chem. Eng. Sci. 11, 1 (1959)) who realized that the segregated flow
model represents one limit on micromixing within the constraints imposed by the RTD of
the system, i.e represents segregation by age and mixing on molecular scale as late as
possible, i.e. only at the exit. The necessary conditions for the segregated flow model can
then be expressed by the requirement that every point in the system has a delta function age
density function I p ( ) =  (   p ) so that mixing of molecules of various ages is
permitted only at the exit.

16

3.2

Maximum Mixedness Model

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)

internal age distribution.


These conditions express the requirement that all the elements, which will leave the system
at the same time and hence will be mixed at the outflow, are mixed already during all the
time that they stay in the system. This means that the elements of the longest residence time
are met and mixed constantly with elements of lesser residence times, but which have the
same life expectancy and with which they will form the outflow. Schematically the
conceptual maximum mixedness model can be described as sketched below.


v
 + 
VI (  +  )C A (  +  ) 

()

VI  rA 

VI (  )C A (  )

Q
C A ( = 0)

()

QC Ao E  

Q
CAo

FIGURE 7:

Schematic of the Maximum Mixidness Model for Reactant A.

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      +  .

()

Then VI ( ) = VI ( +  )  + Q E  

(77)

Taking the limit as    0 yields:


 t lim

  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  

()

()

Taking the limit as   O


I (  +  ) C A (  +  )  I (  ) C A (  )
=
  O

CAo lim E   t lim I  rA  CA 
t lim

  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)

In practical applications of the maximum mixedness model one evaluates C(  = ) by


implementing condition (80b) into eq. (80a) and solving
lim

 

W ( ) CAo  CA ( =  )
=
E ( )
rA  (CA ( =  ))

(81)

Once CA ( = ) = CA is obtained from eq (81), eq (80a) can be integrated backwards


from  =  where CA = CA  to  = 0 where the desired exit reactant concentration
CAexit = CA (  = 0) is obtained. For numerical integration one cannot use an infinite
integration range. Then  = 5t to  = 10t

is often chosen to represent infinity


dCA
( = ) by checking that at such high values of  the value of
is almost zero.
d
Another alternative is to remap  (0, ) to x[0,1] space by a transformation such as
(a)

x = e  ;

 = ln

1
x

(82a)

or
(b)

x =

1
;
1 + 

=

1x
x

(82b)

20

This still introduces a singularity at x = 0 and special precaution has to be taken to


integrate eq (80a) properly.
The maximum mixedness model of Zwietering given by eqs. (80a) and (81) gives us the
d 2r
lower bound on conversion for monotonic rate forms with
> 0 and upper bound on
dCA
d 2r
conversion for monotonic rates with
2 < 0 .
dCA
Now the segregated flow and maximum mixedness model bound the performance of
reactors for all monotonic rate expressions provided the reactor RTD is known. These
models are extremely important conceptually. In practice, their use is limited by our
inability to know reactor RTDs a priori.
These models can be quite useful in assessing the performance of a CSTR for which by a
tracer test it can be proven that there are no stagnancies and no bypassing and that the Ecurve is exponential. For water and liquid hydrocarbons the power input is usually
sufficient to generate small eddies. Kolmogoroff's isotropic turbulence theory predicts that
the microscale of turbulence (i.e the smallest size of turbulent eddies below which turbulent
velocity fluctuations are highly damped by molecular viscosity) is given by:
  3 1/ 4
k 

  

(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

D = O (10 8 cm2 / s) . Then  D =

(10 )
10 8

21

= 10 4 s and a CSTR with a premixed feed

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)

and the maximum mixedness model gives


t=

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)

Then eq (80a) becomes


dCA
= rA 
d

(90)

22

If we substitute t = t   the maximum mixedness PFR equation becomes


dCA
=  rA 
dt

(91)

t = 0

(92)

CA = CAo

and the desired exit concentration is given by CA (  = 0 ) = CA (t = t ) which is exactly


the expression given above for the segregated flow model.
Because the PFR model and its RTD prohibit mixing of elements of different ages, there is
no difference in the prediction of the segregated flow or the maximum mixedness model.
They both yield an identical result.
It remains for the reader to show that a first order reaction in a system of any RTD (not that
of a CSTR or PFR) yields the same result according to the segregated flow model


CA =

C

Ao

e kt E (t )dt

(93)

and the maximum mixedness model


E(  )
dCA
(CAo  CA )
= kCA 
d
W ( )

  

(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.


Fluid passes through the vessel in a regular fashion. The


longer the fluid element has been in the vessel (the larger its
age) the more probable that it will exit in the next time
interval.

23

Bypassing. Some fluid has increased probability of leaving


at a very young age. Rest of the fluid behaves normally.
t

Stagnancy. After main flow leaves the probability of the


fluid caught in dead zones to exit decreases until very large
ages are reached at which point all fluid must ultimately
leave.

For a PFR

For a CSTR


1
t

3.3

Other Micromixing Models

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.

Danckwerts, P.V., Chem. Eng. Science, 2, 1 (1953)


Danckwerts, P.V., Chem. Eng. Science, 8, 93 (1953)
Zwietering, T.W., Chem. Eng. Science, 77, 1 (1959)

24

Mixing Effects in Chemical Reactors


(CHE 471)
M.P. Dudukovic
Chemical Reaction Engineering Laboratory
(CREL),
Washington University, St. Louis, MO

You might also like