Limit Cycles in A CISTR

Universiteit Twente
faculteit der
chemische technologie
Limit cycles
The elimination of limit cycles in a CISTR through
the application of a feedback process controller
Onno Kramer
April 2002
139923
Universiteit Twente
faculteit der
chemische technologie
PO box 217
7500 AE Enschede
The Netherlands
Phone: +31 (0)53 4899111
Internet: http://www.utwente.nl
Limit cycles
Graduate assignment
The elimination of limit cycles in a CISTR through

the application of a feedback process controller
Group:
Process development and design
Faculty: Chemical Technology
Head: Prof. dr. ir. G.F. Versteeg
Commission:
Chairman: Prof. dr. ir. G.F. Versteeg (Twente University)
Member: Prof. dr. ir. J.A.M. Kuipers (Twente University)
Mentor: Dr. ir. D.W.F. Brilman (Twente University)
Company Mentor: Dr. ir. E.P. van Elk (Procede Twente b.v.)
Keywords:
CISTR, cooling, limit cycle, bifurcation, stability, control
Author:
Ing. O.J.I. Kramer
Date:
Friday 7 December 2001
SUMMARY
It is possible to eliminate the occurrence of limit cycles in a continuously ideally stirred tank reactor
(CISTR), in which an irreversible exothermic reaction A(l) P(l) takes place. This can be done through the
application of a traditional feedback process controller, which regulates the extent of cooling. Controlling
the throughput is possible, however nevertheless or the reactor temperature or the conversion can be
preserved, not both.
Problem Analyses
Due to the interaction between the heat withdrawal and the heat generation in the CISTR, dynamical
instability can appear. The resulting oscillating behaviour is generally unwanted and must therefore be
eliminated. The scope of this report is to determine whether a process controller can fulfil this assignment.
In addition, a second objective of this report is to acquire knowledge with respect to the dynamical
behaviour in relation with the controller and phenomenon like the presumed delay.
Solving Method
The key to success is the use of the bifurcation program LOCBIF. This software package has numerical
routines, which can contribute to the construction of distinct stability maps. These stability maps provide
full information with regards to the static and dynamical behaviour of a particular process. Through a
formulated mathematical model, describing the actual physical process and the intervening process
controller, the dynamical behaviour is examined for several distinguished parameters. These parameters
can be both process parameters, like the cooling capacity or throughput, as well as controller parameters
like the proportional gain or the integral time. Three different control methods have been selected: coolant
temperature control, coolant flowrate control and throughput control. All three types are investigated on
their dynamical behaviour for both proportional and for proportional-integral control. For every specific
case, the sensitivity of relevant parameters like the cooling capacity and the presumed delay is
considered.
Results
Changes in continuous process operation are due to internal changes (limit cycles) and external changes
(disturbances). A process, which exhibits distinct dynamical unstable behaviour, can be controlled through
increasing the proportional gain. In many cases, limit cycles shrink and stability is acquired. However,
precaution must always be taken for the processes in which multiplicity exists. Controllers with large
proportional gain values can unmistakably manipulate a process variable (like the coolant flowrate) too far
away from the desired set point. Consequently, transition can be inevitable. Therefore, the choice of the
magnitude of the controller parameters has to be done deliberately. In case of controlling the extent of
cooling, this is certainly possible. Nevertheless, in spite of an appropriate controller configuration, large
disturbances in particular if also the system is dynamical unstable, transition remains apparent. In case of
controlling the throughput, the restriction appears that the deviancy is in view of reactor design is limited.
Therefore, large proportional gain values, which are required to eradicate both limit cycles as
perturbations, are often not feasible. Moreover, in case of flowrate control, the risk of extinction or runaway
is to be concerned with due to multiplicity. The implementation of the integral action in the mathematical
controller model results in the offset being eliminated. To increase the physical realism of the
mathematical model, the presumed delay is introduced. Small delay cannot destabilise a proper controlled
process. Conversely, large delay provokes dynamical instability. Traditional control configuration tuning
methods like the Ziegler and Nichols technique are not appropriate to deal with dynamical instabilities.
Large cooling capacity, in general has a positive effect on the stability of a process in particular in
combination with a suitable proportional gain.
ii
OVERVIEW
Table 1 Suitable base case controller configuration to eliminate limit cycles and preserve stability.
Description
Coolant
Coolant flowrate
Throughput
temperature
control
control
control
Presumed delay
Perturbation
Inlet coolant temperature
Proportional gain
Integral time
Desired base case reactor
temperature
Required conversion
d = 30 [s]
T = 20 [K]
Kc = 5 [-]
I = 600 [s]
T = 468 [K]
d = 30 [s]
T = 20 [K]
Tcool,0 = 303 [K]
Kc = 0.0003 [m3 s-1 K-1]
I = 600 [s]
T = 468 [K]
= 0.68 [-]
= 0.68 [-]
d = 30 [s]
T = 2 [K]
Kc = 0.01 [m3 s-1 K-1]
I = 600 [s]
3 -1
V = 0.012 [m s ] T = 468
[K] with = 0.48 [-]
3 -1
V = 0.002 [m s ] = 0.72
with T = 454 [K]
iii
VOORWOORD
Een voorwoord in een afstudeerverslag is doorgaans bedoeld om o.a. de afstudeercommissie, het
afstudeerbedrijf en de mentor te bedanken. Het bedanken komt straks uitgebreid aan bod, maar eerst zou
ik van de gelegenheid gebruik willen maken om eens te verwoorden waarom ik pas op 34 jarige leeftijd
afstudeer. De primaire oorzaak ligt voor mij zeker bij het feit dat ik een visuele handicap heb. Mensen
vragen altijd: Wat kan je nou wel of niet zien? Dat is lastig. Leg maar eens uit aan een persoon die
gewoon alles kan zien wat 10% zicht is waarbij constant alles heen en weer beweegt zonder perspectief
en ook nog eens wazig. Draai een verrekijker maar eens om want dan krijg je een idee, zeg ik meestal.
Mijn handicap zie ik niet als een excuus voor het feit dat dit studeren zon enorm lang traject werd, wel
een belangrijke oorzaak.
Toen ik een klein jongetje was adviseerden medici mijn ouders om mij naar een school voor slechtzienden
te sturen. Gelukkig zagen mijn vader en moeder in dat ik dat niet wilde en ging ik naar een normale
basisschool. Al vanaf de basisschool werd mij meerdere malen dringend verzocht ander onderwijs te
nemen. Daarom maar naar de HAVO in plaats van het VWO, waar ik koos voor wis-, natuur- en
scheikunde en economie. Er werd mij op het hart gedrukt dat het geen verstandige keuze was om al deze
-vakken te kiezen. Talen, tekenen e.d. was het advies. Ik liet als compromis economie vallen voor
biologie. Zonder doubleren haalde ik de eindstreep. Op de HTS werd mijn handicap steeds meer een
probleem. Ik bleef voor de tweede maal zitten. Het advies was om te stoppen. Maar ik wilde toch verder
en dat was zeer ongebruikelijk. Dankzij de docenten Bronkhorst en van der Meij mocht ik blijven. Door
een nieuw medicijn en in combinatie met een bril met kijker kon ik weer de lessen volgen. Wederom ging
alles daarna op rolletjes. Mijn stages bij de Akzo, de afdeling procesontwikkeling en techniek bij de
Gemeentewaterleidingen Amsterdam o.l.v. Eric Baars en de Heidemij-Arcadis ingenieursbureau waren
uiterst succesvol. Voor mijn afstuderen haalde ik 2 maal een 9 en een 10. Toen ik het beste
afstudeerverslag had van regio Amsterdam
en omstreken en mee deed met de afstudeerprijsvraag van de
e
ingenieursvereniging NIRIA en 4 werd van Nederland begon ik toch in te zien dat je met een handicap en
wilskracht een heel eind kan komen. Iedereen adviseerde mij om door te gaan voor de hoofdprijs te weten
de Universiteit. Deze stap nam ik in 1991 zonder echt goed na te denken of ik dat nou wel wilde. Nu ik dit
voorwoord een week voor mijn afstudeerpraatje schrijf maakt het mij allemaal niet zoveel meer uit. Toch
heb ik eerlijkheidshalve mijzelf wel duizend keer afgevraagd waarom ik dit ben gaan doen. Een heel
belangrijk argument is voor mij dat ik volgens de vele oogartsen de enige ben met zon handicap die zover
is gekomen. Als ik zou stoppen kunnen de ooit toekomstige studenten met dergelijke handicaps te horen
krijgen dat zelfs die hardnekkige Onno Kramer uiteindelijk afhaakte. Nee, dat nooit. Ik wil niet op mijn
geweten hebben dat anderen daarom wellicht negatief advies zouden kunnen krijgen. Een ander
argument is dat ik vind dat je iets moet afmaken waar je aan begint. Wat heb ik dan gedaan al die tijd? In
mijn studieperiode zit een gat van 6 jaar waarvan ik 1 jaar heb gewerkt bij Gelink Adviesbureau als ITr.
Veel gesport: wedstrijdroeien, buitenlandse wedstrijden (topsport), marathons lopen. Zeker 5 jaar heb ik
besteed aan mijn muzikale carrire (veel optredens etc.). De kans is groot dat de beurs die ik heb
aangevraagd wordt toegewezen zodat ik naar Sevilla kan gaan om mijn muzikale kennis te verfijnen. In
deze tussenperiode van weinig studeren, moest ik erachter komen hoe ik in deze wereld sta. De door mijn
ouders met de paplepel ingegoten normen en waarden, klaar staan voor anderen e.d. is voor mij de rode
draad in mijn leven. Toen ik besloot de draad weer op te pakken en begon met de literatuuropdracht bij
professor Kuipers en deze afsloot met een 9 kwam het geloof weer terug. Professor Versteeg liet met het
vak Inleiding procesontwikkeling en ontwerp zien dat je twee uiterste bij elkaar kan brengen: van
fundamenteel wetenschappelijk tot met je schoenen in de mest-praktijk. Chemie is dus niet stoffig en saai.
Niet geheel toevallig zitten deze twee heren in mijn afstudeercommissie. De eerder goedgekeurde
afstudeeropdracht bij de GW-Amsterdam-Kiwa werd niet meer geaccepteerd. Mij werd een alternatieve
opdracht aangeboden bij Edwin van Elk werkzaam bij Procede Twente B.V. Zo gezegd, zo gedaan. De
commissie moest nog een vierde lid bevatten en dat werd Wim Brilman. Het avontuur was lang en ik heb
naast de chemische technologie heel veel geleerd, n, ben benaderd voor een baan bij, DSM, AKZO en
Arcadis. Dus het doorzetten heeft toch nut gehad. Het lijkt er een beetje op dat ik mijzelf aan het
verdedigen ben met dit voorwoord. Ten dele is dat ook waar, omdat ik zo vaak met onbegrip en
kortzichtigheid ben geconfronteerd. Als je probeert om normaal te functioneren vergeten men dat alles
mij meer moeite kost en dat terwijl ik best weinig om hulp vraag.
En nu het bedanken van de mensen die mij hebben begeleid. Ik wil Geert Versteeg bedanken voor het
laten zien dat je onder andere met gezond verstand en een goede planning zelfs de meest saaie en
moeilijke opdrachten leuk en interessant kunt maken. Hans Kuipers wil ik bedanken voor de positieve
ondersteuning die hij aan mij heeft gegeven. Wim Brilman bedank ik voor zijn altijd vriendelijk doch zeer
kritische feedback. Edwin heeft meerdere malen bewezen niet alleen superhandig te zijn met computers
maar ook zijn enorme kennis en inzicht hierin te kunnen vertalen. Een ding vraag ik mij eigenlijk al een
hele tijd af: Waarom houd je niet van sporten maar loop je wel van Oost naar West Engeland? Dat moet
je na 7 december maar eens uitleggen. Edwin heel erg bedankt.
Dit verslag draag ik volledig op aan mijn ouders, omdat zij mij altijd onvoorwaardelijk hebben gesteund en
altijd in mij hebben geloofd.
Onno Kramer
Amsterdam, november 2001
iv
TABLE OF CONTENTS
SUMMARY
OVERVIEW
VOORWOORD
TABLE OF CONTENTS
LIST OF STABILITY MAPS
II
III
IV
V
VIII
PART I
INTRODUCTION
MODEL OF A COOLED CISTR

2.1
INTRODUCTION
2.2
THE MASS AND HEAT BALANCE OF THE REACTION MIXTURE
2.3
MODELLING OF THE COOLING
2.3.1 COOLANT HEAT BALANCE
2.3.2 COOLANT TEMPERATURE EXPRESSION
2
2
2
4
4
5
BEHAVIOUR OF THE MODEL

3.1
INTRODUCTION
3.2
STATIC BEHAVIOUR OF THE COOLED CISTR
3.3
DYNAMIC BEHAVIOUR OF THE COOLED CISTR
3.3.1 STABLE SELF-REGULATING PROCESS
3.3.2 UNSTABLE PROCESS
3.3.3 LIMIT CYCLES
6
6
6
7
7
8
9
THEORETICAL PREDICTION OF LIMIT CYCLES

4.1
INTRODUCTION
4.2
LINEARISED PERTURBATION METHOD
4.3
BIFURCATION THEORY
4.4
LOCBIF
4.5
STABILITY MAP
4.5.1 EQUILIBRIUM CURVE
4.5.2 FOLD BIFURCATION CURVE
4.5.3 HOPF BIFURCATION CURVE
4.5.4 ORBIT CURVE
4.5.5 CONSTRUCTING THE STABILITY MAP
4.5.6 CLASSIFICATION OF THE STABILITY MAPS
4.5.7 INTERPRETING STABILITY MAP
4.6
THEORETICAL PREVENTION AND ELIMINATION OF LIMIT CYCLES
10
10
10
10
11
12
12
12
13
13
13
15
16
18
PROCESS CONTROL
5.1
INTRODUCTION
5.2
FEEDBACK CONTROLLERS
5.2.1 PROPORTIONAL CONTROL
5.2.2 PROPORTIONAL-INTEGRAL CONTROL
5.2.3 OFFSET
5.2.4 COMPARISON P AND PI CONTROLLER
5.3
DISTURBANCES
5.4
CONCEPT OF DEAD TIME
5.5
PROCESS CAPACITY
5.5.1 PROCESS TIME CONSTANT
5.5.2 COOLING TIME CONSTANT
5.6
CONTROLLER TUNING
5.6.1 PERFORMANCE CRITERIA
5.6.2 METHODS OF ADJUSTING FEEDBACK CONTROLLER SETTINGS
5.6.3 ZIEGLER AND NICHOLS METHOD
5.6.4 CONTROLLER CONFIGURATION
19
19
19
20
20
21
22
22
23
23
24
24
25
25
26
26
27
PART II
6
STARTING POINT OF RESEARCH

6.1
INTRODUCTION
6.2
BASE CASE DEFINITION
6.3
STABILITY OF THE BASE CASE
6.3.1 REACTOR DESIGNER STABILITY MAPS
6.3.2 REACTOR CONTROLLER STABILITY MAPS
6.4
OVERVIEW
6.5
BASE CASE LIMIT CYCLES
28
28
28
29
29
31
32
33
STRUCTURE AND OBJECTIVES OF PART II
34
PRELIMINARY CONSIDERATION CONTROLLABILITY

8.1
CONTROLLING THE COOLANT TEMPERATURE
8.2
CONTROLLING THE COOLANT FLOWRATE
8.3
CONTROLLING THE THROUGHPUT
35
35
35
35
PROPORTIONAL CONTROL
9.1
REACTOR DESIGNER STABILITY MAPS
9.1.1 CONTROLLING THE COOLANT TEMPERATURE
9.1.2 CONTROLLING THE COOLANT FLOWRATE
9.1.3 CONTROLLING THE THROUGHPUT
9.2
REACTOR CONTROLLER STABILITY MAP
36
36
36
40
50
57
58
59
60
vi
10
PROPORTIONAL-INTEGRAL CONTROL
10.1 REACTOR DESIGNER STABILITY MAPS
10.2 REACTOR CONTROLLER STABILITY MAPS
10.3 ALTERNATIVES
61
61
61
61
65
68
70
11
CONTROLLED SYSTEM WITH DELAY

11.1 PROPORTIONAL CONTROL
11.2 PROPORTIONAL-INTEGRAL CONTROL
11.3 PROCESS CAPACITY
11.3.1 PROCESS TIME CONSTANT
11.3.2 COOLING TIME CONSTANT
71
71
71
74
77
79
79
81
83
85
85
86
12
CONTROLLER CONFIGURATION TUNING

12.1 CONTROLLING THE COOLANT TEMPERATURE
12.1.1 ZIEGLER NICHOLS TUNING METHOD
12.1.2 CONTROLLER PARAMETER VARIATION
12.2 CONTROLLING THE COOLANT FLOWRATE
12.3 CONTROLLING THE THROUGHPUT
12.4 BASE CASE CONTROLLER CONFIGURATION
90
90
90
91
93
93
93
95
95
96
99
13
THE EFFECT OF FOULING ON THE BASE CASE STABILITY

13.1 INTRODUCTION
13.2 COOLANT TEMPERATURE PROPORTIONAL CONTROL
13.3 COOLANT FLOWRATE PROPORTIONAL CONTROL
13.4 COOLANT FLOWRATE PROPORTIONAL-INTEGRAL CONTROL
CONCLUSIONS
RECOMMENDATIONS
NOMENCLATURE
REFERENCES
APPENDICES
100
100
101
101
102
104
105
106
109
112
APPENDIX 1. DERIVATION CISTR

APPENDIX 2. LITERATURE
APPENDIX 3. VALUE ESTIMATION
APPENDIX 4. SHUTTING DOWN THE INPUT FLOWRATE
APPENDIX 5. COOLING CAPACITY COOLANT FLOWRATE RELATIONSHIP
APPENDIX 6. LOCBIF NOMENCLATURE
APPENDIX 7. LOCBIF SOURCE CODES
vii
LIST OF STABILITY MAPS

Table 2 Base case stability maps for no control
Stability map
Reference
T versus Tcool (varying UA)
T versus Tcool (varying V)
Tcool versus UA (varying V)
Figure 19
Figure 20
Figure 21
p30
p31
p32
Table 3 Reactor designer stability maps (proportional control)

Stability map Reference
Coolant temperature Coolant flowrate
T versus Tcool
T versus UA
Figure 23
Figure 24
p37
p38
T versus V,cool
Figure 43 (Tcool,0 = 303 [K])

Figure 28 (0 0.005 [m3 s-1])
Figure 41 (0 0.02 [m3 s-1])
Figure 42 (0 0.1 [m3 s-1])
Table 4 Reactor controller stability maps

Stability map
Reference
Kc versus I
Controller
Proportional-integral
Coolant temperature control
Coolant flowrate control
Throughput control
Figure 74
Figure 80a (Tcool,0 = 303 [K])
Figure 80b (Tcool,0 = 400 [K])
Figure 80c (Tcool,0 = 435 [K])
Figure 85
Throughput
p45
p46
p47
p49
p41
p43
p44
Figure 67
p56
Figure 50 (KC 0)
Figure 54 (Kc 0)
p51
p53
Kc versus UA
Proportional
Kc versus UA
Proportional-integral
p62
p65
Figure 71
Figure 72
p58
p59
p69
Figure 73
p60
Figure 79
Figure 86
p64
p66
p67
p68
p70
Table 5 Reactor controller stability maps regarding delay

Stability map
Reference
Kc versus d
Proportional control
Proportional-integral control
Coolant temperature control
Coolant flowrate control
Throughput control
Figure 88
Figure 95
Stability map
Process capacity
Reference
Kc versus R
Kc versus cool
Kc versus cool
Figure 108
Figure 110
p73
p76
p76
p76
p77
Figure 103
p80
p82
Figure 107
p84
p85
p87
p88
p89
p89
viii
Chapter I
1 INTRODUCTION
This report is part of the research into the stability of process-operation of gas-liquid reactors performed at
the research group Industrial process development and design of the University of Twente. In gas-liquid
reactors, it appears that at specific conditions sustained temperature oscillations, which are called limit
cycles, occur. In plant operation, these conditions have to be avoided, because they may adversely affect
product quality, catalyst stability and downstream operations and can lead both to serious difficulty in
67
process control and unsafe reactor operations. In some extraordinary cases limit cycles are useful,
namely when they serve as a driving force for a second reactor in series.
It is known from literature that limit cycles do not only occur in gas-liquid reactors, but in many
distinguished processes. The phenomenon of limit cycles is very similar for these processes despite the
different reaction mechanisms, values of system variables, process parameters and conceivably controlmethod. Due to these similarities, it was decided to study the relatively simple case of limit cycle in a
CISTR. The results and obtained knowledge of this study can subsequently be used to study processes
that are more complicated.
To be able to study limit cycles, a mathematical description is required. The mathematical description of
the system of a controlled reactor contains two or more differential equations, which cannot be solved
satisfactory through traditional analytical or graphical methods. With the introduction of the software
package LOCBIF, a new approaching strategy seems to exist.
LOCBIF was already used in other fields than the chemical industry and science. It is capable of
interpreting mathematical problems involving differential equations far more easily than existing methods
could, by supplying stability maps. Thereby it made it possible to study problems which could not be
studied before because of the mathematically complexity. The first results were very promising. The
22
construction of a stability map by van Elk et.al. of the system of a reactor with LOCBIF showed good
36
agreement with the results obtained by method developed by Heiszwolf and Fortuin and Vleeschhouwer
77
and Fortuin . But moreover LOCBIF did this in a much more efficient way.
Although the results were very promising, more investigation had to be done to prove the applicability and
reliability of LOCBIF. Van Elk investigated a system with proportional control, but did not account for
integral control and delay, of which the latter improves the simulation of the real physical process.
In the present report, these aspects are investigated. Thereby, next to the control of the coolant
temperature, also the control of the coolant flowrate and the throughput is investigated with LOCBIF.
As a consequence, the following research objectives can be formulated:
Can limit cycles be predicted for certain reactor design and specific circumstances?
If they can be predicted, can they be prevented?
If they can be prevented, can they be prevented by applying a certain process controller?
If the controller is not configured adequately, what are the consequences for the operating
process?
The resulting outcome of the research is presented in two main parts. The first part discusses the theory
of the reactor system described by mathematical equations. Thereafter, the current status of research is
presented, which is the starting point of the simulations done in the second part.
Chapter II
2 MODEL OF A COOLED CISTR
2.1 Introduction
Although this report is part of the research to the stability of process-operation of gas-liquid reactors, due
to the coincidence of the phenomenon of limit cycles, a relatively simple case of limit cycle in a cooled
CISTR with first order irreversible exothermic reaction is considered. In this chapter, the apparent mass
and heat balance will be postulated.
2.2 The mass and heat balance of the reaction mixture

Consider the reactor shown in Figure 1.
Figure 1 Schematic drawing of a cooled or heated CISTR.

A simple irreversible exothermic reaction 1 takes place in the reactor.
kR
A(l ) P(l )
Reaction: 1
With the following overall reaction rate expression for reaction 1:
r = rA = k R [ A] = k0 e
Eact
RT
(1)
[ A]
3
-1
A liquid enters the reactor with a flow rate of V [m s ] and a temperature T0 [K]. This feed flow contains
component A with concentration [A]0. The tank is considered to be ideally mixed, which implies that the
temperature and concentration of the effluent is equal to the temperature and concentration of the liquid in
the tank V, T, [A] and [P]. The reactor is cooled by a coolant that for example flows through a jacket
around the reactor or flows through a construction of cooling pipes.
The dynamic behaviour of a first-order reaction system in a continuously ideally stirred tank reactor
(CISTR) can be described with a process model consisting of two differential equations based on a
material balance and an energy balance. The mathematical model of the process consists of mass and
energy balances and introducing appropriate constitutive equations. For a complete mathematical
derivation, the reader is referred to Appendix 1.
The general form of the non-stationary transient mass and energy balance for the CISTR is represented
by state equations 2 and 3:
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
= C P V (T0 T ) + ( H R )VR k 0 e RT [ A] UA (T Tcool )
dt
VR
(2)
(3)
Both equations 2 and 3 contain non-linear functions of T.

Equations 2 and 3 are coupled, in the sense, that it not possible to solve one equation
independently of the other.
Commonly, dimensionless numbers have been introduced to simplify the set of equations (2 and 3) and
75
making it more generally applicable. The Uppal et.al. notation is commonly applied in this field of
research in which the following specific numbers and dimensionless variables are essential:
Average residence time:
R =
VR
V
(4)
Conversion:
= 1
[ A]
[ A]0
(5)
The number of transport units which in the literature often is called the Stanton number St:
NTU =
UA
V C p
(6)
Adiabatic temperature rise:
Tad =
H R
[A]0
C P
(7)
The adiabatic temperature rise or adiabatic factor is the extent of adiabatic temperature change if it is
assumed that the reaction goes to completion i.e. =1, The mass balance 2 and heat balance 3 can be
rewritten into respectively equation 8 and 9 implementing equations: 4, 5, 6 and 7,
Eact
d
= (1 )k 0 e RT
R
dt
Eact
dT (T0 T ) NTU (T Tcool )
=
+ Tad k 0 e RT (1 )
R
dt
(8)
(9)
The steady state solution of the CISTR can be found if the left-hand side of equation 2 or 8 and 3 are set
equal to zero i.e. d[A]/dt or d/dt = 0 and dT/dt = 0.
[ A] = [ A]0
E
+ V k e act RT
R 0
V
(10)
=
1+
(11)
1
k 0 R e
Eact
RT
2.3 Modelling of the Cooling

2.3.1
Coolant heat balance
In the literature of process-control, the assumption can be made that the temperature of the cooling fluid
can be described, using the CISTR concept resulting in one extra differential equation (equation 12).
cool C P,coolVcool
dTcool
= coolC P,cool V ,cool (Tcool,0 Tcool ) + UA (T Tcool )
dt
(12)
The left term in equation 12 represents the accumulation term in the cooling apparatus with volume Vcool.
The middle term is the difference in heat per unit of time by coolant convection and the right term is the
transferred heat per unit of time from the reactor. In expression 12 the same assumptions are used for the
heat balance of the reaction mixture (Appendix 1). In the steady state situation equation 12 becomes in
terms of the coolant flowrate:
V ,cool =
UA(T Tcool )
coolC P,cool (Tcool Tcool,0 )
(13)
The coolant flowrate can also be determined using equations 2, 3 and 12 resulting in:
Eact
V ,cool
C P V (T0 T ) + ( H R )k 0VR e
=
cool C P,cool (Tcool Tcool,0 )
RT
[ A]
(14)
Correlations 13 and 14 can henceforth be used to determine the required extent of flowrate in the steady
state situation, or to determine the required extent in case of changed process state.
2
60
Aris et.al. and Ogunnaike et.al. formulated the number of transport units for the coolant, similar to
equation 6 considering the cooled CISTR:
NTU cool =
UA
cool C P,cool V ,cool
(15)
2.3.2
Coolant temperature expression
When modelling the cooling of the reactor one should chose an expression for the coolant temperature,
while the latter is not constant. To acquire a value for Tcool, several methods have been described in the
literature.
1. Average temperature value assumption
The simplest expression for Tcool is the mean of the inlet and outlet cooling fluid temperature:
Tcool = Tcool =
Tcooli,0 + Tcool,1
(16)
2. Logarithmic temperature mean

80
68
According to Westerterp and Roffel an expression for Tcool is the logarithmic temperature mean:
Tlog = T Tcool,log =
(T T
cool ,0
) (T T )
cool,1
T Tcool,0
ln
T
T
cool ,1
(17)
Approximation 1 is preferred because of its mathematical simplicity and is a good approximation when
Tcool,0 T cool,1 . Expression 17 is used when high accuracy is required.
Chapter III
3 BEHAVIOUR OF THE MODEL
3.1 Introduction
In the following chapter, the static and dynamic behaviour of a cooled CISTR will be discussed. Firstly, the
static behaviour is examined, including the concept of multiplicity. Thereafter, various types of dynamic
behaviour including the phenomenon limit cycles will be studied.
3.2 Static behaviour of the cooled CISTR

The static thermal behaviour of a cooled CISTR in which an irreversible exothermic reaction A P takes
place can be presented schematically by Figure 3. The heat of reaction is removed by a coolant medium,
which flows through a jacket around the reactor (Figure 2).
The curve that describes the amount of heat released by the exothermic reaction is a sigmoidal or Sshape function of the temperature T in the reactor (curve A in Figure 3). On the other hand, the heat
removed by the coolant is a linear function of the temperature T (curve B in Figure 3). When more than
one intersection of the energy and mass balance curves appears, there is more than one set of conditions
that satisfies both the energy and mass balance and consequently there will be multiple steady states at
which the reactor can be operated. This is called multiplicity.
When the CISTR is at steady state i.e. nothing is changed, the heat produced by the reaction should be
equal to the heat removed by the coolant. This requirement yields the steady states P1, P2 and P3 at the
34
intersection of curves A and B in Figure 3. Steady states P1 and P3 are called stable , whereas P2 is
unstable. To understand the concept of stability, steady state P2 will be considered. Assume that the
reactor is started at temperature T2 and corresponding concentration [A]2. Consider the temperature of the
feed Ti increases. This will cause an increase in the temperature of the reacting mixture T2. At T2 the heat
released by the reaction Q2 is more than the heat removed by the coolant, Q2 consequently leading to
higher temperatures in the reactor and consequently to increased rates of reaction. Increased rates of
reaction produce larger amounts of heat released by the exothermic reaction, which in turn lead to higher
temperatures, and so on. Therefore, an increase in Ti takes the reactor temperature away from steady
state P2 and the temperature will eventually reach the value of steady state P3. Similarly, if Ti were to
decrease, the temperature of the reactor would take off from P2 and end up at P1. By contrast, if the
reactor is operating at steady state condition P3 or P1 and the operation of the reactor is perturbed, it
would return naturally back to point P1 or P3 from which it started. Sometimes one would like to operate
the CISTR at the middle unstable steady state e.g. the upper-temperature steady state P3 may be very
high, causing unsafe conditions, destroying the catalyst for a catalytic reactor and degrading the product
P. In such a particular case, a controller is required to ensure the stability of the operation at the middle
steady state P2. It is of practical importance to carry out a stability analyse to determine which of the
steady states can be realised i.e. stable and which are unstable; any infinitesimal disturbance will move
the reactor away from such an unstable steady state.
kR
Heat removal
AP
Heat production
AP
Figure 2 CISTR with cooling jacket.
Figure 3 The thermal behaviour of an irreversible, exothermic

reaction A P in a CISTR under non-isothermal
conditions (three steady states).
3.3 Dynamic behaviour of the cooled CISTR

Through creating a chart of a chemical process like Figure 3, one can determine whether a process is
operated in a statically stable region. Disadvantage of Figure 3 is that no information is supplied regarding
the dynamics of the process. Therefore, it is uncertain whether a statically stable steady state is also
dynamically stable. If a steady state is statically and dynamically stable, the process can be operated in
this steady state (preferably at high conversion). It will return to the steady state after a sufficiently small
perturbation. If a steady state is only statically stable, the process can only be operated in the
neighbourhood of the steady state. Then, several situations can be distinct: These will be discussed
below:
3.3.1
Stable self-regulating process
Process variable
Process variable
Consider the behaviour of a process variable, such as temperature, concentration or flowrate, shown in
Figure 4. Notice that at time t = t0 the process variable is disturbed by some external factors, but that as
time progresses its value returns to the initial value and stays there. One can say that the process (Figure
4a and b) is stable or self-regulating and needs no external intervention for its stabilisation. It is clear that
no control mechanism is needed to force the process variable to return to its initial value. If the process of
returning back to its initial value evolves smoothly, it is called asymptotic damping (Figure 4a). On the
other hand, the process variable can, if it returns to its initial value, also exhibit considerable oscillations,
which is often called a spiral point (Figure 4b).
Tim e
Figure 4a Asymptotic damping of a

stable system after a
perturbation. (See also
Table 6). Type 1. Area Ia.
Tim e
Figure 4b Oscillating response (spiral

point) of a stable system
after a disturbance. Type 2.
Area Ib.
3.3.2
Unstable process
Process variable
Process variable
In contrast to the behaviour described above, it is probable that the process variable does not return to its
initial value after it has been disturbed by external influences (Figure 4c-i). Processes whose variables
follow this pattern are called unstable processes and require external control for the stabilisation of their
behaviour. If the response of a process is unstable, the distinction can be made between run away e.g. if
too less superfluous reaction heat is removed (Figure 4f and g), or extinguishes e.g. if the cooling is too
strong (Figure 4h and i).
Tim e
Figure 4f
Response of an unstable
system i.e. Ignition or
runaway. Type 5+6. Area III.
Time
Figure 4g Response of an oscillating

unstable system causing
runaway. Type 5+6. Area III.
Process variable
Process variable
In Figure 4f and h the process variable changes equally until the maximum or minimum value has been
reached. If the response of the process however, is oscillatory (Figure 4g and i) in which the amplitude is
increasing, eventually the process state transverse to another equilibrium state. Therefore, the amplitude
of the peaks is decisive and can cause transition i.e. run away (Figure 4g) or extinction (Figure 4i).
Time
Figure 4h Unstable system. Extinction. Figure 4i

Type 5+6. Area III.
Time
At first the response shows

limit cycles; nevertheless, the
overshoot causes eventually
extinction. Type 5+6. Area III.
Finally, due to the distinguishing characteristic of the kinetic rate term in state equations 2 and 3, and to
the introduction of model extensions i.e. extra differential equations, initially, the response can be in the
opposite direction to where it eventually ends up. This is called inverse response (Figure 4g).
3.3.3
Limit cycles
Process variable
Process variable
Process variable
Next to being stable or unstable as described above, another dynamic mode is possible. Figure 4c-e
shows the phenomenon of self-sustained oscillations. Because this phenomenon occurs in principle
spontaneously (Figure 4d), a constrained perturbation is not inevitable to acquire these limit cycles (Figure
4c). Therefore: The occurrence of limit cycles is not exclusively dependent of a disturbance.
Additionally, the pattern of the limit cycles can vary. The process variable can expose a symmetric
sinusoid curve (Figure 4c) or an asymmetric curve (Figure 4e). The symmetry increases in case the
process state is located towards the transition between dynamic stability and instability.
Tim e
Figure 4c Stable self sustained
Tim e
Tim e
Figure 4d Limit cycles, which occur
Figure 4e Asymptotic oscillations (limit
naturally without a
perturbation. Type 3. Area
II.
oscillations (limit cycles) after

a perturbation. Close to Hopf
(4.5.3). Type 3. Area II.
cycles). Unstable system.

Type 3. Area II.
The process will start oscillating forever in a limit cycle around the statically stable operating point at a
fixed frequency and fixed amplitude (Figure 5). The following example will illustrate the phenomenon limit
cycles more.
1
465
Conversion
Ttemperature
470
2
3
4
0.70
Time
Figure 5a Self-sustained reactor

temperature oscillations
(limit cycles).
0.69
4
0.68
0.67
460
0.71
Conversion
0.71
475
Figure 5b Conversion limit cycles.
0.69
0.68
0.67
460
Time
0.70
4
3
465
470
Temperature
475
Figure 5c Limit cycle around a static

stable but dynamic unstable
steady state. Composition of
Figure 4a and b.
Assume that the CISTR conditions are represented by point 1. Here T and are higher than corresponds
to the stable operating point (dotted in Figure 5), which can be derived if the ODE system 2 and 3 is
solved (steady state P3 in Figure 3). Therefore, the reaction rate is too high; the reaction mixture is heated
up, the temperature increases and simultaneously is reduced. When point 2 is reached, rA has been
reduced already in such an extent that the reaction mixture starts to be cooled by the cold feed, however
rA is still so high that the consumption of reactant A is still higher than the supply of fresh A with the feed.
In point 3 reaches a minimum, rA is now so low that build-up of A in the reaction mixture, which is still
being cooled by the cold feed and the cooling equipment, starts again. When point 4 is reached sufficient
A has been built up to increase the reaction rate in such an extent that the reaction heat evolution is now
so high that the cooling of the reaction mixture stops and the heating up is renewed again. The build-up of
A continues until point 1 is reached, then the whole cycle starts again.
Chapter IV
4 THEORETICAL PREDICTION OF LIMIT
CYCLES
4.1 Introduction
The stability of a process is important for both uncontrolled and controlled processes. Disturbances can
adversely affect the stability of a process and accordingly result in decreased performance or in worstcase scenario imply malfunctioning. Therefore, since long one has attempted to predict the process
stability theoretically. In this chapter, the bifurcation theory which provides the theoretical background to
predict limit cycles, and its application by means of software package LOCBIF, is presented. This yields
the basis to construct and interpret stability maps.
4.2 Linearised perturbation method

One approach to predict the process-stability is the linearised perturbation method. It has been developed
37 68
by Himmelblau and Bischoff
and offers a simple, analytical approach to analyse stability using
linearised versions of differential equations. It can solve up to two differential equations. In this report, the
mathematical description of the problem demands more than just two differential equations. Whereas the
underlying theory is beyond the scope of this report, its outcome is useful for the understanding of the
report and will therefore be presented in 4.5.5 Table 6.
4.3 Bifurcation theory

Although only in simple cases, the analytical method mentioned above could provide some useful
information pertaining to the stability of a process. Advanced methods are required in case more complex
systems are involved. In reactor dynamics, it is particularly important to find if multiple stationary points
exist or if sustained oscillations (limit cycles) can arise. The bifurcation theory contains distinct curves,
which provide useful information about the boundary between stability and instability.
The bifurcation theory of time-depended differential equations have found application in many fields of
research, including chemical reactor engineering, mathematical biochemistry and combustion. The
bifurcation theory provides several methods in which static and dynamical behaviour of processes can be
48
analysed. The exact definition of bifurcation has been given by Kuznetsov :
Bifurcation has been defined as the appearance of topological phase portraits under variation of
parameters. A bifurcation is a change of the topological type of the system as its parameters pass
through a bifurcation.
A bifurcation analysis is aimed at locating the set of parameter values for which multiple steady states will
48
occur. According to Kuznetsov a bifurcation point i.e. critical value, is a point at which two branches of a
curve coalesce as a parameter is varied (point B in Figure 6). Consider function (x,) in which x is a
scalar variable and is a parameter. Figure 6 shows curves AB, BC and BD for which:
f ( x, ) = 0
(18)
is satisfied. In case is increased along AB, there is only one value of x for a given that will satisfy
*
equation 18. However, if is continued to be increased along AB, a bifurcation point is reached beyond
10
which there are two values of x that satisfy equation 18 for a given value of . Consequently, the system
contains multiple solutions i.e. multiplicity.
x
D
x
Figure 6 Bifurcation diagram.

In a bifurcation diagram in fact a curve, which represents the dynamical behaviour of a system (e.g.
reactor temperature), is examined in relation which a distinguished parameter (e.g. cooling capacity). The
bifurcation diagram classifies in a very condensed way all possible modes of behaviour of the system and
transitions between them i.e. bifurcations, under parameter variation. It is desirable to obtain the
bifurcation diagram as a result of the qualitative analyses of a given dynamical system. The bifurcation
diagram, which is a called a stability map, can provide important information about the behaviour of the
system being studied. This stability map will be discussed elaborately in paragraph 4.5.
Suitable methods are available e.g. bifurcation software package LOCBIF (see 4.4). Creating stability
maps using LOCBIF does require neither special scientific skills, nor knowledge of complex non-linear
mathematics or numerical methods. The only condition is that the engineer is capable of describing the
dynamic behaviour of the reaction system using a set of differential and algebraic equations e.g. mass and
energy balances 2 and 3.
4.4 LOCBIF
This chapter reviews the general features of the LOCBIF bifurcation program. The review is limited to the
parts of LOCBIF that are relevant for the analysis of the dynamic and static behaviour of CISTRs. For
45
other LOCBIF features, the reader is referred to the LOCBIF manual . LOCBIF has been developed at
the Institute of Mathematical Problems of Biology of the Russian Academy of Sciences in the Moscow
45
Region in 1992. A. Khibnik and E. Nokolaev developed the numerical algorithms. LOCBIF has no special
system requirements and runs on any modern PC system.
LOCBIF is a very useful tool, which has the numerical routines to explore the existence and stability of
equilibria in dynamic models and to generate stability maps with limited effort. Once familiar with LOCBIF,
it is straightforward to create a LOCBIF input file and analyse the dynamic and static behaviour of a typical
system instead of producing a complete perturbation analyses which is very profoundly and therefore
laboriously. Another advantage is that the linearising of the differential equations, needed for the process
controller, can be skipped. The mathematical equations describing a particular process can directly be
implemented in a LOCBIF source code called RHS (right hand side).
A dynamical system described with one ore more differential equations can be programmed in a LOCBIF
source, in a specific Right Hand Side source code. This subject is discussed comprehensively in the
45
LOCBIF manual . In Appendix 7, LOCBIF right hand side source codes are included and in Appendix 6
the accompanying nomenclature.
11
LOCBIF is a scientific tool for bifurcation analysis of systems of differential equations, which depend upon
parameters. It helps to explore interactively the existence and stability of equilibria (steady states) of
dynamical models. LOCBIF supports systems of up to 10 differential equations.
LOCBIF is based on continuation procedures for relevant local bifurcation curves up to co dimension
three. This means that LOCBIF can calculate bifurcation curves in a multi-dimensional space and can do it
very fast thanks to the efficient continuation algorithms. LOCBIF plots two-dimensional projections of the
bifurcation curves during the computations. The numerical background of the LOCBIF continuation
methods is discussed extensively in the LOCBIF manual.
Besides calculation of bifurcation curves, LOCBIF can also solve the differential equations with respect to
time, like any other ODE solver e.g. Maple. This option is very useful to check the predicted dynamic
behaviour or to find stable steady state solutions. Four different LOCBIF versions exist, but for the
application of interest here, LBEP.EXE is the only LOCBIF version needed.
LOCBIF supports about 15 different curves, however only four of them are of direct interest for our
application i.e: the equilibrium curve, the fold curve, the Hopf curve and the orbit curve.
Some restrictions regarding the LOCBIF program:
LOCBIF cannot deal with integral equations or partial

differential equations.
Maximum number of differential equations, which can be
implemented.
Very limited possibilities to extend the RHS with
mathematical disparities.
Occasionally, LOCBIF has noticeable troubles finding an
initial point.
TIME variable declaration

PHASE variable declaration
PAR parameter declaration
FUN function declaration
COMMON variable declaration
function definition(s)
RHS definition(s)
INIT initial values computations
LOCBIF rhs 1 Structure of ODEs specification in LOCBIF.
4.5 Stability map

A stability map characterises the dynamic behaviour of a system, as a function of important system
parameters like the cooling capacity or the temperature of the coolant. A stability map can be very
convenient to determine quickly whether or not a particular system proceeds to instability. Such a stability
map is composed of equilibrium curves and/or Hopf and/or Fold curves and must be checked by means of
orbit curves. For a good understanding of these stability maps, firstly the main ideas of the above
mentioned curves are discussed after which an example of constructing a stability map will be given.
4.5.1
Equilibrium curve
The equilibrium curve is for chemical reactors also known as the S-curve (Q-T in Figure 3 here T-Tcool)
and represents the solution to:
f ( x, p) = 0
(19)
All parameters in equation 19 but p1 are fixed. This is the steady state solution of the ODE system as a
function of one active parameter. The curve can for example describe the steady state solution
(concentrations and temperature(s)) of a reactor system as a function of the coolant temperature.
4.5.2
Fold bifurcation curve

48
According to Kuznetsov this bifurcation has many names i.e. limit point, saddle node bifurcation and
turning point. The fold bifurcation curve represents the boundary between static stability and static
instability as a function of two active parameters and is defined by:
12
f ( x, p ) = 0
(20)
1 ( x , p ) = 0
In equation 20 all parameters but p1 and p2 are fixed. The fold curve is characterised by one zero
Eigenvalue 1. The projection of a fold bifurcation defines a curve on which a pair of equilibrium points
appears or disappears when parameters are changed in such a way to cross this curve transversally.
The curve can for example define the border between static stable and static unstable steady states of a
reactor system as a function of the coolant temperature and the cooling capacity.
4.5.3
Hopf bifurcation curve
The Hopf bifurcation curve represents the border between dynamic stability and dynamic instability as a
function of two active parameters and is defined by:
f ( x, p ) = 0
1 ( x , p )+ 2 ( x , p ) = 0
(21)
In equation 21 all parameters but p1 and p2 are fixed. The Hopf curve is characterised by two equal
Eigenvalues of opposite sign. This can be achieved in two different ways:
1 ( x, p )= 2 ( x , p )
1 , 2 ( x , p ) = i ( > 0 )
(22)
(23)
Neutral saddle, real Eigenvalues (equation 22) and Andronov-Hopf bifurcation, imaginary Eigenvalues pair
(equation 23). The last situation is of special interest and will result in an oscillating system. The period of
the oscillation being about t = 2/. The Hopf bifurcation can for example define the border between
dynamic stability and dynamic instability of a reactor system as a function of the coolant temperature and
the cooling capacity.
4.5.4
Orbit curve
Whereas a stability map contains the above equilibrium, Fold and the Hopf curves it does not contain orbit
curves, because it represents the solution of the ODE system with respect to time and no parameters are
active except the time, so that all system parameters and an initial condition have to be specified. Orbit
curves are useful to investigate or to verify the dynamical behaviour and stability of a particular chosen
point.
x = f ( x, p)
(24)
All parameters p in equation 24 are fixed. The computation of the curve means numerical integration of
the system in time.
4.5.5
Constructing the stability map
The distinct curves explained above can be drawn in a bifurcation diagram often called stability map,
consisting of several equilibrium curves (dT/dt = 0 and d/dt = 0) (in which an active bifurcation parameter
is varied), Fold curve (border between static stability and static instability) and the Hopf curve (border
between dynamic stability and dynamic instability). The behaviour of the considered system is different
(stable/unstable) at both sides of the bifurcation curves. A bifurcation diagram can provide useful
information about the stability and is therefore called a stability map.
13
Primarily the stability maps have a rather complex structure. Therefore, the following explicating stability
map is constructed. The example provides a clear perceptive of the influence of process parameters like
cooling capacity or coolant temperature and becomes comprehensible but moreover, in a certain
particular case, the actual problem of limit cycles becomes clear. In the stability map, the coolant
temperature Tcool is drawn on the abscissa and the active reactor temperature T on the ordinate. The
22
axes-scale is in accordance with the figures printed in the article by van Elk et.al. resulting in Figure 7 in
which represents a particular process value (in this report the base case).
Ttemperature [K]
510
500
490
480
470
460
450
440
430
420
410
410 420 430 440 450 460
Cooling tem perature [K]
Ttemperature [K]
Figure 7a
Equilibrium curve through

base case value.
represents the base case.
510
500
490
480
470
460
450
440
430
420
410
410 420 430 440 450 460
Figure 7d Addition Hopf curve,

transition between dynamic
stability and unstable.
510
500
490
480
470
460
450
440
430
UA
420
410
410 420 430 440 450 460
Figure 7b Additional equilibrium
Ttemperature [K]
510
500
490
480
470
460
450
440
430
420
410
410 420 430 440 450 460
Figure 7c Addition Fold curve,
curves for varying values of

cooling capacities.
510
500
490
480
470
460
450
No pipes
440
430
420
410
410 420 430 440 450 460
transition between static

stability (right) and unstable
(left).
Ttemperature [K]
Ttemperature [K]
Initially, one equilibrium curve is drawn for one particular active parameter value (cooling capacity)
(Figure 7a).
To examine the influence of the cooling capacity on the process, a number of equilibrium curves is
appended (Figure 7b).
The Fold curve is plotted, in which the boundary between static stability and instability becomes
visible (Figure 7c). The base case is drawn which apparently is located in the static stable
region i.e. no multiplicity in this particular case.
The Hopf curve is plotted, in which the boundary between dynamic stability and instability can be
determined (Figure 7d). The base case is located between the Fold and Hopf curve.
Consequently, its dynamical behaviour is unstable: sustained oscillations i.e. limit cycles.
In Figure 7e the Treactor = Tcool curve has been added. It is obvious that the reactor temperature
should be larger than the coolant temperature. The equilibrium curve corresponding the cooling
capacity for exclusively the reactor wall i.e. no cooling device, is drawn.
Figure 7f portrays clearly the distinct regions. In subsequent stability maps, these coloured
regions will be indicated with I, II or III. Although perhaps more regions can be identified, merely
the most relevant regions will be considered.
Ttemperature [K]
510
500
490
480
470
460
450
440
430
420
410
410 420 430 440 450 460
Figure 7e Addition T = Tcool curve and Figure 7f

equilibrium curve
concerning no cooling
pipes.
Indicating regions. region

I: stable. region II: limit
cycle. region III: instability.
Due to the complex structure of a stability map, it is not clear at once whether a region indicates stability or
instability. Therefore, a particular case has to be verified through producing some curves. In case
multiplicity is involved, a region II can seemingly point towards dynamical instability, which would imply
14
limit cycles. However, sometimes the existence of static instability overrules the existence of a limit cycle,
which means transition and no sustained oscillations. In many articles, this process is a result of a strong
attracter, which forces a process towards transition. Summarising: the labelling procedure of a particular
process regarding dynamical behaviour and stability must be carried out with caution. A check with orbit
curves is indispensable.
The stability map is divided in marked regions:
Region I:
The dynamical behaviour of the system is considered as stable. After a constrained

disturbance, the process variable returns to its steady state value (Figure 4a and b).
Region II: In the area enclosed between the Fold and Hopf curve, the system is dynamically unstable
i.e. limit cycles can exhibit (Figure 4c-e).
Region III: The system is statically unstable due to multiplicity (extinguish or runaway) (Figure 4f and i).
37
Table 6 Various types of behaviour based on linearised perturbation analyses (Himmelblau et.al. ).
Behaviour after perturbation
Area
Figure 4
Asymptotic damping
Spiral point
Stable oscillations, close from Hopf (4.5.3)
Asymptotic oscillations, far from Hopf
Transition
Transition
4.5.6
Ia (point stable)
Ib (point stable)
II (limit cycle)
II (limit cycle)
III (static unstable)
III (static unstable)
a
b
c and e
d
f and g
h and i
Classification of the stability maps

22
In the literature , in general stability maps can be found in which the reactor temperature is plotted
against one or more important process variables. A common plot is the T-Tcool stability map in which the
reactor designer adequately can determine the stable and unstable regions. However, in case of process
control, which involves extra mathematical correlations, this type of stability map is not suitable anymore
to determine stability. This is because too many relevant variables have to be analysed simultaneously.
For proportional-integral control, the process engineer is in particular interested in the relationship
between the proportional gain and the cooling capacity, or the proportional gain and the integral time. If
additionally delay is studied, the T-Tcool stability map is certainly not useful. Therefore, the following
stability map classification is proposed:
1. Reactor designer stability map
2. Reactor controller stability map
1. Reactor designer stability map

In this type of stability map, the reactor temperature is always drawn on the y-axis. At the x-axis one
important process variable, like the coolant temperature is drawn by means of equilibrium curves. The
latter is done, several times in relation with another process parameter, e.g. the cooling capacity, to study
the effect. Through changing the process variables and parameters, the reactor designer can determine
the reactor temperature and the conversion with respect to the behaviour and furthermore the stability.
2. Reactor controller stability map

The reactor control engineer is mainly interested in how to solve and/or prevent unwanted unstable
reactor behaviour. Through creating reactor controller stability maps, one can very easy determine for
instance the new controller configuration. Therefore, the reactor temperature is not drawn on the axes, but
in general one of the crucial process parameters (UA) or the controller parameters (Kc and I). The reactor
controller stability map consist no equilibrium, Fold or orbit curves, merely Hopf curves, in which the stable
and unstable region clearly can be located and indicated.
15
4.5.7
Interpreting stability map
Reactor designer stability map

The same stability map will be considered as in 4.5.5. The regions I, II and II indicating respectively
stability, limit cycles and instability have been marked (Figure 8a). The reactor engineer is for instance
interested in the dynamical behaviour of the process for Tcool = 450 [K] for three different cooling
capacities. These points have additionally been drawn in (Figure 8a) and examined for their dynamical
behaviour through orbit curves.
Table 7
Point
1
1
2
2
3
3
Data regarding the examined points in stability map Figure 8a.

Coolant
temperature
Tcool [K]
450
450
450
445
420
420
Cooling
capacity
UA [kW K-1]
35
35
38
38
25
25
Perturbation
T
Orbit curve
+20
-20
0
0
+10
-10
Figure 9a
Figure 9b
Figure 10a
Figure 10b
Figure 11a
Figure 11b
stst reactor
temperature
T [K]
510
510
505
499
509
401
stst
conversion
[-]
0.94
0.94
0.92
0.90
0.94
004
Region
Behaviour
I
I
II
II
II
I
Stable
Stable
Limit cycles
Limit cycles
Limit cycles
Transition
II
III
Figure 8a Reactor designer stability
700
700
Tem perature [K]
P1
510
500
490
480
470
460
450
440
430
420
410
410 420 430 440 450 460
Temperature [K]
Ttemperature [K]
The steady state values, derived from the equations 2 and 3 (or 8 and 9) have been displayed in Table 7.
As initial point, the steady state value will be taken. If necessary, a perturbation is constrained.
600
500
400
0
400
10800
T = 20 [K] becomes a spiral

point.
3600
7200
Tim e [s]
10800
T = -20 [K]. Point 1 returns

to its steady state (spiral
point).
Point 1 in Figure 8a, is located in region I. From Figure 9a and Figure 9b it can be seen that after
a step disturbance the system returns to P1.
II
III
Figure 8b Reactor designer stability

map. Point 2 is located
beneath the Hopf curve in
region II.
700
700
Tem perature [K]
510
500
P2
490
480
470
460
450
440
430
420
410
410 420 430 440 450 460
Tem perature [K]
Ttemperature [K]
500
Figure 9a Orbit curve for a disturbance Figure 9b Orbit curve for a disturbance
map. Point 1 is located above

the Hopf curve in region I.
P1
3600
7200
Time [s]
600
600
500
400
0
3600
7200
Tim e [s]
10800
Figure 10a Orbit curve. Point 2 in

Figure 8a is located beneath
the Hopf curve in region II
accordingly limit cycles.
600
500
400
0
3600
7200
Tim e [s]
10800
Figure 10b Lowering the coolant

temperature 5 [K] results in
strong limit cycles.
16
Point P2 in Figure 8a is located just beneath the Hopf curve in region II. Figure 10a demonstrate
that limit cycles emerge from the steady state without a perturbation. If the coolant temperature is
slightly decreased Tcool = -5 [K], Figure 10b shows that limit cycles appear strongly.
II
III
Figure 8c Reactor designer stability

map. Point 3 is located above
the Hopf curve in region I,
however with multiplicity.
P3
700
700
Temperature [K]
P3
510
500
490
480
470
460
450
440
430
420
410
410 420 430 440 450 460
Temperature [K]
Ttemperature [K]
P2
600
500
400
0
3600
7200
Time [s]
10800
600
500
400
0
3600
7200
Time [s]
10800
Figure 11a Orbit curve. Point 3 in Figure Figure 11b Point 3 in Figure 8a.
8a is located above the Hopf
curve in region II. T = +10 [K]
limit cycles.
Disturbance T = -10 [K]

results in T = 401 [K], =0.04
reaction extinguishes due to
the existence of multiplicity.
Point P3 in Figure 8a is as point P1 located above the Hopf curve within region I. Therefore,
stability is to be expected. Nonetheless, after a disturbance T = 10 [K], the temperature will
slightly decrease and then pursue the limit cycle which becomes clear in Figure 11a. After a
disturbance T = -10 [K] the temperature will decrease to the lowest steady state i.e. the reaction
-1
-1
extinguishes (Figure 11b). Point P3 with UA = 25 [kJ s K ] exhibits multiplicity and is very
receptive for disturbances. Even a small deviation from the steady state of 10 [K] forces the
system to the lowest steady state.
Although the regions in a reactor designer stability map perhaps point towards a stable dynamical
behaviour, disturbances can force a system towards a dynamical unstable region or even in case of
multiplicity to a lower or higher steady state vale. By means of creating orbit curves, the latter must be
verified.
In the reactor designer stability map in a 2D plot, more than two variables are varied. Therefore, it will be
clear that the reactor designer stability maps are in fact 2D-presemntation of 3D phenomena and is
therefore sometimes difficult to interpret. In Figure 8 T, Tcool, and UA have been varied. Moreover, the
complexity is become worse, if also more ordinary differential equations are added to the mathematical
model.
Reactor controller stability map

In reactor controller stability maps, exclusively the Hopf curves are drawn. Equilibrium curves and Fold
curves like in reactor designer stability maps are excluded. Subsequently, the acquired plots remains well
ordered and clear to interpret. Moreover, in cases like the delay or integral time are considered, no Fold
curves exist. The latter are described mathematically by differential equation in which in the infinity (in the
steady state) the time aspect is eliminated and static instability is then not anymore concerned with. The
Fold curve is then merely a dot (solely one steady state solution of the ODE-system) and does therefore
not make any useful contribution to the stability map. Another explanation is that the equilibrium curves do
not diverge for varying integral times or delay and therefore no Fold curve exists. The advantage of the
reactor controller stability maps is the simplicity, ordered structure and easy to interpret the dynamical
behaviour. Once the Hopf curve has been drawn, the regions with respect to (in)stability have only to be
marked. In most cases, one or two orbit curves for both sides of the Hopf curve provide enough
information about the actual stability.
Consider the following example. The reactor engineer is for instance interested in the dynamical behaviour
of a particular process in relation with the delay. Therefore, a reactor controller stability map is created
with a Hopf curve for parameters Kc and d. In Figure 12, the asymptotic values are visible. Beyond a
certain delay d > 108 [s] apparently the process cannot anymore become stable. If d 0 the minimum
proportional gain can be found, whose value precisely should preserve stability in case delay is not
17
concerned with. It is obvious that delay makes a process less stable because the controller can be too late
with correcting a certain disturbance. Therefore, the right-hand side of the Hopf curve is region II (or III if
multiplicity is involved) and the left-hand side is region I. Large Kc values can eliminate limit cycles and
preserve stability.
II
4
2
Hopf
P1
P2
Temperature [K]
Proportiona l ga in
700
Temperature [K]
700
10
600
500
400
60
120
De la y [s ]
500
400
3600
7200
Time [s]
10800
Figure 12a Reactor designer stability map. Figure 12b Orbit curve. Point 1,
At the left-hand side, the
system behaves stable, at the
right-hand side unstable. Two
orbit curves are needed to
verify the behaviour.
600
behaves stable after a

disturbance T = +20 [K],
spiral point.
3600
7200
Time [s]
10800
Figure 12c Point 2 behaves unstable

(limit cycles) from the initial
steady state. Beyond the
asymptotic delay, instability is
inevitable.
Point 1 shows that after a disturbance, the controller (Kc = 2) is robust enough to reduce the limit cycles,
even with a delay d = 60 [s]. The same Kc = 2 with a delay d = 110 [s], which is slightly larger than the
asymptotic value d = 108 [s] implies definitely an unstable process. Suppose a delay d = 80 [s], the
engineer can effortlessly read from Figure 12a that Kc = 3 signifies no more limit cycles.
4.6 Theoretical prevention and elimination of limit cycles

80
According to Westerterp et.al. the heat capacity of the reactor and in particular in relatively small
reactors, has to be taken into account if limit cycles are to be considered. Westerterp postulated in general
that limit cycles will not occur if, in a high conversion reactor, a combination of the following factors
coincide:
A weak temperature dependence of the reaction rate (E/R low).

A low reactant concentration in the feed (Tad) low.
A low heat effect (Tad) low.
A high reactor temperature (Ts high).
A not too high conversion (1-s relatively large).
A long residence time (R high).
A preferable small reactor (A/VR large).
Good heat transfer coefficient or heat capacity (UA large).
In case more advanced processes are considered, these points are conversely not sufficient enough to
solve a limit cycle appearance and more knowledge about the process is required.
18
Chapter V
5 PROCESS CONTROL
5.1 Introduction
During the process of a chemical reactor, several requirements such as safety, production specifications,
operational constrains and economics must be satisfied, in the presence possible external disturbances.
These requirements dictate the need for continuous monitoring of the operation of a chemical reactor and
external intervention i.e. control, to guarantee achieving the operational objectives. This is accomplished
through a rational arrangement of measuring devices, valves, controllers, computers etc., which
constitutes the control system. There are three general classes of needs, which a control system is
required to satisfy:
1. Suppressing the influence of external disturbances
2. Ensuring the stability of a chemical process
3. Optimising the performance of a chemical process
5.2 Feedback controllers

In every control configuration, the controller is the active element that receives the information from the
measurements and takes appropriate control actions to adjust the values of the manipulated variables.
Consider the generalised process shown in Figure 13. It has an output y, a potential disturbance d also
known as process load and an available manipulated variable m. The disturbance changes in an
unpredictable manner and the control objective is to keep the value of the output y at desired levels.
A feedback control action can take the following steps:
1. Measure the value of the output using the appropriate measuring device. Let ym be the value
indicated by the measuring device.
2. Compare the indicated value ym to the desired value ysp (set point) of the output. Let the deviation
(error) be = ysp ym.
3. The value of the deviation is supplied to the main controller. The controller in turn changes the
value of the manipulated variable m in such a way as to reduce the magnitude of the deviation .
Usually, the controller does not affect the manipulated variable directly but through another device
known as the final control element.
The system in Figure 13 is known as feedback-controlled or closed loop system and pictorially
summarises the foregoing three steps.
d
ysp
c
Controller
Final
control
element
Process
m
Measuring
device
ym
Figure 13 Process and corresponding feedback loop. Ysp = set point or desired value = error d = disturbance m = measured
value c = controlled value.
19
The controller comes between the measuring device and the final control element. Its function is to receive
the measured output signal ym(t) and after comparing it with the set point ysp to produce the actuating
signal c(t) in such a way as to return the output to the desired value ysp. Therefore, the input to the
controller is the error (t) = ysp ym(t), while its output is c(t). The various types of continuous feedback
controllers differ in the way they relate (t) to c(t). Four basic types of feedback controllers can be
distinguished:
1.
2.
3.
4.
Proportional (P controller)
Proportional integral (PI controller)
Proportional derivative (PD controller)
Proportional Integral derivative (PID controller)
Due to the scope of this report, exclusively the proportional and the integral action will be examined.
5.2.1
The principle of the P-controller is the proportional correction of the error. The response has a high
maximum deviation and there is a significant time of oscillation. The period of this oscillation is moderate.
For a sustained change in load, the controller variable is not returned to its original value i.e. the desired
value, but attains a new equilibrium value termed the control point. The difference between desired value
and control point is called the offset or droop. Its actuating output is proportional to the error:
c(t ) = K c (t ) + c s
(25)
In equation 25 parameter Kc is called the proportional gain of the controller and cs = controller bias signal
i.e. its actuating signal when = 0. A proportional controller is described by the value of its proportional
gain Kc or equivalently by its proportional band PB, where PB=100/Kc.. The proportional band
characterises the range over which the error must change in order to drive the actuating signal of the
controller over its full range. Usually (1 PB 50). The larger the gain Kc or equivalently, the smaller the
proportional band, the higher the sensitivity of controllers actuating signal to deviations will be.
Define the deviation c(t) of the actuating signal by
c' (t ) = c(t ) cs
= K c (t )
5.2.2
(26)
Proportional-Integral control
Considerable improvements in the quality of the resulting control can be obtained if a different control law
is used known as proportional-integral control or proportional-plus-reset controller where the proportional
and the integral action are combined. Its actuating signal is related to the error by equation 27:
Kc t
c(t ) = K c (t ) +
(t )dt + c s
I 0
(27)
where I is the integral time constant or reset time often expressed in minutes. The reset time is the time
needed by the controller to repeat the initial proportional acting change in its output. The reset time is an
adjustable parameter and is usually varied in the range (0.1 I 50 [min]). The integral action causes the
controller output c(t) to change as long as an error exists in the process output. Therefore, such a
controller can eliminate even small errors. The maximum deviation of the controlled variable is determined
by the settings of both Kc and I, The integral term of a PI controller causes its output to continue changing
as long as there is a non-zero error. Often the errors cannot be quickly eliminated and given enough time
they produce larger and larger values of the integral term, which in turn keeps increasing the control action
20
until it is saturated (e.g. the valve completely open or closed). This condition is called integral windup.
Then, even if the error returns to zero, the control action will remain saturated. A PI controller needs
special provisions to cope with integral windup.
0
P
Proportional
+ Integral
Kc
P0
Proportional
Integral
Time
Figure 14 Response of PI controller to unit step change in error.
5.2.3
Offset
The question arises: Why does proportional feedback control result in steady state offset unequal to zero?
The proportional controller operates according Equation 26 so that in the transient state, while things are
still changing, the rate at which the process input is being changed is given by:
dc' (t )
d(t )
= Kc
dt
dt
(28)
It is obvious that at steady state, all the derivates will vanish, however is it possible for d(t)/dt to be zero
without (t) itself being zero? The answer, of course, as soon as (t) becomes constant, whether at zero,
or at a nonzero value, steady state is achieved. Observe from equation 26 that for (t) to be zero, c(t)
must be zero. From the process model it is obvious that c(t) will never be zero for non-zero set-point or
non-zero disturbance i.e. there will always be a discrepancy i.e. steady state offset.
Offset = (new set point) (ultimate value of the response)
The offset is the characteristic effect of proportional control. It decreases as Kc becomes larger and
theoretically: offset 0 if Kc .
Another question arises: Why does PI control not result in steady state offset? For the PI controller, control
action is determined according to equation 27. Under transient conditions, the rate at which c(t) changes
is given by equation 29 by differentiating equation 27.
dc' (t )
d(t ) K c
= Kc
+
(t )
I
dt
dt
(29)
Regardless of controller parameters, or the specific process in question, whenever steady state is
achieved, all the derivatives in equation 29 will vanish and (t) will always be zero i.e. no steady state
offset. Therefore, integral action eliminates any offset.
21
5.2.4
Comparison P and PI controller
The advantage of the PI-controller is the offset removing characteristic. However, the disadvantages of PI
control are that it gives rise to a higher maximum deviation, a longer response time and a longer period of
oscillation than with proportional action only. This type of control action is therefore used when the above
can be tolerated and offset is undesirable. It is however a very frequently used combination.
Controlle d va ria ble
No c ontrol
P+I
Off s et
T im e
Figure 15 Response of controlled variable to step disturbance in load.
5.3 Disturbances
Changes in continuous process operation can generally be divided into the following categories:
1. Dynamic operations
2. Internal changes
3. External changes
Dynamic operations
Dynamic operations take place when the process is not operated under constant conditions e.g. during
start-up or switchover to other conditions.
Internal disturbances
Due to changes in one or more state variables like T and or , without any external disturbances, the
process can become dynamical unstable (limit cycles). Drastic internal changes can be caused by failures
in process equipment or control instrumentation in which the state variable change. In addition a process
can become gradually dynamical unstable by the slow change of conditions of a continuous process by
35
fouling, poisoning of catalyst, coke depositing in furnace tubes etc. Henson et.al. .
External disturbances
For this type of change, continuous operation is disturbed by external influences. In many cases, it is a
matter of disturbances, which enter the process by feed, heating, catalyst flows etc. Here disturbances
mean relatively small changes in process conditions without malfunctioning of process apparatus and
control instruments.
22
Process behaviour
Because all disturbances from internal and external causes are undesirable at constant operation, the
process has to be corrected (process control to be examined in chapter 5). These disturbances or
perturbations cause changes in the dynamical behaviour of a process.
5.4 Concept of dead time

Whenever an input variable of a system changes, there is a time interval during which no effect is
observed on the outputs of the system. This time interval is called dead time or delay. This response
poses a difficult problem for a feedback controller, because its next move depends on the response of the
process to its most recent move. The presence of dead time can easily destabilise the dynamic behaviour
of a controlled system to an external perturbation.
f(t)
y(t)
Process
y = t - td
Dead
time
Figure 16 Process with dead time.
yout (t ) = yin (t t d )
(30)
When a time-delay element is present in a feedback control loop, the following obvious control problems
arise:
With a delay, control action will be based on delayed, hence obsolete process information that is
usually not representative of the current situation within the process.
If the process has an input delay, then the effect of the control action will not be immediately felt
by the process, compounding the problem even further.
Virtually all physical processes will involve some dead time between the input and the output. The passing
through the core of a measuring device, the controlling process, the transfer of heat from a heating device
to the reactor etc., will take time.
To avoid the use of complex mathematical descriptions of every process, which contributes to delay and
due to the lack of physical information regarding dead time, the following overall equation can be used:
dyd
= y yd
dt
(31)
In equation 31, all the different kind of delay is allowed for by means of the parameter d.
5.5 Process capacity

Capacity is where a process stores variable amounts of mass or energy. If the flowrate of mass or energy
into and out of a process are not equal, then their difference accumulates within the capacity of the
1
process i.e. the process act as a buffer. If a variable x is a measure of that accumulation, the following
general differential equation can be formulated:
Liquid level is an indication of the variable volume of liquid stored in the capacity of a tank, temperature is a measure of the
energy contained in the heat capacity of a mass of material and stream composition is also indicative of a particular species
accumulating in a mixed vessel such as a reactor or mass transfer process.
23
dx
= in out
dt
5.5.1
(32)
Process time constant
In equation 32 in and out are inflow and outflow respectively and the time constant of the capacity in
the same units as time. The time constant of a process is a measure of the time necessary for the process
to adjust to a change in its input. The smaller the value of , the steeper the initial response of the system.
Consequently, capacity plays a certain role in the dynamical behaviour and stability of a process, which
74
will be examined in next sections. According to Stephanopoulos the CISTR can be observed as a multi
capacity process of mass and energy. By dividing differential equations 2 and 3 with CP V, equations
33 and 34 are obtained.
VR
V
VR
V
Eact
V
d[ A]
= ([ A]0 [ A]) R k 0 e RT [ A]
V
dt
V ( H R ) Eact RT
dT
UA
(T Tcool )
= (T0 T ) + R
k0 e
[ A]
V C P
V C P
dt
(33)
(34)
It is obvious that in this particular case the time constant represents the average residence time R of the
CISTR according to:
VR
V
R =
(4)
If one takes into account also the heat capacity of the equipment present in the reactor (i.e. stirrer etc.) the
equation can be extended to
R =
VR mequipmentC P,equipment
+
+ ...
C P V
V
5.5.2
(35)
Cooling time constant
If equation 12 is divided by cool CP,cool V,cool, than equation 36 is obtained:
cool
dTcool
UA
= (Tcool,0 Tcool )+
(T Tcool )
coolCP,coolVcool
dt
(36)
In which cool is:
cool =
Vcool
V ,cool
(37)
In fact, cool is partially composed of the time constant of the coolant and additionally partially composed of
the physical properties of the equipment. Analogue to equation 35 one can write the coolant time constant
correlation 38 according:
24
cool =
m
C
Vcool
+ equipment P,equipment + ...
V ,cool coolCP,cool V ,cool
(38)
5.6 Controller tuning

Generally, when a set point changes, the response of the process deviates and the controller tries to bring
the output again close to the desired set point (Figure 15). The distinguished controllers have different
effects on the response of the controlled process. The following questions arise:
What type of feedback controller should be used to control the process?

What are the best values for the adjustable parameters of the selected controller?
What performance criteria should be used for the selection and the tuning of the chosen
controller?
5.6.1
Performance criteria
Even though specific details of what is considered acceptable performance can be interpreted throughout
a different perspective, some general principles can be applied universally. An effective close loop system
is expected to be stable and to be capable of causing the system output ultimately to attain its desired set
point value. In addition, the approach of this system output to the desired set point should be neither too
sluggish, nor too oscillatory. A careful examination of criteria by which closed-loop system performance
may be assessed in general:
Stability criteria
Steady state criteria
Dynamic response criteria
Of these, the first two are very easy to specify. The only reasonable specification is that the system must
be stable. Also the key steady state criteria is that there be little or no steady state offset i.e. the error is
brought close to zero at steady state. The third class of criteria specifies how the system responds under
closed-loop control. The evaluation of the dynamic performance of a closed loop system is based on the
following types of commonly used performance criteria:
Overshoot
Rise time
Settling time
Decay ratio
Windup
Keep the maximum deviation (error) of the response as small as possible.

Check the time needed to reach the desired value for the first time.
Return to the desired level of operation as soon a possible (or within 5%).
Verify if the ration between the first and the second peak is approximately 4:1.
Minimize the integral of the errors until the process has settled to its desired set
point and prevent integral windup.
All of the mentioned characteristics above could be used by the designer as the basic criteria for selecting
the controller and the values of its adjusted parameters. It must be emphasised, though, that one simple
characteristic does not suffice to describe the desired dynamic response. Usually, more satisfied
objectives are required i.e. minimum overshoot, minimum settling time etc. Unfortunately, controller
designs based on multiple criteria lead to conflicting response characteristics.
25
5.6.2
Methods of adjusting feedback controller settings
Many procedures exist for estimating optimum settings for controllers:

74
1. Ratio-overshoot method (Stephanopoulos )

One of the usual bases employed is that the system response should have a decay of i.e. the ratio
of overshoot of the first peak to the overshoot of the second peak is 4:1 (Figure 17). There is no direct
mathematical justification for this but it is a compromise between a rapid initial response and a short
response time. The response time is the time required for the absolute value of the system response
to come within a small-specified amount of the final value of the response.
81
2. Ziegler and Nichols method (Ziegler et.al. )
This method is often applied and will be discussed in the next paragraph.
3. Process response methods
25
11
Common methods are e.g. the Loop tuning method by Erickson et.al. , the Cohen Coon method ,
56
the Integral relation method by Murrill and many more which are not suitable for the propose of this
report due to the fact that a process response is practically obtained.
4. Remainder tuning methods
74
35
Other tuning methods are the Nyquist stability criteria and non-linear process control methods .
Many more methods are available in the literature.
Syste m re sponse
a:b = 4:1
Finbal level value

of res pons e
Res pons e tim e
T ime
Figure 17 System response with decay ratio.14 25 74 60.
5.6.3
Ziegler and Nichols method
This method is characterised by finding the gain at which the system is marginally stable and the
frequency of oscillation at this point. From these two parameters, the controller parameters are calculated.
81
The Ziegler Nichols tuning technique is intended to produce a closed-loop damping ration of 1:4.
It goes through the following steps:
1. Bring the system to the desired operational level i.e. design condition.
2. Using proportional control only (or with maximum I) and with the feedback loop closed, introduce
a set point change and vary the proportional gain until the system oscillates continuously i.e. limit
cycle. The frequency of continuous oscillation is the crossover frequency . Let M be the
amplitude ratio of the systems response at the crossover frequency.
3. Compute the following two quantities:
The ultimate gain or proportional gain for sustained oscillations:
26
KU = 1
(39)
The ultimate period of sustained cycling [minutes]:
PU =
(40)
4. Using the values of KU and PU, Ziegler Nichols recommended the following settings for feedback
controllers:
Table 8 Ziegler Nichols controller settings.
Control action Controller settings
Kc
I [min]
P
KU / 2
PI
KU / 2.2
U / 1.2
The settings above reveal the rationale of the Ziegler Nichols methodology.
1
2
For proportional control alone, a gain margin equal to 2 is recommended.

For PI control, a lower proportional gain is advised because the presence of the integral
control mode introduces additional phase lag in all frequencies with destabilizing effects on
the system. Therefore, lower Kc maintains approximately the same gain margin.
5.6.4
Controller configuration
74
The remaining question is, which one to select first: Kc or I? According to Stephanopoulos in general the
proportional gain is selected first, in such way that the controller has the necessary strength to repress
disturbances. Afterwards a suitable integral time is chosen for offset elimination. Together Kc and I must
satisfy the control performance criteria e.g. the decay ratio.
In case a system is not exclusively concerned with external disturbances but additionally is concerned with
dynamical instability, the mentioned tuning methods have to be interpreted with caution and it is not
certain if these methods are still valid. Inevitably, a more robust control action is needed i.e. Kc has to
increase. This might have sincere consequences for the e.g. Ziegler Nichols controller settings. Therefore,
the following strategy is suggested to acquire stability through searching and optimising controller
configuration parameters.
system
parameters
determine Kc
process
Condition:
Stability
Ratio a/b = 4/1
determine I
results stability
Condition:
Fast: offset 0
I >> d
Figure 18 Strategy to acquire stability through searching and optimising controller parameters.
27
Chapter VI
6 STARTING POINT OF RESEARCH
6.1 Introduction
In this chapter, the current level of research of limit cycles is examined. Therefore, two distinctive reactor
designer stability maps and one reactor controller stability map will be considered conjoined with
accompanying orbit curves. To apply the theory, a certain system has to be chosen. In this report, a
22
system previously studied by Van Elk et.al. , is selected. Primarily, in the next paragraph, this system will
be quantified. Subsequently, a method to study the stability will be presented. Then, the resulting stability
maps will elaborately be presented and discussed.
6.2 Base case definition

Despite that the numerical values are taken from the article presented by Van Elk, the system studied in
this report does differ from the system studied by van Elk. He studied the stability and dynamic behaviour
of gas-liquid reactors extensively. Adopting a rigorous gas-liquid reactor model, he demonstrated the
possible existence of dynamic instability (limit cycles) in gas-liquid processes. Since limit cycles can
already emerge in a system that can be described with merely two differential equations i.e. mass and
energy balance, some adaptations can be been made to simplify the model. Primarily, non-relevant
aspects with regards to limit cycles have been removed from the model. The simplest case will be
considered because more detailed issues, through adding more mathematical equations to the model, can
be concerned with afterwards.
The first simplification for the base case model is that the reactor of interest in this report is a
homogeneous liquid phase reactor. Therefore, all the parameters needed for the description of two-phase
reactors and mass transfer (like hold-up and the Hatta number) can be ignored. Because van Elk used a
pseudo-first order reaction constant approximation, therefore the same constant k0 can be used for the
first order reaction of the present case. It is assumed that the liquid is completely saturated with
component A. As a result the concentration of reactant A, can be considered as constant and accordingly
be included in the value of k0 (Appendix 3).
The simplified model based on van Elk will be called the base case. The required values of the
parameters are listed in Table 9.
22
Table 9 Main system parameters (van Elk et.al. ).

Parameter
Symbol
Value
Unit
v
VR
[A]0
[P]0
T0
0.005
5
5000
0
303
[m3 s-1]
[m3]
[mol m-3]
[mol m-3]
[K]
Heat transfer coefficient (transfer) area or cooling capacity

Reactor temperature
Coolant temperature
Fixed parameters
UA
T
Tcool
55
468
441
[kJ s-1 K-1]

[K]
[K]
Density
Heat capacity
Reaction enthalpy
Activation energy
Arrhenius frequency factor
CP
HR
Eact
k0
800
2
-160
90
2.505 107
[Kg m-3]
[kJ kg-1 K-1]
[kJ mol-1]
[kJ mol-1]
[s-1]
System conditions
Liquid flow rate
Reactor volume
Feed concentration
Feed temperature
Important base case variables
28
Additional data
Density coolant
Heat capacity coolant
Feed temperature coolant
cool
CP,cool
Tcool,0
1000
4.2
303, 400, 435
[Kg m-3]
[kJ kg-1 K-1]
[K]
If the mass and heat balance are solved for the steady state situation using the values portrayed in Table
9, a small dissimilarity can be noticed: T = 466 [K] instead of T = 468 [K]. Because the latter should imply
an initial perturbation of T = 2 [K], if necessary the exact numerical solution of the mass and heat
balance is used for the simulations in this report. In Appendix 3, the magnitude concerning the base case
values presented in Table 9 has been evaluated.
6.3 Stability of the base case

The following stability naps will be considered:
Table 10 Base case stability maps.
Stability map
Reference
T versus Tcool varying UA
Figure 19
Figure 20
T versus Tcool varying V
Tcool versus UA varying V Figure 21
Although, in fact Figure 21, which will be discussed later, is sufficient to prove the unstable character of
the base case, likewise, the traditional stability maps Figure 19 (and Figure 20) will be presented. Similar
22
stability maps like Figure 19 have been printed in the article by van Elk et.al. . To be able to compare the
stability maps of this report easily, the axes-scale is as much as possible in accordance with the figures
printed in the article.
6.3.1
Reactor designer stability maps
The following stability maps (theory discussed in 4.5) will be considered: Reactor temperature versus
coolant temperature for initially (UA variation) and afterwards (V variation). The mathematical base case
model discussed in 2.2 states:
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
VR
(2)
(3)
Reactor temperature versus coolant temperature (UA variation)

The first base case designer stability map is created (Figure 19) in which the coolant temperature is drawn
on the abscissa and the active reactor temperature on the ordinate. The stability map consists of
equilibrium curves (dT/dt = 0 and d/dt = 0) as a function of one active parameter, which is in this case the
coolant temperature. Several equilibrium curves have been plotted, for various values of UA with the
intention that the influence of the cooling capability becomes more comprehensible. Additionally the Fold
and Hopf curves are included as a function of two active parameters, which are the coolant temperature
and the cooling capacity UA. Eventually the stability map is divided into marked regions. This marking
process, of the particular regions in the stability map, is sometimes easy and sometimes vague due to the
existence of multiplicity. Then it is a matter of trial and error. However, through creating orbit curves, the
dynamical behaviour can be determined or verified.
29
UA=10
510
Hopf
UA=15
500
490
Reactor temperature [K]
UA=20
II
UA=25
480
UA=30
Hopf
470
Base case
Fold
UA=35
UA=40
UA=45
460
UA=50
450
UA=55
UA=60
440
III
430
UA
UA=65
UA
UA=70
420
410
410
420
430
440
450
Hopf
UA=75
Fold
Fold
460
Coolant temperature [K]
470
Hopf
Tcool = Tr
Figure 19 Reactor designer stability map for the base case model considering the irreversible exothermic reaction A P in a
cooled CISTR. The steady state temperature is plotted as function of the coolant temperature in which the cooling
capacity is varied. The symbol points towards the base case value. The LOCBIF source LOCBIF rhs 2 is included
in Appendix 7. The base case is located in a region II which points towards dynamical instability.
Discussion of the reactor designer stability map:
In Figure 19 a dotted line T=Tcool has been drawn which indicate that the reactor temperature
cannot decrease beneath the coolant temperature curve.
LOCBIF provides the bifurcation maximum and minimum values of the cooling capacity and the
-1 -1
coolant temperature in the Hopf and Fold curve. In case: UA > 45 [kJ s K ] and Tcool = Tcool (base
case) or Tcool > 434 [K] and UA = UA (base case), no multiplicity will occur.
In Figure 19 the transition from the region in which multiplicity occurs and the region with one
operation point is marked with an arrow. .
-1
-1
Absolutely no limit cycles exists in Figure 19 provided that UA > 105 [kJ s K ] (Hopf maximum)
irrespective of the coolant temperature.
It is obvious that the relevant base case, symbolised with , is located in region II in which
undesired limit cycles can emerge.
V variation)
Reactor temperature versus coolant temperature (
In the following designer stability map, instead of the cooling capacity, the throughput is considered.
Firstly, the equilibrium curves are drawn as a function of one active bifurcation parameter Tcool for
increasing flowrate values. Subsequently, the Fold and Hopf curves are added as a function of two active
parameters Tcool and V.
30
510
500
Fold
Fold
II
Hopf
Hopf
Tcool =
Treactor
v
Qv=0.001
490
480
Base case
v
Qv=0.002
470
v
Qv=0.003
460
v
Qv=0.004
II
450
III
v
Qv=0.005
Hopf
v
Qv=0.006
440
Fold
v
Qv=0.007
430
v
Qv=0.008
420
410
410
v
Qv=0.009
Equilibrium
420
430
440
450
460
470
Figure 20 Reactor designer stability map of the base case model considering the irreversible exothermic reaction A P in a
cooled CISTR. The steady state temperature is plotted as function of the coolant temperature in which the flowrate
value is varied. For V > 0.006 [m3 s-1] (Fold maximum) the equilibrium curve exhibit multiplicity which means that
extinction is also possible.
If the dotted line T=Tcool in Figure 20 is taken as the reference point and the flowrate is increased
3 -1
to V = 0.0014 [m s ], the Hopf curve is intercepted at Tcool = 434 [K] which is according to
LOCBIF a Hopf maximum. Up to this flowrate, the system is stable. Larger flowrate values imply
dynamical instability.
If the flowrate is increased, the equilibrium curves increase virtually counter clockwise up to V =
3 -1
0.006 [m s ] from where multiplicity is involved.
3 -1
If V > 0.007 [m s ] the system is unstable because the reaction extinguishes. Apparently too
much heat is withdrawn from the CISTR.
LOCBIF provides the Fold and Hopf maximum and minimum values. According to LOCBIF two
multiplicity transitions exist close to one another (Tcool = 439 [K]) i.e. no multiplicity stands for a
3 -1
flowrate V < 0.006 [m s ]. The existence of multiplicity can be clarified by the fact that in case the
flowrate is considerably increased, this implies that on the one hand, the reaction increases due to
the supply of fresh feed, on the other hand the reactor is also cooled by the cold feed inlet. If for
instance too much cold feed enters the reactor, it can conversely remove too much reaction heat
causing the reaction to slow down and eventually extinguish.
LOCBIF also found the minimum Fold in which certainly no limit cycles occur: Tcool > 464 [K]. The
base case indicated with symbol is obviously also in this stability map located in region II in
which certainly limit cycles exhibit.
6.3.2
Reactor controller stability maps
In reactor controller stability maps, not the reactor temperature but a process controller parameter is
concerned with. Because in this chapter the process controller is not present, the stability map is adapted
and examined for as well UA, Tcool as V.
31
110
vF=0.001
III
100
vF=0.002
Cooling capacity [kJ/s/K]
90
vF=0.003
80
II
70
vF=0.004
Base case
vF=0.005
60
vF=0.006
50
v
40
vF=0.007
vF=0.008
30
vF=0.009
20
vF=0.01
10
410
420
430
440
450
460
470
Figure 21 The stability map concerning the base case. Select one particular flowrate. The inside of the loops (Hopf curves) point
towards limit cycles, the left hand side of the curves indicate multiplicity and the right hand side means stability.
Discussion of the stability map:
In Figure 21 the Hopf curves have been drawn for active LOCBIF parameters UA and Tcool.
The stability map provides in one plot, useful process information with regards to as well the
cooling capacity, coolant temperature as well as the throughput.
The inside of all the loops point towards limit cycles, the left hand side of the curves indicate
multiplicity, often this means instability and the right hand side means stability.
A stable operating point can easily be located. For instance, consider the base case with a fixed
UA and V, Tcool shifts to region I in case Tcool > 464 [K].
6.4 Overview
Table 11 Base case relevant stability values (no control). Static and dynamical behaviour.
System variable Cooling capacity Coolant temperature Throughput
Behaviour
Dimension
variation
variation
variation
Cooling capacity
[kJ s-1 K-1]
Coolant temperature
[K]
Throughput
[m3 s-1]
15 < UA < 45
45 < UA < 78
UA > 78
441
441
441
0.005
0.005
0.005
55
55
55
410 < Tcool < 434
434 < Tcool < 464
Tcool > 464
0.005
0.005
0.005
55
55
55
441
441
441
0 < V < 0.0014
0.0014 < V < 0.006
V > 0.007
Multiplicity
Limit cycles
Stability
Multiplicity
Limit cycles
Stability
Stability
Limit cycles
Instability
32
6.5 Base case limit cycles

Primarily, the orbit curve for the base case reactor temperature is created with LOCBIF, which results in
Figure 22a and the orbit curve pertaining to the conversion in Figure 22b.
The mathematical steady state solution of the mass and heat balance (ODE 2 and 3) has been chosen as
initial point, which is represented by the dotted line. The abscissa of Figure 22a and Figure 22b are
identical. The figures conversely can be combined to phase plot Figure 22c in which the principle of the
phenomenon limit cycle can be elucidated. No perturbation is necessary to initiate the self-sustained
oscillations. The system variables (conversion and the reactor temperature) will cycle around the steady
state point however the steady state itself will never be reached.
1.0
700
Conversion [-]
750
650
600
550
500
450
0.8
0.6
0.4
Initial = steady
state
0.2
0.0
400
0
3600
7200
10800
Time [s]
7200
10800
Time [s]
Figure 22a Orbit curve for the exothermic reaction A P in a
Figure 22b Orbit curve for the base case definition in case limit
CISTR. The initial condition is the steady state

situation. The limit cycles occur consistently with a
temperature deviation T = 284 [K] in a cycle time
t 34 [min].
cycles emerge. The conversion deviation due to

the limit cycles appears to be 0.6 every
t 34 [min].
Table 12
1.0
0.8
Conversion [-]
3600
0.6
The base case value (from Table 9).
Variable
Value
UA
Tcool
V
Tstst
stst
55
441
0.005
466
0.68
Dimension
[kJ s-1 K-1]

[K]
[m3 s-1]
[K]
[-]
Initial =
steady state
0.4
0.2
0.0
350
450
550
650
750
Figure 22c After approximately t = 8 [min], the maximum

conversion has been reached i.e. 0.98 and
decreases to 0.40 and so on. Subsequently, the
system is proverbially dragged into the limit cycle.
33
Chapter VII
7 STRUCTURE AND OBJECTIVES OF PART II
Part II, will have the following elements in every chapter.
It will make the difference in two kinds of stability-maps:
1. Reactor designer stability maps: These are indispensable when designing a new reactor. In this
case reactor- and coolant-temperature, cooling capacity, coolant flowrate and throughput
dependence of the process stability are indispensable parameters for the designer to develop a
new reactor. Due to the number of relevant parameters and their dependence, this kind of stability
maps have a complex nature and are often difficult to interpret.
2. Reactor controller stability maps: These are useful simplifications of the reactor designer stability
maps which focus on a particular case in which reactor- and coolant-temperature, cooling
capacity, coolant flowrate and throughput are already fixed (no degrees of freedom) and stability
only can be attained by applying process control.
The reactor controller stability maps are an extension of the use of data generated with LOCBIF and are a
major improvement in handling stability data and obtaining a profound idea of stability-dependence of the
base case (or any other existing or hypothetical possible reactor).
Therefore the objectives of part two are:
1. Providing indispensable information for the design of a new reactor
2. Providing a practical tool to configure the controller in which reactor stability can be guaranteed
The base case model will initially be examined with respect to its dynamical behaviour in the absence of a
controller. Subsequently the control system is implemented in the model. The control method can be
subdivided into the proportional and the proportional-integral control action. Although the derivative control
action is widely used in the process control, this type of control is not suitable for the use in particular
bifurcation software packages and will therefore not be considered. Three different control methods have
been selected and investigated on its influence on the dynamical behaviour. The system can be controlled
adjusting the temperature of the cooling fluid, adjusting the coolant flowrate or the throughput.
The next step is to analyse every situation for sensitivity of relevant parameters respectively
1.
2.
3.
4.
The cooling capacity UA

The proportional gain Kc
The integral time I
The presumed delay d
For as well the control of cooling temperature as coolant flowrate and throughput, the following will be
presented:
1. Mathematical system will be given
2. Reactor designer stability maps will be given
3. Reactor controller stability maps will be given (to monitor the improvement of the base case due to
process control)
4. Conclusion concerning the being of the base case (this being changes as a result of applied
process control)
34
Chapter VIII
8 PRELIMINARY CONSIDERATION
CONTROLLABILITY
8.1 Controlling the coolant temperature
In paragraph 5.2 the default controller correlation has been given. In case of coolant temperature control
the controller correlation becomes:
Tcool = Tcool, sp + K c (Tsp T )
(41)
This mathematical description is easy to comprehend. If the process is disturbed and the reactor
temperature should for instance rise, it is obvious that the coolant temperature has to be decrease to
return the reactor temperature to its steady state (set point) value. The magnitude of the proportional gain
amounts to Kc 1 [-].
8.2 Controlling the coolant flowrate

In case of flowrate control, a reactor temperature increase has the consequence that the coolant flowrate
has to be enlarged to increase the heat transfer from the reactor. Therefore, to maintain positive Kc
values, the controller equation becomes:
V ,cool = V ,cool,sp K c (Tsp T )
(42)
The proportional gain Kc in equation 42 is not dimensionless as in equation 41. The magnitude of the
3 -1 -1
proportional gain value states Kc 0.0001 [m s K ]
8.3 Controlling the throughput

In case of throughput control, there are two possibilities available: Suppose the reactor temperature
increases, on the one hand increasing the throughput means quenching with the cold feed, but on the
other hand increasing the reactant concentration in which subsequently more heat will be produced due to
the increased conversion.
V = V ,sp + K c (Tsp T )
(43a)
V = V ,sp K c (Tsp T )
(43b)
Both controller equations are possible. Nevertheless, it is uncertain which one can provide stability. If the
proportional gain is too large, the reaction can run away or extinguishes dependent on the Kc sign. This
will be clear in the following chapter.
3 -1
-1
The proportional gain value in equation 43 roughly is Kc 0.001 [m s K ] which is slightly larger than in
case of coolant flowrate control.
35
Chapter IX
9 PROPORTIONAL CONTROL
In this section, the effect of the proportional controller on the stability of the cooled CISTR with irreversible
exothermic reaction is studied.
9.1 Reactor designer stability maps

The following stability maps will be examined:
Table 13 Stability maps P-controlled coolant temperature.
T versus Tcool Figure 23
T versus UA
Figure 24
Table 14 Stability maps P-controlled coolant flowrate.
Stability map
Tcool,0 = 303 [K]
Tcool,0 = 400 [K]
Figure 41
T versus V,cool Figure 28
3 -1
3 -1
0 0.02 [m s ]
0 0.005 [m s ]
T versus Tcool
Figure 43
T versus UA
Figure 44
Figure 45
Tcool,0 = 435 [K]

Figure 42
3 -1
0 0.1 [m s ]
Figure 49
Table 15 Stability maps P-controlled throughput.

Figure 50
T versus V
(Kc 0)
Figure 54
T versus V
(Kc 0)
T versus UA
Figure 67
9.1.1
Controlling the coolant temperature

74 71 55 59
. The main and deciding

In the industry, the concept of master-slave control is often applied
process variable (master) is often the yield of the reaction. Although eventually this quantity is deciding for
the complete process, the determination of the magnitude is often very dilatory and therefore in general
not the active manipulated controller variable, which preserves stability. Mostly, a CISTR with exothermic
reaction is controlled through adjusting the degree of cooling (slave) which response is considerably
faster. The coolant temperature in heat balance 3 expresses together with the cooling capacity the extent
of cooling. The value of the cooling capacity UA is mostly fixed and determined during the design of the
process. Commonly, cooling is realised through the application of cooling pipes or a cooling jacket. The
coolant is pumped into the cooling device in which heat is transferred from the reactor content through the
wall towards the coolant. To describe the temperature of the cooling fluid accurate, information is needed
about the flowrate pattern and mixing behaviour, which in many cases is not available or incomplete.
64
75
According to Ratto et.al. and Uppal et.al. the assumption is made that the coolant temperature is
constant and can be constrained virtually instantaneously. This assumption is permitted in case the
coolant flowrate is considerably large i.e. the average residence time of a coolant particle coolant is rather
short compared to the average residence time R of a particle in the CISTR:
R >> coolant
(44)
36
1. Reactor temperature versus cooling temperature
Mathematical model
Equations 2, 3 and 41 mathematically describe the proportional controlled base case in which the coolant
temperature is manipulated.
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
VR
(2)
(3)
(41)
A designer stability map is created (Figure 23) in which the coolant temperature (set point) and the
proportional gain are the active LOCBIF bifurcation parameters.
510
Fold
Hopf
Fold
II
500
Hopf
490
Kc=
Kc
480
Kc
470
460
Kc=1
Kc=0.9
Kc=0.5
450
440
Kc=10
Base case
II
Kc=0.2
Kc (max)=0.9
Fold
430
Kc=0.1
Hopf
Kc=0
420
Tcool =
Treactor
Equilibrium
410
410
420
430
440
450
460
470
Coolant temperature (set point) [K]
Figure 23 Reactor designer stability map of the coolant temperature proportional controlled base case model. The steady state
reactor temperature is plotted as function of the coolant temperature (set point) in which the gain value Kc is varied.
If Kc > 0.9 no limit cycles will exhibit. The LOCBIF source LOCBIF rhs 3 is included in Appendix 7.
In case the proportional gain is increased, the equilibrium curves rotate clockwise and switch from
region II towards stable region I in case Kc > 0.9.
Considered the base case coolant temperature, multiplicity does not occur over the examined
range (410 Tcool < 470 [K]), which means that static instability will not take place, merely
dynamical instability i.e. the occurrence of limit cycles.
37
In Figure 23 the dotted line T=Tcool has been drawn which indicating the lowest possible reactor
temperature.
Also previous stability map (Figure 19) demonstrated that in case Tcool > 469 [K] no limit cycles
can appear. With regard to stability map (Figure 23) no limit cycles will emerge in case Tcool < 430
[K] and Tcool > 469 [K] which are the acquired Hopf maximum values.
Table 16 Hopf maximum of Figure 23.

Variable Value Dimension
Kc (max)
Tcool (min)
Tcool (max)
0.9
430
469
[-]
[K]
[K]
2. Reactor temperature versus cooling capacity

In this section, it will be demonstrated that increasing cooling capacity improves the stability of the
process.
510
Equilibrium
II
500
490
Hopf
Fold
Fold
Kc=
Hopf
480
Kc=10
470
Base case
Kc
Kc=2
Kc=1
Kc
460
450
Fold
440
II
Kc=0.9
Kc=0.5
Kc=0.2
Tcool
Hopf
430
Kc=0.1
Kc=0
420
Tcool
410
10
30
50
70
90
110

Figure 24 Reactor designer stability map of the proportional controlled base case model. The steady state reactor temperature
is plotted as function of the cooling capacity in which the gain value Kc is varied. The base case shifts from region II to
stable region I for Kc > 0.9.
The subsequent stability map Figure 24, with active LOCBIF bifurcation parameters: UA and Kc,
-1 -1
confirm that if UA > 78 [kJ s K ] in view of the base case coolant temperature Tcool = 441 [K],
absolutely no limit cycles will occur. Compare this fact with the results retrieved from stability map
-1 -1
Figure 19 in which for any arbitrarily Tcool a cooling capacity UA > 105 [kJ s K ] will prevent limit
cycles completely.
A significantly large cooling capacity improves the overall stability of the process. This has been
80
confirmed by Westerterp et.al. (4.6).
Figure 23 and Figure 24 confirm that no limit cycles will occur if Kc > 0.9 (Figure 27).
Nevertheless, too large UA forces the reactor temperature towards the coolant temperature
38
because the heat removal becomes significant. Therefore, a considerable cooling capacity
combined with a suitable proportional gain can maintain sufficient stability. In this report it is not
addressed whether or not a considerable cooling capacity is practically achievable. However, if
the implementation of a suitable controller is conceivable, less cooling area A is needed which can
yield a cost reduction. The choice of UA and Kc is up to the reactor designer and controller.
Appendix 3 provides information about the required cooling pipes in case this cooling device is
applied.
For the reactor designer, the stability map shows for one particular Kc in what extent the process
can be changed and what the consequences are for the stability. Take for instance Kc = 1 which
-1 -1
means stability. If the cooling capacity is decreased to for example 40 [kJ s K ] then the process
is located in a dynamical unstable region.
Orbit curves
Orbit curves have been made for three increasing proportional gain values. Figure 25 shows that a relative
small Kc value hardly has any effect on the appearing limit cycles. Increasing Kc slightly dims the
overshoot (Figure 26) and finally increasing Kc will eradicate the limit cycles (Figure 27). The orbit
simulations demonstrated that actually for a Kc = 0.9 the limit cycle is converted into a spiral point
(consistent with Figure 27) which means stability.
Table 17 Data regarding orbit curves for increasing Kc.
Proportional gain Average reactor temperature Average conversion
Kc [-]
T [K]
[-]
467
467
468
0.680
0.686
0.688
600
500
400
0
3600
7200
Time [s]
10800
700
Temperature [K]
700
Temperature [K]
Temperature [K]
0.5
0.8
0.9
600
500
400
0
3600
7200
Time [s]
10800
Figure
Figure 25
Figure 26
Figure 27
700
600
500
400
0
3600
7200
Time [s]
10800
Figure 25 Small Kc values do not affect Figure 26 The proportional gain is not Figure 27 Kc = 0.9 determined from
the emerging limit cycles
tremendously compared to
not controlled base case
according Figure 22.
robust enough in which limit

cycles still emerge.
Figure 23 can hardly

eliminate limit cycles.
Moreover, the average conversion (Table 17) is slightly increased in case Kc rises. Figure 25 with the
strongest limit cycles has a conversion approximately 1% lesser than the conversion of Figure 27.It is not
clear if this is the result of inaccuracy of the LOCBIF simulations or this is due to the strong oscillations.
Summarising
-1
A considerable large cooling capacity can contribute to the decline of limit cycles. In case UA > 78 [kJ s
-1
K ] no limit cycles will emerge at all for the base case. Nevertheless, very large UA values push the
reactor temperature towards the coolant temperature. Therefore, a considerable cooling capacity
combined with a suitable proportional gain is required to preserve a stable system.
A proportional controller, which affects the coolant temperature, assuming that the coolant temperature
can be instantaneously constrained, can easily eliminate limit cycles. In this particular situation, it counts
that the height of Kc determines the stability. The stability maps provide the information that for the specific
base case situation a proportional gain: Kc > 0.9 is suitable to eliminate the self-sustained oscillations. For
any other cooling capacity a Kc = 5.7 is enough for an absolute stable base case.
39
9.1.2
Controlling the coolant flowrate
In many industrial processes, reactors in which exothermic chemical reactions take place are cooled
continuously. The produced heat is removed from the reactor through for instance a cooling jacket or
through a composition of cooling pipes. This section is concerned with the question if limit cycles can be
prevented by applying a controller, which manipulates the cooling fluid flowrate. Although various
procedures exist in which a reactor can be cooled e.g. cooling water, pressurised steam, reflux vaporiser
etc., the objective of this report is to study the fundamental principles of the cooling effect on the CISTR
stability, not the practical applied method. This implies that matters as Hvap, are not concerned with if for
instance heated water under pressure is assumed as cooling fluid.
The coolant heat balance 12, which describes the coolant flowrate as function of the coolant temperature
will be implemented in the existing mathematical base case model. The use of the cooling differential
equation is justified if the coolant temperature is considered constant in the cooling device. If a liquid with
no vapour is used in combination with a recycle stream, the latter assumption is justified. A major recycle
stream means that a cooling fluid particle, which carries transferred heat, can be transported and replaced
by the cold coolant feed in a very short time. If the recycle stream is large enough, a less significant
temperature profile is then acquired, which has the advantage that the temperature of the coolant can be
74
assumed equal. Stephanopoulos reports that a recycle stream improves the heat transfer in general.
71
Furthermore, according to Shinskey the dynamical system response is improved, after the process
68
controller has intervened. This is due to the faster coolant substitution i.e. less delay. Additionally Roffel
argued that the use of a recycle stream is practically indispensable. This is because if the coolant velocity
becomes too small, heavy fouling can emerge. (In Appendix 5 more about the size of the recycle stream).
The logarithmic temperature mean (expression 17) could be introduced to improve the model accuracy.
Despite this improvement, Tlog will not be considered in this report due to the implicit character of
equation 17 which is not solvable in the bifurcation package LOCBIF.
Mathematical model
The base case mass and heat balance and the coolant temperature balance represented by respectively
equations 2, 3, 12 and 42 describe the proportional controlled base case according:
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
dT
cool C P,coolVcool cool = coolC P,cool V ,cool (Tcool,0 Tcool ) + UA (T Tcool )
dt
VR
(2)
(3)
(12)
(42)
Using the base case values printed in Table 9 the steady state coolant flowrate (set point) can be
calculated with correlations 13 and 14.
Moreover, there is one degree of freedom in the system model which is the inlet coolant temperature
Tcool,0. To acquire knowledge regarding the influence of the inlet coolant temperature, three distinct values
have been chosen: (Table 18). These values have been chosen arbitrary and vary from cold to average
and high in which for the latter no multiplicity is possible.
Several reactor designer stability maps can provide useful information about the stability of this case. The
active LOCBIF parameters are: the coolant flowrate (set point) and the proportional gain Kc.
It is however, possible that the controller calculates negative coolant flowrate values, which in fact means
that heat is generated (perpetual mobile). Therefore, restrictions have been implemented In the LOCBIF
models (V,cool 0). Practically this should mean that the coolant throughput supply is halted temporarily.
40
Table 18 Coolant flowrate.

Tcool,0 [K] V,cool [m3 s-1]
303
400
435
0.00241
0.00811
0.05540
510
500
480
Base case
Kc=1
Kc=0.1
Kc=0.01
Kc
Kc=0.005
Kc=0.001
Kc
III
460
Tcool
III
470
Hopf
Equilibrium
II
490
Fold
Fold
Hopf
Kc=0.0005
III
450
Kc=0.0002
Kc=0.0001
Kc=0.000075
440
Tcool
430
Fold
Kc=0.00005
Kc=0.00001
Hopf
Kc=0.000025
420
Kc=0.000001
410
0.000
Kc=0
0.001
0.002
0.003
0.004
0.005
Coolant flowrate (set point) [m3/s]
Figure 28 Reactor designer stability map. Proportional controlled coolant flowrate. Tcool,0 = 303 [K]. Active LOCBIF bifurcation
parameters: V,cool (set point). Kc > 0.00025 [m3 s-1 K-1] (Hopf maximum) means more stability. Fold maximum Kc >
0.000067 [m3 s-1 K-1].
The equilibrium curves in stability map (Figure 28) rotate counter clockwise if the proportional gain
is increased. The existence of multiplicity in the stability map (Figure 28) is conspicuous.
3 -1
-1
In case Kc > 0.000067 [m s K ] i.e. Fold maximum, the hazard for multiplicity near the base
case situation is reduced. Nevertheless, still extinction is possible although Figure 28 does not
provide that information.
3 -1
-1
If Kc > 0.00025 [m s K ] i.e. Hopf maximum, the equilibrium curves traverse region I instead of
region II which implies that stability is possible. Through the following orbit curves, this will be
demonstrated.
The marking process i.e. region I, II and III in Figure 28 must be interpreted with caution. Because
extinction is always possible in the system is large disturbed. Then, region I indicates actually a
region III. In case a region is indicated as region I this implies that stability is possible, even if the
process it slightly disturbed, however it must be considered with caution.
Orbit curves
Orbit curves have been created for increasing Kc values for both the reactor temperature and the coolant
flowrate. Concerning the first orbit curve (Figure 29) no controller action is applied and a strong limit cycle
system should be expected due to the fast temperature rise (like Figure 22). Nevertheless, due to the
existence of multiplicity, the reactor temperature is in fact attracted towards a higher steady state value
after the first peak and the limit cycles disappear. Figure 32 is the accompanying coolant flow flowrate
curve, which shows that the coolant flowrate remains at the same value, which is obvious because the
controller does not intervene.
41
600
500
400
0
3600
7200
Time [s]
10800
Figure 29 Kc = 0 [m3 s-1 K-1].
T = 0 [K] transition to
higher steady state.
600
500
400
0
3600
7200
Time [s]
derived from equation 13.

T = 0 [K] extinction.
0.0025
0.0000
0
3600
0.0050
0.0025
0.0000
0
Time [s]
400
10800
T = 0 [K]. Instability and

extinction.
3600
0.0075
0.0050
0.0025
0.0000
7200 10800
Time [s]
3600
7200 10800
Time [s]
Figure 32 Kc = 0 [m3 s-1 K-1]. T = 0 Figure 33 Kc = 0.00006 [m3 s-1 K-1]

[K] Coolant flow is not
adjusted.
3600
7200
Time [s]
0.0100
0.0075
7200 10800
500
Figure 31 Kc = 0.00019 [m3 s-1 K-1].
Coolflow [m 3/s]
Coolflow [m 3/s]
0.0050
600
0.0100
0.0075
700
10800
Figure 30 Kc = 0.00006 [m3 s-1 K-1]
0.0100
Coolflow [m 3/s]
Temperature [K]
700
Temperature [K]
Temperature [K]
700
derived from equation 13.

T = 0 [K] extinction.
Figure 34 Kc = 0.00019 [m3 s-1 K-1].
T = 0 [K]. Instability and

extinction.
3
Using equation 13 a first approximation can be made to determine the proportional gain: Kc = 0.00006 [m
-1
-1
s K ]. Figure 30 clearly shows that with this Kc the proportional controller is not robust enough to restrict
the evolving limit cycles. After the first peak, which in fact is a limit cycle, transition towards another steady
state is inevitable causing i.e. extinction. In the stability map Figure 28, for this particular Kc value, the
multiple steady states are visible. The latter can be explained because the flowrate manipulation is not
accurate. The initial temperature rise has to be repressed through increasing the coolant flowrate. This
has the consequence that the reactor temperature drops distinguished. To halt the temperature decline,
the proportional gain of Figure 30 and Figure 33 is not robust enough resulting in extinction i.e. too much
heat has removed from the system. If Kc is increased, Figure 31 shows that still the proportional gain is not
powerful enough to suppress the self-sustained oscillations (Figure 34). After approximately two hours, the
overshoot is considerably enlarged and the temperature is again attracted to a lower steady state i.e.
extinction towards the feed temperature. Again increasing Kc means stability, nonetheless small
temperature limit cycles emerge and considerable for the manipulated coolant flowrate (Figure 35). If for
the same Kc value a disturbance is constrained, due to the large overshoot extinction is inevitable (Figure
36 and Figure 39).
600
500
400
0
3600
7200
Time [s]
10800
Figure 35 Kc = 0.0002 [m3 s-1 K-1].
T = 0 [K]. Limit cycles.
700
700
Temperatur [K]
Tem perature [K]
Temperature [K]
700
600
500
400
0
3600
7200
Tim e [s]
10800
Figure 36 Kc = 0.0002 [m3 s-1 K-1].
T = +10 [K]. Extinction T

T0.
600
500
400
0
3600
7200
Time [s]
10800
Figure 37 Kc = 0.00025 [m3 s-1 K-1].
T = +10 [K]. Spiral point.
42
0.0075
0.0050
0.0025
0.0000
0
3600
0.0100
Coolflow [m 3/s]
0.0100
Coolflow [m 3/s]
Coolflow [m 3/s]
0.0100
0.0075
0.0050
0.0025
0.0000
7200 10800
3600
Time [s]
7200 10800
0.0075
0.0050
0.0025
0.0000
0
Time [s]
Figure 38 Kc = 0.0002 [m3 s-1 K-1]
7200 10800
Time [s]
Figure 39 Kc = 0.0002 [m3 s-1 K-1].
3600
Figure 40 Kc = 0.00025 [m3 s-1 K-1].
T = +10 [K]. Extinction T

T0.
T = +10 [K]. Spiral point.
A proportional gain larger than the Hopf maximum accomplishes a stable system after a step disturbance
(Figure 37). The controlled coolant flowrate varies however extraordinary (Figure 40) and besides has to
be halted for a particular period. Proportional controlling the coolant flowrate can provide a stable process.
However, the severe coolant flowrate variation is disputed. Consider a recycle stream which is applied
with a recycle flow which is 10 times the throughput (Appendix 5), the maximum fluid velocity after t 4
-1
[min] (Figure 40) becomes vcool =1.2 [m s ] which is permitted.
The next stability map (Figure 41) with a less cooler inlet coolant temperature Tcool,0 = 400 [K] is basically
similar to the Figure 28. The main difference is the range on the abscissa, which is obvious because in
case the coolant temperature is higher, the temperature difference is smaller; consequently, more coolant
throughput is required to transfer the surplus heat.
510
Fold
Hopf
500
II
490
480
Tcool
Kc
Kc=1
Kc=0.1
III
470
Hopf
Fold
Equilibrium
Kc=0.01
Base case
Kc=0.005
Kc=0.001
Kc
460
Kc=0.0005
III
450
III
Kc=0.0002
Kc=0.0001
Kc=0.000075
440
Tcool
430
Fold
Kc=0.00005
Kc=0.000025
Hopf
Kc=0.00001
420
Kc=0.000001
410
0.000
Kc=0
0.005
0.010
0.015
0.020
Figure 41 Reactor designer stability map. Proportional control on coolant flowrate. Tcool,0 = 400 [K]. Active parameters: V,cool
(set point) and Kc. Hopf maximum Kc > 0.00078 [m3 s-1 K-1]. Fold maximum Kc > 0.00026 [m3 s-1 K-1].
43
The equilibrium curves in Figure 41 rotate likewise counter clockwise if the proportional gain is
increased. The existence of multiplicity in Figure 28 is evident.
3 -1
-1
In case Kc > 0.00026 [m s K ] i.e. Fold maximum, multiplicity is less dangerous because the
controller cooling is stronger and prevent too large overshoot caused by the limit dynamical
3 -1
-1
instability. If Kc > 0.00078 [m s K ] i.e. Hopf maximum, region II in which limit cycles emerge, is
excluded. The gain is apparently large enough to restrict the limit cycles. Nevertheless, large
external disturbances can always cause the transition towards the lower steady state.
The equilibrium curves intersect regions I and III. Once again, region I implies that the controller
can stabilise the process accurately, although extinction remains thinkable possibility in case of
large perturbations. Region III refers to extinction.
510
500
II
Hopf
Fold
Fold
Hopf
Tcool
490
Kc=1
Equilibrium
Kc=0.1
480
Hopf
470
Base case
Kc=0.005
Hopf
Kc=0.001
Kc
460
Kc=0.0005
II
Fold
450
440
Kc=0.01
Kc
Kc=0.0002
Kc=0.0001
Kc=0.000075
Tcool
Kc=0.00005
Kc=0.000025
430
Kc=0.00001
420
Kc=0.000001
410
0.00
Kc=0
0.02
0.04
0.06
0.08
0.10
Figure 42 Reactor designer stability map. Proportional control on coolant flowrate. Tcool,0 = 435 [K]. Active bifurcation
parameters: V,cool (set point) and Kc. Hopf maximum Kc > 0.012 [m3 s-1 K-1]. No multiplicity is involed.
If the inlet coolant temperature is enlarged to Tcool,0 = 435 [K], multiplicity near the base case value
which has formerly been found in (Figure 19), is not an issue according to Figure 42.
Nevertheless, due to the relatively small temperature difference between the reactor temperature
and the coolant temperature, a considerable coolant throughput is needed to cool the reactor i.e.
about 20 times larger than Figure 28 Tcool,0 = 303 [K]. Despite the advantage that multiplicity is
less prominent, too large fluid velocity in the cooling apparatus can be an impediment.
Differential equation 12 describes the dynamical behaviour of the coolant temperature in which
the assumption has been made that the temperature difference between the coolant inlet and exit
is negligible. This assumption is merely acceptable if a recycle stream is involved in which the
recycle flow is several times the throughput. If due to the minor temperature difference a large
coolant throughput is needed and additionally a substantial recycle stream is required, the coolant
velocity can practically be too large.
3 -1 -1
Finally, LOCBIF noticed the Hopf maximum Kc = 0.012 [m s K ], which means that larger values
involve stability and no hazard of extinction is possible because no multiplicity is involved.
44
Table 19 Proportional gain values P-control coolant flowrate.

Description Tcool,0 = 303 [K]
Tcool,0 = 400 [K]
Tcool,0 = 435 [K]
Kc > 0.00067 [m3 s-1 K-1]
Kc > 0.00025 [m3 s-1 K-1]
Fold maximum
Hopf maximum
Kc > 0.00026 [m3 s-1 K-1]

Kc > 0.00078 [m3 s-1 K-1]
Not existing: no multiplicity

Kc > 0.012 [m3 s-1 K-1]
Summarising
On the one hand, a low inlet coolant temperature has the advantage that superfluous heat can be
withdrawn from the reactor very easy and that coolant fluid velocities are not too large. The disadvantage
is that multiplicity in case of improper proportional control can be inevitable, in particular, if the process is
considerable disturbed. On the other hand, a rather high inlet coolant temperature has the advantage that
multiplicity is less dangerous, but in combination with a recycle stream, too large coolant fluid velocities
can be a hindrance. Therefore, the choice of reactor designer for Tcool,0, Depends greatly on the safety of
the process and in second place on economic aspects e.g. reproduction of heat.
Comparison coolant temperature and coolant flowrate control (T0=303 K)

In previous paragraph, the stability map (Figure 23) for coolant temperature control has been drawn with
active parameters Tcool and Kc. In this chapter, in which coolant flowrate is concerned with, the same
stability map is drawn to investigate the effect of the introduction of the extra differential equation 36. The
implementation implies that multiplicity is probable. The latter is not the case in case the coolant
temperature proportional is controlled.
Equilibrium
510
Fold
500
Fold
III
490
Hopf
Hopf
T=Tcool
Kc=1
Kc=0.1
480
Base case
Kc
470
460
Kc=0.01
Kc=0.005
Kc=0.001
Kc
III
450
Kc=0.0005
Kc=0.0002
III
Kc=0.0001
Kc=0.000075
440
Kc=0.00005
430
Hopf
III
420
410
410
Kc=0.00001
Fold
Kc=0.000025
Kc=0.000001
Kc=0
420
430
440
450
460
470
Coolant temperature (set point) [K]
Figure 43
Reactor designer stability map. Proportional control on coolant flowrate. Tcool,0 = 303 [K]. Active LOCBIF bifurcation
parameters: Tcool (set point) and Kc. Hopf maximum Kc > 0.00025 [m3 s-1 K-1]. Fold maximum Kc > 0.000079 [m3 s-1
K-1]. Compare this stability map with Figure 19. The LOCBIF source LOCBIF rhs 5 is included in Appendix 7.
45
-1
-1
Figure 43 show clearly that for Kc < 0.000079 [m s K ] more steady state solutions are possible.
The main difference between Figure 23 and Figure 43 is the shape or trajectory of the drawn
equilibrium curves. In case of coolant temperature control, increasing the proportional gain causes
the equilibrium curves to move towards the stable region and prevents the occurrence of limit
cycles. In case of coolant flowrate control, the risk for multiplicity is always a point of concern,
despite a dynamically stable behaviour. The stable and unstable regions in Figure 23 are very
clear in contrast to Figure 43 in which, due to multiplicity, these regions are more complicated to
interpret.
According to Figure 23 the (Hopf) maximum Tcool = 469 [K]. Beyond this coolant temperature,
instability is impossible. According Figure 43 this maximum is considerable higher: Tcool = 477 [K].
Again, this can be explained through the fact that, in case of coolant temperature control
multiplicity can cause instability, which is not possible if the coolant temperature is controlled i.e.
only limit cycles.
Reactor temperature versus coolant capacity (T0=303 K)

In previous chapters, it has been found that the cooling capacity, in combination with a suitable
proportional gain, has a stabilising effect on the dynamical behaviour of the concerned base case process.
The following stability maps, in which the coolant inlet temperature is varied, respectively Figure 44,
Figure 45 and Figure 49 will be examined.
510
Equilibrium
500
Fold
II
Fold
Hopf
No pipes
490
Tcool
Kc=1
480
Hopf
470
Kc=0.1
Base case
Kc=0.01
Kc
Kc=0.005
Kc
Kc=0.001
460
Kc=0.0005
III
450
III
Kc=0.0002
Kc=0.0001
Kc=0.000075
440
Tcool
430
Fold
Kc=0.00005
Kc=0.000025
Hopf
Kc=0.00001
420
Kc=0.000001
No pipes
Kc=0
410
10
30
50
70
90
110
Figure 44 Reactor designertability map. Proportional control on coolant flowrate. Tcool,0 = 303 [K]. Active parameters: UA and Kc.
Hopf maximum Kc > 0.00025 [m3 s-1 K-1]. Fold maximum Kc > 0.000079 [m3 s-1 K-1]. UAmax = 78 [kJ s-1 K-1].
46
In Figure 44 the equilibrium curves rotate counter clockwise if the proportional gain is increased.
In this stability map the existence of multiplicity is visible.
3 -1
-1
In case Kc > 0.00025 [m s K ] a dynamically stable system is possible. Nevertheless, extinction
is still probable if the appropriate proportional gain is not used.
Suppose that the reactor designer decides to change the set point reactor temperature, through
following the equilibrium curves the consequences regarding the stability can be determined.
-1
-1
In case UA > 79 [kJ s K ] stability is more certain, although in case of large disturbances,
extinction still is a risk, besides the disadvantage of the economical aspects of the application of
large cooling capacity. A suitable combination of UA and Kc is required.
Table 20 Survey stability P-control coolant flowrate. Tcool,0 = 303 [K].

-1
-1
Cooling capacity [kJ s K ] Dynamical behaviour Explanation
0 UA < 15
15 UA < 65
65 UA < 78
UA > 78
Not existing
Dangerous region
Moderate safe region
Stable region, no limit cycles
UA of the reactor > 0

Multiplicity can disturb the stability
Stability if Kc > Kc (Hopf)
Reactor temperature versus cooling capacity (T0=400 K)
510
II
500
490
Fold
Fold
Hopf
No pipes
Tcool
Equilibrium
Kc=1
480
Base case
Hopf
470
Kc=0.1
Kc=0.01
Kc
Kc=0.005
Kc
460
III
450
Kc=0.001
I
III
Kc=0.0005
Kc=0.0002
Kc=0.0001
Tcool
Kc=0.000075
440
Kc=0.00005
Hopf
Fold
430
Kc=0.000025
Kc=0.00001
420
Kc=0.000001
No pipes
Kc=0
410
10
30
50
70
90
110
Figure 45 Reactor designer stability map. Proportional control on coolant flowrate. Tcool,0 = 400 [K]. Active parameters: UA and
Kc. Hopf maximum Kc > 0.00078 [m3 s-1 K-1]. Fold maximum Kc > 0.00026 [m3 s-1 K-1]. UAmax = 79 [kJ s-1 K-1].
47
Like, Figure 44 the equilibrium curves in Figure 45 rotate counter clockwise if the proportional gain
3 -1 -1
is increased. In case Kc > 0.00078 [m s K ] a dynamically stable system is possible.
3 -1 -1
In Figure 45 in case Kc > 0.00026 [m s K ] (Fold maximum) it seems that the cooling capacity is
-1
-1
limited UAmin 22 [kJ s K ]. For the base case temperature T = 468 [K] this is the minimum
attainable cooling capacity.
Apparently, as the inlet coolant temperature is increased the minimum cooling capacity shifts
upwards.
-1
-1
In case UA > 79 [kJ s K ] stability is more certain, although in case of large disturbances,
extinction still can emerge.
Survey stability P-control coolant flowrate. Tcool,0 = 400 [K].

Table 21
-1
-1
0 UA < 15
15 UA < 22
22 UA < 56
56 UA < 78
UA > 78
Not existing
Not existing
Dangerous region

Minimum Hopf UA for base case
Multiplicity can disturb the stability
Orbit curves
600
500
400
0
3600
7200
Time [s]
10800
700
Temperature [K]
700
Temperature [K]
Temperature [K]
The influence of UA on the stability is demonstrated through of the following orbit curves with equal Kc
-1
-1
value and varying UA values. In Figure 46, UA = 45 [kJ s K ] the proportional gain is not powerful
enough resulting in exhibiting limit cycles. If the cooling capacity is increased to the base case value, the
cooling capacity is considerably large i.e. much more heat can be transferred, however the overshoot
cannot be restrained, which causes the reaction to extinguish. (Figure 47). If the cooling capacity is even
larger, the overshoot can be shortened, in which the transition towards the lower steady state is excluded
(Figure 48).
600
500
400
0
3600
7200
Time [s]
10800
700
600
500
400
0
3600
7200
Time [s]
10800
Figure 46 Tcool,0 = 400 [K], Kc = 0.001 Figure 47 Tcool,0 = 400 [K], Kc = 0.001 Figure 48 Tcool,0 = 400 [K], Kc = 0.001
[m3 s-1 K-1], UA = 45 [kJ s-1
K-1], T = +10 [K]. Limit
cycles.
[m3 s-1 K-1], UA = 55 [kJ s-1

K-1], T = +10 [K]. Extinction
T T0.
[m3 s-1 K-1], UA = 95 [kJ s-1

K-1], T = +10 [K]. Stable
system.
48
Reactor temperature versus cooling capacity (T0=435 K)

Equilibrium
510
Fold
500
II
490
480
No pipes
Tcool
Kc=1
Base case
Kc
470
Hopf
Fold
Kc=0.1
Kc=0.01
Kc
Hopf
Kc=0.005
450
Fold
440
Kc=0.001
460
II
Kc=0.0005
Kc=0.0002
Kc=0.0001
Hopf
Kc=0.000075
Tcool
Kc=0.00005
430
Kc=0.000025
Kc=0.00001
420
Kc=0.000001
No pipes
Kc=0
410
10
30
50
70
90
110
Figure 49 Reactor designer stability map. Proportional control on coolant flowrate. Tcool,0 = 435 [K]. Active parameters: UA and
Kc. Hopf maximum Kc > 0.012 [m3 s-1 K-1]. UAmax = 83 [kJ s-1 K-1].
-1
-1
In Figure 49 for the base case temperature the minimum cooling capacity is UAmin 46 [kJ s K ].
This means that for lower UA values, the reactor temperature increases distinctive, no matter the
proportional gain and strong limit cycles occur.
In Figure 49 it is evident that multiplicity near the base case value is not anymore concerned.
Increasing the proportional gain again rotates counter clockwise the equilibrium curves. In case Kc
3 -1 -1
> 0.012 [m s K ] the system behaves dynamically stable.
-1 -1
In case UA > 83 [kJ s K ] stability is guaranteed.
Table 22 Survey stability P-control coolant flowrate. Tcool,0 = 435 [K].

-1
-1
0 UA < 15
15 UA < 46
46 UA < 82
UA > 82
Not existing
Not existing

Minimum Hopf UA for base case
Table 20, Table 21 and Table 22 provide a survey of the dynamical behaviour for several cooling inlet
temperatures. A system that appears to be stable can become unstable after a considerable perturbation.
Therefore, these tables have to be examined with caution.
49
9.1.3
Controlling the throughput
In the process industry where one aims to operate at maximum yield, where primarily the process
conditions based on process design are in general fixed, it is not likely that tremendous changes in
flowrate are permitted. In case slight changes (i.e. in between a range of a few per cent) are allowable, it
can be interesting to study if the residence time or throughput control can contribute to the possibilities of
eliminating limit cycles.
Mathematical model
Equations 2 and 3 are similar regarding coolant temperature control, except the controller correlation 43:
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
VR
(2)
(3)
(43)
Both the possibilities will be examined i.e.: Kc > 0 and Kc < 0.

Comparison coolant temperature and flowrate control
The main difference between coolant temperature control and flowrate control is explained by the fact that
changing the coolant temperature affects directly the heat balance, which has, due to the fact that the
balances are coupled, indirectly consequences for the conversion. Flowrate alteration affects both directly
the heat and mass balance. It changes the existing situation in the CISTR whereas the coolant
temperature can be seen more as an external parameter. Cooling a reactor by decreasing the coolant
temperature mainly implies a larger driving force to heat transfer. Controlling a reactor for instance by
increasing the flowrate will on the one hand decrease the reactor temperature due to the supply of cold
feed flow, but on the other hand it will increase the reactor temperature due to the increasing heat
production by the chemical reaction. Thereby, flowrate control is more complicated due to the existence of
multiplicity, which is not the case for coolant temperature control regarding the base case. This role of
multiplicity will be demonstrated in a later stadium in this chapter. Although the mathematical description
does not differ, the impact on the reactor temperature and conversion is considerably different.
Positive Kc value
The controller equation becomes:
V = V ,sp + K c (Tsp T )
(43a)
A reactor designer stability map is created (Figure 50) applying positive proportional gain values (Kc > 0).
The active bifurcation parameters are: the flowrate (set point) and the proportional gain.
50
Fold
510
III
500
Fold
Tcool
490
Hopf
Kc=0.00001
Kc=0.00005
Base case
Kc
470
Kc=0
Equilibrium
Kc
II
480
Hopf
Kc=0.0001
Kc=0.00012
460
450
440
Kc
III
Max Hopf
Kc=0.000125
Kc=0.0002
III
Kc=0.0003
Tcool
Kc=0.0004
430
Kc=0.0005
Kc=0.001
420
410
0.000
Kc=0.01
0.002
0.004
0.006
Flowrate (set point) [m3/s]
Figure 50
0.008
0.010
Kc=0.1
Kc=oo
Reactor designer stability map proportional controlled flowrate (Kc > 0). The considered base case model exhibiting
limit cycles cannot be controlled satisfactorily due to unstable regions in the vicinity of the base case. The region II
in fact is a III region because too large Kc values can cause the reaction to extinct, however for low Kc values limit
cycles exhibit. Hopf maximum Kc = 0.000125 [m3 s-1 K-1].
In case the proportional gain is increased, the equilibrium curves rotate a little clockwise.
In view of the controlled flowrate, multiplicity is of significance over the entire range (0 V < 0.01
3 -1
[m s ]), which is clearly visible in the stability map because of the equilibrium curves.
The equilibrium curves alter due to the increased proportional gain nevertheless cannot prevent
the process to behave unstable. No region I can be observed.
51
Orbit curves
700
600
500
400
0
3600
7200
Tim e [s]
10800
Figure 51 Kc = 0.00007 [m3 s-1 K-1].

The proportional controller
cant eradicate the evolving
limit cycles.
Temperature [K]
0.010
Flow rate [m 3/s]
Tem perature [K]
Considered positive Kc values the following orbit curves have been drawn to demonstrate that for Kc > 0
throughput control is not applicable.
0.008
0.006
0.004
0.002
0.000
0
3600
7200
Time [s]
600
500
400
0
10800
Figure 52 Kc = 0.00007 [m3 s-1 K-1].
700
3600
7200
Time [s]
10800
Figure 53 In case the controller gain is
The flowrate varies

exceedingly which is
presumable unwanted.
slightly increased towards Kc

= 0.00008 [m3 s-1 K-1] the
reaction extinguishes.
Using the steady state solution of the mass and heat balance 2 and 3, a proportional gain can be
3 -1
-1
estimated: Kc = 0.00007 [m s K ]. Apparently, this value is too small to elucidate the limit cycles, which
is confirmed through Figure 51. A larger Kc however cannot achieve the desired steady state; actually, a
steady state with a much lower conversion is reached (Tcool) (Figure 53). Figure 52 displays the eminent
variation of the throughput, which is in industrial process unacceptable.
The question arises: Why cant the process be controlled? In case a disturbance occurs, due to external
factors or due to the self-sustained oscillations, the proportional controller determines the corrected
flowrate, which should adapt the process towards the set point. On the one hand, a large Kc value is
required to repress the evolving limit cycles causing overshoot. On the other hand, small Kc values are
inevitable to regulate the process mildly and not to disrupt the actual process. Too rigorous flowrate
variations, however causes the process to jump (attract) to other steady state situations. In case too less
heat is withdrawn runaway occurs and if the throughput is too large the reaction extinguishes called the
26
blow-out velocity (Fogler ).
Based on stability map (Figure 50) and displayed orbit curves one can conclude that it is dissuadable to
control the base case applying a proportional controller with positive values for Kc.
Negative Kc value
Using the steady state solution of the mass and heat balance 2 and 3, a proportional gain can be
3 -1 -1
estimated: Kc = -0.0007 [m s K ]. Therefore, due to the negative sign equation 43 changes into:
(43b)
A new stability map (like Figure 50 for Kc > 0) is created for Kc < 0 (Figure 67).
52
510
III
500
Hopf
Kc=0
Hopf
490
Fold
Fold
Kc=0.00001
Kc=0.00005
480
II
470
Base case
Kc
III
Kc
Equilibrium
Hopf
440
Kc=0.00012
Kc=0.000125
460
450
Kc=0.0001
Kc=0.0002
Kc=0.0003
Fold
Kc=0.0004
Kc=0.0005
Tcool
Kc=0.001
430
Kc=0.01
420
Kc=0.1
410
0.000
Tcool
0.002
0.004
0.006
0.008
0.010
Flowrate (set point) [m3/s]
Figure 54 Reactor designer stability map. Proportional controlled flowrate (Kc < 0). Hopf maximum Kc = 0.000125 [m3 s-1 K-1].
The region marked II is actually a region III because the limit cycles can cause extinction or runaway but is marked
with II because the base case exhibits limit cycles in case no controller is interfering.
The Hopf and Fold curves are similar like Figure 50.
The equilibrium curves rotate counter clockwise in case the proportional gain is increased.
A particular small region marked I is found which implies that the process may behave stable. The
consequences are that the desired set point T = 468 [K] cannot be maintained, through following
the equilibrium curves. If however, it is not allowed to change this temperature, one has to accept
a lower conversion. = 0.47 instead of = 0.69.
Orbit curves
3
-1
-1
For region I in Figure 54 the dynamical behaviour is studied. A proportional gain Kc = 0.000124 [m s K ]
3 -1
-1
slightly smaller than the Hopf maximum value Kc = 0.000125 [m s K ] clearly shows limit cycles (Figure
55). The maximum Hopf Kc value (Figure 56) shows that the process behaves stable after a step
disturbance of T = 10 [K] through a spiral point, nevertheless with higher temperature and conversion.
Conversely, Kc is not large enough to maintain the set point steady state. Increasing Kc results in the first
place runaway (Figure 57) and subsequently realises stability (Figure 58). Very large Kc values with a
positive temperature disturbance preserve the stability and the desired set point (Figure 59). However, in
case the disturbance is negative, it is apparent that Kc causes the temperature to drop to the coolant
temperature (Figure 60).
53
600
500
400
0
3600
7200
700
Temperature [K]
Temperature [K]
Tem perature [K]
700
600
500
400
10800
3600
Tim e [s]
400
0
3600
500
400
10800
Time [s]
Figure 58 Kc = 0.00140 [m3 s-1 K-1].

Step disturbance T = 10 [K].
The process behaves stable.
However with a lower
conversion = 0.47 instead
of = 0.69.
700
600
500
400
0
3600
7200
10800
Figure 57 Kc = 0.00139 [m3 s-1 K-1]

and T = 10 [K]. Slightly
higher Kc values cause
runaway. T 774 [K].
Temperature [K]
600
7200
Time [K]
and T = 10 [K]. The system

is stable (spiral point) but
jumps to a higher steady state
i.e. T 524 [K]. = 0.92.
Temperature [K]
Temperature [K]
700
7200
500
10800
Figure 56 Kc = 0.000125 [m3 s-1 K-1]
and T = 10 [K]. A limit

cycle evolves around a
different steady state.
3600
600
Time [s]
Figure 55 Kc = 0.000124 [m3 s-1 K-1]
7200
700
700
600
500
400
10800
Time [s]
Figure 59 Kc = 0.01 [m3 s-1 K-1].
3600
7200
10800
Time [s]
Figure 60 Kc = 0.01 [m3 s-1 K-1].
If T = 10 [K] stable system,

nevertheless the set point T
can be maintained, not the
conversion.
Step disturbance T = -10

[K]. The process behaves
unstable. The temperature
drops to the coolant
temperature due to the too
rigorous control.
Figure 60 demonstrated another complication. It is remarkable that a disturbance T = +10 [K] (Figure 59)
can be controlled, however it makes the system to extinguish in case T = -10 [K] as is illustrated in
Figure 60. The sign of the offset is evidently decisive. In case the offset > 0 the reaction appears to
extinguish (T Tcool) i.e. the controllers manipulates with the wrong sign, the same phenomenon as in
case Kc > 0. If the sign is reverse, the controller will decrease the flowrate causing temperature drop
eventually until the coolant temperature is reached (extinction). Because V 0 the controller cannot
anymore manipulate (multiply) the flowrate and the reactor is actually halted.
Modification:
Downwards approaching temperature disturbance i.e. offset < 0: sign > 0.
Upwards approaching temperature disturbance i.e. offset > 0: sign < 0.
Disadvantage of this trick is the virtually endless converging process around the set point. Nonetheless,
an adaptation in the LOCBIF source achieves Figure 60 to become like Figure 59.
The problems caused by the sign of the offset can be solved using a gimmick but is not a very
scientifically solution. A controller with knowledge of the process in which the Kc value is fitted to the
particular situations e.g. (Kc < 0 / Kc > 0) or (Kc robust / Kc gentle) can help to improve the controlling
procedure. In this case, non-linear process control could be an option.
54
Conversion versus reactor temperature

If not the reactor temperature but the base case conversion has to be remained, the proportional controller
equation becomes:
V = V ,sp K c ( sp )
(45)
The controller changes the flowrate to acquire a base case conversion = 0.68 which has the corollary
that the reactor temperature decreases and finally another steady state is reached.
Orbit curves
Orbit curves will be created applying controller equation 45. Due to proportional control, steady state offset
however is inevitable.
Table 23 Stability maps P-control throughput.
Proportional Final
Final reactor Final
gain
conversion temperature flowrate
3 -1
3 -1
Kc [m s ]
T [K]
[-]
V [m s ]
5 10-6
6 10-6
10 10-6
average 0.84
0.72
0.91
Taverage 464
454
444
0.0027
0.0021
0.0003
Dynamical
behaviour
Reference
Limit cycles
Spiral point
Asymptotic damping
Figure 61, Figure 62

0.8
0.6
0.4
0.2
0.0
0
3600
7200
1.0
1.0
0.8
0.8
Conversion
1.0
Conversion [-]
Conversion [-]
Figure 61 demonstrates that the proportional controller is not powerful enough to eradicate emerging limit
cycles. Both conversion and temperature (Figure 62) still exhibits sustained oscillations. Increasing Kc
implies stability (Figure 63 and Figure 64) nevertheless with an offset. The conversion finally becomes =
0.72 [-] with an accompanying temperature T = 454 [K].
0.6
0.4
0.2
0.0
10800
Figure 62
3600
7200
10800
600
500
400
10800
600
500
400
Tim e [s]
Figure 64 Kc = 5 10-6 [m3 s-1].
Kc = 6 10-6 [m3 s-1].
700
3600
7200
10800
Time [K]
Figure 65
3600
7200
10800
Time [s]
Figure 63 Kc = 1 10-5 [m3 s-1].

Temperature [K]
Temperature [K]
Tem perature [K]
700
7200
0.2
Time [s]
Figure 61 Kc = 5 10-6 [m3 s-1].
3600
0.4
0.0
Time [s]
0.6
Kc = 6 10-6 [m3 s-1].
700
600
500
400
0
3600
7200
10800
Time [s]
Figure 66 Kc = 1 10-5 [m3 s-1].
If throughput control is applied, the reactor designer and controller have to decide which process variable
(temperature or conversion) is decisive.
55
Reactor temperature versus cooling capacity

A significant large cooling capacity can have a positive effect on the stability of the process. The study to
the not controlled base case resulted in a stable region, in which no limit cycles exhibit, in case UA > 78
-1
-1
[kJ s K ] viewed in Figure 24. Hence the effect of the cooling capacity is examined for the base case +
proportional controller on the flowrate. A stability map is created regarding active parameters, which are
the proportional gain and the cooling capacity resulting in Figure 67.
Fold
510
Hopf
500
Hopf
490
Fold
Kc=0.00838
(Hopf max)
Kc=0
Kc=0.00001
Hopf
480
Kc=0.0001
Kc=0.0002
Equilibrium
470
Kc=0.0003
Kc
Kc=0.0004
Base case
460
Kc=0.0005
Hopf
Kc=0.000838
Fold
Kc=0.001
Kc
450
Kc=0.0015
440
Tcool
Kc=0.00175
Kc=0.002
430
Kc=0.003
Hopf
420
Kc=0.1
Tcool
410
10
30
50
70
90
110
Figure 67 Reactor designer stability map. Proportional controlled flowrate. Hopf maximum UA > 84 [kJ s-1 K-1]. Fold maximum Kc
= 1.05 10-4 and Fold minimum Kc = 2.2 10-4 [m3 s-1 K-1].
The 3D effect is here clearly visible, therefore,

In case the proportional gain is increased, the equilibrium curves at first rotate clockwise until Kc =
3 -1 -1
0.00838 [m s K ] then the equilibrium curves are indented and resemble a parabolic shape. It is
apparent to perceive in Figure 67 that more than one steady state solution is possible and
definitely with respect to the base case.
Multiplicity is of significance over the entire examined range. Therefore, the region marked with I
must be considered with the knowledge that runaway or extinction is probable in case an external
disturbance is eminent. In case the simulation is started from the base case steady state situation,
the system behaves stable.
In the contrary, to the stability maps formerly created in, (e.g. Tcool control) in which the stable
region could be marked clearly, is in case of flowrate control due to the prominent existing of
multiplicity much more complicated. A specific chosen Kc could mean that a point is stable
throughout the perspective of the base case, however can make a system to become unstable in
case a external disturbance is imposed.
-1
-1
LOCBIF found the Hopf maximum for UA = 73 [kJ s K ], which is almost the same value UA =
-1 -1
75 [kJ s K ] as for coolant temperature control (Figure 24) which has formerly been found.
56
Orbit curves
-1
-1
600
500
400
0
3600
7200
10800
700
Temperature [K]
700
Tem perature [K]
Temperature [K]
If the cooling capacity UA = 73 [kJ s K ] derived from the Hopf maximum value. Limit cycles appear in
case no controller is correcting (Figure 68). In case Kc is increased the limit cycles are visibly shrinking
(Figure 69). For slightly larger Kc values limit cycles, vanish (Figure 70).
600
500
400
0
Time [s]
-1
7200
10800
-1
Figure 69 UA = 73 [kJ s-1 K-1]
Kc = 0 [m s K ].
= 0.50, T = 452 [K].
-1
-1
Kc = 3 10 [m s K ].
T = 10 [K]. Fading limit
cycles. = 0.52, T = 452 [K].
-5
-1
600
500
400
0
3600
7200
10800
Time [s]
Tim e [s]
Figure 68 UA = 73 [kJ s-1 K-1]

3
3600
700
-1
Figure 70 UA = 73 [kJ s-1 K-1].
3 -1
-1
Kc = 3.1 10-5 [m s K ].
T = 10 [K]. Spiral point.
= 0.53, T = 451 [K].
-1
-1
With the base case cooling capacity UA = 55 [kJ s K ], stability can be acquired if Kc > 0.002 [m s K ]
although with lower conversion. The choice of the proportional gain must be done with caution in
consideration of instability.
Summarising
Based on the stability map and orbit curves it can be noted that proportional control is possible.
Nevertheless, suitable choice of the proportional gain is necessary to avoid unstable situations. The
acquired Kc values through the bifurcation software program LOCBIF can merely indicate that Kc must be
strong enough to exclude the self-sustained limit cycles but cannot be too large to prevent instability due
to the existence of multiplicity. The choice has to be made to, or to preserve the reactor temperature or the
conversion. In case the proportional controller is applied which regulated the throughput, both cannot
simultaneously be maintained. An advanced controller with knowledge of the process device is therefore
recommended. With regards to the cooling capacity, larger UA values have a stabilising effect. The
LOCBIF source LOCBIF rhs 7 is included in Appendix 7.
9.2 Reactor controller stability map

The following stability maps will be examined:
Table 24 Reactor controller stability maps P-control.
Stability map
Reference
Kc versus UA Controlling Tcool
Figure 71
Kc versus UA Controlling V,cool Figure 72
Figure 73
Kc versus UA Controlling V
57
9.2.1
Mathematical model
The mathematical model with regards to the coolant temperature P-controlled base case states:
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
VR
(2)
(3)
(41)
Proportional gain [-]
A very distinct stability map can be created using active LOCBIF bifurcation parameters UA and Kc
resulting in Figure 71.
10
9
8
7
6
5
4
3
2
1
0
Hopf
Base case
II
10
30
50
70
I
90
110
Figure 71 The reactor controller stability map of the proportional controlled base case model. To obtain a stable process, a
larger proportional gain has to be chosen in case UA decreases. UAmax (Hopf) = 78 [kJ s-1 K-1] and Kc, max = 5.7 [-].
Discussion of the reactor controller stability map:
The same Hopf maximum can be found values as for the reactor designer stability map Figure 24.
The stability map shows that Kc decreases when the cool capacity increases.
The dotted line represents the UA pertaining to the reactor with no cooling pipes. Through Figure
71, the maximum required proportional controller gain could accordingly be derived. If Kc = 5.7,
the base case will never exhibit limit cycles no matter the extent of cooling capacity, provided that
the coolant temperature can be controlled instantaneously, which in practise is not completely the
case. Nevertheless, for this particular dealing scenario it is wise to choose a Kc as large as
possible.
By means of increasing the cooling capacity, less proportional gain is required to prohibit limit
cycles.
The reactor controller stability map is a very useful tool to determine an appropriate Kc value in
case for instance the cooling capacity is changed. Therefore, a particular UA value is chosen, e.g.
-1 -1
UA = 40 [kJ s K ], to acquire a stable process, Kc = 1.8 to shift from region II to stable region I.
58
9.2.2
Mathematical model
The mathematical model with regards to the coolant flowrate P-controlled base case states:
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
dT
dt
VR
(2)
(3)
(12)
(42)
For coolant flowrate control, the following reactor controller stability map is created, similar as for coolant
temperature control (Figure 71).
Proportional gain [m3/s/K]
Tcool,0 = 303 [K]
0.010
0.009
0.008
0.007
0.006
0.005
0.004
0.003
0.002
0.001
0.000
III
Tcool,0 = 435 [K]
I
Hopf
maximum
Base
case
10
No pipes
Figure 72
Tcool,0 = 400 [K]
30
50
70
90
110
The reactor controller stability map. Hopf curves using active bifurcation parameters UA and Kc for various coolant
inlet temperatures. For each Hopf curve, the process is stable at the right-hand side and unstable at the left-hand
side.
Although In the stability map Figure 72 different axes are used, compared to Figure 71, The
stability map shows again that Kc decreases when the cool capacity increases.
The stable region I is clearly appearing, besides for every Tcool,0 the minimum and the maximum
cooling capacity.
The introduction of the coolant differential equation, which implies multiplicity, is clearly visible. For
Tcool,0 = 435 [K], one can see that multiplicity is not anymore existing with regards to the base
case.
The stability map demonstrates that the maximum UA Hopf value is dependent of the inlet coolant
temperature.
59
9.2.3
Mathematical model
The mathematical model with regards to the proportional controlled throughput base case states:
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
VR
(2)
(3)
(43b)
The subsequent reactor controller stability map is drawn to investigate the influence the relationship
between the cooling capacity and the proportional gain with respect to the stability.
0.0025
0.0020
0.0015
Hopf
0.0010
0.0005
III
0.0000
10
30
50
70
90
110
Figure 73 Reactor controller stability map. Proportional controlled throughput. Hopf maximum UA > 73 [kJ s-1 K-1] and Kc =
0.00084 [m3 s-1 K-1]. In case the cooling capacity UA > 78 [kJ s-1 K-1] or Kc > 0.002 [m3 s-1 K-1] stability is certain,
nevertheless with lower conversion.
Analogous to previous stability maps (Figure 71 and Figure 72), Figure 73 the clearly
demonstrates the problem with multiplicity.
The same Hopf maximum values have been acquired as for reactor designer stability map Figure
67 in which the latter is rather complicated to interpret.
A clear stable region I is however located on the right hand side of the Hopf curve, however with
lower conversion.
Summarising
Proportional controlling the throughput is possible. Disadvantage is that or the reactor temperature or the
conversion must be chosen as a set point. Both contemporaneous is not possible because due to the
existence of multiplicity the process is forced to one steady state situation. Due to the dynamical unstable
situation of the base case, a certain proportional gain has to be chosen to resist the limit cycles. And here
the problem arises, that too large Kc values cause inevitable transition towards other steady states and in
particular if the process is also disturbed. Therefore, if only one reactant is involved, like the base case,
controlling the throughput is not a suitable controlling method if both the reactor temperature as well the
conversion must be maintained.
60
Chapter X
10 PROPORTIONAL-INTEGRAL CONTROL
In this section, it is considered to what extent the proportional + integral controller affects the stability of
the cooled CISTR with exothermic reaction.
10.1
Reactor designer stability maps
The integral action merely eliminates the offset in which eventually a particular process variable coincides
with the controller set point. Therefore, the equilibrium, Fold and Hopf curves remain unchanged.
Therefore, the reactor designer stability maps of proportional control and proportional-integral control are
similar.
10.2
Reactor controller stability maps
The following reactor controller stability maps will be examined:

Table 25 Stability maps PI-control coolant temperature.
Figure 74
KC versus I
KC versus UA Figure 79
Table 26 Stability maps PI-control coolant flowrate.
Stability map Tcool,0 = 303 [K] Tcool,0 = 400 [K]
Figure 80a
Figure 80b
Kc versus I
KC versus UA Figure 82
Figure 83
Figure 114
Kc versus cool Figure 113
Tcool,0 = 435 [K]

Figure 80c
Figure 84
Figure 115
Table 27 Stability maps PI-control throughput.

Figure 85
Kc versus I
Kc versus UA Figure 86
10.2.1
The mathematical notation by means of equations 8, 9 and 46 is:

Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
K 1
Tcool = Tcool,sp + K c (Tsp T ) + c (Tsp T )dt
I 0
VR
(2)
(3)
(46)
61
The first interesting issue is to find out in what way the integral time I is related with the proportional gain
Kc. Therefore a stability map is created with the proportional gain at the ordinate and the integral time at
the abscissa (Figure 74).
5
4
3
2
1
Hopf
II
0
0
600
1200
Intergral time constant [s]
Figure 74 Reactor controller stability map with Hopf curve in

which the active parameters are Kc and I. Kc(I 0)
= 1.35 and Kc(I ) = 0.9.
Kc is apparently not very dependent on the magnitude of I. In case I 0, Kc = 1.35 and when I
, Kc = 0.9 which is in fact the Kc for the solely proportional controlled base case.
In Figure 74 the Hopf curve is plotted indicating the transition between dynamically stable and
unstable.
Figure 74 demonstrates that increasing the integral time, in fact slightly improves the stability.
Orbit curves
600
500
400
0
Figure 75a
3600
7200
Time [s]
Kc = 1 I = 60 [s].
10800
700
Temperature [K]
700
Temperature [K]
Temperature [K]
The following step is to examine the dynamic behaviour for one selected proportional gain Kc = 1 and
varying integral time. Figure 75a with an integral time I = 60 [s] displays that limit cycles emerge. In case
a proportional controlled has exclusively been applied Kc = 0.9 was sufficient to eliminate limit cycles.
Evidently, small integral times imply larger Kc values to provide stability. Figure 75a-e show that increasing
the integral time is causing the originating of the limit cycles is to be postponed. If I = 600 [s] no limit
cycles will arise according to Figure 75f.
600
500
400
0
3600
7200
Time [s]
Figure 75b Kc = 1 I = 120 [s].
10800
700
600
500
400
0
Figure 75c
3600
7200
Time [s]
10800
Kc = 1 I = 180 [s].
62
600
500
400
0
Figure 75d
3600
7200
Time [s]
700
600
500
400
10800
Kc = 1 I = 240 [s].
Temperature [K]
Temperature [K]
Temperature [K]
700
Figure 75e
3600
7200
Time [s]
700
600
500
400
10800
Kc = 1 I = 300 [s].
3600
7200
Time [s]
10800
Kc = 1 I = 600 [s].
Figure 75f
In previous sections it has been found that Kc = 0.9 could narrowly preserve stability. Figure 76a-c
demonstrate that I have to be increased to acquire a stable system.
K c (PI) .
I
600
500
400
0
3600
7200
Time [s]
700
600
500
400
10800
Figure 76a Kc = 0.9 I = 60 [s].
Temperature [K]
700
Temperature [K]
Temperature [K]
In fact: K c (P) =
3600
7200
Time [s]
700
600
500
400
10800
Figure 76b Kc = 0.9 I = 600 [s].
3600
7200
Time [s]
10800
Kc = 0.9 I = 3600 [s].
Figure 76c
600
500
400
0
Figure 77a
3600
7200
Time [s]
10800
Kc = 1.35, I = 60 [s]
derived from Figure 74.
Constrained temperature
disturbance T = 20 [K].
700
Temperature [K]
700
Temperature [K]
Temperature [K]
For the next case the proportional gain is selected which provides stability for every chosen integral time
Kc = 1.35 (value derived from Figure 74). Figure 77a demonstrates that the limit cycles vanish after a
temperature disturbance T = 20 [K], although despite the long lasting spiral point. Decreasing the integral
time will not change the stability in general (Figure 77b), primarily the oscillation time will be decreased. A
large Kc accomplishes the desired steady state value within a few minutes and not too large overshoot
(Figure 78).
600
500
400
0
3600
7200
Time [s]
10800
Figure 77b Kc = 1.35, I = 5 [s]
derived from Figure 74.

500
490
480
470
460
450
0
300
Time [s]
600
Figure 78 Kc = 4 I = 60 [s].
Ratio a/b = 4/1 discussed
in 5.6.2 (Figure 17).
The LOCBIF source LOCBIF rhs 4 is included in Appendix 7.
63
Proportional gain versus cooling capacity

In previous section in which exclusively the proportional gain has been considered, the cooling capacity, in
combination with a suitable proportional gain, had a stabilising effect on the process stability (Figure 71).
Regarding the PI-controller, the cooling capacity still has this positive consequence alomst irrespective of
the integral action. Figure 79 represents the stability map with applied active parameters UA and Kc in
relation with several values of the integral time.
10
9
no pipes
8
7
tauI=0.1
6
tauI=5
5
4
tauI=60
II
1
No pipes
tauI=600
Hopf
0
10
30
50
70
90
110
Figure 79 Reactor controller stability map of the proportional-integral controlled base case model. If the cooling capacity
decreases, a larger Kc value has to be chosen. The integral time constant has a minor influence. If I the stability
map matches with the proportional controlled base case showed in Figure 71.
The Hop curves are similar to the curves for proportional control only (Figure 71). In case the
cooling capacity decreases, a larger proportional gain must be selected to preserve stability.
The integral time constant has little influence on the stability. Kc is mainly decisive.
The reactor controller can easily determine or a suitable combination of UA and Kc in which
stability is guaranteed.
Summarising
In case the base case is proportional-integral controlled through regulating the coolant temperature, the
integral action has the advantage that the offset is eliminated, however short integral times has to be
avoided because then the process behaves less stable. The controller stability maps Kc-UA can be used
to effortless determine a stable process.
64
10.2.2
Equations 2, 3, 36 and 47 mathematically describe the proportional + integral controlled base case.
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
dT
dt
K 1
V ,cool = V ,cool,sp K c (Tsp T ) c (Tsp T )dt
I 0
VR
(2)
(3)
(12)
(47)
0.0020
0.025
Hopf
Hopf
In Table 19 the proportional gain values regarding the P-control have been printed. These values are
again obtained if (I ) is taken in the Hopf curves of Figure 80a,b and Figure 80c.
0.0015
Tcool,0 = 400 [K]
0.0010
0.0005
Tcool,0 = 303 [K]
0.0000
0.020
0.015
0.010
Tcool,0 = 435 [K]
III
0.005
0.000
300
600
900
Figure 80a,b Reactor controller stability map with Hopf curves

for Tcool,0 = 303 [K] and Tcool,0 = 400 [K].
Proportional gain versus integral time. For each
Hopf curve the upper area I refers to stability and
the lowest region II to instability.
300
600
900
Figure 80c
Reactor controller stability map with Hopf curves for

Tcool,0 = 435 [K]. Proportional gain versus integral
time.
Discussion of the reactor controller stability maps:
For each Hopf curve the upper region points towards a stable region I and the lower region means
instability. Moreover, due to the existence of multiplicity, transition to another steady state remains
present, in particular if the process is considerably disturbed. Therefore, the stable regions always
have to be interpreted with caution.
The stability maps show that for short integral times, for instance if the process response has to
be fast, the proportional gain must be selected considerably larger to maintain stability.
If I > 300 [s] the problem of integral windup almost not anymore exists.
If I the same values are obtained:
3 -1 -1
For Tcool,0 = 303 [K]: Kc > 0.00025 [m s K ] retrieved from Figure 28,
3 -1 -1
For Tcool,0 = 400 [K]: Kc > 0.00078 [m s K ] retrieved from Figure 41,
3 -1 -1
For Tcool,0 = 435 [K]: Kc > 0.012 [m s K ] retrieved from Figure 42.
65
Orbit curves
600
500
400
0
3600
7200
Time [s]
600
500
400
10800
Tcool,0 = 303 [K], Kc =

0.00025 [m3 s-1 K-1],
I = 60 [s] T = +10 [K].
Limit cycles.
Figure 81a
700
Temperature [K]
700
Temperature [K]
Temperature [K]
The following orbit curves demonstrate the influence of the integral time on the stability for a particular
proportional gain (Hopf maximum value).
3600
7200
Time [s]
700
600
500
400
10800
Tcool,0 = 303 [K], Kc =

Figure
0.00025 [m3 s-1 K-1],
I = 600 [s] T = +10 [K].
Spiral point.
Figure 81b
81c
3600
7200
Time [s]
10800
Tcool,0 = 303 [K], Kc =

0.00025 [m3 s-1 K-1],
I = 3600 [s] T = +10
[K]. Stable spiral point.
That a large integral time improves the stability can be confirmed through the orbit curves Figure 81a-c.
Small I make the system unstable in which substantial limit cycles emerge. A larger integral time makes
66
the process more stable (Ratto et.al. ) (Figure 81b) and a very large I (Figure 81c) is similar to Figure 36
in which merely proportional control is concerned. The LOCBIF source LOCBIF rhs 6 is included in
Appendix 7.
In previous chapters, it has been found that the cooling capacity, in combination with a suitable
proportional gain, has a stabilising effect on the dynamical behaviour of the concerned base case process.
In the following stability maps, the cooling capacity has been drawn at the abscissa and the proportional
gain at the ordinate for various integral times.
0.0020
tauI=60
tauI=120
0.0015
tauI=300
0.0010
tauI=600
tauI=3600
0.0005
no pipes
Hopf
Kc(min)
0.0000
10
No pipes
30
50
70
90
110
Figure 82 Reactor controller stability map. Hopf curves with bifurcation parameters: Proportional gain and the cooling capacity
for Tcool,0 = 303 [K]. For each Hopf curve the stable region I is located at the right hand side and the unstable at the
left hand side.
66
Figure 82 (Tcool,0 = 303 [K]) proves that the cooling capacity has a stabilising effect on the
dynamical behaviour of the process. In case a larger integral time is chosen, the stability is even
increased. This is also the case for Figure 83 (Tcool,0 = 400 [K]) and Figure 84 (Tcool,0 = 435 [K]).
Typical for the considered stability maps is the fact that beneath a specific cooling capacity, the
proportional gain increases severely. This can be explained through the fact that the controller
ought to increase the coolant flowrate in case less heat can be transferred through the wall.
-1 -1
-1 -1
In Figure 82 (Tcool,0 = 303 [K]) Kc if UA < 40 [kJ s K ]. UAmin UAvessel = 15 [kJ s K ].
-1 -1
-1 -1
In Figure 83 (Tcool,0 = 400 [K]) Kc if UA < 50 [kJ s K ]. UAmin 30 [kJ s K ].
-1 -1
-1 -1
In Figure 84 (Tcool,0 = 435 [K]) Kc if UA < 60 [kJ s K ]. UAmin 50 [kJ s K ].
The stability maps as function of the inlet coolant temperature confirm that the stable region
increases for a low inlet coolant temperature. This is due to the larger temperature driving force in
the coolant differential equation. In situations in which the cooling capacity could decline, the
lowest inlet coolant temperature provides the best stability. Although still large perturbations can
cause extinction.
0.0020
Hopf
Ho
tauI=60
0.0015
tauI=120
tauI=300
0.0010
tauI=600
0.0005
tauI=3600
No pipes
no pipes
0.0000
10
30
50
70
90
110
left hand side.
67
0.025
Hopf
tauI=60
0.020
tauI=120
0.015
tauI=300
0.010
tauI=600
0.005
tauI=3600
no pipes
No pipes
0.000
10
30
50
70
90
110
left hand side.
Summarising
In case the base case is proportional-integral controlled through regulating the coolant flowrate, like with
coolant temperature control, the same integral action has the advantage that the offset is eliminated,
however short integral times has to be avoided because the process behaves then less stable. The
controller stability maps Kc-UA can be used to determine effortless a stable process. A low inlet coolant
temperature implies that instability is less prominent. However again, large external disturbances can
force the process to a lower steady state.
10.2.3
Mathematical model
Equations 2, 3 and 48 describe the proportional + integral controlled base case.
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
K 1
V = V ,sp K c (Tsp T ) c (Tsp T )dt
I 0
VR
(2)
(3)
(48)
In the previous section, the influence of a proportional controller is examined. It was found that controlling
through regulating the flowrate is possible, however, a suitable Kc is necessary and a lower conversion is
inevitable.
An integral action causes the controller output to change as long as an error exists in the process output.
In the vicinity of the constrained set point, the integral action eliminates the offset. However, in case Kc
has not been chosen properly or the offset becomes substantial, the integral action can make the process
even more unstable i.e. integral wind-up.
In previous section, it was also found that the way a considered disturbance is introduced in the model is
of concern. In case a higher temperature perturbation is imposed i.e. offset < 0 the P-controller could
68
stabilise the system. However, if the temperature was lower than the set point, the controller could
regulate in the wrong direction. The applied trick in which the sign is reversed caused the process variable
the virtually endless revolve around the set point. It could be expected that, due to the offset removing
ability off the I-action, the wobbling effect would decline. Unfortunately, this is not the case according to
several simulations.
The next step is to examine the relation between the proportional gain and the integral time by creating a
suitable stability map in which Kc at the ordinate was viewed versus I at the abscissa.
0.0025
Hopf
0.0020
III
0.0015
0.0010
0.0005
0.0000
600
III
900
1200
Integral time [s]
Figure 85 Process controller stability map with Hopf curves. The proportional gain versus integral time is asymptotic shaped in
which Kc 0.0002 [m3 s-1 K-1] and I.740 [s].
The same pattern arises for this particular stability map: smaller integral time means larger
proportional gain values.
At both sides of the Hopf curve, the system behaves unstable. Compared to Figure 73 in which Kc
3 -1 -1
has been drawn against UA, it became clear that for Kc > 0.000 [m s K ] stability is possible,
although with lower conversion. LOCBIF did not generate such curve or maximum value.
LOCBIF could not find any solutions for small integral time i.e. I.< 740 [s].
3 -1 -1
According to the attained Hopf curve, Kc = 0.0002 [m s K ] if I. which also will be
confirmed in case Kc is drawn against d (Figure 95). Because both sides of the curve are
indicated as III, this value is not relevant anymore.
Figure 85 demonstrates once again that a larger integral time has an increasing stabilising effect.
Nevertheless, merely the integral action cannot provide absolute stability i.e. region I.
Still the same disadvantage stands for that or the reactor temperature or the conversion has to be
chosen.
The proportional gain must be robust enough to be able to adjust the flowrate into the right direction i.e.
set point after both a external disturbance and both to restrain the evolving limit cycles. This has the
disadvantageous effect that even a slight correction has a tremendous effect on the system. The latter is
with respect to the obligation that the throughput variation must be in between a few percent not any more
valid.
The stabilising effect of a large integral time can be explained as follows: apparently, large integral times
make this process executing more smoothly. Consider a small integral time and a large proportional gain.
After a disturbance, the controller determines the new e.g. higher flowrate and because the integral time is
small, in a very short time the flowrate is adjusted. The fresh supported cold feed lowers the reactor
temperature (decreased reaction rate) resulting in fast blow out and subsequently lower conversion. If the
progress of this process is too hasty, extinction is the unwanted and inevitable result. The LOCBIF source
LOCBIF rhs 8 is included in Appendix 7.
69
In the previous section in which merely proportional gain was concerned, LOCBIF had severe problems
finding a solution (Figure 73). Adding the integral action does not improve this solving process.
Based on previous chapters an increasing cooling capacity can have a stabilising effect on the process,
however too large UA can cause extinction T towards Tcool or T0 which can be explained by the fact that a
tremendous amount of heat quickly can be removed from the system which implies fast temperature drop
and consequently reaction extinction.
0.0025
0.0020
0.0015
Hopf
0.0010
0.0005
III
0.0000
10
30
50
70
90
110
Figure 86 Reactor controller stability map. Hopf curve: proportional gain versus cooling capacity.
In view of the base case, once more the existence of multiplicity is confirmed. A large Kc value is
needed to stay away from region III.
3 -1
-1
The same maximum Hopf value is found as for P-control (Figure 73) Kc = 0.00084 [m s K ].
Nevertheless the maximum Hopf value with regards to the cooling capacity is considerably lower
-1 -1
-1 -1
UA = 60 [kJ s K ] compared to the previously found UA = 73 [kJ s K ].
Although region I means stability. The same disadvantage of the lower conversion for one
particular reactor temperature is involved.
Summarising
If the base case is proportional-integral controlled through regulating the throughput, no stable system can
be obtained for small integral time values. Still the same disadvantage stand for that or the reactor
temperature or the conversion has to be chosen. In general, controlling the throughput is compared to
regulating the extent of cooling not preferable.
10.3
Alternatives
In previous section, the influence of the cooling capacity is examined. Additionally, the study to the effects
on the stability of the feed temperature and the coolant temperature can perhaps contribute to the search
for stability. In case the coolant temperature is decreased from Tcool = 441 [K] toward Tcool = 438 [K] the
stability map like Figure 54 changes sincere in which the stable region I is slightly increased. Decreasing
Tcool = 441 [K] beneath Tcool = 435 [K] which represents the transition to the region exhibiting multiplicity
(Figure 19), deteriorate the stability. Increasing Tcool > 441 [K] diminish the stable region. The effect of the
feed temperature on the process is minor. Decreasing T0 to environmental temperature slightly increases
the stable region. In fact the complete stability map are translated upwards in view of the abscissa. Finally,
the introduction of a recycle stream can be useful, because the cooling effect of the fresh feed supply can
26 76
be repressed which has its effect on the controller equation. (Fogler
).
70
Chapter XI
11 CONTROLLED SYSTEM WITH DELAY
In this chapter, the effect of delay on the overall stability (explained in paragraph 5.4) is to be considered.
Due to the ease of LOCBIF, a presumed delay correlation can effortless be implemented to an existing
mathematical model.
The following reactor controller stability maps will be examined:
Table 28 Stability maps delay
Stability map
Reference
Kc versus d
Coolant temperature Figure 88
Coolant flowrate
Throughput
Figure 95
Stability map
Reference
Process capacity
Figure 108
Kc versus R
Figure 110
Kc versus cool
Kc versus cool
Figure 103
Figure 107
11.1
11.1.1
Mathematical model
The mathematical system consist of the following equations:
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
Tcool = Tcool,sp + K c (Tsp Td )
VR
dTd
= T Td
dt
(2)
(3)
(41)
(49)
The overall value of d will be estimated, considering the following partial d s:

68
17
Measuring device and transmitter. According to Roffel and Ding et.al. , the delay caused by the
measuring device varies commonly between 5-10 [s] for liquid reactors and 10-30 [s] for gas-liquid
53
reactors. According to Marlin the typical dynamical response is generally about 0-5 [s].
Heat transfer through pipes. The delay caused by heat transfer process is relatively small due to
the assumption that the heat transfer coefficient for the steel cooling pipes is high and is probably
5-10 [s]. The delay by the controller device is assumed to be 0-5 [s].
71
71
Control valve and activator response time. The valve response time according to Shinskey is
53
practically 0-5 [s]. According to Marlin the final control element response time is 1-4 [s] and for
signal conversion about 0.5 1.0 [s].
Coolant residence time. The cooling fluid is rapidly pumped through the cooling pipes or jacket
and is estimated 0-5 [s] based on 2.3.2 and calculation made in Appendix 3.
Dead zone. Additionally, industrial reactors often are concerned with the existence of a dead zone
in the reactor, in which apparently nothing happens. For this case i.e. the CISTR, it is assumed
that no dead zone exists.
Consequently, the value of the presumed delay can be composed from the points mentioned above.
Roughly 10 < d < 60 [s]. In the succeeding paragraphs d 30 [s] is presupposed. It is not necessary to
determine the exact magnitude of d, primarily due to the lack in information. Furthermore, the purpose of
the introduction of delay is to obtain qualitative comprehension about the effect on the dynamical
behaviour. To determine the effect of the delay on the dynamical behaviour several orbit curves have
been created.
Orbit curves
Primarily, the proportional gain is taken Kc = 0. Figure 87 demonstrates that in case the delay is increased,
the appearance of the oscillations change i.e. the peaks of the measured temperature Td are dimmed and
stretched in which the maximum amplitude emerges later in time than the maximum of the actual reactor
temperature (Figure 89). For very large delay (Figure 87f) the measured reactor temperature is almost
completely diminished i.e. Td (t ) = T .
600
500
400
0
3600 7200
Tim e [s]
700
600
500
400
10800
Figure 87a The measured (delayed)
Temperature [K]
700
Temperature [K]
Tem perature [K]
Figure 87b
600
500
400
0
3600 7200
Tim e [s]
Figure 87d Kc = 0, d = 120 [s].
10800
Kc = 0, d = 30 [s].
600
500
400
Figure 87e
3600 7200
Time [s]
500
400
0
Figure 87c
700
600
10800
Temperature [K]
700
Temperature [K]
Tem perature [K]
reactor temperature.
Kc = 0, d = 5 [s].
3600 7200
Time [s]
700
10800
Kc = 0, d = 600 [s].
3600 7200
Time [s]
10800
Kc = 0, d = 60 [s].
700
600
500
400
0
Figure 87f
3600 7200
Time [s]
10800
Kc = 0, d = 3600 [s].
Td468 [K].
According Equation 49 it is presumed that the reactor temperature is delayed in time. If instead of the
reactor temperature, conversely the coolant temperature is delayed in time according equation 50 the
same results have been obtained. Moreover, the fundamental principle of delay remains the same.
72
cool,d
dTcool,d
dt
= Tcool Tcool,d
(50)
The following reactor controller stability map is created using active LOCBIF bifurcation parameters d and
the proportional gain Kc.
Pro p o rtio n al g ain [-]
M e a sure d te m pe ra ture [K]
700
10
9
Hopf
8
7
6
5
4
3
2
II
1
0
tdelay = 0
tdelay = 30
tdelay = 60
600
tdelay = 120
tdelay = 600
500
400
60
1800
120
2400
3000
T im e [s]
D e la y [s ]
Figure 88 Reactor controller stability map Hopf curves active
3600
Figure 89 Orbit curves of the delayed (measured) reactor
parameters Kc and d. If d > 108 [s] limit cycles will

certainly appear.
temperature base case Kc = 0. The LOCBIF rhs 13 is

included in Appendix 7.
The stability map (Figure 88) unmistakably confirms that in case the apparent delay is increasing
the proportional gain has to be taken considerably larger to preserve stability i.e. no limit cycles.
If d > 108 [s] which is apparently an asymptotic value, limit cycles will certainly exhibit, no matter
the magnitude of Kc, unless the initial point is exactly the steady state, than the temperature
29
remains unaffected. The critical value has been confirmed by Giona and Paladino .
From Figure 88 the minimum proportional gain can be derived Kc(d 0) = 0.9, the same value
which has been found in preceding stability maps (e.g. Figure 23 and Figure 24).
Orbit curves
600
500
400
0
3600
700
Treactor
600
Tdelay
500
400
0
Time [s]
Figure 90a
Kc = 0.9. d = 10 [s].
Temperature [K]
700
Temperature [K]
Temperature [K]
For the subsequent orbit curves, the minimum proportional gain is taken which would provide precise
stability in case the base case is considered without presumed delay.
3600
700
600
500
400
0
Figure 90b
Kc = 0.9. d = 60 [s].
3600
Time [s]
Time [s]
Figure 90c
Kc = 0.9. d = 120 [s].
73
A relatively short delay cannot provoke a stable, minimal controlled system to become unstable, according
to Figure 90a. Increased delay makes the process less stable i.e. limit cycles will exhibit (Figure 90b)
which are originated earlier in case delay is further increased (Figure 90c). In case a very high
proportional gain is selected Kc > 100 and a long delay is presumed d = 120 [s] the system will be
absolutely unstable i.e. limit cycles will occur. This is because this particular base case coolant
temperature multiplicity is not in question. Therefore, extinction or runaway will not take place. The most
unstable situation to fall into is the occurrence of limit cycles.
Stability map (Figure 88) shows that an estimated delay d < 60 [s] has no dramatically effect on the
stability of the process i.e. the proportional gain has only to be slightly increased. For a significantly longer
delay d, the proportional gain has to be considerably larger to avoid limit cycles. If the non-ideal aspect
can be translated into a delay d 60 [s] a proportional gain Kc = 1.3 would be sufficient to eliminate limit
cycles. In this particular case, the larger Kc, the more stability is acquired. In next chapter on the other
hand, it will be shown that in other cases too large Kc causes conversely instability, especially if multiplicity
is considered. Another noticeably effect is the fact that the oscillation frequency increases in case Kc is
increased.
Summarising
Consequently, delay makes a process less stable. In case the non-ideal presumed delay equation is
implemented in the model, the proportional gain has to be increased to remove limit cycles. However, too
large delay cannot be controlled properly to prevent instability.
11.1.2
Mathematical model
The delay is mathematically described with equation 49 and appended to the mathematical models of the
proportional controlled coolant flowrate.
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
dT
dt
V ,cool = V ,cool,sp K c (Tsp Td )
VR
dTd
= T Td
dt
(2)
(3)
(12)
(42)
(49)
Orbit curves
Lets consider the proportionally controlled base case, which behaves narrowly dynamically stable. If the
delay is implemented, primarily limit cycles emerge (Figure 91a and b), which eventually causes extinction
for a slightly larger delay (Figure 91c). A very small different initial condition can drastically affect the
21
behaviour of a system. This has been confirmed by Doherty and Ottino .
74
600
500
400
0
700
Temperature [K]
Temperature [K]
Temperature [K]
700
600
500
400
3600
600
500
400
3600
Time [s]
Figure 91a
700
Time [s]
Tcool,0 = 303 [K], Kc =

Figure
0.0002 [m3 s-1 K-1], d = 1
[s]. Small limit cycles.
91b
3600
Time [s]
Tcool,0 = 303 [K], Kc =

Figure
0.0002 [m3 s-1 K-1], d = 9
[s] limit cycles.
91c
Tcool,0 = 303 [K], Kc =

0.0002 [m3 s-1 K-1], d = 10
[s]. Extinction.
600
500
400
0
Figure 92a
3600
7200
Time [s]
700
Temperature [K]
700
Temperature [K]
Temperature [K]
For the following orbit curves, the non-ideal behaviour has been arbitrarily chosen and presumed to be d
30 [s]. The choice of the proportional gain is rather complicated. In Figure 92a, Kc is evidently not large
enough causing extinction. In Figure 92b Kc is however, more robust though still limit cycles emerge. In
Figure 92c Kc is too strong, in which the combination with delay eventually leads to the lower steady state
i.e. extinction. This is due to the too late intervention of the controller.
600
500
400
10800
Tcool,0 = 303 [K], d = 30 Figure

[s], Kc = 0.0001 [m3 s-1 K1
],. Extinction
3600
7200
Time [s]
92b
700
600
500
400
10800
Tcool,0 = 303 [K], d = 30

[s], Kc = 0.001 [m3 s-1 K1
], limit cycles.
Figure 92c
3600
7200
Time [s]
10800
Tcool,0 = 303 [K], d = 30

[s], Kc = 0.01 [m3 s-1 K-1],
Extinction.
600
500
400
0
0.010
Cool speed [m/s]
700
Coolflow [m 3/s]
Tem perature [K]
In Figure 93a the Kc value is sufficient enough to preserve stability after a disturbance. The disadvantage
is that the stability is very unsteady and additionally the settling time is considerably long and the coolant
flowrate varies enormously to eventually obtain a stable system (Figure 93b and c).
0.008
0.006
0.004
0.002
0.000
3600
Tim e [s]
Figure 93a
Tcool,0 = 303 [K], d = 30 Figure

[s], Kc = 0.0006 [m3 s-1 K1
], Stable system. T =
+10 [K].
3600
0
0
93b
The Coolant flowrate

fluctuates fairly. The
valve actually has to be
halted for several
periods.
3600
Time [s]
Time [s]
Figure 93c
Coolant fluid velocity

supposed the recycle
stream = 10
throughput (equation 84)
75

The following reactor controller stability maps have been constructed using active parameters Kc and the
delay d. These stability maps provide useful information with regard to the stability of the system. In
Figure 94a,b the Hopf curves have been plotted for respectively Tcool,0 = 303 [K] and Tcool,0 = 400 [K].
0.25
0.010
0.009
0.008
0.007
0.006
0.005
0.004
0.003
0.002
0.001
0.000
Tcool,0 = 400
Tcool,0 = 303
Tcool,0 = 435
0.20
0.15
III
0.10
0.05
Hopf
Hopf
0.00
0
30
60
Delay [s]
30
60
90
Delay [s]
Figure 94a,b Reactor controller stability map. Hopf curves for
Figure 94c
Tcool,0 = 303 [K] and Tcool,0 = 400 [K]. Proportional

gain versus delay. The process is stable at the
left-hand side of each Hopf curve and unstable
at the right-hand side.
Reactor controller stability map. Hopf curves for

Tcool,0 = 435 [K]. Proportional gain versus delay. The
process is stable at the left-hand side of each Hopf
curve and unstable at the right-hand side. The
LOCBIF source LOCBIF rhs 15 is included in
Appendix 7.
In case delay is implemented in the model with respect to the proportional controlled coolant
flowrate, a proportional gain can be found in between a particular range, in which a stable process
without limit cycles can be acquired.
The area at the left hand side of each Hopf curve encloses the stable region I, at the right hand
side the unstable region III. For Tcool,0 = 303 [K] this results to Figure 94c.
In Table 29 the range is printed in which stability is concerned. The Kc regarding Figure 93, in
which the system appears to be stable, matches with the results in Table 29.
The delay is less dangerous for larger inlet coolant temperatures.
Table 29 P-control coolant flowrate. Presumed delay d = 30 [s].

Inlet coolant temperature Tcool,0 [K] Kc (min) Kc (max)
303
400
435
0.00036
0.00088
0.011
0.00065
0.0046
0.25
Summarising
In case of proportional controlled coolant flowrate, delay makes a process less stable. In case the nonideal presumed delay added to the mathematical model, the proportional gain has to be increased to
remove limit cycles. However, too large delay means extinction and subsequently the lower steady state.
76
11.1.3
Mathematical model
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
V = V ,sp K c (Tsp Td )
VR
dTd
= T Td
dt
(2)
(3)
(43b)
(49)
The following reactor controller stability map has been constructed using active parameters Kc and the
delay d.
0.0025
Hopf
III
0.0020
0.0015
III
0.0010
0.0005
Hopf
II
0.0000
0
150
300
450
600
Delay [s]
Figure 95
Reactor controller stability map. Proportional controlled throughput with delay. The maximum delay is d = 458 [s].
Figure 95 consists more than one Hopf curve, which confirm the appearance of multiplicity.
3 -1 -1
Kc > 0.00014 [m s K ] indicates the transition from region II to region III.
Larger delay implies larger Kc values. Asymptotic value: d = 458 [s]. Nonetheless, because
throughput control is difficult and transition is always to be concerned with, this value is not as
useful as in case of coolant temperature control, because at the left hand side no region I can be
found.
3 -1
-1
In the previous section through Figure 73, a Kc > 0.002 [m s K ] is pointed out to acquire a
stable process. Figure 95 does not confirm this. This can be due to the fact that LOCBIF
sometime has problems with finding an initial point. This is probable because of the occurrence of
multiplicity.
77
Orbit curves
600
500
400
0
3600
7200
Tim e [s]
0.08
0.06
0.04
0.02
0.00
0
10800
Figure 96 Kc = 0.002 [m3 s-1 K-1]
0.10
Temperature [K]
700
Temperature [K]
Tem perature [K]
To study the effect of delay, several orbits curves are produced with LOCBIF. Again, the proportional gain
3 -1
-1
has been chosen Kc = 0.002 [m s K ] for all the to examine situations. In Figure 96 a delay d = 30 [s]
cannot provoke a system become unstable after an external temperature disturbance. In Figure 97 the
flowrate fluctuates in the first minute and becomes quickly stabilised, however the flowrate fluctuations are
too large. The initial error (constrained temperature disturbance) causes a flowrate deviation of 75%. The
controller corrects primarily in the right direction but oscillates sincerely resulting in a overshoot of +18%
and accordingly 82%, 56% etc and is finally after about 1 hour at the desired steady state, which
adversely is not the base case steady state situation i.e. a steady state with higher temperature although
with lower conversion. In Figure 98 the delay is slightly increased to d = 32 [s] and causes the system to
runaway.
3600
7200
Time [s]
600
500
400
0
10800
Figure 97 Kc = 0.002 [m3 s-1 K-1]
T = 10 [K], d = 30 [s],
Stability. Higher steady
state.
700
3600
7200
Time [s]
10800
Figure 98 Kc = 0.002 [m3 s-1 K-1]
T = 10 [K], d = 30 [s],
A higher stable steady state is
reached R = 350 [s] = 0.47
and T = 472 [K].
= 10 [K], d = 32 [s].
Runaway.
700
Throughput [m 3/s]
Temperature [K]
In Figure 99 the delay is again increased to d = 60 [s] but no disturbance is constrained. The system
remains scarcely stable. The flowrate causes the jump to a higher steady state with lower conversion
(Figure 100).
600
500
400
0
3600
7200
Time [s]
10800
Figure 99 Kc = 0.002 [m3 s-1 K-1]

T = 0 [K] delay = 60 [s]. The
system is considered stable.
0.10
0.08
0.06
0.04
0.02
0.00
0
3600
7200
Time [s]
10800
Figure 100 Kc = 0.002 [m3 s-1 K-1]

T = 0 [K] delay = 60 [s].
A higher steady state however
is reached R = 350 [s] = 0.47
and T = 472 [K].
In case in Figure 99 is simulated with a disturbance, runaway pursues (Figure 101) caused by the flowrate
fluctuations (Figure 102).
78
Throughput [m 3/s]
Temperature [K]
700
600
500
400
0
3600 7200
Time [s]
10800
Figure 101 Kc = 0.002 [m3 s-1 K-1]

T = 10 [K] delay = 60 [s].
Runaway. The LOCBIF
source LOCBIF rhs 17 is
0.10
0.08
0.06
0.04
0.02
0.00
0
3600
7200
10800
Time [s]
Figure 102 Kc = 0.002 [m3 s-1 K-1]

T = 10 [K] delay = 60 [s].
Flowrate varies suchlike that
the process cannot be
controlled reliable.
Summarising
In previous sections is became clear that controlling the throughput had the disadvantage that it is not
possible to preserve the base case situation. Due to the dynamical instability, on the one hand, a large Kc
is necessary, however on the other hand, large Kc cause in particular after a disturbance a considerable
overshoot causing often transition. Even if transition does not occur, the flowrate fluctuation is sincere and
one cannot preserve the commitment to keep the flowrate deviation between only a few percent. The
introduction of delay principally deteriorates the stability in general.
11.2
11.2.1
Mathematical model
The mathematical system consist of the following equations:
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
Kc 1
(T T )dt
Tcool = Tcool,sp + K c (Tsp Td )+
I 0 sp d
VR
dTd
= T Td
dt
(2)
(3)
(46)
(49)
79

A stability map is created (Figure 103) similar to the proportional controlled base case model (Figure 88).
10
tauI=10
9
8
Hopf
curves
II
tauI=30
tauI=60
tauI=120
tauI=180
II
tauI=300
tauI=600
tauI=3600
3
I
2
P-controller
Kc (min)
1
Kc (max)
0
0
30
60
90
120
Delay [s]
Figure 103 Reactor controller stability map with Hopf curves for increasing integral time vales. The dotted line (I ) is similar
to the proportional controlled base case model Figure 88. If d > 108 [s] limit cycles will certainly appear regardless
the extent of delay. At the left-hand side of one particular Hopf curve the behaviour is dynamically stable and ate the
right-hand side ergo unstable. The LOCBIF source LOCBIF rhs 14 is included in Appendix 7.
In the actual stability map, Hopf curves are drawn for increasing values of the integral time. The
dotted line (I ) is similar to the Hopf curve in Figure 88 and coincides with the ordinate for the
formerly found value Kc = 0.9.
To determine the stable and unstable regions in Figure 103, one has to select a particular Hopf
curve in accordance with a specific integral time. The region on the left hand side is the stable
region and on the right hand side the region in which limit cycles occur.
In case the delay is increased, the Kc value rises severely to preserve stability. For (I ) the
same delay d =108 [s] is found (Figure 88) if merely P-controller is considered. For smaller
integral times, the asymptotic value decreases.
Figure 103 proves that the integral action n combination with delay can adversely affect the
stability. The shorter the integral time, the higher the proportional gain. An arbitrary integral time
constant is chosen e.g. I = 60 [s].
Orbit curves
Through Figure 74 it was determined that Kc = 1.35 is robust enough to preserve stability. Introducing an
estimated delay d = 10 [s] disrupts the stability and causes limit cycles (Figure 104a). Decreasing the
integral time close to the assumed delay, nevertheless larger, causes limit cycles (Figure 104b). In case
the delay exceeds the integral time, the system becomes seriously unstable (Figure 104c). Therefore, it is
80
600
500
400
0
3600
7200
Time [s]
10800
700
Te m pe r atur e [K]
700
Tem perature [K]
Temperature [K]
essential to acquire an integral time considerably larger than the assumed delay. Several simulations
indicate that approximately I 10 d has good results, because the integral time amply overlaps the
delay. Subsequently a proportional gain can be selected in which the overshoot is optimal i.e. ratio a/b =
4/1, (5.6.2). The main conclusion is that large Kc values are needed to eliminate the self-sustained
oscillations. However, too large Kc values in combination with delay can provoke instability. The
advantage of the integral action is the fact that the offset is removed. It is desired that the offset is
decreased as soon as possible i.e. small I. However, if the assumed delay exceeds the integral time, a
very unstable process is definitely resulting.
600
500
400
0
700
600
Run aw ay
500
400
3600
Tim e [s]
3600
T im e [s ]
Figure 104a Orbit curve. Kc = 1.35, I Figure 104b Orbit curve. Kc = 1.35, I = Figure 104c Orbit curve. Kc = 1.35, I =
= 60 [s], d = 10 [s].
Dynamical behaviour:
limit cycles.
10 [s], d = 5 [s].
Dynamical behaviour: limit
cycles.
10 [s], d = 20 [s].
Dynamical behaviour: run
away i.e. integral windup.
Summarising
Delay can be the cause of instability during a process, which operates at stable steady state. The
proportional gain is decisive for its stability. In case of coolant temperature control in which no delay is
involved, the integral time has a minor effect on the stability of the process. If delay however is taken into
account, it is crucial that the integral time is considerably larger than the delay. Otherwise, instability is
inevitable.
11.2.2
Mathematical model
The PI-controlled coolant flowrate with delay is mathematical described with:
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
dT
dt
1
K
V ,cool = V ,cool,sp K c (Tsp Td ) c (Tsp Td )dt
I 0
VR
dTd
= T Td
dt
(2)
(3)
(12)
(47)
(49)
81
Orbit curves
600
500
400
0
3600
7200
Time [s]
700
600
500
400
10800
Figure 105a Tcool,0 = 303 [K], Kc =
0.00025 [m3 s-1 K-1], I =

600 [s], T = +10 [K],
d = 1 [s]
Temperature [K]
700
Temperature [K]
Temperature [K]
The destabilising effect of delay will be demonstrated through the following orbit curves. Lets consider the
proportional-integral controlled base case, which primarily behaves dynamical stable. After an external
disturbance, a small delay d = 1 [s] cannot disturb the stable situation (Figure 105a). A larger delay d =
10 [s] disorders the process which exhibits limit cycles (Figure 105b). Larger delay eventually causes
extinction (Figure 105c).
3600
7200
Time [s]
10800
Figure 105b Tcool,0 = 303 [K], Kc =
0.00025 [m3 s-1 K-1], I =

600 [s], T = +10 [K],
d = 10 [s].
700
600
500
400
0
3600
7200
Time [s]
10800
Figure 105c Tcool,0 = 303 [K], Kc =
0.00025 [m3 s-1 K-1], I =

600 [s], T = +10 [K],
d = 30 [s].

For the stability map (Figure 106), active LOCBIF bifurcation parameters Kc and the delay are used. Due
to the fact that LOCBIF couldnt find a (initial) solution, the stability map is exclusively created for Tcool,0 =
303 [K]. The principles remain the same for Tcool,0 = 400 [K] and Tcool,0 = 435 [K].
Hopf
Tcool,0 = 303
0.0025
I= 60 [s]
I= 600 [s]
0.0020
III
I= 3600 [s]
0.0015
0.0010
P-control
0.0005
0.0000
0
Figure 106
15
30
Delay [s]
45
60
Reactor controller stability map with Hopf curves for Tcool,0 = 303 [K]. Proportional gain versus delay. The LOCBIF
source LOCBIF rhs 16 is included in Appendix 7. The left-hand side of each Hopf curve for one particular inlet
coolant temperature indicates stability and vice versa the right hand side instability.
82
The area at the left-hand side of each Hopf curve encloses the stable region I, at the right-hand
side the unstable region III.
The figure shows that for smaller integral time, a considerably larger proportional gain is needed.
The dotted line matches Figure 94 for merely proportional control.
For particular small integral times, a considerable large proportional gain is required to maintain
stability.
Beyond a particular value, instability is inevitable.
Summarising
In case of the integral action is considered combined with delay instability can be the consequence. Large
integral time values increase the stable region.
11.2.3
Mathematical model
The mathematical notation with regards to the PI-controlled throughput with delay:
Eact
d[ A]
= V ([ A]0 [ A]) VR k 0 e RT [ A]
dt
Eact
dT
C PVR
dt
1
K
V = V ,sp K c (Tsp Td ) c (Tsp Td )dt
I 0
VR
dTd
= T Td
dt
(2)
(3)
(48)
(49)
83

For the stability map, Kc has been drawn versus d.
0.0025
tauI=60
tauI=90
0.0020
tauI=300
tauI=600
0.0015
tauI=900
tauI=1200
0.0010
tauI=3600
0.0005
tauI= Pcontroller
Kc (min)
0.0000
0
60
120
180
240
300
Delay [s]
Figure 107
Reactor controller stability map. Proportional + integral control flowrate with delay. The LOCBIF source LOCBIF rhs
18 is included in Appendix 7.
The implementation of delay in the model resulting in the stability map Figure 107 shows that for
small integral time the proportional gain increases dramatically from which can be concluded that
the integral time should be larger than the delay with respect to stability.
Increasing the integral time results in shifting the Hopf curve to the right. In fact I. the curves
3 -1 -1
matches the lower curve of Figure 95 (Kc = 0.00014 [m s K ]). Therefore, the upper Hopf curve
in Figure 95, for varying integral times, cannot be found with LOCBIF. This is due to the existence
of multiplicity. LOCBIF requires suitable initial numerical values, which are sometime not
obtainable.
The asymptotic delay value for Kc states d = 192 [s]. In case d > 192 [s] LOCBIF provides
no solutions for Figure 107. One can expect, despite the lack in Figure 107, that at the right-hand
side for considerably large delay, the process behaves unstable.
In exclusively P-control (Figure 95), an asymptotic value: d = 458 [s] has been found. For PIcontrol, the asymptotic value is d = 192 [s] has been found, whereas for I. the asymptotic d
= 458 [s] should be expected.
Apparently, the combination of delay and integral time together with multiplicity can sometimes
give some simulation problems. This is perhaps due to the appearing of more asymptotic Hopf
curves (maximum values), which might interfere with other Hopf, curves. It could be so that the
asymptotic value of one predominate the other. Therefore, another asymptotic value will be found.
Summarising
Controlling the throughput, using PI-action, combined with presumed delay confirms the formerly
presented conclusion: which is that throughput control compared to cooling control is not a suitable
method to eliminate limit cycles and preserve stability with the same reactor temperature and conversion.
84
11.3
Process capacity
The process capacity is explained in paragraph 5.5.
11.3.1
Process time constant
Mathematical model
The mathematical model for the CISTR states:
Eact
V
d[ A]
= ([ A]0 [ A]) R k0 e RT [ A]
V
dt
V ( H R ) Eact RT
dT
UA
(T Tcool )
R
= (T0 T ) + R
k0 e
[ A]
V CP
V CP
dt
(33)
(34)
In which the process time constant becomes:
R =
VR mequipmentC P,equipment
+
+ ...
V
C P V
(51)

If the reactor controller is interested in the influence of process equipment on the dynamical behaviour of a
process, a plot like Figure 108 can be constructed.
1.2
1.0
0.8
0.6
I
Hopf
II
0.4
0.2
0.0
1000
Average
residence
time CISTR
1100
1200
1300
Process time constant [s]

Figure 108
Reactor controller stability map with Hopf curve. The proportional gain has to be increased slightly to preserve
stability, if the influence of the CISTR equipment is taken into account.
If exclusively the base case is considered, the particular time constant is equal to the average
residence time of the CISTR. R = VR/V = 1000 [s]
If the external elements are concerned with through equation 35, a slower responding process is
to be expected e.g. heating up of the CISTR wall etc. Consequently, a larger proportional gain is
then required to maintain stability.
85
If the process time is increased as Hopf bifurcation parameter, Figure 108 confirms the increase
of Kc. Nevertheless, if the time constant is composed of reactor wall, cooling equipment, stirrer
etc., its maximum value states approximately R 1200 [s] compared with R = 1000 [s], which
means according to Figure 108 that Kc should be increased with 9%. Therefore, the influence of
process equipment moderate affects the value of the proportional gain.
11.3.2
Cooling time constant
In 5.5.2 the coolant differential equation 36 has been postulated in which cool represents the cooling time
constant.
Mathematical model
The mathematical model states:
cool
dTcool
UA
= (Tcool,0 Tcool ) +
(T Tcool )
dt
coolC P,coolVcool
(36)
In which the cooling time constant becomes:
cool =
mequipmentCP,equipment
Vcool
+
+ ...
V ,cool coolCP,cool V ,cool
(38)
Equation 38 clearly demonstrate that cool increases if the physical properties of the equipment is
concerned with. A larger cool implies that the response of the process becomes slower. The latter can
68
perhaps affect the behaviour of the system. According to Roffel through the lack of information, the
physical meaning of equation 38 is difficult to determine. Therefore, analysing exclusively the parameter
cool can provide some information with respect to the dynamical behaviour of the process.
Based on the data printed in Appendix 3 (Table 38) and equation 38 the value of the cooling time constant
can be estimated. It is presumed that the cooling equipment has been constructed from stainless steel.
The results from Table 30 confirm the capability of the cooling equipment to store energy which is obvious
considering that if the reactor temperature increases it takes time to bring the cooling equipment to the
same temperature as the reactor contents resulting in a slower response to changes.
Coolant time constant [s]
300
Tcool
240
With pipes
180
120
60
No pipes
300
350
400
450
Inlet coolant temperature [K]
Figure 109
Coolant time constant versus inlet coolant temperature.
86
Table 30 Time constants coolant. Results (dotted lines in Figure 110).

Time constant coolant Tcool,0 = 303 [K] Tcool,0 = 400 [K] Tcool,0 = 435 [K]
cool(no pipes)
cool(with pipes)
117
205
35
61
5
9
Reactor controller stability map proportional control

If coolant flowrate proportional control is considered, the following stability map can be drawn. In Figure
110 the proportional gain is drawn against the cooling time constant.
0.025
Tcool,0 = 303
Tcool,0 = 400
0.020
Tcool,0 = 435
No pipes 303
0.015
No pipes 400
0.010
No pipes 435
With pipes 303
0.005
With pipes 400

With pipes 435
0.000
0
300
600
Cooling time constant [s]
Figure 110
Reactor controller stability map. Influence cooling time constant for various inlet coolant temperatures.
The magnitude of the cooling time constant becomes considerably relevant in case
For Tcool,0 = 303 [K] cool > 650 [s]
For Tcool,0 = 400 [K] cool > 250 [s]
For Tcool,0 = 435 [K] cool > 0 [s].
The latter can be explained by the fact that due to the minor temperature difference between the
cooling set point and the coolant inlet temperature even the slightest disturbance has a
tremendous effect.
87
Orbit curves
700
Temperature [K]
Temperature [K]
The cooling time constant in fact causes a particular delay. Figure 111 and Figure 112 prove that the
cooling device can cause instability.
600
500
400
0
3600
7200
Time [s]
700
600
500
400
0
10800
Figure 111 Orbit curve. Kc = 0.00025
[m3 s-1 K-1], Tcool,0 = 303 [K],

T = 10 [K], V,cool (set
point) = 5.6 10-3 [m3 s-1].
cool = 56 [s].
3600
7200
Time [s]
10800
Figure 112 Orbit curve. Kc = 0.00025
[m3 s-1 K-1], Tcool,0 = 303 [K],

T = 10 [K], V,cool (set
point) = 2.6 10-3 [m3 s-1].
cool = 98 [s]
Proportional gain [m 3/s/K]
Reactor controller stability map proportional-integral control

0.10
tauI = 60
0.08
tauI = 600
0.06
tauI = 3600
0.04
No pipes
0.02
With pipes
0.00
0
No pipes
Figure 113
300
With pipes
600
900
1200
Reactor controller stability map. The influence of the coolant time constant.Tcool,0 = 303 [K].
The values have been retrieved from Table 30 results (dotted lines in Figure 113, Figure 114 and
Figure 115).
If the I the curves from Figure 110 are again acquired (proportional control only).
The figures demonstrate that in case a small integral time is chosen, the coolant time constant
becomes more crucial.
88
0.10
tauI = 60
0.08
tauI = 600
0.06
tauI = 3600
0.04
No pipes
0.02
With pipes
0.00
0
60
No pipes
Figure 114
120
180
240
300
With pipes
0.10
tauI = 60
0.08
tauI = 600
0.06
tauI = 3600
0.04
No pipes
0.02
With pipes
0.00
0
No pipes
Figure 115
60
With pipes
120
180
240
300
Summarising
In fact, the process time constant acts as a certain delay. If of a particular process the mathematical
description is available, the process engineer can easily modify the time constant of the actual differential
equations, in which the influence of process equipment can be examined. With respect to the base case, if
the base case is adequately controlled i.e. not too low Kc and I, in which a low inlet coolant temperature is
chosen, the process capacity does not affect the overall stability.
89
Chapter XII
12 CONTROLLER CONFIGURATION TUNING
In this chapter, a suitable controller configuration is searched for the base case to eliminate limit cycles
and preserve stability. Therefore the reactor temperature is disturbed T = 20 [K] and an arbitrary chosen
delay d = 30 [s] is concerned with to improve the physical realism. Firstly, the Ziegler Nichols controller
tuning method (explained in 5.6.3) is examined. Accordingly, the proportional gain and the integral time
are varied to study the effect on the process stability. Finally, the important issues and strategy regarding
stability control is to be considered (presented in (5.6.4).
12.1
12.1.1
Ziegler Nichols tuning method
If the Ziegler Nichols methodology is applied for proportional control only, the self-sustained oscillations
have been acquired (Figure 116a). According to the method the integral time is set to infinite I . A
proportional gain KU = 0.96 is obtained and an ultimate period of sustained cycling PU = 16 [min]. The
contradiction appears due to the fact that limit cycles system exhibit even if Kc = 0. Therefore, a
proportional gain is searched for the transition between limit cycles and spiral point. The resolved Ziegler
Nichols controller-tuning configuration becomes:
Kc = 0.43 [-]
and
I = 760 [s].
510
Temperature [K]
Tem perature [K]
According to Figure 116b this configuration is not robust enough to eliminate limit cycles. In this particular
case the Ziegler Nichols tuning method is not suitable to provide a stable system. Apparently, this tuning
method does not take into account dynamically instabilities like limit cycles and is developed for static
instability problems only.
485
460
435
410
60
75
90
105
Tim e [min]
120
700
600
500
400
0
3600
7200
Time [s]
10800
Figure 116a KU = 0.98 [-], I [s], T Figure 116b Kc = 0.43 [-] I = 760 [s].,
= +20 [K], d = 30 [s], PU =
15.2 [min].
T = +20 [K], d = 30 [s].
90
12.1.2
Controller parameter variation
The following directives can provide a reasonably appropriate controller configuration, which preserves
stability:
The generally accepted -decay criteria (5.6.2) is qualitatively applied to determine a

proportional gain, which minimises the settling time i.e. time needed for the response to settle
within 5% of the desired value.
Because delay can be the occasion of instability during a process, the integral time must be
considerably larger than the delay.
Overshoot cannot be avoided, however peaks, which are too strong, are undesirable.
The frequency of oscillation cannot be too high due the practical limitations of the actual process
controller.
Table 31 Results various controller parameters. Step disturbance T =20 [K], delay d = 30 [s].
Proportional gain Integral time
I = 60 [s]
I = 600 [s] I = 3600 [s]
Kc = 1 [-]
Figure 117 Figure 118 Figure 119
Kc = 5 [-]
Kc = 10 [-]
Kc = 100 [-]
1200
2400
Tim e [s]
510
500
490
480
470
460
450
440
430
420
410
3600
Figure 117 Kc = 1, I = 60.
Temperature [K]
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
Tem perature [K]
Figure 117-Figure 119 show visibly that the proportional gain Kc = 1 [-] is not robust enough to reduce the
self-sustained oscillations.
1200
2400
Time [s]
510
500
490
480
470
460
450
440
430
420
410
3600
Figure 118 Kc = 1, I = 600.
1200
2400
Time [s]
3600
Figure 119 Kc = 1, I = 3600 [s].
1200
2400
Tim e [s]
Figure 120 Kc = 5, I = 60 [s].
3600
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
Tem perature [K]
In case the proportional gain is increased to Kc = 5 [-] the limit cycles can be repressed. However, a small
integral time results in a long settling time (Figure 120). An integral time I = 600 [s] conversely offers
rapidly a stable system with a small settling time and without too excessive overshoot (Figure 121). Larger
integral time has no improvement effect on the stability (Figure 122).
1200
2400
Time [s]
Figure 121 Kc = 5, I = 600 [s].
3600
510
500
490
480
470
460
450
440
430
420
410
0
1200
2400
Time [s]
3600
Figure 122 Kc = 5, I = 3600 [s].
91
1200
2400
Tim e [s]
510
500
490
480
470
460
450
440
430
420
410
3600
Figure 123 Kc = 10, I = 60 [s].
Temperature [K]
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
Tem perature [K]
A larger proportional gain Kc = 10 [-] Figure 123- Figure 125, decreased the settling time. Although this
improvement, it is not recommended tot chose a too large proportional gain due to the fact that the in case
the assumed delay d = 30 [s] has been underestimated, the stability of the process could be at risk.
1200
2400
Time [s]
510
500
490
480
470
460
450
440
430
420
410
3600
Figure 124 Kc = 10, I = 600 [s].
1200
2400
Time [s]
3600
Figure 125 Kc = 10, I = 3600 [s].
Further increasing the proportional gain (Figure 126-Figure 128) however decreases the settling time
nevertheless raises the oscillation frequency. A process controller has a practical limited response time
after an exhibiting disturbance, therefore too large Kc are unlikely.
Literature
The literature raises several directives to acquire controller parameters.
72
According to Shinskey I > 4 d. to achieve a self-regulating process, hence I > 120 [s].
I =
mR C P
UA
(52)
59
1200
2400
Tim e [s]
Figure 126 Kc = 100, I = 60 [s].
3600
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
Tem perature [K]
According to Perry the magnitude of the integral time can be estimated using correlation 52,
consequently: I 160 [s]. These values are clearly too small to acquire stability. An integral time I = 600
[s] however is more appropriate value for this particular base case.
1200
2400
Time [s]
Figure 127 Kc = 100, I = 600 [s].
3600
510
500
490
480
470
460
450
440
430
420
410
0
1200
2400
Time [s]
3600
Figure 128 Kc = 100, I = 3600 [s].
92
12.2
12.2.1

3
-1
-1
A proportional gain KU = 0.00052 [m s K ] is acquired and an ultimate period of sustained cycling PU = 9

[min] (Figure 129). Therefore, the Ziegler Nichols controller-tuning configuration becomes:
3
-1
-1
Kc = 0.00023 [m s K ]
I = 444 [s].
and
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
Temperature [K]
This configuration is not robust enough to eliminate limit cycles (Figure 130), although the provoked
disturbance can be removed. Apparently, the Ziegler Nichols tuning method is suitable to provide a stable
system in which external disturbances can be eliminated. Nevertheless does not take into account the
internal disturbances like exhibiting limit cycles. Therefore, the Ziegler Nichols tuning method has to be
applied with caution.
10
20
Time [m in]
30
700
600
500
400
0
3600
7200
Time [s]
10800
Figure 129 KU = 0.00051908 [m3 s-1 K-1], Figure 130 Kc = 0.000234 [m3 s-1 K-1]
I [s], PU = 9 [min].
12.2.2
I = 444 [s].
To acquire a suitable controller configuration, the distinct parameters are varied and evaluated on basis of
its dynamical behaviour.
Table 32 Results various controller parameters. Presumed delay d = 30 [s].
Proportional gain, perturbation
Integral time
I = 60 [s]
I = 600 [s] I = 3600 [s]
3 -1 -1
Figure
131
Figure
132 Figure 133
Kc = 0.00020 [m s K ] T =20 [K]
3 -1 -1
Kc = 0.00025 [m s K ] T =20 [K] Figure 134 Figure 135 Figure 136
3 -1 -1
3 -1 -1
3
-1
-1
Figure 131 and Figure 132 show clearly that the proportional gain Kc = 0.0002 [m s K ] is not robust
enough to reduce the self-sustained oscillations. In case of Figure 133 the integral time is too large
resulting in extinction.
93
1200
2400
Tim e [s]
510
500
490
480
470
460
450
440
430
420
410
3600
Figure 131 Kc = 0.0002 [m3 s-1 K-1],
Temperature [K]
Temperature [K]
Tem perature [K]
510
500
490
480
470
460
450
440
430
420
410
1200
2400
Time [s]
3600
Figure 132 Kc = 0.0002 [m3 s-1 K-1],
I = 60 [s].
510
500
490
480
470
460
450
440
430
420
410
0
1200
2400
Time [s]
3600
Figure 133 Kc = 0.0002 [m3 s-1 K-1],
I = 600 [s].
I = 3600 [s].
3
-1
-1
1200
2400
Tim e [s]
510
500
490
480
470
460
450
440
430
420
410
3600
Figure 134 Kc = 0.00025 [m3 s-1 K-1],
Temperature [K]
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
Tem perature [K]
In case the proportional gain is increased to Kc = 0.00025 [m s K ], the overshoot is reduced. However,
limit cycles still emerge (Figure 134 and Figure 135). An integral time I = 3600 [s] again is too large
(Figure 136).
1200
2400
Time [s]
3600
Figure 135 Kc = 0.00025 [m3 s-1 K-1],
I = 60 [s].
510
500
490
480
470
460
450
440
430
420
410
0
-1
3600
Figure 136 Kc = 0.00025 [m3 s-1 K-1],
I = 600 [s].
1200
2400
Time [s]
I = 3600 [s].
-1
1200
2400
Tim e [s]
3600
Figure 137 Kc = 0.0003 [m3 s-1 K-1],

I = 60 [s].
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
Tem perature [K]
A larger proportional gain Kc = 0.0003 [m s K ], (Figure 137 and Figure 138) is robust enough to
eliminate the constrained disturbance and the self-sustained oscillations i.e. spiral point. However if the
integral; time becomes too long instability is inevitable (Figure 139).
1200
2400
Time [s]
3600
Figure 138 Kc = 0.0003 [m3 s-1 K-1],
510
500
490
480
470
460
450
440
430
420
410
0
1200
2400
Time [s]
3600
Figure 139 Kc = 0.0003 [m3 s-1 K-1],
I = 600 [s],
I = 3600 [s].
3
-1
-1
Further increasing the proportional gain to Kc = 0.0005 [m s K ], can preserve stability for a small
integral time (Figure 140) however for larger integral times (Figure 141 and Figure 142) the reaction
irrevocably extinguishes.
94
1200
2400
Time [s]
510
500
490
480
470
460
450
440
430
420
410
3600
Figure 140 Kc = 0.0005 [m3 s-1 K-1],
Temperature [K]
Temperature [K]
Temperature [K]
510
500
490
480
470
460
450
440
430
420
410
1200
2400
Time [s]
3600
Figure 141 Kc = 0.0005 [m3 s-1 K-1],
I = 60 [s].
510
500
490
480
470
460
450
440
430
420
410
0
1200
2400
Time [s]
3600
Figure 142 Kc = 0.0005 [m3 s-1 K-1],
I = 600 [s].
I = 3600 [s].
In case of coolant temperature control, a larger integral time implies more stability. Due to the extra
cooling equation 36 multiplicity is concerned, which can cause the reaction to extinguish.
The best controller configuration regarding the PI-control coolant flowrate has been given in Figure 138
3 -1 -1
i.e. Kc = 0.0003 [m s K ], I = 600 [s].
The constrained disturbance has a tremendous effect on the controlled coolant flowrate. However due to
the considerable cold coolant inlet temperature, the maximum fluid velocity will not become critical. In
-1
case of a recycle stream (10) vcool < 2 [m s ].
3
-1
-1
Considered the best controller configuration Kc = 0.0003 [m s K ], I = 600 [s] (Figure 138) the
dynamical behaviour of the coolant flowrate is portrayed in Figure 143. Despite the large overshoot, in
case of steady state and minor disturbances, the coolant flowrate will vary much less.
Coolflow [m 3/s]
0.025
0.020
0.015
0.010
0.005
0.000
0
Figure 143
1200
2400
Time [s]
3600
Kc = 0.0003 [m3 s-1 K-1], I = 600 [s], T = +20 [K], d = 30 [s].
12.3
12.3.1

3
-1
-1
A proportional gain KU = 0.00012 [m s K ] and an ultimate period of sustained cycling PU = 52 [min] is

acquired (Figure 144). According to the Ziegler Nichols method the controller-tuning configuration
becomes:
3
-1
-1
Kc = 0.000054 [m s K ]
and
I = 2600 [s].
In previous chapters, it has been obvious that the Ziegler Nichols method was not suitable to determine
appropriate controller setting. In addition, for throughput control, this method is not applicable (Figure
145).
95
Temperature [K]
Temperature [K]
700
600
500
400
0
60
120
Time [m in]
700
600
500
400
180
Figure 144 KU = 0.00012 [m3 s-1 K-1],
10800
Figure 145 Kc = 0.000054 [m3 s-1 K-1]
I [s], T = +20 [K],

d = 30 [s], PU = 52 [min].
12.3.2
3600
7200
Time [s]
I = 2600 [s], T = +20 [K],

d = 30 [s].
Consider a constrained disturbance in the reactor temperature T and a presupposed delay of d = 30 [s].
The proportional gain and the integral time are varied to study the effect on the stability.
Table 33 Results various controller parameters. Presumed delay d = 30 [s].
Proportional gain, perturbation
Integral time
I = 60 [s]
I = 600 [s] I = 3600 [s]
3 -1 -1
Figure
146
Figure
147 Figure 148
Kc = 0.0001 [m s K ] T =20 [K]
3 -1 -1
3 -1 -1
Kc = 0.01 [m s K ]
T =20 [K] Figure 152 Figure 153 Figure 154
3 -1 -1
Kc = 0.1 [m s K ]
T =20 [K] Figure 155 Figure 156 Figure 157
3 -1 -1
Kc = 0.01 [m s K ]
T =2 [K]
3
-1
-1
1200
2400
Tim e [s]
3600
Figure 146 Kc = 0.0001 [m3 s-1 K-1],

I = 60 [s].
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
Tem perature [K]
Figure 146, Figure 147 and Figure 148 show visibly that the proportional gain Kc = 0.0001 [m s K ] is
not robust enough to reduce the self-sustained oscillations. Actually, the temperature perturbation causes
the reaction immediately to runaway.
1200
2400
Time [s]
510
500
490
480
470
460
450
440
430
420
410
3600
Figure 147 Kc = 0.0001 [m3 s-1 K-1],
1200
2400
Time [s]
3600
Figure 148 Kc = 0.0001 [m3 s-1 K-1],
I = 600 [s].
I = 3600 [s].
3
-1
-1
In case the proportional gain is increased with a factor 10 to Kc = 0.001 [m s K ] (Figure 149) initially the
controller corrects the flowrate in the right direction towards the base case value, however cannot avoid
runaway directly after the first peak. Increasing the integral time has a considerable stabilising effect
(Figure 150). Actually, after a reasonable short time the reaction is stabilised with the desired reactor
temperature T = 468 [K] however with a poorer conversion = 0.48 [-] which is 25% lower than the
required base case conversion = 0.68 [-]. The lower conversion is due to the fact easily another steady
3 -1
state is reached corresponding with a throughput V = 0.001 [m s ] which is 5 times lower than the base
case value. An integral I = 3600 [s] is conversely too large and causes once more runaway (Figure 151).
96
1200
2400
Tim e [s]
510
500
490
480
470
460
450
440
430
420
410
3600
Figure 149 Kc = 0.001 [m3 s-1 K-1],
Temperature [K]
Temperature [K]
Tem perature [K]
510
500
490
480
470
460
450
440
430
420
410
1200
2400
Time [s]
510
500
490
480
470
460
450
440
430
420
410
3600
Figure 150 Kc = 0.001 [m3 s-1 K-1],
-1
3600
Figure 151 Kc = 0.001 [m3 s-1 K-1],
I = 600 [s]. = 0.48,

3 -1
V = 0.012 [m s ].
I = 60 [s].
1200
2400
Time [s]
I = 3600 [s].
-1
1200
2400
Tim e [s]
510
500
490
480
470
460
450
440
430
420
410
3600
Figure 152 Kc = 0.01 [m3 s-1 K-1],
Temperature [K]
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
Tem perature [K]
Increasing Kc 10 times to Kc = 0.01 [m s K ] however with a short integral time implies substantial
instability (Figure 152). A larger integral time conversely preserves stability, unfortunately again with a
lower conversion (Figure 153). Larger integral times do not change anything in view of the stability (Figure
154). These figures prove that the integral time should not be too small i.e. the influence of the delay can
be reduced by choosing a considerably larger integral time, which overlaps abundantly the delay. Because
a change in the flowrate, caused by the proportional gain, affects directly the dynamical behaviour of the
system, the integral time spreads the large peaks in which risk for transition to other steady states is
decreased.
1200
2400
Time [s]
510
500
490
480
470
460
450
440
430
420
410
3600
Figure 153 Kc = 0.01 [m3 s-1 K-1],
3600
Figure 154 Kc = 0.01 [m3 s-1 K-1],
I = 600 [s], = 0.48.
I = 60 [s].
1200
2400
Time [s]
I = 3600 [s], = 0.48.
1200
2400
Tim e [s]
Figure 155 Kc = 0.1 [m3 s-1 K-1],

I = 60 [s].
3600
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
Tem perature [K]
Further increasing the proportional gain (Figure 155) once more causes runaway, although it takes more
time to reach that critical situation. Increasing the integral time provides stability although with high
oscillation frequencies, which are practically unlikely to realise i.e. opening and closing the throughput
valve (Figure 156). Further increasing the integral time has no effect on the stability (Figure 157).
1200
2400
Time [s]
Figure 156 Kc = 0.1 [m3 s-1 K-1],
I = 600 [s, = 0.48.
3600
510
500
490
480
470
460
450
440
430
420
410
0
1200
2400
Time [s]
3600
Figure 157 Kc = 0.1 [m3 s-1 K-1],
I = 3600 [s], = 0.48.
97
1200
2400
Tim e [s]
3600
Figure 158 Kc = 0.01 [m3 s-1 K-1],

I = 60 [s] T = +2 [K].
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
510
500
490
480
470
460
450
440
430
420
410
Temperature [K]
Tem perature [K]
The previous cases, demonstrate that a considerable disturbance on the one hand can cause instability in
the reactor and on the other hand imply stability however corresponding with a lower conversion i.e. the
base case situation can not be preserved.
If the actual disturbance is decreased T = +2 [K], the steady state temperature can be remained. (Figure
158, Figure 159 and Figure 160). Unfortunately, due to the unstable operating point the conversion is only
temporarily close to the base case steady state value and is irrevocable attracted to the closest steady
state, i.e. the steady state corresponding with lower throughput and consequently lower conversion.
1200
2400
Time [s]
3600
Figure 159 Kc = 0.01 [m3 s-1 K-1],

I = 600 [s] T = +2 [K].
510
500
490
480
470
460
450
440
430
420
410
0
1200
2400
Time [s]
3600
Figure 160 Kc = 0.01 [m3 s-1 K-1],

I = 3600 [s] T = +2 [K].
Figure 159 and Figure 160 show that the reactor temperature can be remained at base case temperature
however with a lower conversion.
Finally, if the two possible steady state situation in which or the reactor temperature or the conversion is
maintained are compared:
3
-1
Production = V [A]0 = 0.005 5000 0.68 = 17 [m s ]

3 -1
Production = V [A]0 = 0.012 5000 0.48 = 28 [m s ]
It is apparent that the controller, which uses the reactor temperature as set point, has the largest
production, although it is not certain if the larger flowrate is allowed. Besides, in case of conversion as set
point the reactor temperature is higher. In view of safety, the latter is perhaps not allowed. Nevertheless, it
is not the scope of this report to decide which choice must be made.
98
12.4
Base case controller configuration
To obtain a suitable controller configuration, on the one hand a robust proportional controller gain is
needed to restrain evolving limit cycles, however in addition, too large Kc can cause a certain overshoot,
which can cause unwanted or even dangerous situations. If multiplicity is concerned with, transition is
inevitable like if the throughput control.
The integral time has to be chosen in such a way that the assumed delay is inferior. A too large integral
time has hardly any improvement towards the stability.
The conventional controller tuning methods like the Ziegler Nichols method are not adequate to adopt for
processes with distinct dynamically unstable behaviour like the base case.
Consider the base case situation, with an initial temperature disturbance of T = 20 [K] and a presumed
delay d = 30 [s] to improve the physical realism, the following controller configuration preserved a stable
process:
Table 34 Suitable PI-controller settings coolant temperature control
Description
Value
Presumed delay
Perturbation
Proportional gain
Integral time
Desired base case reactor temperature
Required conversion
d = 30 [s]
T = 20 [K]
Kc = 5 [-]
I = 600 [s]
T = 468 [K]
= 0.68 [-]
Table 35 Suitable PI-controller settings coolant flowrate control

Description
Value
Inlet coolant temperature
Presumed delay
Perturbation
Proportional gain
Integral time
Required conversion
Tcool,0 = 303 [K]

d = 30 [s]
T = 20 [K]
Kc = 0.0003 [m3 s-1 K-1]
I = 600 [s]
T = 468 [K]
= 0.68 [-]
Table 36 The best achievable PI-controller settings throughput control

Description
Value
Presumed delay
Perturbation
Proportional gain
Integral time
Required conversion
Throughput
d = 30 [s]
T = 2 [K]
Kc = 0.01 [m3 s-1 K-1]
I = 600 [s]
T = 468 [K]
= 0.48 [-]
3 -1
V = 0.012 [m s ]
It is obvious that these values have no practical significant importance, merely give a qualitative
impression of the magnitude.
99
Chapter XIII
13 THE EFFECT OF FOULING ON THE BASE
CASE STABILITY
13.1
Introduction
In existing industrial processes, reactor vessels with exothermic chemical reaction are cooled
27
continuously. Simultaneously, the cooling apparatus is contaminated called fouling (Foust ). As a
consequence of fouling the cooling capacity can decrease during the process. Consider the UA reduction
per time.
dUA
= UAf
dt
(53)
In correlation 53, f stands for the fouling factor i.e. the amount of reduction of the cooling capacity per unit
of time. An arbitrary chosen f is to be assumed.
55
50
45
40
35
30
0
12
18
24
30
36
Time [h]
Figure 161 Presuppose that every hour 1% of the cooling

capacity is reduced. According equation 53
f = -(1-0.99)/3600 = -2.8 10-6 [s-1].
Consider the controlled base case process, operating at a stable point. From the moment t = 0 the
process of fouling depletion according equation 53 is to be considered. The simulation starts using the
default base case values, however with an initial temperature disturbance of T = +10 [K]. The following
systems will be considered:
1.
2.
3.
4.
Coolant temperature proportional control

Coolant temperature proportional-integral control
Coolant flowrate proportional control
Coolant flowrate proportional-integral control
100
13.2
Coolant temperature proportional control
The process is proportional controlled (Kc = 0.95). It has been previously found through Figure 27 that the
minimum required proportional gain to provide stability is Kc = 0.9. During the process, the cooling
capacity declines as has been portrayed in Figure 161. Initially, Figure 162 confirms that the controller is
suitable to eliminate the disturbance. Nevertheless, after approximately 6 hours, the system becomes
unstable and exhibits undesired limit cycles. The example demonstrates that a controlled real process,
which is marginal stable, can definitely become unstable in case the process variables change.
1 Coolant temperature proportional-integral control
540
Through Figure 74 it was determined that Kc = 1.35 is sufficient to preserve stability no matter the value of
the integral time. Nonetheless, an arbitrary integral time constant is chosen I = 60 [s]. The controller
initially removes the constrained disturbance, but after approximately 7 hours, the system becomes
unstable which has been shown in Figure 163.
490
440
540
490
440
0
9 10
Time [h]
9 10
Time [h]
Figure 162 Tcool P-control. Initial step disturbance T = 10 [K].

After approximately t = 6 [h], the process becomes
unstable and exhibits large limit cycles. The LOCBIF
source LOCBIF rhs 9 is included in Appendix 7.
Figure 163 Tcool PI-control. Initial step disturbance T = 10 [K].

After approximately t = 7 [h], the process becomes
unstable. The LOCBIF source LOCBIF rhs 9 is
Because the offset is not removed in case a P-controller is applied which is clearly visible in Figure 162,
the process becomes unstable more early according Figure 163. Correspondingly, in case of coolant
temperature control, the advantage of the integral action is demonstrated.
13.3
Coolant flowrate proportional control

3
-1
-1
A proportional gain Kc = 0.00025 [m s K ] is selected (derived from Figure 37) at which narrowly a
stable system is achieved. The proportional controller does not remove the offset which is confirmed in
Figure 164a. After approximately 50 hours, the process becomes unstable and exhibits limit cycles. Due to
the reduced cooling capacity, more heat has to be removed through the coolant flowrate, which will
consequently increase (Figure 164b) whereas the coolant temperature decreases (Figure 164c).
101
13.4
Coolant flowrate proportional-integral control

3
-1
-1
540
The same Kc = 0.00025 [m s K ] is selected at which a stable system is acquired. The integral action
clearly removes the offset (Figure 165a). Due to UA decline, the coolant flowrate has to be increased to
remove the surplus heat Figure 165b although less compared to P-control only. After approximately 32
hours, the system becomes unstable. In case of proportional control, only it took 50 hours before the
system became unstable. The proportional + integral action (I 600 [s]) makes the system more quickly
unstable compared the proportional control action exclusively (I ). This can be derived from Figure 80
in which Kc has been plotted against I.
490
440
540
490
440
0
12 18 24 30 36 42 48
12
Time [h]
Figure 164a Tcool,0 = 303 [K]. Initial step disturbance T = 10
0.010
0.009
0.008
0.007
0.006
0.005
0.004
0.003
0.002
0.001
0.000
0
12 18 24 30 36 42 48
Time [h]
Figure 164b V,cool P-control. Initial step disturbance T = 10

[K]. The coolant flowrate increases due to UA
decline.
24
30
36
Time [h]
Figure 165a Tcool,0 = 303 [K], I = 600 [s]. Initial step disturbance
T = 10 [K]. The process becomes unstable after
approximately after t = 32 [h].
Coolant flowrate [m3/s]
Coolant flowrate [m3/s]
[K]. The process becomes unstable after

approximately t = 50 [h].
18
0.010
0.009
0.008
0.007
0.006
0.005
0.004
0.003
0.002
0.001
0.000
0
12
18
24
30
36
Time [h]
Figure 165b Tcool,0 = 303 [K], I = 600 [s]. Initial step disturbance
T = 10 [K]. The coolant flowrate increases less
than P-control.
102
510
500
490
480
470
460
450
440
430
420
410
0
12 18 24 30 36 42 48
Time [h]
Figure 164c Tcool,0 = 303 [K]. Due to the increasing coolant

flowrate the coolant temperature declines. The
Appendix 7. The influence of the heat capacity of
the tube wall is discussed by Roffel68.
510
500
490
480
470
460
450
440
430
420
410
0
12
18
24
30
36
Time [h]
Figure 165c Tcool,0 = 303 [K]. Due to the increasing coolant

flowrate the coolant temperature declines. The
Appendix 7.
103
CONCLUSIONS
Limit cycles which exhibit in a cooled CISTR with irreversible exothermic reaction A P can completely
be eliminated through the application of a proportional-integral controller. This can be achieved
successfully if the PI-controller regulates the extent of cooling.
The key to success is the software program LOCBIF in which distinct stability maps and orbit curves can
be produced based on the mathematical descriptions. A stability map can give valuable information about
the static and dynamical behaviour of a considered process.
Throughput control is not a suitable control method due to the fact that adjusting the flowrate decreases
the reactor temperature by the supply of cold feed but on the other hand increases the reaction rate due to
the fresh reactant. To match both internal disturbances (self-sustained oscillations) and external
disturbances (irregularities) the controller must be robust enough. Nevertheless, too drastic manipulation
will on the contrary disrupt the dynamics of a process. If a dynamical unstable process is considered, like
the base case, the proportional gain is certainly too large and can cause transition. An inevitable lower
proportional gain means either the reactor temperature is chosen as set point with lower conversion, or
the conversion is chosen with a higher reactor temperature. With respect to the base case with throughput
control, no appropriate control configuration could be obtained.
Cooling control implies practically the regulation of the temperature of a cooling fluid. In this report, two
distinguished models have been examined. At first the instantaneously manipulation of the coolant
temperature and secondly the implementation of a cooling differential equation which the latter describe
the coolant temperature as function of the coolant flow. Considered the first model, it is relatively simple to
elucidate limit cycles and preserve stability. In the worst case scenario, limit cycles emerge. Regarding the
second model, suitable control can be achieved, although the configuration of the controller parameters
must be tuned adequately. Inaccurate control can imply for instance too much heat withdrawal from the
reactor causing extinction.
To enlarge the physical realism, the apparent delay is introduced. Small delay cannot destabilise a proper
controlled process. Conversely, large delay provokes dynamical instability and if multiplicity is involved
even static instability (extinction or runaway). Beyond a certain delay, instability is inevitable in which no
controller can provide a stable process. The simulated delay is apparently less dangerous for higher inlet
coolant temperatures. In the contrary, the influence of process equipment on the stability of a process,
which in fact represents a certain delay, is less affecting the stability for low coolant inlet temperatures.
The integral time removes any offset and is therefore a welcome supplement. The conditions are that the
integral time constant is considerable larger than the apparent delay.
Traditional control configuration tuning methods like the Ziegler and Nichols technique are not appropriate
to deal with dynamical instabilities.
80
Large cooling capacity, in general has a positive effect on the stability of a process . However too large
UA values, which are practically not thinkable due to economic aspects, can contribute to a process in
which too much heat is removed from the system if the controller manipulates too drastic.
If a controller is accurate and can safely preserve stability, a larger proportional gain may imply less UA
i.e. installed heat transfer area (cost reduction).
In the literature in general, stability maps are composed in which one focuses on process variables such
as the reactor temperature or conversion. In this report, such stability maps are called reactor designer
stability maps. This is because the influence of important process parameters on the process stability can
be analysed in the stability map. Due to the complex structure of the mathematical description and the fact
that more than 2 variables are varied in a 2D plot, the curves in these stability maps become rather
entangled and are therefore often difficult to interpret, in which the advantage of such maps in fact
perishes. The directly study of particular process parameters such as proportional gain or cooling capacity
through Hopf curves is less elaborate and often much more easy to interpret. Besides asymptotic limits
can also be obtained. Additionally, stable and unstable regions can be indicated easier. Therefore,
another stability map is introduced in this report, called the reactor controller stability map, in which easily
one can determine the controller settings in which stability can be maintained or acquired, in particular
after the process has been changed.
104
RECOMMENDATIONS
Validation
Although in this report it has been demonstrated that the base case, which exhibits limit cycles, can be
controlled adequately, it is of crucial importance that the developed models will be validated for real
physical existing processes.
This report is part of the research to the stability of process-operation of gas-liquid reactors. The next and
23 22 78 26
logical step in research is the extension of the model with respect to a two-phase reactor.
Model enhancement and extension

17 19 38 76
In this report, throughput control for one reactant is not recommended. From literature
, it is know
that in case more reactants are concerned with, which the latter is practically mostly the case, flowrate
control of one or more inlet flowrate streams might contribute to stabilising a process.
The application of n CISTRs in series with the same total volume n VR,i may be in advantage compared to
46 55 26 68 67
the use of 1 CISTR. This has to be investigated.
The volume in the CISTR can be modified by changing the liquid height in the reactor. Consequently, the
heat transferred through the cooling pipes changes. Through regulating the height, flowrate control is
25
perhaps possible .
In case of cooling control, the assumption is made that the coolant temperature is constant. In practise,
this is not always the case. If the flow pattern though the cooling device (e.g. pipes) is taken into account,
60 59 4 37 55 26 68
one can make the outgrowth probably more realistic.
The accuracy regarding the coolant temperature can be improved in case the logarithmic mean
80 68
temperature according equation 17, is implemented in the mathematical model.
37 4 49
In the models, the assumption is made that the reactor is perfectly mixed. According to
the stirrer
affects both the conversion as the heat capacity. A more advanced description of this influence on the
conversion and UA value can be an improvement to the model.
Process control improvement

Derivative control can make a control much more stable than a controller with merely proportional or
proportional-integral action. Derivative control involves an extra integral equation. Unfortunately,
implementation of derivative control action is not possible in LOCBIF, because this package is not adopted
for the use of integral equations; conversely, the UNIX package AUTO has special features, which can
handle integral equations. The PID controller can then be examined.
In case of substantial deviation of essential system variables, feed forward control can be a suitable
74 9 10 14 75 35 55 59 72 53
solution, or conceivably the combination of feed forward and feedback.
74
All the considered controller tuning methods have been developed to deal with static instability. Methods
have to be acquired or developed to determine suitable controller parameters which preserve a both static
as dynamical stable process. The new and advanced kind of process control is the non-linear process
35 35 74 53 10 51 16
control which provides extraordinarily new features.
The control objective of the controller is to keep the temperature T of the reacting mixture constant at a
desired value. Possible disturbances to the reactor include the feed temperature T0 and the coolant
temperature Tcool. The temperature in the reactor responds much faster to changes in T0 than to changes
74
in Tcool, in case of a coolant flow control method.
105
NOMENCLATURE
Symbol
Description
Dimension
A
As
c
c(t)
CP
CP,cool
Dimp
Dpipe
Dvessel
E
Eact
f
(x)
g
Gc(s)
~
H subscript
Heat transfer or exchanger area

Surface area available for heat transfer per unit length
Controlled variable (equation 25)
Output signal (equation 26)
Heat capacity
Heat capacity of cooling fluid
Impeller diameter
Pipe diameter / thickness
Vessel diameter
Total energy of reacting mixture
Activation energy (equation 1)
Fouling factor, rate of UA depletion (equation 53)
Mathematical function description
Gravitational acceleration
Transfer function (for a controller)
Partial molar enthalpies of component subscript (equations 61)
HR
k0
Kc
Ke
Ki
kR,m,n
KU
Kx
Lpipe
L2pipes
M
nsubscript
Nimp
NTU
NTUcool
NuR
Nucool
pi
P
P
PrR
Prcool
PU
Q
Qp
Qfeed
Qreaction
Qtransfer
Qcool
R
ReR
Recool
Ri
ri
S
Reaction heat / enthalpy of reaction based on R

[J mol ]
-1
Arrhenius pre-exponential factor or frequency factor (equation 1)
[s ]
3 -1 -1
3 -1
Proportional control gain
[-] or [m s K ] or [m s ]
-1 -1
Estimated gain conversion parameter (equations 42 and 47)
[m s K ]
Integral control gain
[-]
3(m+n-1)
(m+n-1)
Reaction rate constant
[m
/ mol
s]
Ultimate gain of final control element
[-]
-1
Mass transfer coefficient of phase x
[m s ]
Pipe length in reactor
[m]
Distance between two pipes in square vessel
[m]
Ziegler Nichols amplitude ratio
[-]
Number of moles (equation 55)
[mole]
-1
Stirring speed
[# s ]
Number of transport units (equation 6)
[-]
Number of transport units coolant (equation 15)
[-]
Nusselt number reaction mixture
[-]
Nusselt number coolant
[-]
Bifurcation parameter (equation 24)
[-]
Potential energy
[J]
-2
Pressure
[N m ] [Pa]
Prandtl number reaction mixture
[-]
Prandtl number coolant
[-]
Ziegler Nichols ultimate period of sustained cycling
[min]
-1
Amount of removed heat
[J s ]
59
Heat generation potential Perry et.al.
[-]
-1
Required heat to bring the feed to operation temperature (equation 65)
[J s ]
-1
Heat produced by reaction (equation 66)
[J s ]
-1
Heat transferred through cooling/heating area (equation 67)
[J s ]
-1
Heat, which is removed by coolant (equation 68)
[J s ]
-1 -1
Ideal gas constant
[J mol K ]
Reynolds number reaction mixture
[-]
Reynolds number coolant
[-]
2
-1
Scale resistance
[m K s J ])
3(m+n-1)
(m+n-1)
Reaction rate
[m
/ mol
s]
2
Cross-section area
[m ]
[m ]
-1
[m m ]
[.]
[.]
-1 -1
[J Kg K ]
-1 -1
[J Kg K ]
[m]
[m]
[m]
[J]
-1
[J mol ]
-1
[s ]
[-]
-2
[m s ]
[-]
[-]
2
-1
106
St
t
T
Tcool
Tdelay
T
Tad
Tlog
td
tU
UA
U
U
uG
uL
Vcool
VR
Vvessel
x
Xi
Stanton number = NTU

Simulation time variable
Temperature
Coolant temperature
Temperature after an apparent delay
Temperature of the fluid around the wall
Adiabatic temperature difference
Logarithmic temperature mean (equation 17)
Time (delay or dead)
Ultimate period
Cooling capacity
Internal energy
Overall heat transfer coefficient
Gas velocity
Average liquid velocity based on cross section
Cooling fluid volume
Reactor volume
Vessel volume
Mathematical variable (equation 32, 24)
Conversion rate
Greek
Description
R
cool
pipe
vessel
cool
wall
R
pipe
cool
T
coolant
cool
R
dead
d
D
I
m
U
Angle
Heat transfer coefficient reactor
Heat transfer coefficient coolant
Thickness of cooling pipe (equation 84)
Thickness of pipe in a cooled vessel (equation 84)
Thickness vessel wall
Hold-up
Error or deviation
Time transformation
Parameters used for ecological-bifurcation model
Liquid viscosity
Coolant viscosity
Viscosity at pipe wall
Dimensionless temperature
Dimensionless time
Bifurcation active parameter (equation 18)
Thermal conductivity
Thermal conductivity reaction mixture
Thermal conductivity cooling coil / pipes
Thermal conductivity coolant
2
59
Temperature sensibility E/RT Perry et.al.
Time constant in a process (equation 32)
Average residence time of a coolant particle in cooling device
Time constant of the cooling fluid (equation 38)
Average resident time (equation 4)
Dead time constant
Presumed dead time
Derivative time constant
Integral time constant / reset time
Measuring time constant
Apparent time constant
Integral control state variable
Ziegler Nichols crossover frequency
Production rate
[-]
[s]
[K]
[K]
[K]
[K]
[K]
[K]
[s]
[s]
-1 -1
[J s K ]
[J]
-2 -1 -1
[J m s K ]
-1
[m s ]
-1
[m s ]
3
[m ]
3
[m ]
3
[m ]
[-]
-1
[s ]
Dimension
[-]
-2 -1 -1
[J m s K ]
-2 -1 -1
[J m s K ]
[m]
[m]
[m]
3
-3
[m m ]
[-]
[-]
[-]
-1 -1
[kg m s ]
-1 -1
[kg m s ]
-1 -1
[kg m s ]
[-]
[-]
[-]
-1 -1 -1
[J m s K ]
-1 -1 -1
[J m s K ]
-1 -1 -1
[J m s K ]
-1 -1 -1
[J m s K ]
-1
[K ]
[s]
[s]
[s]
[s]
[s]
[s]
[s] [min]
[s] [min]
[s]
[s]
[K]
-1
[s ]
-1
[mol s ]
107
v
v,cool
cool
subscript
-1
[m s ]
3 -1
[m s ]
[-]
[-]
-3
[kg m ]
-3
[kg m ]
[-]
[mol]
Abbreviations
Subscripts
0
1, 2,
a, b
ana
act
bc
c
cool
conv
crit
d
e
(g)
i
(l)
lm
log
m
meas
min
max
num
p
pipes
s
sp
stst
t
u
Volumetric flow rate

Flowrate cooling fluid
Place transformation
Ratio height and diameter vessel (equation 75)
Fluid density
Density of cooling fluid
Relative conversion
Absolute conversion
Inlet or fresh feed

Reactor 1, 2,
Component A, B
Analytic solution
Activation
Base case
Controller
Coolant / cooling (media)
Converted
Critical
Dead time
External
Gas phase
Initial / species i / inlet
Liquid phase
Linear median / mean
Logarithmic mean
Measured
Measuring device
Minimum
Maximum
Numerical
Process
Cooling pipes
Specific / steady state
Set point
Steady state
Total
Ultimate
ART
BR
CISTR
CORR
D
G
HPR
HWR
I
IAE
ITAE
L
ODE
PDE
PFR
P
PI
PID
PB
RCSM
RDSM
RHS
RT
Average residence time

Batch reactor
Continuously ideally stirred tank reactor
Corrective
Derivative
Gas
Heat production rate
Heat withdrawal rate
Integral
Integral of the absolute error
Integral time absolute error
Liquid
Differential equations
Partial differential equations
Plug flow reactor
Proportional
Proportional integral
Proportional integral derivative
Proportional band
Reactor designer stability map
Right hand side
Residence time
Symbols
Base case
LOCBIF List of symbols

Printed in Appendix 6 at p123.
108
REFERENCES
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
[18]
[19]
[20]
[21]
[22]
[23]
[24]
[25]
[26]
[27]
[28]
[29]
[30]
Aris, R. and Amundson, N.R., 1958, An analysis of chemical reactor stability and control - part I
to III, Chem. Eng. Sci. 7, p121-155.
Aris, R, 1969, Elementary chemical reactor analyses, Prentice-Hall, Englewood Cliffs, NJ.
Bazykin, A.D., 1985, Mathematical biophysics of interacting populations, Nauka Moscow Russia.
Beek, W.J. and Muttzall, K.M.K.,1983, Transport phenomena, Wiley Intersscience Publication,
ISBN: 0-471-06173-5.
Bush, S.F., 1969, Proc. Royal. Soc., 309A, p1-26.
Bykov, V.I., Yablonskii, G.S. and Kim, V.S., 1978, On the simple model of kinetic self-oscillations
in catalytic reaction of CO oxidation, Dokl. Sov. Math., p637-639.
Bilous, O. and Amundson, N.R., 1955, Chemical reactor stability and sensitivity, A.I.Ch.E. J. 1,
p513-521.
Bykov, V.I., Yablonskii, G.S., Kim, V.F., 1978, On the simple model of kinetic self-oscillations in
catalytic reaction of CO oxidation, Dokl. Sov. Math., p637-639.
Calvet, J.P. and Arkun, Y., 1988, Feedforward and feedback linearisation of non-linear systems
and its implementation using internal model control, Ind. Eng. Chem. Res. 27 p1822-1831.
Calvet, J.P. and Arkun, Y., 1990, Design of P and PI stabilizing controllers for quassi non-linear
systems., Comp. Chem. Eng. 14, p415-426.
Cohen, G.H. and Coon, G.A., 1953, Theoretical consideration of retarded control, Trans.
A.S.M.E., 75, p827.
Coughanowr, D.R., and Koppel, L.B., 1965, Process systems analysis and control, McGraw-Hill,
th
64 ed.
Coulson, J.M., Richardson, J.F., Backhurst, J.R. and Harker, J.H., 1993, Chemical Engineering,
th
Pergamon Press, Volume 1, 4 edition, ISBN: 0-08-037948-06.
Coulson, J.M., and Richardson, J.F., 1979, Chemical Engineering, Pergamon Press, Volume 3,
nd
2 edition, ISBN: 0-08-023818-1.
Coulson, J.M., Richardson, J.F. and Sinnott, R.K., 1991, Chemical Engineering, Pergamon
th
Press, Volume 6, 1 edition, ISBN: 0-08-022969-0.
Daoutidis, P. and Kravaris, C., 1991, Dynamic output feedback control of minimum-phase nonlinear processes, Chem. Eng. Sci., 47-4, p837-849.
Ding, J.S.Y., Sharma, S. and Luss, D., 1973, Steady state multiplicity and control of the
chlorination of liquid n-decane in an adiabatic continuously stirred tank reactor, Ind. Chem. Eng.
13, p76-82.
Doedel, E.J. and Heinemann, R.F., 1983, Numerical computations of periodic solutions branches
and oscillatory dynamics of the stirred tank reactor, Chem. Eng Sci., 38 ,p1493-1499.
Dolnik, M., Banks, A.S. and Epstein, I.R., 1997, Oscillatory chemical reaction in a CISTR with
feedback control of flowrate., J. Phys. Chem., 101, p5148-5154.
Dougles, J., 1972, Process dynamics and control, volume 1, Analyses of dynamic systems,
Prentice-Hall, Englewood Cliffs.
Doherty, M.F. and Ottino, J.M., 1988, Chaos in deterministic systems: strange attractors,
turbulence and applications in chemical reactors , Chem. Eng. Sci., 43-2, p139-183.
Elk, E.P., van, Borman, P.C., Kuipers, J.A.M. and Versteeg, G.F., 1999, Modelling of gas-liquid
reactors Stability and dynamic behaviour of gas-liquid mass transfer accompanied by
irreversible reaction, Chem. Eng. Sci., 54, p4869-4879.
reactors Implementation of the penetration model in dynamic modelling of gas-liquid processes
with the presence of a liquid bulk, Chem. Eng. J., 76, p223-237.
reactors Stability and dynamical behaviour of a hydroformylation reactor, Chem. Eng. Sci., 56,
p1491-1500.
Erickson, K.T., and Hedrick, J.L., 1999, Plantwide process control, John Wiley & Sons, ISBN: 0471-17835-7.
nd
Fogler, H.S., 1992, Elements of chemical reactor engineering, Prentice Hall Int., 2 ed., ISBN 013-263534-8.
Foust, A.L., Wenzel, L.A., Clump, C.W., Maus, L. and Andersen, L.B., 1980, Principles of unit
operations, John Wiley & Sons, New York, ISBN 0-471-26897-6.
Franks, G.E.F., 1966, Mathematical modelling in chemical engineering, John Wiley & Sons, Inc.
Giona, M, and Paladino, O, 1994, Bifurcation analysis and stability of a controlled CSTR,
Computers and chem. eng. J., 18, p877-887.
Hancock, M.D., Kenney, C.N., 1976, The stability and dynamics of a gas-liquid reactor, Chem.
Eng. Sci. 32, p629-636.
109
[31]
[32]
[33]
[34]
[35]
[36]
[37]
[38]
[39]
[40]
[41]
[42]
[43]
[44]
[45]
[46]
[47]
[48]
[49]
[50]
[51]
[52]
[53]
[54]
[55]
[56]
[57]
[58]
[59]
[60]
[61]
[62]
[63]
Harold, M.P. and Luss, D., 1984, An experimental study of steady state multiplicity features of
two parallel catalytic reactors, Chem. Eng. Sci., 40, p35-52.
Harold, M.P., Ostermaier, J.J., Drew, D.W., Lerou, J.J. and Luss, D., 1996, The continouslystirred dacanting reactor: Steady state and dynamic features, Chem. Eng. Sci. 51, p1777-1786.
Heemskerk, A.H., Dammers, W.R., Fortuin, J.M.H., 1980, Limit cycles measured in a liquidphase reaction system, Chem. Eng. Sci. 35, p439-445.
Heerden, van, C., 1953, Ind. Eng. Chem., p1242.
Henson, M.A., and Seborg, D.E., 1997, Nonlinear process control, Prentice Hall PTR., ISBN: 013-625179-X.
Heiszwolf, J.J. and Fortuin, M.H., 1997, Design procedure for stable operations of first-order
reaction systems in a CSTR, A.I.Ch.E. J. 43, p1060-1068.
Himmelblau, D.M., & Bischoff, K.B., 1968, Process analyses and simulation, John Wiley & Sons,
Inc.
Hjelmfelt, A. and Ross, J., 1994, Experimental stabilization of unstable steady states in oscillatory
and excitable reaction systems, J. Phys. Chem. 98, p1176-1179.
Hoffman, L.A., Sharma, S. and Luss, D., 1975, Steady state multiplicity of adiabatic gas-liquid
reactors: I. The single reactor case, A.I.Ch.E. J. 21, p318-326.
Hopf, E., 1942, Abzweigung einer periodischen Losung eines Differentialsystems. Aus den
Berichten der Mathematisch-Physikalischen Klasse de Schsischen Akademie der Wissenhaften
zn Leipzig XCIV, p1-22.
Hork, J., Belohlav, Z., Rosol, P. and Madron, F., 1987, Analysis of the oscillatory behaviour of
an industial reactor for oxination of propene; combined models, Coll. Czechoslovak Chem.
Cummun. 52, 2865-2875.
Huang, D. T.-J. and Varma, A., 1981, Steady-state and dynamic behavior of fast gas-liquid
reactions in non-adiabatic continuous stirred tank reactors, Chem. Eng. J. 21, p47-57.
Huang, D. T.-J. and Varma, A., 1981, Steady-state uniqueness and multiplicity of nonadiabatic
gas-liquid CSTRs, A.I.Ch.E. J. 27, p481-489.
Khibnik, A.I., Bykov, V.I. and Yablonskii, G.S., 1987, 23 phase portraits of the simplest catalytic
oscillator, J. Phys. Khim., p1388-1390.
Khibnik, A.I, Kuznetsov, Y.A., Levitin, V.V. and Nikolaev, E.V., 1992, Interactive LOCal
BIFurcation Analyzer.
Kuntha, A. and van Elk, E.P., 1999, Modelling of gas-liquid reactors: stability and dynamic
behaviour of gas-liquid mass transfer accompanied by irreversible reaction in CISTRs in series,
Procede Twente b.v..
Kravaris, C. and Palanki, S., 1988, Robust non-linear state feedback under structured
uncertainty, AIChE J., 34-7, p1119-1127.
Kuznetsov, Y.A., 1995, Elements of applied bifurcation theory, Springer-Verlag, ISBN: 0-38794418-4.
Levenspiel, O.,1972, Chemical reaction engineering, John Wiley & sons, ISBN: 0-471-53019-0.
Li, R.S. and Horsthemke, W., 1991, Effects of product occupancy on self-organization in
heterogeneous catalysis. I. Multiple steady states and temporal oscillations, J. Chem. Phys. 95,
5785-5789.
Limquerco, L.C. and Kantor, J.C., 1989, Non-linear output feedback control of an exothermic
reaction, Comp. Chem. Eng., 14, p427-437.
Lopez, A.M., Smith, C.L. and Murill, P.W., 1969, An advanced tuning method, Brit. Chem. Eng.,
14, p1533.
Marlin, T.E., 2000, Process control, designing processes and control systems for dynamic
nd
performance, McGraw-Hill Chem. Eng. Series., 2 ed., ISBN 0-07-039362-1.
McAuley, K.B., Macdonald, D.A. and McLellan, P.J., 1995, Effects of operating conditions on
stability of gas-phase polyethylene reactors, A.I.Ch.E. J. 41, 868-879.
th
Mohilla, R. and Ferencz, B., 1982, Chemical process dynamics, Elsevier 4 ed., ISBN 0-44499730-x.
Murrill, P.W., 1967, Automatic control of processes, International Textbook Co.
Narendra, K.S. and Monopoli, R.V., 1980, Applications of adaptive control, Academic Press New
York.
Pellegrini, L. and Biardi, G., 1990, Chaotic behaviour of a controlled CSTR, Computers and
chem. eng. J., 18, p877-887.
th
Perry, R.H., Green, D., 1984, Perrys chemical engineers handbook, McGraww-Hill Int. Ed., 50
ed., ISBN 0-07-Y66482-X.
Ogunnaike, B.A., and Ray, W.H., 1994, Process dynamics modelling and control, Oxford
Unversity Press, ISBN: 0-19-509119-1.
Olsen, R.J. and Epstein, I.R., 1993, Bifurcation analysis of chemical reaction mechanisms. II.
Hopf bifurcation analysis, J. Chem. Phys. 98, p2805-2822.
Pinto, J.C., 1995, The dynamic behavior of continuous solution polymerization reactors - a full
bifurcation analysis of a full scale copolymerization reactor, Chem. Eng. Sci. 21, 3455-3475.
Raghuram, S., Shah, Y.T. and Tierney, J.W., 1979, Multiple steady states in a gas-liquid reactor,
Chem. Eng. J. 17, p63-75.
110
[64]
[65]
[66]
[67]
[68]
[69]
[70]
[71]
[72]
[73]
[74]
[75]
[76]
[77]
[78]
[79]
[80]
[81]
Ratto, M. and Paladino, O., 2000, Controllability of start-ups in CISTRs under PI control,
Dipartimento di Ingegneria Ambientale, Universit di Genova, Italy.
Ratto, M. and Paladino, O., 1999, Robust stability and sensitivity of real controlled CSTRs
Chem. Eng. Sci., 55, p321-330.
Ratto, M., 1998, A theoretical approach to the analyses of PI-controlled CSTRs with noise,
Computers and chem. eng. J., 22, p1581-1593.
Ray, A.K., 1995, Performance improvement of a chemical reactor by non-linear natural
oscillations, Chem. Eng. J., 59, p169-175.
Roffel, B. and Rijnsdorp, J.E., 1982, Process dynamics, control and protection, Ann Arbor
Science, ISBN 0-250-40483-4.
Scott, S.K., 1987, Oscillations in simple models of chemical systems, Acc. Chem. Res., 20,
p186-191.
Sharma , S., Hoffman, L.A. and Luss, D., 1976, Steady state multiplicity of adiabatic gas-liquid
reactors: II. The two consecutive reactions case, A.I.Ch.E. J. 22, p324-331.
Shinskey, F.G., 1994, Feedback controllers for the process industries, McGraw Hill, ISBN 0-07056905-3.
th
Shinskey, F.G., 1996, Process control systems, McGraw Hill, 4 ed., ISBN 0-07-057101-5.
Singh, C.P.P., Carr, N.L. and Shah, Y.T., 1982, The effect of gas feed temperature on the steady
state multiplicity of an adiabatic CSTR with a fast pseudo-first-order reaction. Chem. Eng. J. 23,
p101-104.
Stephanopoulos, G., 1984, Chemical process control: an introduction to theory and practice,
Englewood Cliffs, N.J.: Prentice Hall Inc., ISBN: 0-13-128629-3.
Uppal, A., Ray, W.H. and Poore, A.B., 1974, On the dynamic behavior of continuous stirred tank
reactors, Chem. Eng. Sci. 29, p967-985.
Uppal, A., Ray, W.H. and Poore, A.B., 1976, The classification of the dynamic behavior of
continuous stirred tank reactors - Influence of reactor residence time, Chem. Eng. Sci. 31, p205214.
Vleeschhouwer, P.H.M. and Fortuin, J.M.H., 1990, Theory and experiments concerning the
stabilty of a reacting system in a CSTR. A.I.Ch.E. J. 36, p961-965.
Vleeschhouwer, P.H.M., Garton, R.D. and Fortuin, J.M.H., 1992, Analysis of limit cycles in an
industrial oxo reactor. Chem. Eng. Sci. 47, p2547-2552.
th
Weast, R.C., 1984, Handbook of Chemistry and Physics, CRC Press, 64 ed.
Westerterp, K.R., van Swaaij, W.P.M. and Beenackers, A.A.C.M., 1990, Chemical Reactor
Design, John Wiley & Sons, ISBN: 0-471-90138-0.
Ziegler, J.G. and Nichols, N.B., 1942, Optimum settings for automatic controllers, Trans.
A.S.M.E., 64, p759.
Table 37 References list and categories.

Subject
Limit cycles
Dynamic instability
Bifurcation
Process control
CISTR
Performance
Dead time / delay
Dimensionless notation
Experiments
General mathematics
References
18 19 33 35 38 48 51 54 61 62 67 69 71 75 76 77 78
1 7 16 17 18 21 30 31 32 33 34 35 36 37 38 40 42 47 48 51 54 60 61 62 64 67 68 71 72 75 76 77 78
1 7 18 19 21 32 36 38 43 48 50 54 61 62 64 67 75 76 77 78
1 10 11 12 14 16 17 19 25 28 29 35 37 47 51 52 53 56 55 57 58 60 64 65 68 71 72 74
1 4 10 14 17 18 19 21 26 30 31 32 33 34 35 38 41 42 43 47 51 54 55 59 60 61 62 63 64 67 68 70 73 75 76 80
14 60 61 67 50
14 25 37 53 60 71 72 74
1 10 32 35 36 42 43 51 60 62 70 73 75 76 77
5 14 19 27 30 31 33 38 54 62 69 77 78
3 28 45 48 59 79
Bifurcation theory
68 25 37 74 60
Linearisation
25 60 53 68 74
Laplace transforms
7 14 60 74
Routh criterium
55 68
Partial differential equations
68
Discretisation and continuisation
111
APPENDICES
Appendix 1. Derivation CISTR
The continuously ideally stirred tank reactor
A simple exothermic reaction 1 takes place in the Continuously Ideally Stirred Tank Reactor which has
3
-1
been displayed in Figure 1 at p2. A liquid enters the reactor with a flow rate of V [m s ] and a
temperature T0 [K]. This feed flow contains component A with concentration [A]0. The tank is considered to
be perfectly mixed, which implies that the temperature and concentration of the effluent is equal to the
temperature and concentration of the liquid in the tank V, T, [A] and [P]. The reactor is cooled by a
coolant that for example flows through a jacket around the reactor or flows through a system of cooling
pipes (see also Appendix 3).
The fundamental dependent quantities for a reactor are:
1
2
3
Total mass of the reacting mixture in tank

Mass of chemical substance A in the reacting mixture
Total energy of the reacting mixture in the tank
Remarks:
1
2
The mass of component P can be found from the total mass and the mass of component A.
Therefore, it is not an independent fundamental quantity and therefore is superfluous.
The momentum of the CISTR does not change under any operating conditions for the reactor
and will be neglected.
80
74
68
53
Assumptions (Westerterp et.al. , Stephanopoulos , Roffel and Marlin ):

1.
2.
3.
4.
5.
6.
7.
8.
9.
The liquid phase is ideally mixed and can be described as a CISTR.

The irreversible, elementary first order chemical reaction states A P.
The potential and kinetic energies of the inlet and outlet streams are equal.
The liquid in and out flowrate are assumed equal.
The physical parameters are considered constant i.e. = 0 and CP = 0.
The CISTR is well insulated so that negligible heat is transferred to the surroundings
The accumulation of energy in the reactor walls, agitator and cooling coil is negligible
80
compared with the accumulation in the liquid. According to Westerterp et.al. this is
permissible for large industrial reactor.
60
According to Ogunnaike et.al. The coolant is at quasi-steady state.
The shaft work done by the impeller of the stirring mechanism, has been neglected.
The conservation principles on the three fundamental quantities will be applied:
Total mass balance
accumulation
of total mass
=
time
input
output
total mass
of total mass of total mass generated or consumed
time
time
time
or
Appendices
112
d VR
= (V )0 (V ) = 0
dt
(54)
where: 0, = densities of the inlet and outlet streams.

V,0, V, = volumetric flow rates of the inlet and outlet streams.
VR = volume of the reacting mixture in the tank.
Mass balance on component A
accumulation
of A
time
input
of A
time
output
of A
time
disappearance of A
due to reaction
time
or
d (n A ) d (VR [ A])
=
= ([ A]V )0 ([ A]V ) rVR
dt
dt
(55)
where: r = rate of reaction per unit volume.

[A], [A]0 molar concentrations of vA in the inlet and outlet streams
nA = number of moles of A in the reacting mixture.
In expression 55 r represents the overall chemical reaction rate (reaction 1) defined as:
r = rA = k R [ A] = k 0 e
Eact
RT
(1)
[ A]
In equation 1, Eact is the activation energy and k0 the Arrhenius frequency factor. Expression 1 was
59 14 26
considered by Arrhenius
. It represents the effect of temperature on the rate constants (of simple
reactions) i.e. it symbolises the amount of energy in excess of the average energy level, which the
reactants must have in order for the reaction to proceed.
The specific state variables are [A] and VR. Algebraic manipulation on equation 55 leads to:
Eact
d (VR [ A])
dV
d[ A]
= [ A] R + VR
= (V [ A])0 (V [ A]) k0 e RT [ A]VR
dt
dt
dt
(56)
this expression becomes:

Eact
d[ A] V
([ A]0 [ A]) k0e RT [ A]
=
dt
VR
(57)
Total energy balance
accumulation
of total energy
=
time
input of total output of total energy removed

energy of feed energy of feed by coolant
time
time
time
74
Extensive thermodynamic derivations can be found in Stephanopoulos .
Appendices
113
dH
= 0 V ,0 h0 (T0 ) V h(T ) Q
dt
(58)
Where Q is the amount of heat removed. The amount of heat supplied by e.g. steam or removed by the
80
coolant per unit time is given by following expressions:
Qheat = UA (Tsteam T )
Qcool = UA (T Tcool )
(59)
(60)
In equation 60 U is the overall heat transfer coefficient and A represents the total area of heat transfer.
Equations 54, 55 and 58 are not in their final and most convenient form for process control design studies.
To bring them to such form the appropriate state variables will be identified.
From the thermodynamics, it is known that the enthalpy of a liquid system is a function of the temperature
and its composition:
VR C P
dT
~
~
= V ,0 0 C P,0 (T0 T ) + H A H P rVR Q
dt
(61)
CP is the specific heat capacity of the reacting mixture and in equation 61 H A and H P are the partial
molar enthalpies of A and P. Since H A
~
H P = H R
which represents the heat of reaction at
temperature T and = 0, CP = CP,0.
VR
H R
Q
dT
= V ,0 (T0 T ) +
rVR
C P
C P
dt
(62)
From equation 62 the temperature T is the state variable that characterise the total energy of the system.
Summarizing all the steps above in the mathematical modelling of a CISTR results in state variables [A]
and T and state equations 63 and 64 according
Eact
d[ A] V
=
([ A]0 [ A]) k0 e RT [ A]
dt
VR
Eact
UA (T Tcool )
dT V
(T0 T ) + H R k 0 e RT [ A]
=
C P
C PVR
dt VR
(63)
(64)
Heat classification
Heat can be categorised in respectively:
heat, which is needed to bring the feed to operation temperature
Q feed = C P V (T T0 )
(65)
heat produced by reaction
Qreaction = H RVR k 0 e
Eact
RT
[ A]
(66)
transferred heat
Qtransfer = UA (T Tcool )
(67)
heat, which is removed by coolant according:
Qcool = cool C P,cool V ,cool (Tcool Tcool,0 )
(68)
heat production rate:
Appendices
114
HPR = Tad
k R R
1 + k R R
(69)
and heat withdrawal rate (derived from equations 2, 3, 10, 4, 6 and 7)
HWR = (T0 T ) NTU (T Tcool )
(70)
Linearisation of the non-isothermal CISTR

In process control applications, generally differential equations are linearised to examine its dynamic
25 37 59 60 72 74
behaviour. The literature
deals often with the linearisation of the non-linear CISTR model
74
according the state equations 2 and 3. According to Stephanopoulos the linearised mass and energy
balances could be perceived as a capacity in which energy can be stored. The response of [A] and T to
inlet changes is second order. A CISTR is distinguished by the kinetic rate term, which denotes
disappearance of component A. Such terms may produce not only overdamped but also underdamped
and inverse response.
Appendix 2. Literature
Although many articles can be found regarding CISTRs, stability and process control, little information can
be found with respect to limit cycles and process control. The information that can be gathered is
concerned with proportional-integral control in which the coolant temperature is adjusted. Not taken into
account is, instead of the coolant temperature, the coolant flowrate nor the derivative control action,
assumed delay and the non-ideal process behaviour, which can seriously affect the process stability.
Besides, the literature often focuses on complex mathematical models, which requires very profound
background knowledge. Nevertheless, some articles are refereed which may contribute to understand the
limit cycle and stability problem in CISTRs.
Multiplicity and stability in CISTRs

68
Roffel examined state equations 2 and 3 applying an information flow diagram in which the stable
equilibrium points can be derived from specific trajectories. This method however is rather laborious in
7
case model improvements are concerned. Bilous et.al. presented the first modern stability analyses of the
equilibrium states in a CISTR with a single exothermic reaction. It has been shown that instabilities may
exist in the three steady states. Analytical criteria have been for the determination of stability. The article
provides much information about the dynamical behaviour of a CISTR. However, in case the model is
21
extended, the methods become very complex. Doherty and Ottino explained in their article through
chaos theory the sensitive dependence of solution of differential equations on initial conditions. Hoffman
39
et.al. developed a model to determine steady state multiplicity in a gas-liquid CISTR with exothermic
58
reaction. Pellegrini et.al. determined the region of asymptotic stability for, and the chaotic behaviour of a
17
CISTR. Ding et.al. demonstrated through chlorination of liquid n-decane experimentally in an adiabatic
continuously stirred tank reactor the existence of steady state multiplicity. They investigated the influences
31
of a control scheme. Harold et.al. investigated experimentally the steady state multiplicity features of two
32
parallel catalytic reactors. They found methods to determine all possible steady states. Ostermaier et.al.
43
compared the CISTR with the decanting reactor for its dynamical behaviour. Huang et.al.
studied the
78
second order reaction with regards to multiplicity. Vleeschhouwer et.al. proposed an analytical
perturbation analysis of a cobalt catalysed oxo reactor. They described the process with a pseudo
homogeneous pseudo first order system, which can be described by merely two differential equations.
36
Heiszwolf and Fortuin give a very clear example of a design procedure for stable operation of a first67
order system in a CISTR using an analytical perturbation method. Ray discussed the performance
50
improvement of a chemical reaction by non-linear natural oscillations. Li and Horsthemke studied the
product occupancy in heterogeneous catalyses in relation with multiple steady states and temporal
oscillations.
Appendices
115
Limit cycles theoretical results

1
34
75
33
Aris and Amundson , van Heerden , Uppal et.al. , Heemskerk and Dammers , published the theoretical
treatment of instabilities in continuously ideally stirred tank reactors based on the fundamental mass and
energy balance. In these articles the appearance of self-sustained oscillations i.e. limit cycles have been
69
7
theoretically demonstrated. Scott and Bilous et.al. considered the oscillations in simple models in
18
chemical reactors and represented some experimental cases. Doedel and Heinemann used the
bifurcation theory to investigate the oscillatory behaviour in a CISTR with consecutive ABC reaction.
They used a continuation technique to compute response diagrams exhibiting stable and unstable
80
periodic branches that contain multiple limit points. Finally, Westerterp et.al. mentioned the work on
laboratory scale in which the existence of limit cycles have been demonstrated.
Limit cycles experimental results

33
Heemskerk et.al. described experimental investigations concerning limit cycle behaviour in a one-phase
reaction system i.e. the acid-catalysed hydrolyses of 2,3-epoxypropanol-1b and showed a good
30
agreement with theory based on the fundamental mass and energy balance. Hanckock et.al. analysed
the instabilities in a CISTR with temperature fluctuations in one-phase and two-phase reactors reactor
theoretically and carried out experiments e.g. the formation of methyl-chloride from methanol and
5
hydrogen-chloride. They found the principle features of oscillatory behaviour. Bush discovered the
32
consternation of practical engineering in commercial situations. Harold et.al. constructed maps with
parameter regions with qualitatively different steady state and dynamical bifurcation diagrams. These
maps could clearly describe the desirable regions of operation and point out the potential stability. They
suggested to maintain a sufficiently small difference between the reactor temperature and the coolant
42
temperature to avoid oscillatory behaviour. Huang et.al. did research on predicting the steady state and
dynamical behaviour of fast gas-liquid reactions in a non-adiabatic CISTR. The models they developed
77
were experimentally identified i.e. multiplicity in the chlorination of n-decane. Vleeschhouwer and Fortuin
have been concerned with the theoretical and experimental aspects with respect to stability in a CISTR.
76
Uppal and Ray investigated theoretically the influence of the average residence time on the dynamical
38
behaviour of a CISTR. Hjelmfelt and Ross used a proportional flow feedback method to stabilise
78
unstable stationary states in experiments with chlorite-iodide reactions. Vleeschhouwer et al.
demonstrated that sustained temperature oscillations can occur in a commercial scale gas-liquid oxoreactor.
Process control
47
Kravaris et.al. proposed a robust non-linear feedback control scheme. This article is an improvement of
the global input/output linearisation approach. The article is difficult to expand with differential equations.
64
Ratto et.al. analysed the PI-controlled CISTR with fluctuating and uncertain parameters. In the article,
66
they applied the Hopf bifurcation plot to determine the stability. Ratto made also a theoretical approach
10
to the analyse of PI-controlled CISTRs with noice. Calvet et.al. presented a method to compute gains to
design a PI-controller based on perturbations in non-linear systems. In the article, the influence of limit
16
cycles was not considered. Daoutidis et.al. are concerned with a different approach of non-linear
systems i.e. the state-space approach. Although the author proclaimed stability, the applied methods are
29
rather laborious. Giona and Paladino developed a comprehensive bifurcation analyses and stability
consideration of a controlled CISTR in the presence of a th order exothermic reaction. They studied the
65
PI-control and introduced controller time delay. Paladino and Ratto developed a procedure for stability,
64
robustness and sensitivity of real controlled CISTRs through bifurcation analyses. Paladino and Ratto
also studied the controllability of start-ups in CISTRs under PI-control.
Limit cycle control

51
Limquerco et.al. proposed a non-linear output feedback control of an exothermic reaction in a CISTR. In
case the system was perturbed, they found both ignition and extinction, however more interesting they
measured limit cycles under certain conditions. The disadvantage of this article like the article of Uppal
75
et.al. is the fact that instantaneous change of the coolant temperature is assumed. Interesting point of
Appendices
116
19
the article is the non-linear approach. Dolnik et.al. studied numerical and experimental the chlorine
dioxide-iodide reaction in a CISTR with feedback regulation of flow rate. The difference with this report is
the existence of more than one reactant.
Appendix 3. Value estimation

Despite the fact that this report has a qualitative analytical character and the chosen base case, definition
22
has been extracted from the article by van Elk et.al. it is indispensable to investigate the magnitudes of
the base case parameters; consequently, the chosen values will be compared with regular values and
validated. To examine parameters: UA, -HR, E, k0, , VR, CP.
Dpipe
Hvessel
Dvessel
Figure 166 Theoretical cylindrical vessel containing cooled through n pipes.
UA value estimation
In this section, the value of product U and A will be considered.
-1
-1
In the base case definition, the value UA = 55 [kJ s K ] has been chosen. To validate this magnitude,
inevitably several calculations will be made.
Assume a chemical reactor has to be cooled through n pipes. Consider a cylindrical vessel with volume
Vvessel. The vessel contains a reacting substance with volume VR. Heat is transferred by the pipes and the
-1 -1
reactor wall. The base case UA-value states UA = 55 [kJ s K ].
The reactor volume is the sum of the liquid volume and the volume of the cooling pipes:
Vvessel = V pipes + VR
(71)
In the first place, the number of pipes n is estimated.

The total area regarding heat transfer states:
Apipes = nDpipe H vessel
2
2
+ Dvessel H vessel 2n D pipe
Avessel = 2 Dvessel
4
4
Atotal = A pipes + Avessel
(72)
(73)
(74)
The ratio vessel height Hvessel and vessel diameter Dvessel is introduced:
Appendices
117
H vessel = Dvessel
(75)
The volume of the vessel states:
Vvessel =
2
Dvessel H vessel
4
(76)
The number of pipes in the vessel can be estimated combining equations 72, 73, 74, 75 and 76 under the
assumption that Vvessel VR.
2
Atotal 4VR 3 1
( + )
n=
1
4VR 3
2
Dpipe 12 D pipe

(77)
Subsequently the vessel height and diameter is derived.
The volume of n pipes states:
2
V pipes = n D pipe
H vessel
4
(78)
Combining expressions 71, 75, 76 and 78 provides the following cubic equation:
3
+ a1 D pipe + a2 = 0
Dvessel
(79)
Where:
2
a1 = nDpipe
a2 =
(80)
4VR
The vessel diameter can be calculated with:

1
3
2a1
3
2
Dvessel = 16 108a2 + 12 12a1 + 81a2 +
108a + 12 12a 3 + 81a 2

2
1
2
(81)
Due to the assumption Vvessel VR in equation 77, a slightly divergent value is obtained, using equation 73.
However, merely a few recalculations is sufficient enough to eliminate this discrepancy. Better is to apply
a solver e.g. Excel or Maple.
Using equation 82 the volumetric percentage of all the pipes in the reactor can be acquired:
2
D pipe
100%
Vol% = n
D
vessel
(82)
The specific area of the pipes becomes:
Appendices
118
D pipe
a pipes = 4n
(83)
2
Dvessel
Suppose the cooling fluid flows through pipes with thickness pipe and with an average fluid velocity uL,
V ,cool = n
2
2
(
)uL
D pipe pipe
4
(84)
14 15
summarised typical values for the overall heat transfer coefficient U for various types of
Coulson et.al.
80
heat exchanger and specific media. According to Westerterp et.al. , a general applicable value for the
-2 -2 -1
overall heat transfer coefficient for a cold liquid is U 900 [J m s K ].
The overall heat transfer coefficient value can additionally be verified using several calculation methods
59
36
14
68
according to Perry et.al. , Heiszwolf et.al. , Coulson et.al. and Roffel .
Table 38 represents the essential calculations and results for the UA value estimation.
Table 38 Example UA calculation value in cylindrical vessel.
Description
Calculation
Reactor with volume
Desired UA value
Thickness of 1 pipe
Thickness of pipe wall
Velocity cooling fluid
Reynolds number reaction mixture
Reynolds number cooling fluid
Prandtl number reaction mixture
Prandtl number cooling fluid
Nusselt number reaction mixture
Nusselt number cooling fluid
Heat transfer coefficient reactor
Heat transfer coefficient coolant
Overall Heat transfer coefficient
Ratio vessel height and diameter
Total area for heat transfer
Vessel diameter
Vessel height
Number of pipes in vessel
Vessel volume
Volume of n pipes
Mass of n pipes
Pipe volume percentage
Specific area pipes
Total area n pipes
Total area vessel
UA contribution n pipes
UA contribution vessel
Total pipe length
Residence time coolant
Appendices
Base case value[22] ( = 0)

Base case value22
Assumed value
10% pipe thickness
Based on equation 84
Default equation14 36 59 68
Default equation14 36 59 68
Specific correlation14 36 59 68
Specific correlation14 15 36 59 68 80
Assumed value (arbitrary chosen)
UA / U
Equation 75
Equation 76 (first approximation Vvessel VR = 5 [m3])
Equation 77
Equation 71
Equation 78
m = Vpipes
Equation 82
Equation 83
Equation 72
Equation 73
U Apipes
U Avessel
n Hvessel
Equation 84
V cool / V,cool
Result
VR = 5 [m3]
UA = 55 103 [J s-1 K-1]
Dpipe = 1 = 0.0254 [m]
pipe = 0.00254 [m]
uL = 1.7 [m s-1]
ReR = 2 105 [-]
Recool = 8 104 [-] turbulent
PrR = 4 [-]
Prcool = 3.4 [-]
NuR = 106 [-]
Nucool = 300 [-]
-2 -1
-1
R = 1740 [J m s K ]
-2 -1
-1
cool = 5920 [J m s K ]
U = 900 [J m-2 s-1 K-1]
= 1 [-]
Atotal = 61.1 [m2]
Dvessel = 1.75 [m]
Hvessel = 2.19 [m]
n = 255 [-]
Vvessel = 5.28 [m3]
Vpipes = 0.23 [m3]
m = 1770 [kg]
Vol% = 5.3 [%]
8.4 [m2 m-3]
Apipes = 44.5 [m2]
Apipes = 16.6 [m2]
UA = 40 103 [J s-1 K-1]
UA = 15 103 [J s-1 K-1]
Lpipes,total = 558 [m]
3
-1
V,cool = 13 [m min ]
1.3 [s]
119
Number of pipes
The first calculation gives n = 263 for Vvessel = 5.29 [m ]. The solver results in n = 255 for Vvessel = 5.28 [m ].
Therefore, the number of pipes has no decisive effect on the vessel volume. The vessel requires
approximately 250 pipes to cool the exothermic reaction. Consequently, 5% of the vessel is filled with
-1 -1
pipes. Concisely, the base case value UA = 55 [kJ s K ] is allowable.
300
200
100
0
10
30
50
UA [kJ/s/K]
Figure 167 The number of pipes in the

vessel (Table 38) rises if the
required UA value increases. No
cooling pipes are needed if
required UA < 15 [kJ s-1 K-1] due
to the heat transfer through the
vessel wall.
-
HR value estimation
80
According to Westerterp et.al. HR is defined as the heat absorbed by the system when the reaction
proceeds completely in the direction indicated by the arrow, at constant temperature and pressure.
79
In Weast et.al. tables with values of the heat formation at H (T = 25 [C]) for various species have been
5
-1
given. On account of these tables the reaction enthalpy HR = -1.6 10 [J mol ] given in the base case
22
definition van Elk et.al. , can be deduced for instance organic-water mixtures.
Eact value estimation

80
According to Westerterp et.al. typical values of Arrhenius activation energy vary globally between 4 10
5
-1
4
-1
Eact < 3 10 [J mol ]. The value Eact = 9 10 [J mol ] which has been chosen in the base case
definition is a common quantity for organic or organic-water mixtures compounds and therefore
acceptable.
k0 value estimation
22
In the article by van Elk et.al. the value of the Arrhenius pre-exponential factor or frequency factor
5
(equation 1) has been chosen based on a gas-liquid system in a kinetic controlled regime: k0 = 5 10 [s
1
]. According to the base case definition, no gas is present. Therefore the A-concentration [A]L,bulk = 50.1
-3
5
7
-1
[mol m ] will be integrated in the base case k0 value i.e.: k0 = 50.1 5 10 = 2.505 10 [s ] to obtain
the same results published in the articles published by van Elk et.al.
value estimation
-3
14
In the base case definition, the density states = 800 [Kg m ]. According to Coulson et.al. , Fogler
26
68
et.al. and Roffel this value is a common quantity for regular organic compound or organic-water
14
-3
mixtures. The density for stainless steel is assumed to be = 7830 [Kg m ].
Appendices
120
CP value estimation
-1
-1
[59
In the base case definition, a heat capacity CP = 2000 [J Kg K ] has been chosen. According to Perry
2
3
regular value for liquids are on the order of 5 10 CP < 2 10 , Therefore the chosen CP value is
59
-1 -1
acceptable for the base case. The heat capacity for steel is assumed to be CP = 500 [J Kg K ].
VR value estimation
22
The reactor volume discussed in van Elk et.al. , has been set to VR = 10 [m ]. Additionally 50% of the
reactor is filled with gas i.e. the hold-up = 0.5. According to the base case definition, no gas is present in
3
the reactor i.e. = 0, consequently the current reactor volume becomes VR = 10 (1 - 0.5) = 5 [m ].
Appendix 4. Shutting down the input flowrate

In case the stability of the process cannot be maintained, due to e.g. control failure, the possibility of
shutdown the feed can than be considered. If has been decided to close the feed flow i.e. not the drain,
the state equations 2 and 3 of the CISTR change into respectively equations 85 and 86:
Eact
d[ A]
= V [ A] VR k0 e RT [ A]
dt
Eact
dT
CPVR
= CP V T + ( H R )VR k0 e RT [ A] UA(T Tcool )
dt
VR
(85)
(86)
750
Maximum overshoot
650
550
Steady state
450
350
0
300
600
900
Tim e [s]
Reactor tem perature [K]
Reactor tem perature [K]
The reaction will drain (tR) and the reaction will expire (Figure 168). Loss of production is not the only
disadvantage of shutdown the reactor. The shutdown operation itself will waste energy and material stored
in the process, which must be removed, and the subsequent start-up will require a similar amount of
energy and material to be added to reach operating levels again. However, in case the temperature
71
overshoot. More about this subject can be found in Shinskey .
1
0.8
0.6
0.4
0.2
0
0
300
600
900
Tim e [s]
Figure 168 The reactor temperature drops to Figure 169 The conversion (equation 86)
T=385 [K] in case the input
flowrate is halted (equation 86)
and is mainly stabilised after
t 500 [s].
Appendices
rises until the maximum value

has been reached, which is
obviously (t = R) 1.
121
Appendix 5. Cooling capacity coolant flowrate relationship

27
60
Foust et.al. and Ogunnaike et.al. proposed an empirical relationship 87 in which the cooling capacity
can be substituted by the cooling flowrate.
UA = a1 (V ,cool ) 2
a
(87)
Table 39 Values empirical heat capacity relation.

Parameter Value
Dimension
[J s-1 K-1]
[-]
1.2 105
0.5
A1
A2
Combination of equations 13 and 87 can be used to determine the recycle stream ratio. The base case
-1 -1
cooling capacity UA = 55 [kJ s K ] corresponds with a recycle ration of approximately 10.
Foust
100
80
60
40
20
0
0.0
0.1
0.2
0.3
0.4
0.5
Coolant flow rate [m3/s]
Figure 170 Empirical relationship 87 between the cooling capacity and the coolant flowrate.
2
More about this subject Aris et.al. .
Appendices
122
Appendix 6. LOCBIF Nomenclature

LOCBIF
Symbol
Description
A
cA
cA0
Cp
Cpe
Delay
Dpipe
Eact
fouling
g
dH
k
Kc
kr
Lpipe
NTU
NTUe
A
[A]
[A]0
CP
CP,cool
d
Dpipe
Eact
f
g
HR
k0
Kc
kR,m,n
Lpipe
NTU
NTUcool
Heat transfer or exchanger area

Concentration component A
Concentration component A
Heat capacity
Heat capacity of cooling fluid
Apparent delay
Pipe diameter
Activation energy (equation 1)
Fouling factor, rate of UA depletion
Gravitational acceleration
Reaction heat / enthalpy of reaction based on R
Arrhenius pre-exponential factor or frequency factor
Proportional control gain
Reaction rate constant
Pipe length in reactor
Number of transport units
Number of transport units coolant
phi
(T
Qfeed
Qreaction
Qtransfer
Qcool
R
Spiral
t
T
Te
Tdelay
thickness
dTad
td
UA
ve
Vcool
Vr
Dpipe
e
taue
taue
tau
tauI
F
Fe
rho
rhoe
Qfeed
Qreaction
Qtransfer
Qcool
R
mpipesCp,pipes
t
T
Tcool
Tdelay
Dpipe
Tad
td
UA
uL
Vcool
VR
coolant
cool
R
I
v
v,cool
cool
sp
T )dt
Appendices
Dimension
2
[m ]
-3
[mol m ]
-3
[mol m ]
-1 -1
[J Kg K ]
-1 -1
[J Kg K ]
[s]
[m]
-1
[J mol ]
-1
[s ]
-2
[m s ]
-1
[J mol ]
-1
[s ]
3 -1 -1
[-] or [m s K ]
3(m+n-1)
(m+n-1)
[m
/ mol
s]
[m]
[-]
[-]
[K s]
-1
Required heat to bring the feed to operation temperature

[J s ]
-1
Heat produced by reaction
[J s ]
-1
Heat transferred through cooling/heating area
[J s ]
-1
Heat, which is removed by coolant
[J s ]
-1 -1
Ideal gas constant
[J mol K ]
-1
Additional (cooling device) term in heat balance coolant flowrate
[J K ]
Simulation time variable
[s]
Temperature
[K]
Coolant temperature
[K]
Temperature after an apparent delay
[K]
Pipe thickness
[m]
Adiabatic temperature difference
[K]
Time (delay or dead)
[s]
-1 -1
Cooling capacity
[J s K ]
-1
Average liquid velocity based on cross section
[m s ]
3
Volume cooling fluid
[m ]
3
Reactor volume
[m ]
Thickness of cooling pipe
[m]
3
-3
Hold-up
[m m ]
Average residence time of a coolant particle in cooling device
[s]
Time constant of the cooling fluid
[s]
Resident time
[s]
Integral time constant / reset time
[s] [min]
3 -1
Volumetric flow rate
[m s ]
3 -1
[m s ]
-3
Fluid density
[kg m ]
-3
Density of cooling fluid
[kg m ]
123
Appendix 7. LOCBIF Source codes

PHASE zeta, T
PAR UA,Tcool,F
COMMON Vr,rho,cA0,Cp,dH,k,E,R,T0,ret
FUN tau,NTU,dTad
tau=Vr/F
NTU=UA/(rho*Cp*F)
dTad=-dH/(rho*Cp)*cA0
zeta'=(1-zeta)*k*exp(-E/R/T)-zeta/tau
T'=((T0-T)-NTU*(T-Tcool))/tau+dTad*(1-zeta)*k*exp(-E/R/T)
INIT={F=0.005 Vr=5 rho=800 Cp=2000 cA0=5000 dH=-1.6E5 k=2.505E7 E=9E4 R=8.31441 T0=303} ret
LOCBIF rhs 2 Bifurcation sample base case.

PHASE zeta, T
PAR UA,Tset,Teset,Kc
COMMON F,Vr,rho,cA0,Cp,dH,k,E,R,T0,ret
FUN tau,NTU,dTad,Te,offset
tau=Vr/F
NTU=UA/(rho*Cp*F)
offset=Tset-T
Te=Teset+Kc*offset
T'=((T0-T)-NTU*(T-Te))/tau+dTad*(1-zeta)*k*exp(-E/R/T)
LOCBIF rhs 3 Base case + coolant temperature proportional control.

PHASE zeta, T, phi
PAR UA,Tset,Teset,Kc,taui
FUN tau,NTU,dTad,Te
tau=Vr/F
NTU=UA/(rho*Cp*F)
Te=Teset+Kc*((Tset-T)+phi/taui)
phi'=Tset-T
LOCBIF rhs 4 Base case + coolant temperature proportional-integral control.

PHASE cA, T, Te
PAR Tset,UA,Kc,Te0,Teset,recycle
COMMON F,Vr,rho,cA0,Cp,dH,k,E,R,T0,Vcool,rhoe,Cpe,Dpipe,thickness,n,pi,spiraal,ret
FUN Qfeed,Qreaction,Qtransfer,Qcool,zeta,Feset,offset,Fe,ve,NTUe
Qfeed=rho*Cp*F*(T0-T)
Qreaction=(-dH)*k*exp(-E/R/T)*cA*Vr
Qtransfer=UA*(T-Te)
Qcool=rhoe*Cpe*Fe*(Te0-Te)
zeta=1-cA/cA0
Feset=-UA*(Tset-Teset)/(rhoe*cpe*(Te0-Teset))
offset=Tset-T
Fe=(sgn(Feset-Kc*offset)+1)*(Feset-Kc*offset)/2
ve=(1+recycle)*(Fe/n/(pi/4*(Dpipe^2-thickness^2)))
NTUe=UA/(rhoe*Cpe*Fe)
cA'=F/Vr*(cA0-cA)-k*exp(-E/R/T)*cA
T'=(Qfeed+Qreaction-Qtransfer)/(rho*Cp*Vr)
Te'=Fe/Vcool*(Te0-Te)+UA/(spiraal+rhoe*Cpe*Vcool)*(T-Te)
INIT={F=0.005 Vr=5 rho=800 Cp=2000 cA0=5000 dH=-1.6E5 k=2.505E7 E=9E4 R=8.31441 T0=303
Vcool=0.283 rhoe=1E3 Cpe=4200 spiraal=0 Dpipe=0.0254 thickness=0.00254 n=255
i 3 1415926536}
t
LOCBIF rhs 5 Base case + coolant flowrate proportional control.

PHASE cA, T, Te, phi
PAR Tset,UA,Kc,Te0,Teset,recycle,taui
COMMON F,Vr,rho,cA0,Cp,dH,k,E,R,T0,Vcool,rhoe,Cpe,Dpipe,thickness,n,pi,spiraal,Kp,ret
FUN Feset,Fe
Fe=Feset-Kc*((Tset-T)+phi/taui)
T'=((rho*Cp*F*(T0-T))+((-dH)*k*exp(-E/R/T)*cA*Vr)-(UA*(T-Te)))/(rho*Cp*Vr)
Te'=((UA*(T-Te))+(rhoe*Cpe*Fe*(Te0-Te)))/(spiraal+rhoe*Cpe*Vcool)
phi'=(Tset-T)
pi=3 1415926536} ret
LOCBIF rhs 6 Base case + coolant flowrate proportional-integral control.
Appendices
124
PHASE cA, T
PAR UA,Te,Tset,Fset,Kc
FUN tau,Qfeed,Qreaction,Qtransfer,zeta,offset,F,deviation
tau=Vr/F
Qtransfer=UA*(T-Te)
zeta=1-cA/cA0
offset=Tset-T
F=(sgn(Fset-Kc*(Tset-T))+1)*(Fset-Kc*(Tset-T))/2
deviation=(Fset-F)/Fset*100
INIT={Vr=5 rho=800 Cp=2000 cA0=5000 dH=-1.6E5 k=2.505E7 E=9E4 R=8.31441 T0=303} ret
LOCBIF rhs 7 Base case + flowrate proportional control.

PHASE cA, T, phi
PAR UA,Te,Tset,Fset,Kc,taui
tau=Vr/F
Qtransfer=UA*(T-Te)
zeta=1-cA/cA0
offset=Tset-T
F=(sgn(Fset-Kc*((Tset-T)+phi/taui))+1)*(Fset-Kc*((Tset-T)+phi/taui))/2
phi'=offset
LOCBIF rhs 8 Base case + flowrate proportional-integral control.

PHASE zeta, T, UA
PAR Tset,Teset,Kc
FUN tau,NTU,dTad,Te,fouling
tau=Vr/F
NTU=UA/(rho*Cp*F)
Te=Teset+Kc*(Tset-T)
fouling=-UA*(1-0.99)/3600
UA'=fouling
LOCBIF rhs 9 Fouling effect at a coolant temperature proportional controlled base.

PHASE zeta, T, phi, UA
PAR Tset,Teset,Kc,taui
FUN tau,NTU,dTad,Te,fouling
tau=Vr/F
NTU=UA/(rho*Cp*F)
Te=Teset+Kc*((Tset-T)+phi/taui)
fouling=-UA*(1-0.99)/3600
phi'=Tset-T
UA'=fouling
LOCBIF rhs 10 Fouling effect at a coolant temperature proportional-integral controlled base.
Appendices
125
PHASE cA, T, Te, UA

PAR Tset,Kc,Te0,Teset,recycle
FUN Qfeed,Qreaction,Qtransfer,Qcool,zeta,Feset,offset,Fe,ve,fouling
Qtransfer=UA*(T-Te)
zeta=1-cA/cA0
offset=Tset-T
fouling=-UA*(1-0.99)/3600
UA'=fouling
pi=3.1415926536} ret
LOCBIF rhs 11 Fouling effect at a coolant flowrate proportional controlled base.

PHASE cA, T, Te, phi, UA
PAR Tset,Kc,Te0,Teset,recycle,taui
FUN Feset,Fe,fouling,ve
Fe=Feset-Kc*((Tset-T)+phi/taui)
fouling=-UA*(1-0.99)/3600
phi'=(Tset-T)
UA'=fouling
pi=3.1415926536} ret
LOCBIF rhs 12 Fouling effect at a coolant flowrate proportional-integral controlled base.

PHASE zeta, T, Td
PAR UA,Tset,Teset,Kc,delay
FUN tau,NTU,dTad,Te
tau=Vr/F
NTU=UA/(rho*Cp*F)
Te=Teset+Kc*(Tset-Td)
Td'=(T-Td)/delay
LOCBIF rhs 13 Base case + coolant temperature proportional control with delay.
PHASE zeta, T, phi, Td
PAR UA,Tset,Teset,Kc,taui,delay
FUN tau,NTU,dTad,Te
tau=Vr/F
NTU=UA/(rho*Cp*F)
Te=Teset+Kc*((Tset-Td)+phi/taui)
phi'=Tset-Td
Td'=(T-Td)/delay
LOCBIF rhs 14 Base case + coolant temperature proportional-integral control with delay.
Appendices
126
PHASE cA, T, Te, Td

PAR Tset,UA,Kc,Te0,Teset,recycle,delay
FUN Qfeed,Qreaction,Qtransfer,Qcool,zeta,Feset,offset,Fe,ve
Qtransfer=UA*(T-Te)
zeta=1-cA/cA0
offset=Tset-Td
Td'=(T-Td)/delay
pi=3 1415926536} ret
LOCBIF rhs 15 Base case + coolant flowrate proportional control with delay.
PHASE cA, T, Te, phi, Td
PAR Tset,UA,Kc,Te0,Teset,recycle,taui,delay
FUN offset,Feset,Fe
offset=Tset-Td
Fe=(sgn(Feset-Kc*(offset+phi/taui))+1)*(Feset-Kc*offset)/2
phi'=(Tset-Td)
Td'=(T-Td)/delay
pi=3.1415926536} ret
LOCBIF rhs 16 Base case + coolant flowrate proportional-integral control with delay.
PHASE cA, T, Td
PAR UA,Te,Tset,Fset,Kc,delay
tau=Vr/F
Qtransfer=UA*(T-Te)
zeta=1-cA/cA0
offset=Tset-T
F=(sgn(Fset-Kc*(Tset-Td))+1)*(Fset-Kc*(Tset-Td))/2
Td'=(T-Td)/delay
LOCBIF rhs 17 Base case + flowrate proportional control with delay.

PHASE cA, T, phi, Td
PAR UA,Te,Tset,Fset,Kc,taui,delay
FUN tau,Qfeed,Qreaction,Qtransfer,zeta,offset,F,deviation,dcA,dT
tau=Vr/F
Qtransfer=UA*(T-Te)
zeta=1-cA/cA0
offset=Tset-Td
F=(sgn(Fset-Kc*((Tset-Td)+phi/taui))+1)*(Fset-Kc*((Tset-Td)+phi/taui))/2
dcA=F/Vr*(cA0-cA)-k*exp(-E/R/T)*cA
dT=(Qfeed+Qreaction-Qtransfer)/(rho*Cp*Vr)
phi'=offset
Td'=(T-Td)/delay
INIT={Vr=5 rho=800 p=2000 A0=5000 dH=-1 6E5 k=2 505E7 E=9E4 R=8 31441 T0=303} ret
LOCBIF rhs 18 Base case + flowrate proportional-integral control with delay
Appendices
127
Part I
THEORY
Part II
APPLICATION
Base case model:
Source:
Mechanism:
Reactor:
Dynamical behaviour:
Controllers:
Manipulated variables:
Active parameters:
Limit cycle system

Van Elk et.al. [22, p4874]
Irreversible exothermic A(l) P(l)
Cooled CISTR
Limit cycles
None, P or PI
Tcool, V,cool or V
UA, Kc, I and d
Base case
Van Elk et.al. [22, p4874]
No control
Control
Chapter 6
RDSM
6.3.1
RCSM
6.3.2
Proportional
control
Chapter 9
Proportional
Integral control
Chapter 10
Delay
Chapter 11
UA
UA
UA
UA
Kc
Kc
Kc
Kc
RDSM
9.1
RCSM
9.2
RDSM
10.1
Tcool
Tcool
Tcool
Tcool
V,cool
V,cool
V,cool
V,cool
RCSM
10.2
P
11.1
PI
11.2
process
11.3
Part III
APPENDICES
FAQ
(Frequently Asked Questions)
What is a CISTR?
!"This is the abbreviation of continuously ideally stirred tank reactor i.e. a chemical reactor
frequently applied in the process industry.
What is a limit cycle?
!"A limit cycle is the phenomenon of self-sustained oscillations of the temperature and
concentration in a chemical reactor, permanently around an equilibrium point, which never
reaches the particular steady state value.
What happens if a limit cycle occurs?
!"In the worst case, the temperature can rise to extreme heights; consequently, the reaction can run
away from the desired steady state, which can lead to dangerous situations.
Do limit cycles exhibit in the process industry?
!"Yes, in the literature the existence of limit cycles has been demonstrated and described. It
predominantly occurs in cooled CISTRs.
What is a stability map?
!"In a stability map, the stable and unstable regions, in which limit cycles and multiplicity exhibit, can
be indicated, based on a mathematical description.
When can a limit cycle appear?
!"A limit cycle appears when the system conditions are such that they can be classified as unstable.
This can be derived from a so-called stability map, which is divided into stable and unstable
regions.
How can a limit cycle exhibit?
!"In the worst case when the controller device fails and the temperature changes towards an
unstable region in which limit cycles can arise. Conversely, when due to long time changes of
operation conditions e.g. fouling of the heat transfer unit, the working point has been modified.
How long does a limit cycle oscillation take?
!"The oscillations time varies from minutes to hours. In this report, in the most uncomplicated case,
the time between two peaks is on the order of 30 minutes.
How can a limit cycle be prevented?
!"Of course, when a reactor is perfectly designed and nothing goes wrong
What does perturbation and bifurcation mean?
!"Literally, perturbation means excitement or commotion. In fact, it means an alteration of a specific
quantity. A bifurcation point is a point at which two branches of a curve coalesce as a parameter is
varied i.e. a critical value.
How can a limit cycle be controlled?
!"With a suitable control device and tremendous appropriate strategy plan.
Can a limit cycle be predicted?
!"If the engineer is capable of making the preferred stability maps. The stable and unstable regions
in such a stability map can provide the desired information considering the occurrence of limit
cycles.
What is an orbit curve?
!"In an orbit curve an important system variable (temperature, conversion) is plotted against the
time, in most cases to investigate the dynamical behaviour of this particular system variable.
Waarom heb je geen motto voorin je verslag?
!"Een motto is een punt op een lijn. (Onno Kramer 2001).
END OF REPORT
Printed at: Thursday, 25 April 2002 at 1:11 PM:

Notes:
Notes
GLOSSARY
Automatic feedback control
Autothermal reaction
Base case
Bifurcation
Blowout velocity
Cascade control
Closed loop
CISTR
Continuous process
Controlled variable
Control state
Control strategy
Control type
Damping ratio
Dead time
Derivative action
Distance velocity lag
Disturbance
Droop
Dynamic model
Equilibrium state
Equilibrium curve
Endothermic reaction
Exothermic reaction
Feedback control
Feed-forward control
Fold curve
Fundamental model
Gain
Heat capacity
Heat production rate
Heat withdrawal rate
Hopf curve
Input signal
Integral action
Integral time constant
Integral windup
Limit cycle
Limit point
Line-out time
Linear system
Local stability
Loop
Manipulated variable
Mathematical model
Non-linear
Non-self regulating
Notes
Feedback control in which a controller device monitors the controlled variable of interest and
commands a manipulated variable in order to maintain a desired value of the controlled variable.
A system, which is completely self-supporting in its thermal energy requirements. The essential
feature of an autothermal reactor system is the feedback of reaction heat to raise the temperature and
hence the reaction rate of the incoming reactant stream.
The mathematical base case model is a collection of mathematical relationships between process
variables, which purports to describe the behaviour of the actual physical process.
This means an alteration of a specific quantity. Bifurcation has been defined as the appearance of
topologically phase portraits under variation of parameters. In fact, a bifurcation is a change of the
topological type of the system as its parameters pass through a bifurcation i.e. critical value.
In case the throughput is too large, the cooling effect of the feed flow is causing the chemical reaction
to extinguish because the average reactor temperature drops rapidly.
Control scheme where the manipulated variable output of one controller becomes the set-point of
another controller.
See feedback control.
Continuously ideally stirred tank reactor.
A process in which material passes in a continuous stream through the processing equipment. Once
the process has established in a steady state operating state, the nature of the process does not
depend on the length of time the process is operating.
In a control loop, the variable of which is sensed to originate a feedback signal.
The current condition of the control entity e.g. process management, unit supervision or process
control. The control state defines how the control entity will operate and how it will respond to
command.
The control schemes are composed of descriptive information that identifies the process control types
of objects i.e. equipment module control, loop, device, indicator and/or status.
Any of the control that directs, initiates and/or modifies the execution of procedural control and the
utilisation of equipment entities.
The ratio of the second peak to the first peak in the process variable response, only if the process
variable response is underdamped.
The interval of time between initiation of a input change or stimulus and the start of the resulting
observable response.
Calculation of the controller manipulated variable is based on the rate of change of the process
variable, also called rate action.
See dead time.
Process influence that affects the process variable but is not manipulated by the controller.
See offset.
Mathematic model that attempts to capture the output response of a system as it changes as a
function of time.
This state refers to a system at rest.
This curve is the steady state solution of the ODE system as a function of one active parameter
A reaction which can maintain exclusively if heat continuously is supplied.
A reaction which produces heat, increasing the reactor temperature, which in turn increases the rate
of reaction. This positive feedback loop can result in a thermal runaway.
Control scheme that uses knowledge of the output to take corrective action.
Control scheme that eliminates or reduces the disturbance effect on the process variable by using a
measurement of the disturbance to modify the manipulated variable.
The Fold bifurcation curve represents the border between static stability and static instability as a
function of two active parameters
A mathematic process model; based on fundamental concepts e.g. the conservation of mass and
energy.
The ratio of the change in the steady state output to a step change in the input provided the output
does not saturate.
The product of the overall heat transfer coefficient U and the area A in which the heat is transferred
The heat produced in a reactor by a chemical reaction. The HPR is a function of the properties of the
reaction mixture.
The total heat removal from a reactor is the sum of the heat removed by the cooling medium and the
cold feed supply. The HWR depends on the residence time in the reactor
The Hopf bifurcation curve represents the border between dynamic stability and dynamic instability as
a function of two active parameters
A signal applied to a device, element or system.
The controller manipulated variable is based on the integrated error i.e. the set-point minus the
process variable.
The reciprocal of the integral gain.
When a controller has integral action, a persistent error will cause the integral term to increase or
decrease to a value of larger magnitude.
Self sustained oscillation e.g. reactor temperature or conversion.
Fold bifurcation.
See settling time.
A system in which the time response to several simultaneous inputs is the sum of their independent
time response. A linear system is generally represented by a set of linear differential equations.
Using linearised models the stability at the steady state is determined.
A single-input, single output controller that monitors a transmitter and manipulates a physical quantity,
usually a control valve, in order to force a process variable to the desired set-point.
A quantity varied by the controller in order to affect the controlled variable.
A system of equations whose solution, given specific input data, is representative of the response of a
process to a corresponding set of input.
A system that is not linear, see linear.
A process in which both inflow and outflow are independent of the controlled variable.
Offset
Operating point
Operating state
Orbit curve
Output signal
Overshoot
Peak time
Perturbation
PI-control
PID-control
Process
Process control
Process model
Process operation
Process variable
Proportional action
Proportional gain
Pure delay
Residence time
Response time
Ride hand side
Rise time
Runaway
Saddle node
Self regulating process
Set-point
Settling time
Stable system
Stability map
State
Steady state
Time constant
Transfer function
Transportation lag
Underdamped
The difference between the set-point and the actual value of the process variable when the system
has reached steady state.
The normal operating value for the system variable.
The currents condition of the equipment entity e.g. process cell, unit or equipment module. The
operating state defines how the entity will operate and how it will respond to command.
In a orbit curve an important system variable is plotted against the time. In most cases to investigate
the dynamical behaviour of this particular system variable.
A signal delivered by a device, element or system.
The maximum execution beyond the final steady state value of outputs as result of an input change.
The time from the set-point step change to the time of the first peak of an underdamped process
variable response.
Translated: excitement or commotion. See bifurcation.
Abbreviation for proportional plus integral control. A controller whose output is the sum of proportional
action + integral action.
Abbreviation for proportional plus integral plus derivative control.
Physical or chemical change of matter or conversion of energy e.g. change in temperature or
concentration.
The control entity that encompasses the basic discrete regulatory and equipment module procedural
control elements.
An overall model of the process that describes the processing actions required to convert raw material
into finished product.
Process operations represent major processing activities, which usually result in a chemical or
physical change in the material being processed.
The measured variable of the controller variable that, along with the set-point, is used by the controller
to calculate a value of the manipulated variable.
The controller manipulated variable is calculated as a constant, called the proportional gain, multiplied
by the error, the set-point minus the process variable.
The ratio of the change in the manipulated variable due to proportional action to the change in the
error.
See dead time.
The quotient of the reactor volume and the throughput i.e. the average time the species remain in the
reactor.
See settling time.
LOCBIF source code.
The time required for the process variable to go from 0% to 90% of the steady state change.
See exothermic reaction.
Fold bifurcation.
This process has a steady state relationship between the controlled variable and either inflow or
outflow. The output variables tend to a steady state after the input variables have reached constant
values.
An input variable which sets the desired value of the controller variable.
The time from the set-point change to the time that the process variable response has settled within a
certain percentage band of the final value, usually 2% - 5%.
A system is considered to be stable if a bounded input signal always returns in an output signal that is
also bounded.
In a stability map, the stable and unstable regions, in which limit cycles exhibit, can be indicated
The current condition of the physical or control entity. The state also defines how the entity will
operate and how it will respond to commands.
The long-term output response of a system after it has been disturbed.
In an expression for linear system time response, the time constant is the value in the response term
e-s. In a transfer function, the time constant is the value in the denominator 1+s. For the output of a
first order system whose input is a step signal, the time constant is the time to complete 63.2% of the
total output change.
For a continuous time system, the transfer function is the ratio of the Laplace transform of the output
variable to the Laplace form of the input variable, with all initial conditions assumed to be zero.
See dead time.
The time response of the system to a step signal input has overshoot.
I have a small
limit cycle problem,
can you help me?
Notes

Limit Cycles in A CISTR

Uploaded by

Copyright:

Available Formats

Limit Cycles in A CISTR

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Limit Cycles in A CISTR

Uploaded by

Copyright:

Available Formats

Universiteit Twente

The elimination of limit cycles in a CISTR through

MODEL OF A COOLED CISTR

BEHAVIOUR OF THE MODEL

THEORETICAL PREDICTION OF LIMIT CYCLES

STARTING POINT OF RESEARCH

STRUCTURE AND OBJECTIVES OF PART II

PRELIMINARY CONSIDERATION CONTROLLABILITY

CONTROLLED SYSTEM WITH DELAY

CONTROLLER CONFIGURATION TUNING

THE EFFECT OF FOULING ON THE BASE CASE STABILITY

APPENDIX 1. DERIVATION CISTR

LIST OF STABILITY MAPS

Table 3 Reactor designer stability maps (proportional control)

Figure 43 (Tcool,0 = 303 [K])

Table 4 Reactor controller stability maps

Table 5 Reactor controller stability maps regarding delay

2.2 The mass and heat balance of the reaction mixture

Figure 1 Schematic drawing of a cooled or heated CISTR.

With the following overall reaction rate expression for reaction 1:

Both equations 2 and 3 contain non-linear functions of T.

Adiabatic temperature rise:

2.3 Modelling of the Cooling

Coolant heat balance

cool C P,cool V ,cool

Coolant temperature expression

2. Logarithmic temperature mean

3.2 Static behaviour of the cooled CISTR

Figure 2 CISTR with cooling jacket.

Figure 3 The thermal behaviour of an irreversible, exothermic

3.3 Dynamic behaviour of the cooled CISTR

Stable self-regulating process

Figure 4a Asymptotic damping of a

Figure 4b Oscillating response (spiral

Figure 4g Response of an oscillating

Figure 4h Unstable system. Extinction. Figure 4i

At first the response shows

Figure 4c Stable self sustained

Figure 4d Limit cycles, which occur

Figure 4e Asymptotic oscillations (limit

oscillations (limit cycles) after

cycles). Unstable system.

Figure 5a Self-sustained reactor

Figure 5b Conversion limit cycles.

Figure 5c Limit cycle around a static

4.2 Linearised perturbation method

4.3 Bifurcation theory

Figure 6 Bifurcation diagram.

LOCBIF cannot deal with integral equations or partial

TIME variable declaration

LOCBIF rhs 1 Structure of ODEs specification in LOCBIF.

4.5 Stability map

Fold bifurcation curve

Hopf bifurcation curve

Constructing the stability map

Equilibrium curve through

Figure 7d Addition Hopf curve,

Figure 7b Additional equilibrium

Figure 7c Addition Fold curve,