Chapter 1

Basic Modelling of Dynamic Engineering

Systems

This chapter1 presents an introduction to the basic issues of the modelling of dynamic engineering
systems. We will consider the role and use of models from a point of view that is common in the
engineering sciences.

1.1 Introduction
The use of models has a long tradition in science. Yet, there is a difference in attitude towards the use
and role of models and modelling between the physicists and the engineers. Consider for example the
elementary properties of matter such as specific heat, viscosity, or thermal conductivity. A physicist
will not be entirely satisfied with these concepts because they are not elementary enough. He will try to
explain these concepts in terms of the properties of more elementary concepts such as molecules, atoms,
and electrons. The engineer is also dissatisfied because these concepts in general are too elementary. For
most problems in engineering it is desirable to use these concepts merely as building blocks for more
elaborate concepts, by which the behaviour of larger processes can be explained. Thus the engineer is
interested in discovering the behaviour of macroscopic processes, whereas the physicist will naturally
concentrate upon the microscopic behaviour or will consider nature at even smaller scale.
In engineering science, models constitute a main formalism for describing our knowledge about
engineering processes and systems. A model basically describes a limited number of properties of the
real world in terms of verbal, graphical, or mathematical language. With the assistance of models we
describe that part of the real world that is thought to be relevant to the solution of a particular problem
or to a class of related problems.
In the sequel, a systematic approach will be developed that enables to formulate dynamic mathe-
matical models for physical systems and processes. Essentially, there exist two different approaches for
obtaining models.

Experimental modelling The use of experimental data, obtained by measuring and collecting the
values of variables of the real system or process during well-planned experiments, and through
subsequent data processing, should discover systematic relationships amongst the variables.
These relationships are considered as a model. We further may wish or require this model to be
1 4CM40 Physical and data-driven modelling chapter 1 c O. H. Bosgra† TU/e 2008-2023

1
2 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

able to explain and predict the behaviour of the studied system or process. The joint actions of ex-
periment design, experimental data acquisition, data processing, and model validation are known
under the headings of experimental modelling, empirical modelling or system identification.
Theoretical modelling Using the accepted theory of the underlying sciences, equations are derived
that describe the basic behaviour of our system or process. In general the equations involved will
be conservation laws or other relations that have been derived from first principles. Additionally,
we need to describe the physical properties and other relevant descriptions of the materials in
the system or process. Thus, the underlying theories that we accept to build the model are our
basic axioms and assumptions in the logical process of model construction. This approach will
be called theoretical or fundamental model building. It is also known under the name physical
modelling.
The steps to be followed in experimental model formulation in general require a considerable amount
of a priori knowledge. This a priori knowledge must be brought in a format such that certain basic
decisions regarding variables, experiments, and responses can be taken with the required confidence.
This implies that we must execute certain steps in the process of theoretical modelling in order to have
a sufficient basis for experimental model formulation. Thus the second approach can be considered
as the natural starting point for any model building activity. Additionally, in many cases the first
approach is the natural sequel to the second approach in the sense that it allows to verify or validate the
structure and parameters in a formulated model, or it may enhance and amplify certain insufficiently
understood or incompletely quantifiable parts of a fundamental model. In these notes an approach to
theoretical model formulation will be covered. In several cases the connections with certain aspects of
experimental modelling will be discussed.
The process of modelling is widely applied in many fields in science and engineering. However,
the amount of knowledge available for the formulation of theoretical models is different in the various
fields. In many fields of engineering there is sufficient knowledge available to formulate simple over-all
models that describe the behaviour of systems and processes. However, for example in chemical process
engineering the operation of most chemical process equipment is so complicated that the detailed
behaviour of chemical processes can only in part be described accurately by mathematical models
based upon our present knowledge. Thus there are severe limits in our ability to formulate theoretical
models for engineering systems. In many cases the use of experimental modelling techniques is a
requirement that by necessity complements the theoretical approach.
In these notes the formulation of models is aimed at obtaining a description of the dynamic
behaviour of processes under transient conditions. This implies that we will formulate the equations of
motion of the process variables that describe the evolution of the process as a function of time. Our
models will formulate the process dynamics in a form as required for the understanding of process
operations such as startup and shutdown, or for studying the transitions from one operating condition to
another one as, e.g., required by grade changes in a production plant or by changes in the composition
of the feedstock. Process dynamic models also are of great importance for providing control engineers
with qualitative and quantitative descriptions of the transient behaviour of processes that are to be used
in model based control system design.
By its proper nature, a model is limited in the sense that the real world is approximated by definition.
Only those properties and phenomena that are relevant to the purpose of the model should be included
in the model. The real world is so highly detailed, and possesses so rich a complexity, that even a small
part of the real world can never be exactly represented in an isomorphic mapping sense by a single
mathematical model, notwithstanding its high number of equations and its great fidelity. Thus the
process of model building is fundamentally a process of making appropriate and clever approximations.
1.2. THE FORMAL STEPS OF THE MODELLING PROCESS 3

This fact forces the modeller to consider the object of his activity from an over-all point of view and
from a certain distance. He must have an overview over the most important physical phenomena in
this object, evaluated from the point of view of the intended use or the foreseen purpose of the model.
The modeller must have an attitude that enables him to make decisions regarding the most important
phenomena to be embodied in the model. For this reason, we must consider the process of modelling
as a decision process.
A model can be defined as a “small representation of a planned or existing object.” The adjective
“small” connotes the central idea that a model is something less than reality. Thus, in the sciences,
models are sought which illuminate natural phenomena. The objective is to strip away all phenomena
that are not thought to be essential in view of the intended use of the model. Thus a model must be able
to explain those phenomenological patterns of interest in terms of a set of easily understood elements.
In this context a model is a theory which constitutes a set of propositions or laws from which may
be deduced those facts which are exhibited in the real system or process under consideration. More
fundamentally, a theory must explain the less comprehensive laws, established by observation of reality,
in the sense of introducing ideas which are more familiar or, in some way, more acceptable. Moreover,
theory must predict new laws and these laws must turn out to be true.
This last notion illustrates sharply the approach that has come to be called the scientific method. The
scientific method of establishing an understanding of any physical phenomena is generally identified as
consisting of the following three phases:
• Initial observation

• Formulation of a theory (in our language, a model)

• Experimentation (model (in)validation) and prediction of new observations (model use)

The completion of the last stage will in many cases suggest refinements to the theory (model), and
the process is repeated. The emphasis on observations and experiments has its roots in the empiricist
philosophy which has been at the heart of modern science.

1.2 The Formal Steps of the Modelling Process

The scientific formalization of the approach to modelling sketched in the preceding section can now be
formulated and involves the following main features.
1. Defining the problem. The purpose of the model must be formulated and clarified. This requires
a certain mature familiarity of the modeller with the object to be modelled. Also an understanding
of the future model use will help in deciding which properties of reality should be represented in
the model. In general the modeller must make decisions regarding the desired degree of detail in
the model. Higher granularity in general implies a higher amount of complexity, and in many
cases the dominant phenomena in the model can be represented with higher accuracy than the
minor details. Representing these details may be counter-productive in terms of the intended use
of the model. Thus unnecessary high detailing should be avoided.

2. Define system boundaries and interconnecting variables. The process of modelling is aimed
at describing only a part of the real world. The system boundary defines in its interior the
system as opposed to its exterior which are the surroundings. The modeller himself is part of the
surroundings. The system and its surroundings are coupled through signals which are process
variables that originate in the surroundings or in the system. Those variables that originate in the
4 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

surroundings and influence the system by exerting an action upon the system by passing across
the system boundary are called input variables contained in a vector u(t). We consider these
influences in a causal fashion, i.e., such that the cause precedes in time the effect. A system is
called causal if past and present internal and/or output variables do not depend on future input
values.
From a modelling point of view there is no distinction between the input variables that can
be manipulated at will, and those input variables that occur in the surroundings as a result
of unmodelled processes or of which the origin or cause is not exactly known (disturbances).
Finally, those variables that are the result of the relationships within the system and which pass
outwards across the system boundary such that they become available to the surroundings are
called the output variables of the system, contained in the vector y(t). All variables that occur in
the system and that are not included in the set of input variables are called the internal variables,
represented by the vector z(t). See the representation in the Fig. 1.1. In general, the total set
of variables in a model setup is the sum of the input variables and the internal variables. The
output variables are to be considered as a subset of this total set of variables. Note that the
formulation of a mathematical model by its very nature yields a set of variables that constitute
the total set of variables mentioned before. If this total set of variables is available as a result of
the initial modelling decisions, then it is a decision of the modeller to assume which subset of
these variables is specified by the surroundings and thus constitutes the set of input variables.
The remaining variables are the internal variables. As the input variables are specified by the
surroundings (i.e., by the modeller), these variables must be considered as known variables in
our problem setup. Then our model must be able to (uniquely) determine the internal variables,
given the input variables. This implies that in general we must have a number of (independent)
model relations that equals the number of internal variables.

z(t)
u(t) r internal y(t)
variables
m input output
variables system variables

surroundings

Figure 1.1: The system concept requiring r independent model equations

3. Development of a conceptual model. Once the model variables have been preliminarily iden-
tified, the next step is the development of a conceptual physical model, containing only those
phenomena which are considered as relevant to the intended use of the model. A thorough knowl-
edge of the physics of the process plus intuition gained from practical experience is valuable
and important in this step. In its most simple fashion, the conceptual model only consists in the
mind of the modeller. The model can be verbally described, and the relevance of incorporating
or neglecting certain phenomena can be discussed.
1.2. THE FORMAL STEPS OF THE MODELLING PROCESS 5

4. Formulation of model hypotheses. In this step the real physics in the process are replaced by
a simplified hypothetical physical representation. The precise assumptions needed to arrive at
this simplified representation are to be formulated in a number of model hypotheses or model
assumptions. In a certain sense the predictive power of the resulting model is reduced by adapting
its properties to the specific situation at hand. This will be greatly counterbalanced in general by
the increase in model simplicity. The chosen degree of simplification determines the character of
the mathematical relationships that will finally represent the physical phenomena.

5. Decomposition into subsystems. In many cases, a modeller’s intuition can only keep track of
a limited number of phenomena simultaneously. For this reason it can be beneficial to split
up the complete process into a number of parts by formulating spatial boundaries by which
subsystems are created. The interconnection of all subsystems gives rise to the original system.
The subsystems thus defined must satisfy individually the requirements of a well-posed system.
This implies amongst other things, that the interconnection of the various subsystems should
be consistent, the number of internal variables should be equal to the number of (independent)
model relationships available for each subsystem, and the inputs and outputs of the original
system should be contained in the various inputs and outputs of the subsystems.

6. Formulation of conservation laws and other mathematical equations. The hypotheses formu-
lated in the previous steps must enable now to write down (if appropriate) the basic conservation
laws for energy, momentum and mass. These equations will be complemented with rate equations
for mass and heat transfer, with descriptions of the thermodynamic and other physical properties
of the materials in the process, and with all other fundamental or correlations-based relationships
that are needed in order to arrive at a physically and mathematically consistent model. The
result is a set of algebraic and differential equations having a certain mathematical structure, and
containing a set of parameters that fix the model behaviour quantitatively.

7. Mathematical simplifications and approximations. The mathematical model can be examined

for its properties. Certain basic a priori requirements of a model can be investigated, such as
the properties of causality, continuity, and conditioning with respect to possible numerical
implementation for simulation purposes. Also further properties such as stability, time scales,
and independence of the set of model equations can be established. On the basis of these findings
it can be decided to simplify or to approximate the mathematical representation of the system.

8. Numerical solution of the model equations, simulation. Further confidence in the model can
be gained by solving the model equations for well considered input signals (simulation), and by
observing the behaviour of the internal and output variables. By exercising with the model in
this way, the model can inspire confidence or the modeller can experience certain weaknesses of
the model formulation.

9. Validation of the model and its hypotheses. Comparison of predicted model responses and
similar responses obtained at the real system may provide the final proof of confidence needed to
validate a model. It must be stressed that the experiments must be representative for the intended
purpose of the model.

The main steps of the approach are listed in Table 1.1.

6 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

Table 1.1: Basic decisions in the modelling procedure

reality
⇓ ⇐H intended model purpose
process under consideration
⇓ ⇐H select system boundary, inputs, outputs
system and surroundings defined
determine relevant phenomena
⇓ ⇐H
and formulate model hypotheses
conceptually simplified physical model
formulate conservation laws
⇓ ⇐H
and other model relations
mathematical model having
structure and parameters
⇓ ⇐H mathematical approximations
simplified mathematical model
⇓ ⇐H experimental model validation
model satisfies requirements?
⇓ H⇒ yes, model ready for use
no, return to previous steps

1.3 Formulation of Conservation Laws for Process Models

The equations of conservation of some quantity (mass, momentum, energy) are fundamental laws of
nature. These equations apply both in the steady state and under unsteady state conditions. We first will
consider macroscopic conservation laws for mass and energy. These are integral balance equations that
are applied to a volume inside a macroscopic control surface, as opposed to microscopic conservation
laws that are formulated for an infinitesimal or elemental volume.

1.3.1 Conservation of mass

The density and concentration of a chemical species that possibly reacts are often measured in moles
per unit volume. A mole of a species contains Avogadro’s number (6.022 × 1023 ) of molecules. Moles
can be converted to mass units by multiplying by the molar mass which has units of [kg/mol].
Consider a system defined by its (fixed) macroscopic control surface or system boundary S, enclosing
a (macroscopic) volume V as shown in Fig. 1.2. Stated verbally, a macroscopic conservation law for
mass for a chemical species A within the volume V reads, in conservation form:
     
d Net accumulation Net transport Net generation
= + (1.1)
dt of A in V of A through S of A in V
Note that the left-hand side of (1.1) specifies the rate of change of mass accumulation within the control
volume V , and not the amount of accumulated mass itself.
Let ρ A denote the density of species A, then the accumulation term in the left-hand side of (1.1)
describes the rate of change per unit time of the accumulated mass of A within the volume V :
Z
d
ρ A dV (1.2)
dt V
1.3. FORMULATION OF CONSERVATION LAWS FOR PROCESS MODELS 7

system boundary S

mass V
ρ Ai φi
flow in

mass perfectly
ρ A φo
flow out
mixed

surroundings
Figure 1.2: System boundary for macroscopic mass balance

This form of the left-hand side allows the density to have varying values over the macroscopic volume V .
In fact the integral expression performs a summation over all elementary volume parts d V of the
contribution of each ρ A d V to the total accumulation of mass of A. This formulation in fact does
not reflect the overall character of the macroscopic balance and is not practical. In terms of cause
and effect, we must consider the right-hand side of (1.1) as the cause, creating an effect in terms
of a variation of overall accumulated mass in the left-hand side. Our macroscopic conservation law
determines the overall accumulation and thus does not provide a mechanism to deal with spatially
varying accumulation phenomena. Consequently, the density ρ A should be considered as having the
same value over the whole volume V . This can be achieved by assuming that there exists a physical
mechanism that creates so much interchange of material over the volume V that ρ A has virtually the
same value over V . One possible assumption is that V behaves as if it were perfectly mixed. Thus,
under the assumption of V being perfectly mixed, our macroscopic conservation law for mass reads
dρ A V X
= ρ A k φk + G (1.3)
dt k

where the φk denote the volumetric flows across the boundary S each having density ρ Ak , and G denotes
the generation of species A per time unit within V , e.g., as a result of chemical reaction. In this
expression, φk are positive and negative for in- and outgoing flow, respectively. The k is supposed to
run over all flows that cross the system boundary. Note that under the assumption that the volume
V behaves as a perfectly or ideally mixed fluid, the density ρ A has the same value over the whole
volume V , whereas this value may vary in time. The assumption of an ideally mixed fluid implies that
it is assumed that each element in the volume V will vary in time in exactly the same fashion. It also
implies that any outward flow across the boundary S has a density that equals ρ A . Thus for one flow in
and one flow out, the macroscopic conservation law (1.3) for mass reads
dρ A V
= ρ Ai φi − ρ A φo + G (1.4)
dt
where now the negative sign of the outward-bound flow has been explicitly formulated in the expression,
and thus we consider positive values for the flows in the direction given by the arrows in Fig. 1.2.
8 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

While the general principle of the conservation law (1.1) is maintained, the actual equations may have
other forms than (1.3) or (1.4) depending upon the structure of the system under consideration. Under
the assumption that the density ρ of the fluid is constant, we still can have variations of accumulated
mass within the system, as the following example shows.

1.3.2 Example: Mass accumulation in vessel

Consider the simple flow system shown in Fig. 1.3. The steps to be executed for modelling this process
are as follows.
system boundary S

h(t)

φi (t) φo (t)

Figure 1.3: Fluid buffering vessel

1. Purpose of the model: To describe the behaviour in time of the fluid level h(t) varying with the
purpose of buffering variations in input mass flow

2. System boundary: Around the vessel

3. Variables across boundary: φi , φo , h. As there apparently is no mechanism inside the sys-

tem boundary to determine φi and φo , these variables are supposed to be determined by the
surroundings, i.e., φi , φo are inputs and h is output.

4. Model hypotheses: the fluid is incompressible; the vessel cross sectional area A is constant; the
velocity of the fluid inside the vessel is slow so that local variations in fluid level can be ignored.

5. Conservation law: Mass balance

d
(ρ Ah(t)) = ρφi (t) − ρφo (t) (1.5)
dt
which can be written in differential equation (state space) form:
dh(t) 1 1
= φi (t) − φo (t), h(0) = h 0 (1.6)
dt A A
Note that the model (1.5) can be depicted by the block diagram in Fig. 1.4 and by the structural
representation in Fig. 1.5. Also observe that the system is not asymptotically stable, as the (open)
1.3. FORMULATION OF CONSERVATION LAWS FOR PROCESS MODELS 9

ρ Ai φi (t) vessel
h(t)
ρ A φo (t) model

Figure 1.4: Block diagram of inputs and outputs of fluid buffering vessel

ρ Ai φi (t) initial m(0)

+ Z
1
h(t)
ρA A
−

ρ A φo (t)

Figure 1.5: Internal structure of model for fluid buffering vessel

integrator may give an unbounded response for bounded inputs. This indicates, that the model
represented by (1.6) is only valid for h(t) ≥ 0 and h(t) ≤ h max , where h max indicates the height
of the vessel. A simulation of this process should include these bounds. If one of the bounds is
reached, the model (1.6) should be replaced by a different one, one for the empty vessel case,
and one for the overflow case.

6. Model satisfies requirements? In retrospect, at this stage it is useful to reconsider the assump-
tions and starting points of the modelling exercise, and to make modifications if desired. For
instance, if in reality the flow φo depends on the value h, then one might wish to include this
relationship in the model. Figure 1.6 represents the situation where the outgoing flow is deter-

system boundary S

h(t)
X (t)

φi (t) φo (t)

pi (t) po (t)

Figure 1.6: Fluid buffering vessel with outward flow controlled by valve
10 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

mined by the valve opening and by the pressure difference across the valve. We assume that the
fluid pressure behind the valve po is atmospheric, as is the pressure above the fluid level. We
assume further, that the flow through the valve is turbulent, proportional to the square root of the
pressure drop across the valve. This pressure drop is proportional to the fluid level. We have the
following 5 variables: φi , φo , pi , X , h where X represents the input command determining valve
position. The model now consists of (1.6) and the following two relations:
φo (t) = K Av (X (t)) pi (t)
p
(1.7)
pi (t) = ρgh(t) (1.8)
where Av (X (t)) represents the valve cross-sectional area as a function of valve position X (t),
K is the valve coefficient, and g is the gravitational acceleration. Equation (1.7) represents the
flow through the valve and (1.8) determines the actual pressure drop across the valve. We have 5
internal variables and 3 model relations, thus two variables should act as inputs, i.e., are assumed
to be determined by the surroundings of the process. As φo , pi , and h are determined by the
model relations, we can assume φi and X to be input variables. This is in accordance with
our physical imagination. The block diagram with input and output variables, and an internal
representation of the model structure are given in Figs. 1.7 and 1.8, respectively. The following

h(t)
φi (t) vessel
+ pi (t)
X (t) valve
φo (t)

Figure 1.7: Inputs and outputs of fluid buffering vessel with outward flow determined by valve

φi (t) initial m(0)

dm(t) h(t)
+ dt
Z
1
ρ
− ρA

K
φ0 (t) r pi (t)
ρg


∗

Av (t) pi (t)

Av

X (t)

Figure 1.8: Internal structure of model of vessel and valve

observations can be made:

1.3. FORMULATION OF CONSERVATION LAWS FOR PROCESS MODELS 11

• Comparing Fig. 1.5 with Fig. 1.8, we see that the representation of the conservation law
represented in Fig. 1.5 is contained in Fig. 1.8. However, now there is a feedback loop
around the integrator.
• The sign of the feedback path is such that an increase in accumulation of mass will increase
the outflow of mass, i.e., there is negative feedback.
• Thus we may expect that the system will be asymptotically stabilized by the feedback
relationship. This implies that a permanent change in inflow φi or valve position X (t) will
lead to new equilibrium values for h(t), φo , and pi (t).
• The process is bilaterally coupled to the neighbouring flow system to the left and to the
right. This means that there is a mutual information exchange between the process and
the adjacent systems at both sides. The model expresses this coupling by assuming the
pressure at the right interconnection po to be atmospheric, i.e., not to change in relation to
the variables considered in the process. Thus φo (t) can vary without causing a response
in po (t). At the left interconnection, φi (t) is assumed to be determined by the upstream
flow process, and our model determines the reaction in pi (t) which acts as an input variable
for the upstream process.
• In Fig. 1.8, the block indicated by (∗) represents a multiplication block, i.e., its output
is the multiplication of all (three) input variables. Similarly, the block indicated by the
square root sign determines its output as the square root of its input. These blocks represent
nonlinear signal operations and show that the represented system is nonlinear, i.e., does
not satisfy the superposition principle.

7. Numerical solution of the equations: The block diagram of Fig. 1.8 can directly be used in a
simulation tool for continuous systems (e.g., Simulink). The following parameters are used:

A 1.0 [m2 ] ρ 980 [kg/m3 ]

g 9.8 [m/s2 ] K 0.01 [m3/2 /kg1/2 ]
A∗v 0.001 [m2 ] φo∗ 0.001 [m3 /s]

where the asterix values indicate the used steady state. The resulting value for h ∗ is 1/(0.98)2 =
1.0412. Figure 1.9 shows responses of flow out and fluid level on step disturbances in the valve
position. The nonlinear behaviour causes different responses for the positive and negative step
inputs. Figure 1.10 shows for fixed valve position the response to inlet flow step disturbances.
Figure 1.11 shows the buffering capabilities of the vessel for fast stochastic variations in input
flow, which are considerably reduced in the outlet flow.
12 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

−3 valve area level

x 10
1.1
1.6 1

0.9
1.4
0.8
m2

m
1.2 0.7

0.6
1
0.5

0.8 0.4
0 5000 10000 15000 0 5000 10000 15000

−3 flow out valve inlet pressure

x 10
1.6 11000

10000
1.4
9000
1.2 N/m2 8000
m3/s

1 7000

6000
0.8
5000

0.6 4000
0 5000 10000 15000 0 5000 10000 15000
time (s) −−> time (s) −−>

Figure 1.9: Response of fluid buffering vessel on step in valve opening

−4 flow in level
x 10
12
1

10
0.8
m3/s

8
0.6

6 0.4

4 0.2
0 0.5 1 1.5 2 0 0.5 1 1.5 2
4 4
x 10 x 10
−4 flow out valve inlet pressure
x 10
12
10000

10
8000
N/m2
m3/s

8 6000

4000
6

2000
4
0 0.5 1 1.5 2 0 0.5 1 1.5 2
time (s) −−> 4
x 10 time (s) −−> 4
x 10

Figure 1.10: Response of fluid buffering vessel on step in inlet flow, valve opening fixed
1.3. FORMULATION OF CONSERVATION LAWS FOR PROCESS MODELS 13

−3 flow in
x 10

2
m3/s
1

0
0 5000 10000 15000
−3 flow out
x 10

2
m3/s

0
0 5000 10000 15000
level
1.5

1
m

0.5

0 5000 10000 15000

time (s)

Figure 1.11: Buffering vessel: response on stochastic inlet flow variations, valve opening fixed

1.3.3 Conservation laws for energy

The formulation of correct energy balances is not a trivial matter, see, e.g., Denn [13, Ch. 5]. In general,
the total energy of a system consists of the internal, kinetic and potential energy terms. In many cases
where the energy balance is dominated by thermal effects, the kinetic and potential energy terms can
be neglected. A simple energy balance without generation phenomena for a one-phase incompressible
fluid of constant density ρ and constant specific heat C p in a constant volume V reads

dρV T (t) X
Cp = ρφk (t)C p Tk (t) (1.9)
dt k

where V is the volume within a fixed control surface S, T (t) is the time-varying temperature within
the control surface, and Tk (t) are the temperatures of the volumetric flows φk (t) through the surface S,
having positive and negative values for inflow and outflow, respectively. Note that the above conserva-
tion law is a macroscopic balance, i.e., we have explicitly or implicitly assumed that all fluid particles
inside V behave identically in time, as if V were perfectly mixed.

1.3.4 Example: Energy accumulation in vessel

For the case of one flow in and one out, the process is shown in Fig. 1.12 and the energy balance is
dρV T (t)
Cp = ρφC p [Ti (t) − To (t)] (1.10)
dt
Note that this is a macroscopic overall balance equation in conservation form, summing all heat
balance contributions in the total volume V . Due to this overall summation, each particle within V has
14 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

S
φ V

Ti (t)

φ
To (t)
T (t)

Figure 1.12: Process vessel under ideal mixing assumption

temperature T (t), possibly varying in time but spatially uniform. This corresponds to the assumption
of an ideally mixed internal flow pattern, and we may replace To (t) by T (t). The block diagram of this
model is shown in Fig. 1.13. The resulting energy balance can be written in state-space form, i.e., in

initial E(0)
dE(t)
Ti (t) + dt Z
1 T (t)
ρφC p
ρV C p
−

Figure 1.13: Internal structure of energy balance model for mixing vessel

first-order differential equation form with only a time derivative of a state variable in the left-hand side,
because the parameters ρ, V, C p are constant:

dT (t) φ
= [Ti (t) − T (t)] (1.11)
dt V
Note that this represents only one of the possible choices for a state variable. The effect of the feedback
loop shows a negative sign, which in this case implies that the system is asymptotically stable. The
equation is a linear first-order differential equation, having a time constant
V
τ= (1.12)
φ
which in physical terms equals the average residence time of particles flowing through the vessel.
Figure 1.14 shows the temperature responses for a vessel with time constant (average residence time)
τ = 1 hour.

1.4 Microscopic Conservation Laws

In the case that a concentration, density, temperature or any particles-related quantity can not be assumed
to be uniform over the volume V inside a control surface S, it is not useful to formulate a macroscopic
(integral) balance over the total volume V . In other words, if a volume V cannot be assumed to
1.4. MICROSCOPIC CONSERVATION LAWS 15

inlet temperature inlet temperature

51.2
100
51
90 0.02 width
50.8
degrees

degrees
80 50.6

70 50.4

50.2
60
50
50
49.8
0 1 2 3 4 5 0 1 2 3 4 5

outlet temperature outlet temperature

51.2 51.2

51 51

50.8 50.8
degrees

degrees
50.6 50.6

50.4 50.4

50.2 50.2

50 50

49.8 49.8
0 1 2 3 4 5 0 1 2 3 4 5
time (hrs) time (hrs)

Figure 1.14: Temperature response of mixing vessel on inlet temperature impulse and step

behave like ideally mixed, one could resort to the formulation of microscopic conservation laws. The
formulation of a conservation law for the elemental volume (δx, δy, δz) in V , and subsequently letting
the δx, δy, δz approach zero, yields a microscopic conservation law having the form of a partial
differential equation. For energy, a possible form is

∂T ∂T ∂T ∂T ∂ T ∂2T ∂2T
   2 
ρC p + vx + vy + vz =k + 2 + 2 +G (1.13)
∂t ∂x ∂y ∂z ∂x2 ∂y ∂z
Here the first term represents the variation of accumulation of energy per time units, and the next terms
in the left hand side describe the net bulk transport of energy through the surface S. The first three terms
in the right hand side describe transport through the surface S by thermal diffusion; and G denotes
a net generation term. The temperature T (x, y, z, t) now is a function of both time t and the spatial
coordinates x, y, z. The parameter k represents the thermal conductivity [kgm/s3 K], and defines the
thermal diffusion coefficient
k
DT := (1.14)
ρC p

where DT has units [m2 /s]. The three velocities vx , v y , vz represent the local bulk velocity components
of mass particles in the three coordinate directions. In deriving (1.13) we have assumed that ρ, C p
and k are constant over V . As a preliminary observation, we see:
• A microscopic balance such as (1.13) contains three different types of terms:

1. Variations per time unit of accumulated energy

2. Transport of energy through convective transport, i.e., as a consequence of the fact that
particles move with a velocity (vx , v y , vz ) whereas the coordinate frame is fixed
16 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

3. Transport of energy through diffusion, i.e., through a physical mechanism of molecular

interaction that creates an energy flux proportional to (the negative of) the actual spatial
gradient of the temperature

• The flow pattern in a vessel must be known or must be computed in some way if one wishes to
formulate and to solve a microscopic conservation law

• For the derivation of microscopic balance equations such as (1.13), consult, e.g., Himmelblau
and Bischoff [23] and Bird, Stewart and Lightfoot [9]. A derivation for the one-dimensional case
is given in the next section.

Along similar lines one obtains the following microscopic mass balance in three rectangular spatial
coordinates for a component A in a dilute binary system:

∂c A ∂c A ∂c A ∂c A ∂ cA ∂ 2c A ∂ 2c A
 2 
+ vx + vy + vz = D AB + + + GA (1.15)
∂t ∂x ∂y ∂z ∂x2 ∂ y2 ∂z 2

where c A is the molar concentration [mol/m3 ] of component A, and D AB [m2 /s] is the diffusion
coefficient.
Note that the microscopic balances (1.13) and (1.15) require the availability of the initial values
T (x, y, z, t0 ) and c A (x, y, z, t0 ), respectively, if we wish to determine their evolutionary solution
in time for t ≥ t0 . These initial values require the specification of T and c A at t = t0 over the
whole spatial domain. Additionally, at the boundary of this spatial domain the boundary values for
T (x, y, z, t) and c A (x, y, z, t) and their partial space derivatives have to be specified as a spatial
function of time. Moreover, the formulation of (1.13) and (1.15) assumes that the velocity field
vx (x, y, z, t), v y (x, y, z, t), vz (x, y, z, t) is known as a function of the three spatial coordinates and
time. These facts indicate that in general a microscopic formulation either requires the parallel computa-
tion of the actual possibly time-varying flow pattern, or requires considerable simplifying assumptions
in order to be of practical use. The first option can be realized by extensive numerical computations
using specialised numerical software packages, and will be utilised in cases where a rigorous model is
required. Simplifications may involve the observation that the flow pattern can be idealized, e.g., in
cases where the main flow direction is in only one coordinate direction. In that event a one-dimensional
spatial representation may suffice. Also note that the occurrence of velocity and diffusion coefficients
in (1.13) or (1.15) implies that when these coefficients become large, then the behaviour of the balance
equations approaches the behaviour of the ideally mixed macroscopic case.

1.4.1 One-dimensional microscopic balance equations

We will present the derivation of microscopic balances for the one-dimensional case. Consider the
flow in a pipe having essentially a uniform velocity v(t) over its cross-sectional area A. We call this
type of (idealised) flow pattern plug flow. Applying a (macroscopic) energy balance to an elementary
volume A1x as shown in Fig. 1.15 and assuming constant physical material properties, yields

dT (x + 1x, t)
ρC p A1x =
dt
∂ T (x, t) ∂ T (x + 1x, t)
 
v(t)ρC p A[T (x, t) − T (x + 1x, t)] + Ak − + (1.16)
∂x ∂x

The following issues explain the equation:

1.4. MICROSCOPIC CONSERVATION LAWS 17

x =0 x=L
T (x) T (x + 1x)

v v
A

x x + 1x
x
S S

Figure 1.15: Pipe flow with uniform velocity v(t) (plug flow)

• We have assumed the volume A1x to behave as if it were ideally mixed, i.e., its bulk temperature
equals the local temperature T (x + 1x, t) occurring at the outflow location

• In the right hand side we have two terms. The first one describes convective transport due to
particle velocity v(t)

• The second term describes diffusion. The occurrence of heat conduction or molecular diffusion
introduces a heat flux assumed to be proportional to the spatial temperature gradient. The
representation of this phenomenon can also be used to represent turbulent eddy diffusion.
For the diffusion phenomena, we assume a transport term across the boundary to be proportional to the
negative of the local spatial temperature gradient:
∂ T (x, t)
− (1.17)
∂x
with a proportionality factor given by the thermal conductivity k [kgm/s3 K]. Division of (1.16) by
ρC p A1x, and letting 1x → 0, leads to

∂ T (x, t) ∂ T (x, t) ∂ 2 T (x, t)

+ v(t) = DT (1.18)
∂t ∂x ∂x2
where the thermal diffusion coefficient DT is defined by (1.14). The derived one-dimensional plug-flow
based energy balance (1.18) combines the phenomena of convection and diffusion and is called a
convection-diffusion equation. This can be shown by first investigating the case DT = 0. In that case
we have the one-dimensional microscopic balance for convective transport
∂ T (x, t) ∂ T (x, t)
+ v(t) =0 (1.19)
∂t ∂x
with initial- and boundary conditions

T (x, t0 ) = Tinit (x) (1.20)

T (x0 , t) = Ti (t) (1.21)

Equation (1.20) specifies the initial temperature profile to be given for all x ≥ x0 in the spatial domain,
e.g., x ∈ [0, L] if x0 = 0 is at the inlet and x = L is at the outlet position of the considered pipe flow
system. Equation (1.21) specifies the inlet temperature as a function of time for all t ≥ t0 .
18 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

A different simplification occurs if we assume the velocity to be zero, v(t) = 0, in (1.18). In that
case we obtain the one-dimensional diffusion equation
∂ T (x, t) ∂ 2 T (x, t)
= DT (1.22)
∂t ∂x2
for which we need to specify initial and boundary conditions, which now will involve T (x, t), ∂ T∂(x,t)x
or both. The issues of how to specify initial- and boundary conditions for these partial differential
equations is important. In a future chapter, the notion of bilaterally coupled systems will be considered.
With respect to boundary conditions for the diffusion equation this implies that at each side of the
spatial domain (i.e., at x = 0 and at x = L), one can specify either a temperature or a heat flux as
boundary condition (or a single linear combination of these two), but not both. Specifying a heat flux at
a boundary is equivalent to specifying ∂ T∂(x,t)
x
at that boundary.
The terminology for the models defined in this section is as follows. Microscopic balances are
written in terms of partial differential equations. A system described by a partial differential equation is
called a distributed parameter system or an infinite-dimensional system. A system described by a finite
number of integral balances, i.e., macroscopic balances, is called a lumped parameter system, and its
mathematical representation is in terms of ordinary differential equations.

1.5 Linearity and Linearization of Nonlinear Models

Suppose that a lumped parameter system is described by the possibly nonlinear vector-differential
equation in state-space form
d
x(t) = f (x(t), u(t)) x(t0 ) = x0 (1.23)
dt
where x, f ∈ Rn , u ∈ Rm . The m variables in the vector u(t) are the input variables, to which the
system is assumed to be subjected by the surroundings. The n variables in x(t) are called the state
variables. The n model relations expressed by f (x, u) are assumed to be constant or time-invariant.
It can be observed that most processes behave globally as nonlinear systems, whereas many processes
behave locally as approximately linear systems. This observation provides the motivation for deriving
linearized models.

1.5.1 Definition of linearity

Linearity is a property which is identical to the fact that the superposition principle applies. The model
(1.23) is linear if it satisfies the following three basic superposition properties:
1. Zero-state linearity:

x(φx , t0 , au 1 + bu 2 ) ≡ ax(φx , t0 , u 1 ) + bx(φx , t0 , u 2 ) (1.24)

for all inputs u 1 , u 2 ∈ Rm and for all constants a, b ∈ R, where φx denotes the zero state of the
system.

2. Zero-input linearity:

x(ax01 + bx02 , t0 , φu ) ≡ ax(x01 , t0 , φu ) + bx(x02 , t0 , φu ) (1.25)

for all states x01 , x02 ∈ Rn and constants a, b ∈ R, where φu denotes the zero input function.
1.5. LINEARITY AND LINEARIZATION OF NONLINEAR MODELS 19

3. Decomposition property:

x(x0 , t0 , u) = x(φx , t0 , u) + x(x0 , t0 , φu ) (1.26)

for all x0 ∈ Rn and all inputs u ∈ Rm .

The three conditions in the definition can be given an interpretation as follows. Zero-state linearity
requires superposition to hold for the system when the initial state is zero and a state response occurs
only due to the effect of an input. Similarly, zero-input linearity requires the system to obey superpo-
sition in the response to nonzero initial states, when the inputs are zero. Finally, the decomposition
property requires a state response to be the sum of the effects of a response due to an input for initial
state zero, and a response due to a nonzero initial state for inputs equal to zero.
Now a nonlinear system is a system for which we only specify that the superposition principle
does not hold. This is a rather void characterization in the sense that it states the fact that an idealised
property such as linearity does not hold, without stating any further characterization or property of the
system. Thus we may expect that a further characterization of the properties of a nonlinear system is
required before we can make any assessment about its behaviour.
We may think of linear systems as locally linearized approximations of smooth nonlinear processes.
There are some strong arguments in favour of the study of linear systems :
• For linear systems, analytic solutions and analytically derived system properties are available.
It is often useful to determine analytically the parameter groups that are responsible for such
properties as steady-state gain, time constants, coupling coefficients, stability behaviour, or
sensitivity coefficients with respect to parameter variations.

• For nonlinear systems, such an analysis is in general not possible. We have to resort to numerical
simulations and numerical computations for nonlinear systems. A computer can provide a numer-
ical solution, but it in general does not provide background information about the relationship of
this solution to the parameters of the model.

• For many processes the function f (x, u) in the right hand side of the model (1.23) is a smooth
function of x and u, except possibly in a finite number of values of its arguments. This implies
that local linearization might be feasible in many cases.

• Many of the presently available techniques for stability analysis and for feedback design for
nonlinear systems are direct extensions of corresponding linear techniques and they require
smoothness of f (x, u) in (1.23). Thus a linearized analysis provides a result that also is relevant
to the understanding of the local behaviour of the nonlinear techniques.
On the other hand, we must keep in mind proper respect for the limitations of a linearized analysis.
For example, there exist certain stability issues for the global behaviour of nonlinear systems for which
a linearized analysis does not present any valuable information.

1.5.2 Linearization
Consider

x(t) = x ∗ (t) + 1x(t) u(t) = u ∗ (t) + 1u(t) (1.27)

where both the actual trajectory x(t), u(t) and the nominal trajectory x ∗ (t), u ∗ (t) satisfy the dynam-
ics (1.23). The terms 1x(t), 1u(t) are perturbations about the nominal trajectory. As an example,
20 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

such a nominal trajectory can be the trajectory that is foreseen in a startup or shutdown operation, or
during the transfer of a process from a set of existing operating conditions to a different one, e.g., a
grade change. The nominal trajectory may correspond to a set of stationary operating conditions for
the nonlinear state-space model (1.23). A necessary and sufficient condition for x ∗ , u ∗ to qualify in
defining stationary operating conditions is that they satisfy

0 = f (x ∗ , u ∗ ) (1.28)

Often there are infinitely many x ∗ , u ∗ satisfying (1.28).

In general
d(x ∗ + 1x) d d d ∗
= x ∗ + 1x = f (x, u) and x = f (x ∗ , u ∗ ) (1.29)
dt dt dt dt
so
d
1x = f (x, u) − f (x ∗ , u ∗ ) (1.30)
dt
Now expand the function f (x, u) in a Taylor series about the nominal conditions x ∗ , u ∗ :
∂f ∂f
f (x, u) = f (x ∗ , u ∗ ) + 1x + 1u + higher-order terms (1.31)
∂x ∗ ∂u ∗

By substituting (1.31) in (1.30) we have

d ∂f ∂f
1x(t) = 1x(t) + 1u(t) + higher-order terms (1.32)
dt ∂x ∗ ∂u ∗

For sufficiently smooth f and small 1x(t), 1u(t), the higher-order terms are small compared to the
linear terms and can be neglected. By doing this we approximate the nonlinear equations by linear ones
and this leads to the approximating linear state-space model
d
1x(t) = A1x(t) + B1u(t) 1x(t0 ) = 1x0 (1.33)
dt
where
∂ fi ∂ fi
Ai j = Bik = (1.34)
∂x j ∗ ∂u k ∗
i, j ∈ {1, 2, . . . , n} k ∈ {1, 2, . . . , m}

define the matrices

∂f ∂f
A= ∈ Rn×n B= ∈ Rn×m (1.35)
∂x ∗ ∂u ∗

If our nominal conditions x ∗ , u ∗ are not stationary, but rather constitute a nominal trajectory for the
system, the same concept of linearization still applies, be it now linearization along a (given) nominal
trajectory. This leads to the linearized model
d
1x(t) = A(t)1x(t) + B(t)1u(t) 1x(t0 ) = 1x0 (1.36)
dt
1.5. LINEARITY AND LINEARIZATION OF NONLINEAR MODELS 21

We assume that a nominal trajectory x ∗ (t), u ∗ (t) is available in quantitative form. Due to its inherent
variation in time, the partial derivatives
∂f ∂f
, (1.37)
∂x ∗ ∂u ∗
no longer are constant but vary when passing along the trajectory. Consequently, these matrices can
be expressed as a function of time and lead to time-varying matrices A(t), B(t) in (1.35). Note that
a time-varying system is the result of a linearization about a nonconstant nominal trajectory. If the
nominal trajectory is not varying in time, then the resulting matrices A, B in (1.35) are time-invariant.

1.5.3 Example: Process dynamics of mixing process

Consider the mixing process shown in Fig. 1.16 where two flows φ1 (t) and φ2 (t) having concentration
c1 (t) and c2 (t), respectively, are input flows for a vessel and are mixed. These concentrations describe
a single component occurring in both flows. The incoming flows are controlled by a valve in each of
the flow lines. Outflow occurs under the influence of gravity through an orifice. The vessel is agitated

1 2
φ1 (t) φ2 (t)

c1 (t) c2 (t)

S h(t)
c(t)

φ(t) c(t)

Figure 1.16: Mixing vessel having natural outflow through orifice

by a stirrer. The basic approach of modelling leads to the following steps.

1. Defining the problem. The purpose of the model is to describe how the outgoing flow φ(t)
and its concentration c(t) depend upon the flows φ1 (t) and φ2 (t) and their concentrations c1 (t)
and c2 (t).
2. System boundary. The mixing process takes place in the interior of the vessel; take system
boundary around the vessel.
3. Variables through system boundary. The flows φ1 (t) and φ2 (t) are the result of pressure drop
and valve openings in the surroundings, and thus are inputs to our system. Their concentrations
c1 (t) and c2 (t) describe properties of mass flow particles that pass the system boundary inwards,
and thus also have to be considered as input variables. By similar reasoning, φ(t) and c(t) are to
be considered as output variables.
22 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

4. Model hypotheses.

• The vessel is assumed to be ideally mixed. The concentration in the vessel is not space-
dependent and thus can be represented by a single variable c(t).
• We assume that no leakage or evaporation of the fluid occurs.
• The density ρ is assumed to be independent of the concentration c(t).
• The cross-sectional area A is constant.
• The level in the vessel can be described by a single variable h(t).

5. Development of a conceptual model. The total volumetric incoming and outgoing flows will
determine the level. The outgoing flow will be related to the pressure drop over the orifice which
has a fixed opening. This pressure drop is proportional to the level in the vessel.

6. Decomposition into subsystems. As the valves are left outside the system boundary, there is no
direct reason to partition the system into subsystems.

7. Formulation of conservation laws. In the problem we have the 7 variables φ1 (t), φ2 (t), c1 (t),
c2 (t), φ(t), c(t) and h(t). Of these 7 variables, 4 are given by the surroundings as inputs. Thus
we need 3 independent equations to describe the remaining 3 variables. We have the following
relations:
mass balance (total):
dρ Ah(t)
= ρφ1 (t) + ρφ2 (t) − ρφ(t) (1.38)
dt
mass balance (component):
dc(t)Ah(t)
= c1 (t)φ1 (t) + c2 (t)φ2 (t) − c(t)φ(t) (1.39)
dt
flow through orifice:
p
φ(t) = k h(t) (1.40)

where k is a constant, determined by the opening of the orifice and by material properties of
the fluid. The three model relations (1.38), (1.39) and (1.40) describe the model. There are two
differential equations and one algebraic equation.

8. Mathematical simplifications. The three relations (1.38), (1.39) and (1.40) can be simplified
by eliminating φ(t) from (1.39) by using (1.40), and by inserting (1.38) into (1.39):

dh(t) p
A = φ1 (t) + φ2 (t) − k h(t) (1.41)
dt
dc(t)
Ah(t) = [c1 (t) − c(t)]φ1 (t) + [c2 (t) − c(t)]φ2 (t) (1.42)
dt
We now have two model equations in the form of two differential equations, with h(t) and c(t)
as unknowns. The equation for the output flow (1.40) can be part of the model in the sense that it
does not describe the internal variables of the system, but, given the internal variables, it can be
considered as defining an output variable, i.e., a variable that is transmitted through the system
1.5. LINEARITY AND LINEARIZATION OF NONLINEAR MODELS 23

inlet flow phi2 outlet concentration

1.5
0.5

0.45
1

kmol/m3
m3/hr

0.4

0.5
0.35

0.3
0
0 10 20 30 0 10 20 30

outlet flow phi level

2.5
1.6
1.4 2
1.2
1.5
m3/hr

m
0.8 1
0.6
0.5
0.4
0.2 0
0 10 20 30 0 10 20 30
time (hrs) −−> time (hrs) −−>

Figure 1.17: Responses of mixing process on inlet flow step

boundary to the surroundings. Simulation results in Figs. 1.17 and 1.18 show typical responses
of the process. The following parameter values are used:
A 1.5 [m2 ] k 1.0
c1∗ 0.8 [kmol/m3 ] φ1∗ 0.25 [m3 /hr]
c2∗ 0.2 [kmol/m3 ] φ2∗ 1.25 [m3 /hr]
Figure 1.17 shows process responses for a large step disturbance in inlet flow φ2 . Note the
difference in concentration response speed for the high level and low level situation. The
responses show nonlinear behaviour, as superposition clearly does not hold. Figure 1.18 shows
small signal behaviour in the close vicinity of a stationary operating condition. The process
variables respond almost as if the process were linear. The stationary initial process condition in
both figures is identical.
9. Mathematical approximation. If desired, the model relations (1.41) and (1.42) can be lin-
earized about a set of stationary operating conditions. Write each variable as
h(t) = h ∗ + 1h(t)
c(t) = c∗ + 1c(t)
and similarly for other variables. Here, the variables h ∗ , c∗ indicate the value of the variable
in the stationary operating condition, and 1h(t), 1c(t) are the deviations from these values.
Inserting these variables into (1.42) yields:
d[c∗ + 1c(t)]
A[h ∗ + 1h(t)] = [c1∗ + 1c1 (t) − c∗ − 1c(t)][φ1∗ + 1φ1 (t)]+
dt
[c2∗ + 1c2 (t) − c∗ − 1c(t)][φ2∗ + 1φ2 (t)] (1.43)
24 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

inlet flow phi2 outlet concentration

1.26 0.3006
b
a
0.3004
1.255
0.3002

kmol/m3
m3/hr

1.25 0.3
0.2998
1.245
0.2996
b
1.24 a
0.2994
0.2992
0 5 10 15 0 5 10 15

outlet flow phi level

1.51 2.28
a a
1.505
2.26
m3/hr

m
1.5
2.24
1.495
b b
1.49 2.22

0 5 10 15 0 5 10 15
time (hrs) −−> time (hrs) −−>

Figure 1.18: Small signal responses of mixing process on inlet flow step

In steady-state the relation is:

0 = [c1∗ − c∗ ]φ1∗ + [c2∗ − c∗ ]φ2∗ (1.44)

Subtracting (1.44) from (1.43):

d1c(t) d1c(t)
Ah ∗ + A1h(t) =
dt dt
[c1∗ − c∗ ]1φ1 (t) + φ1∗ [1c1 (t) − 1c(t)] + [1c1 (t) − 1c(t)]1φ1 (t)+
(1.45)
[c2∗ − c∗ ]1φ2 (t) + φ2∗ [1c2 (t) − 1c(t)] + [1c2 (t) − 1c(t)]1φ2 (t)

Linearization of (1.45) amounts to deleting all terms that contain a product of two 1-variables:

d1c(t)
Ah ∗ = [c1∗ − c∗ ]1φ1 (t) + φ1∗ [1c1 (t) − 1c(t)]+
dt
[c2∗ − c∗ ]1φ2 (t) + φ2∗ [1c2 (t) − 1c(t)] (1.46)

A similar operation can be performed by inserting the perturbed variables into (1.41), yielding:

d[h ∗ + 1h(t)] p
A = φ1∗ + 1φ1 (t) + φ2∗ + 1φ2 (t) − k h ∗ + 1h(t) (1.47)
dt
In the steady state the relation is:
√
0 = φ1∗ + φ2∗ − k h ∗ (1.48)
1.6. DIFFERENTIAL-ALGEBRAIC MODEL EQUATIONS 25

The term containing the square root can be linearized by the following observation:

1h(t)
r
p √
h ∗ + 1h(t) = h ∗ 1 +
h∗ 
√ 1h(t)

= h∗ 1 + + higher order terms in 1h(t) (1.49)
2h ∗
The last step can be verified by computing the square of the left and right hand sides of (1.49).
Inserting (1.48) and (1.49) into (1.47) yields the second linearized model relation:
d1h(t) φ ∗ 1h(t)
A = 1φ1 (t) + 1φ2 (t) − (1.50)
dt 2h ∗
where φ ∗ = φ1∗ + φ2∗ . The linearized version of (1.40) is:
√ 1h(t)
 
φ ∗ + 1φ(t) = k h ∗ 1 + (1.51)
2h ∗
so that
φ∗
1φ(t) = 1h(t) (1.52)
2h ∗
Define the parameter θ as the average residence time of mass particles in the vessel:
h∗ A
θ= (1.53)
φ∗
Then (1.46), (1.50) and (1.52) can be formulated in the standard state space form for linear sets
of differential equations:
1φ1 (t)
 
 " 1 1
#
d 1h(t) 1h(t) 0 0 1φ2 (t)
   1 
− 2θ 0 A A
= + ∗ ∗ ∗
c2 −c∗ φ1∗ φ2∗
dt 1c(t) − θ 1c(t)  1c1 (t) 
1
 
c1 −c
0
Ah ∗ Ah ∗ Ah ∗ Ah ∗
1c2 (t)
1φ(t) 0 1h(t)
  A  
= 2θ (1.54)
1c(t) 0 1 1c(t)
The model is in the standard form
d
x(t) = Ax(t) + Bu(t)
dt
y(t) = C x(t) (1.55)
where x(t) is the state vector, u(t) is the input vector, and y(t) is the output vector.

1.6 Differential-algebraic model equations

The model (1.23) is not the most general one encountered in dynamic modelling. In fact, in general we
have coupled differential equations and algebraic equations of the form
d
x1 (t) = f (x1 (t), x2 (t), u(t))
dt
0 = g(x1 (t), x2 (t), u(t)) (1.56)
26 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

The algebraic equations (1.56) are in implicit form. This may be experienced by the occurrence of
algebraic loops, i.e., implicit algebraic relations amongst the variables in x1 , x2 and u. We can try to
replace (1.56) by an equivalent explicit set of equations in which the variables x2 (t) are specified as
explicit functions of x1 (t) and u(t):

x2 (t) = h(x1 (t), u(t)) (1.57)

If indeed the model (1.56) can be brought into the form (1.57), we can insert (1.57) into the differential
equations in (1.56) and obtain a model in state-space form:

d
x1 (t) = f (x1 (t), h(x1 (t), u(t)), u(t)) (1.58)
dt
or

d
x1 (t) = f ∗ (x1 (t), u(t)) (1.59)
dt
The structure of these relations is shown in Fig. 1.19. Note that all signals pass through the integrator
in their route from inputs to outputs.

f∗
x1 (0)
u(t)
u(t) dx1 (t)
x2 (t) x1 (t)
Z
f dt
h

x1 (t)

Figure 1.19: Differential-algebraic equations, transformed to state-space form

For smooth functions f, g in (1.56), the solutions to the linearized model are contained in the set of
solutions to the nonlinear model, be it in an approximate sense and under the conditions of the validity
of these approximations. For the approach shown above for bringing a model formulated in terms of
algebraic and differential equations into an explicit state-space form, further insight into the existence
question of this model transformation can be derived from a linearized analysis.
A linearized version of f and g with respect to x1 , x2 and u about a stationary solution x1∗ , x2∗ , u ∗
leads to the linearized model
d
1x1 (t) = A11 1x1 (t) + A12 1x2 (t) + B1 1u(t)
dt
0 = A21 1x1 (t) + A22 1x2 (t) + B2 1u(t) (1.60)

where
∂f ∂f ∂f
A11 = ∂ x1 ∗
A12 = ∂ x2 ∗
B1 = ∂u ∗
(1.61)
∂g ∂g ∂g
A21 = ∂ x1 ∗
A22 = ∂ x2 ∗
B2 = ∂u ∗
1.7. CONTINUOUS CHEMICAL REACTOR PROCESS DYNAMICS 27

Now a sufficient condition for the existence of an explicit expression that is equivalent to the algebraic
equations (1.60) is that the matrix A22 is nonsingular, i.e., it is invertible. In that case, the explicit
expression is

1x2 (t) = −A−1

22 A21 1x 1 (t) − A22 B2 1u(t)
−1
(1.62)

Inserting this explicit expression for 1x2 (t) into the differential equations in (1.60) leads to the linear
state-space model
d
1x1 (t) = A1x1 (t) + B1u(t) 1x1 (t0 ) = x10 (1.63)
dt
with

A = A11 − A12 A−1

22 A21 (1.64)
B = B1 − A12 A−1
22 B2 (1.65)

If A22 turns out to be singular, then a further analysis is necessary in order to assess the properties of
the system under consideration. In that case we may have to investigate the index of the system, or
it may turn out that our system model is internally or externally nonproper. The discussion of these
issues will be given in a following chapter.
The perturbation variables 1x, 1u will often be replaced by x, u if the context of a linear model is
clear. Then (1.63) is written as
d
x1 (t) = Ax1 (t) + Bu(t) x1 (t0 ) = x10 (1.66)
dt
Assuming that A, B are time-invariant, the solution is:
Z t
A(t−t0 )
x1 (t) = e x10 + e A(t−τ ) Bu(τ ) dτ (1.67)
t0

where the matrix exponential is defined by its series expansion

1 2 2
e At = I + At + A t + ··· (1.68)
2!

1.7 Continuous chemical reactor process dynamics

As an example of a more involved process modelling exercise, a continuous flow system will be
considered in which a nonisothermal chemical reaction takes place. The reaction is assumed to be
exothermic and irreversible. Reactant A forms the product B at a specific reaction rate k, which is
k
indicated by A −→ B. The reaction vessel will be assumed to be well mixed, and such a process is
known under the name CSTR (Continuous Stirred-Tank Reactor). There is a continuous input feed
flow to the process, and a continuous outflow. Density variations are assumed to be neglected, and the
vessel is completely filled. Thus input volumetric flow φ(t) is equal to the output volumetric flow. The
flow can vary in time. The vessel is surrounded by a cooling jacket to remove the heat of reaction. The
amount of cooling water is controlled by a valve, the opening of which is determined by a PI-controller.
The reaction temperature in the vessel is the controlled variable. The process is shown in Fig. 1.20. The
various decisions underlying the modelling procedure closely follow the scheme given in Table 1.1:
28 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

Tset (t)
feed flow
control
TC
Ci (t) Ti (t)

cooling water out

T j (t) C(t)
X (t)
φ j (t) T (t)
jacket
T jo (t)
cooling water in
vessel
product flow

C(t) T (t)

Figure 1.20: Exothermal chemical reactor

1. Intended purpose of the model: To describe the dynamics of the process such that conclusions
regarding stability of operation can be drawn.

2. System boundary: The exothermic reaction needs cooling, and the amount of cooling must
be adapted to the actual heat production if constant operational conditions are required in the
face of disturbances. Thus the system consists of vessel, cooling jacket, cooling flow valve and
temperature controller. The amount of flow through the vessel is assumed to be determined
outside the system.

3. Variables across system boundary: The following input variables must be considered:

• volumetric flow φ(t) [m3 /s] through the vessel

• temperature Ti (t) [K] of the input flow to the vessel
• concentration Ci (t) [mol/m3 ] of component A in the input flow
• coolant flow inlet temperature T jo (t) [K]
• setpoint value Tset [K] for the temperature controller.

As internal variables we have:

• volumetric flow φ j (t) [m3 /s] through the jacket

• temperature T j (t) [K] of the jacket
• concentration C(t) [mol/m3 ] of component A in the vessel
• temperature T (t) [K] of the fluid in the vessel.

This list must be completed on the basis of a further study of the model hypotheses and conser-
vation laws.
1.7. CONTINUOUS CHEMICAL REACTOR PROCESS DYNAMICS 29

4. Model hypotheses:

• The vessel is perfectly mixed

• There are no heat losses to the environment from coolant flow, jacket and vessel
• The jacket is assumed to be perfectly mixed
• The heat capacity of the vessel wall between vessel and jacket is neglected; its heat
resistance is assumed to be incorporated in the heat transfer coefficients between coolant
and wall, and between wall and fluid in the vessel
• The heat transfer between coolant and wall, and between wall and fluid in the vessel is
assumed to be proportional to the difference in temperatures
• The vessel fluid volume V is constant
• The density of the fluid ρ [kg/m3 ] is independent of temperature T (t) and concentra-
tion C(t), so that input volumetric flow φi (t) [m3 /s] instantaneously equals output volu-
metric flow φ(t) [m3 /s]. The input flow is assumed to be determined outside the system.
• All material properties of the fluid in the vessel are constant, independent of concentra-
tion C(t)
• The pressure drop over the coolant valve is assumed to be constant.

5. Subdivision into subsystems: The four components of the process determine four subsystems
which will be discussed regarding variables and interconnection structure. The variables and
interconnection structure are shown in the block diagram of Fig. 1.21.

T jo (t) Ti (t) Ci (t) φ(t)

T j (t) ? ? ? C(t)
Tset (t) - - -
- jacket reactor
controller - valve -  -
- X (t) φ j (t) Q(t) T (t)

Figure 1.21: Block diagram of CSTR subsystems and interconnections

(a) Reactor vessel

• input variables: Ci (t), Ti (t), φ(t), T j (t)
• internal variables: C(t), T (t), heat flow Q(t) [J/s] between vessel and jacket, rate of
reaction R(t) [mol/m3 s]
Thus 4 equations required. The variables C(t), T (t) and Q(t) act as output variables.
(b) Jacket
• input variables: T jo (t), φ j (t), Q(t)
• internal variable: T j (t)
30 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

Thus 1 equation required. The variable T j (t) acts as output variable.

(c) Temperature controller
• input variables: Tset (t), T (t)
• internal variables: Temperature error (t) [K], valve position X (t) [m]
Thus 2 equations required. The variable X (t) acts as output variable.
(d) Coolant valve
• input variable: X (t)
• internal variable: φ j (t)
Thus 1 equation required. The variable φ j (t) acts as output variable.

Note that the jacket and reactor are thermally bilaterally coupled.

6. Conservation laws and other model relations: These will be discussed for each subsystem
separately.

(a) Reactor vessel

Mass balance for component A:

d
{V C(t)} = φ(t)[Ci (t) − C(t)] − V R(t) (1.69)
dt
Energy balance for the vessel

d
Cp {ρV T (t)} = φ(t)ρC p [Ti (t) − T (t)] − Q(t) + (−1i)V R(t) (1.70)
dt
Here −1i [J/mol] is the heat of reaction, C p [J/kgK] is the specific heat. The equation for
the heat flow is:

Q(t) = α Ac [T (t) − T j (t)] (1.71)

where α [J/m2 Ks] is the total heat transfer coefficient between vessel and coolant in the
jacket, and Ac [m2 ] is the total area of heat transfer between vessel fluid and jacket.
The reaction rate for an n th -order reaction can be described by
E
R = C n ke− RT (1.72)

where k is the frequency factor, E [J/kg] is the activation energy of the reaction, and R [J/kgK]
is the specific gas constant. In our case where one species A reacts with a component as-
sumed to be available in excess, the reaction is of first order, i.e., n = 1.
(b) Jacket
Energy balance for the jacket volume:

d
C pj {ρ j V j T j (t)} = φ j (t)ρ j C pj [T jo (t) − T j (t)] + Q(t) (1.73)
dt
where V j [m3 ] is the jacket volume, C pj [J/kgK] is the specific heat of the coolant, ρ j
[kg/m3 ] is the coolant density.
1.7. CONTINUOUS CHEMICAL REACTOR PROCESS DYNAMICS 31

(c) Temperature controller

Suppose that this controller has the characteristics of a PI-controller:
(t) = Tset (t) − T (t) (1.74)
Z t
X (t) = K p (t) + K i (τ ) dτ (1.75)
−∞

Here, K p and K i are the controller parameters for proportional action and for integral
action, respectively.
(d) Coolant valve
If we assume a linear relationship between valve position and flow area, we have:

0
 if X < 0
φ j (t) = K v X (t) if 0 ≤ X ≤ X max (1.76)
K v X max if X > X max



where X max is the value for X (t) for which the valve is maximally open, and K v is a
normalization constant allowing to take X max = 1.
7. Model structure investigation: First the structure of the set of differential and algebraic equa-
tions comprising the model will be analysed. The open-loop model of the reactor and jacket
(without controller and coolant valve) will be considered. The dynamic balance equations contain
a time derivative of the accumulation terms in the left-hand side. These terms will be brought to
a form containing only time derivatives of individual variables. These variables will turn out to
become the state variables in the nonlinear dynamic system.
d φ(t)
C(t) = [Ci (t) − C(t)] − R(t) (1.77)
dt V
d φ(t) Q(t) (−1i)
T (t) = [Ti (t) − T (t)] − + R(t) (1.78)
dt V ρC p V ρC p
d φ j (t) Q(t)
T j (t) = [T jo (t) − T j (t)] + (1.79)
dt Vj ρ j C pj V j
0 = Q(t) − α Ac [T (t) − T j (t)] (1.80)
E
0 = R(t) − C(t)ke− RT (t) (1.81)
Equations (1.77)–(1.81) are in the form of the general set of differential-algebraic equations
discussed earlier in Section 1.5, equations (1.56):
d
x1 (t) = f (x1 (t), x2 (t), u(t))
dt
0 = g(x1 (t), x2 (t), u(t)) (1.82)
In the present case, the vectors x1 (t), x2 (t) and u(t) are:
T jo (t)
 
 φ j (t) 
 
C(t)  
Q(t)
x1 (t) = T (t)  x2 (t) = (t)
 
 u(t) =  T i
 (1.83)
R(t)
T j (t)  Ci (t) 
 

φ(t)
32 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

The set of equations (1.80), (1.81) determining x2 (t) can be written explicitly as

Q(t) = α Ac [T (t) − T j (t)]

E
R(t) = C(t)ke− RT (t) (1.84)

and inserted in (1.77)–(1.79), leading to the state space equations:

d φ(t) E
C(t) = [Ci (t) − C(t)] − C(t)ke− RT (t) (1.85)
dt V
d φ(t) α Ac [T (t) − T j (t)] (−1i) E
T (t) = [Ti (t) − T (t)] − + C(t)ke− RT (t) (1.86)
dt V ρC p V ρC p
d φ j (t) α Ac [T (t) − T j (t)]
T j (t) = [T jo (t) − T j (t)] + (1.87)
dt Vj ρ j C pj V j
The elimination of two internal variables using the two algebraic equations of the model exactly
corresponds to the elimination procedure described in Section 1.5, leading to (1.58), (1.59).

8. Steady-state nonlinearity and stability analysis: The state-space model (1.85)–(1.87) is non-
linear for several reasons. One is the occurrence of products of input and state variables in
the right-hand sides of the state space equations, e.g., the product φ(t)C(t) in (1.85). Another
reason is the appearance of the nonlinear expression for the reaction rate. The character of this
nonlinearity will be investigated by considering the steady-state behaviour of the reactor vessel.
The reactor vessel has as inputs the variables Ci (t), Ti (t), φ(t) and T j (t). In steady state (1.85)
and (1.86) become
E
φ ∗ [Ci∗ − C ∗ ] − V C ∗ ke− RT ∗ = 0 (1.88)
E
φ ∗ ρC p [Ti∗ − T ∗ ] − α Ac [T ∗ − T j∗ ] + V (−1i)C ∗ ke− RT ∗ = 0 (1.89)
where the variables (·)∗ indicate steady-state values. From (1.88) follows an explicit expression
for C ∗ :
φ ∗ Ci∗
C∗ = E
(1.90)
φ ∗ + V ke− RT ∗
The left-hand side of the energy balance (1.89) can be split into two parts; one part describes the
heat generation G(T ∗ ), the other part is the heat removal from the vessel −F(T ∗ ):
E
G(T ∗ ) := AH (−1i)C ∗ ke− RT ∗ (1.91)

Inserting (1.90):
E
φ ∗ Ci∗ (−1i)ke− RT ∗
G(T ) :=
∗
(1.92)
φ∗ E
AH
+ ke− RT ∗

For given input variables φ ∗ , Ci∗ and given material properties and reactor volume this indeed is
only a function of T ∗ . The function G(T ∗ ) is S-shaped and is shown in Fig. 1.22 as a function
of T ∗ . The other terms in (1.89) define the heat removal −F(T ∗ ):

−F(T ∗ ) = φ ∗ ρC p [T ∗ − Ti∗ ] + α Ac [T ∗ − T j∗ ] (1.93)

1.7. CONTINUOUS CHEMICAL REACTOR PROCESS DYNAMICS 33

−F

G, -F

T∗

Figure 1.22: Heat generation G(T ∗ ) and removal −F(T ∗ ) in CSTR

For given input variables φ ∗ , Ti∗ , T j∗ this is a linear function of T ∗ :

−F(T ∗ ) = [φ ∗ ρC p + α Ac ]T ∗ − φ ∗ ρC p Ti∗ − α Ac T j∗ (1.94)

The function −F(T ∗ ) is also shown in Fig. 1.22 as a function of T ∗ . In steady state the heat
generation equals the heat removal. The figure shows three possible steady state equilibrium
solutions: one for low temperature, one for intermediate temperature and one for high temper-
ature. Suppose the process is in the intermediate point of intersection. Let a small increase in
heat generation lead to a slightly higher temperature, i.e., to a temperature to the right of the
intersection point. Then the heat removal is smaller than the heat generation, i.e., the process
tends to be removed away from the intersection conditions. The steady state operational condition
in the intermediate intersection point is called statically unstable, the other intersections are
statically stable. The following remarks apply.

• A dynamic analysis shows that the process operated in the intermediate intersection point
is also asymptotically unstable from a dynamic point of view, i.e., regarding its natural
behaviour under disturbances as a function of time.

• The stability result holds for a fixed jacket temperature T j∗ . If the jacket temperature is not
fixed, but varies as a function of T ∗ , a different heat removal function results, exhibiting
different equilibrium and stability properties.

• In particular if stabilizing proportional temperature feedback control is applied, then the

actual steady-state heat removal function −F(T ∗ ) can be given a greater slope. If stabilizing
integral feedback control of the temperature T ∗ is applied, then the slope of the heat removal
function −F(T ∗ ) will be infinite, so that only one equilibrium intersection point occurs.

9. Dynamic simulation of the CSTR: A dynamic simulation of the complete process shows the
behaviour of the process under disturbances of the various inputs, for different initial steady-state
34 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

operating conditions. The following numerical values are used.

φ 1.13 [m3 /hr] k 7.08 × 1010 [hr−1 ]

AH 1.36 [m3 ] E 6.96 × 107 [J/kmol]
Ci 8.0 [kmol/m3 ] R 8314.3 [J/kmolK]
Ti 294 [K] α 3.07 × 106 [J/m2 hrK]
T jo 294 [K] Ac 23.2 [m2 ]
(−1i) 6.96 × 107 [J/kmol] Vj 0.109 [m3 ]
ρ 801 [kg/m3 ] ρj 998 [kg/m3 ]
Cp 3140.1 [J/kgK] Cj 4186.8 [J/kgK]
Kp −0.5 X max 1
Ki −5 Kv 2.5

The results are shown in the Figs. 1.23, . . . , 1.29.

Steady state reactor temperature vs coolant flow

420

400

380
T (degrees)

360

340

320

300

0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5 1.6
coolant flow (m3/hr)

Figure 1.23: Steady state process operational curve of CSTR

1.7. CONTINUOUS CHEMICAL REACTOR PROCESS DYNAMICS 35

setpoint Tset concentration C

4.4
335
4.2

kmol/m3
degrees

330 4

3.8
325
3.6
0 2 4 6 8 0 2 4 6 8
temperature T coolant flow

335 2
degrees

m3/hr
330 1

325 0
0 2 4 6 8 0 2 4 6 8
coolant temperature Tj reaction rate script−R

5
335
kmol/m3hr
degrees

4
330
3

325
2
0 2 4 6 8 0 2 4 6 8
time (hrs) time (hrs)

Figure 1.24: Process response on step in Tset under closed-loop control

coolant flow coolant flow

0.75 0.65

a a
0.7 0.645

0.64
m3/hr

m3/hr

0.65

0.6 0.635

b 0.63 b
0.55

0 1 2 3 0 1 2 3

temperature T temperature T

330 312

325 b 311.5
degrees

degrees

320 311 b

315 310.5

310 310
a
a
305 309.5
0 1 2 3 0 1 2 3
time (hrs) time (hrs)

Figure 1.25: Responses on step in coolant flow, open loop T ∗ = 310 K

36 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

coolant flow coolant flow

1.5 a 1.418 a

1.45 1.4175
m3/hr

m3/hr
1.4 1.417

1.35 1.4165
b b

1.3 1.416
0 1 2 3 4 5 0 1 2 3 4 5

temperature T temperature T
360.4
370 b
360
360.2 b
350
degrees

degrees
340
360
330
320
359.8 a
310
a
300
359.6
0 1 2 3 4 5 0 1 2 3 4 5
time (hrs) time (hrs)

Figure 1.26: Responses on step in coolant flow, open loop T ∗ = 360 K

coolant flow coolant flow coolant flow

0.6385 0.9365
0.8935

0.638 0.936
0.893
m3/hr

0.6375 0.9355
0.8925

0.637 0.935
0.892

0.6365 0.9345
0 5 10 0 5 10 0 5 10

temperature T temperature T temperature T

300.3 310.3 320.3

300.2 310.2 320.2

300.1 310.1 320.1

degrees

300 310 320

299.9 309.9 319.9

299.8 309.8 319.8

0 5 10 0 5 10 0 5 10
time (hrs) time (hrs) time (hrs)

Figure 1.27: Responses on step in coolant flow, open loop T ∗ = 300, 310, 320 K
1.7. CONTINUOUS CHEMICAL REACTOR PROCESS DYNAMICS 37

coolant flow coolant flow coolant flow

1.538
1.296
1.518
1.5375
1.2955
1.5175
m3/hr

1.537
1.295
1.517
1.5365
1.2945
1.5165
1.536
0 5 10 0 5 10 0 5 10

temperature T temperature T temperature T

330.3 340.3 350.3

330.2 340.2 350.2

330.1 340.1 350.1

degrees

330 340 350

329.9 339.9 349.9

329.8 339.8 349.8

0 5 10 0 5 10 0 5 10
time (hrs) time (hrs) time (hrs)

Figure 1.28: Responses on step in coolant flow, open loop T ∗ = 330, 340, 350 K

coolant flow coolant flow coolant flow

1.2415
1.0605
1.4175

1.241
1.06
1.417
m3/hr

1.2405
1.0595
1.4165

1.24
1.059
1.416

1.2395
0 5 10 0 5 10 0 5 10

temperature T temperature T temperature T

360.3 370.3 380.3

360.2 370.2 380.2

360.1 370.1 380.1

degrees

360 370 380

359.9 369.9 379.9

359.8 369.8 379.8

0 5 10 0 5 10 0 5 10
time (hrs) time (hrs) time (hrs)

Figure 1.29: Responses on step in coolant flow, open loop T ∗ = 360, 370, 380 K
38 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

10. Relaxing the ideally mixing assumptions: The simulation model has been extended by as-
suming a subdivision into two compartments for both the reaction vessel and the jacket. Each
compartment is assumed to be ideally mixed, and both compartments are connected in series
without backflow. The heat exchanging area is subdivided into four equal parts, according to the
block scheme in Fig. 1.30. The model equations have been formulated for each compartment

Ti , Ci , φ

T j1

Q1
T1 , C1 Q2

T j1 T j2
Q3
T jo T2 , C2 T j2
Q4

T2 , C2

Figure 1.30: Block diagram of 2 × 2 compartments CSTR model

separately, taking into account the partial volumes and heat exchange areas. The steady state
process operational curve for the compartmental model is compared with the original model in
Fig. 1.31. Relaxing the ideally mixing assumption has some influence on quantitative steady state
and dynamic characteristics, although the qualitative behaviour is retained. Figure 1.32 shows
the responses for the compartmental model under closed-loop control with step disturbance on
temperature setpoint. The responses show considerable agreement in their form with those of
Fig. 1.24, although the steady state values differ.
1.7. CONTINUOUS CHEMICAL REACTOR PROCESS DYNAMICS 39

Steady state reactor temperature vs coolant flow

420

400

380 2 x 2 compartments
T (degrees)

360
perfect mixing

340

320

300

0.6 0.8 1 1.2 1.4 1.6

coolant flow (m3/hr)

Figure 1.31: Steady state process operational curve of compartmental CSTR

setpoint Tset concentration C1

335 5.6
kmol/m3
degrees

5.4
330
5.2
325
5
0 2 4 6 8 0 2 4 6 8
temperature T1 coolant flow

335 2
degrees

m3/hr

330
1
325
0
0 2 4 6 8 0 2 4 6 8
coolant temperature Tj1 reaction rate script−R1

335 6
kmol/m3hr
degrees

330 5

4
325
3
0 2 4 6 8 0 2 4 6 8
time (hrs) time (hrs)

Figure 1.32: Compartmental CSTR response on step in Tset , closed-loop control

40 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

1.8 Literature for this Chapter

Basic technical and philosophical ideas of modelling as given in Section 1.1 can be found in Kalman
[27], Kwatny and Mablekos [29], [30], Golomb [19], Profos [42], or Bub and Lugner [10]. A system
theory that defines physical models without a priori distinction between inputs and outputs has been
proposed by Willems [52]; see also the earlier work [51], and the full treatment in Polderman and
Willems [39].
The formulation of conservation laws for the general microscopic and macroscopic situations and
for a variety of boundary conditions and coordinate frame descriptions are discussed in Himmelblau
and Bischoff [23] and in Bird et al. [9]. The two volumes (in Dutch) on transport phenomena by Van
den Akker and Mudde [46] and by Hoogendoorn and Van der Meer [25] provide background on the
formulation of basic relations for process models.
The details of the modelling of process heat and mass transfer can be found in books like Welty, Wicks
and Wilson [49] or Hewitt, Shires and Bott [22].
Introductory treatments of process modelling are Guy [20], Roffel and Rijnsdorp [44], and Kecman [28].
The role of dynamic modelling in control engineering is considered in Profos [43]. Ljung and Glad [31]
consider some elementary dynamic modelling issues in their relationship with system identification.
Modelling from a fundamental chemical engineering point of view is considered in Eigenberger
[14], Marquardt [34], and Marquardt and Zeitz [36]. Systematic approaches for model development
can be found in Lohmann and Marquardt [32] and Marquardt [35]. Specifics of phase equilibrium
process modelling are addressed in Ponton and Gawthrop [40], and an approach towards hybrid models
is provided in Barton and Pantelides [8]. A systematic way to formulate process models has been
developed in Preisig, Kimmich and Rippin [41], Weiss and Preisig [48], and Westerweele [50].
The basic techniques of chemical process modelling are treated in Hartmann and Kaplick [21], Luyben
and Wenzel [33], and in Ingham, Dunn and Heinzle [26]. Process modelling for process control is
treated in Ogunnaike and Ray [37]. Classical references on process modelling are Franks [15], [16],
Friedly [17], Aris [1], [2], and Denn [12], [13].
Further discussions can be found in Paynter [38], Asbjornsen [5], [6], Sargent [45], and Georgakis [18].
The continuous stirred tank reactor has been studied in the classical work of Aris and Amundson [3],
[4]. The multiplicity of steady states has been studied by many researchers, e.g., Hlavacek, Kubicek
and Visnak [24], Varma and Amundson [47], Balakotaiah and Luss [7], and Cicarelli and Aris [11].
1.8. LITERATURE FOR THIS CHAPTER 41

References
[1] R. Aris. Mathematical Modelling Techniques. Pitman Publishing Ltd., London, UK, 1978.

[2] R. Aris. Method in the modeling of chemical engineering systems. In C. T. Leondes, editor,
Control and Dynamic Systems, volume 15, pages 41–98. Academic Press, New York, NY, 1979.

[3] R. Aris and N. R. Amundson. An analysis of chemical reactor stability and control. I. The
possibility of local control, with perfect or imperfect control mechanisms. Chemical Engineering
Science, 7:121–131, 1958a.

[4] R. Aris and N. R. Amundson. An analysis of chemical reactor stability and control. II. The
evolution of proportional control. Chemical Engineering Science, 7:132–147, 1958b.

[5] O. A. Asbjornsen. Modeling techniques: theory and practice. Modeling Identification and Control,
6:105–125, 1985.

[6] O. A. Asbjornsen. A systems engineering approach to process modeling. In D. M. Prett and

M. Morari, editors, The Shell Process Control Workshop. Process Control Research: Industrial and
Academic Perspectives, Houston TX, dec. 15-18, 1986, pages 139–182. Butterworth Publ., Boston, MA,
1987.

[7] V. Balakotaiah and D. Luss. Analysis of the multiplicity patterns of a CSTR. Chemical Engineering
Communications, 13:111–132, 1981.

[8] P. I. Barton and C. C. Pantelides. Modeling of combined discrete/continuous processes. AIChE

Journal, 40:966–979, 1994.

[9] R. B. Bird, W. E. Stewart, and E. N. Lightfoot. Transport Phenomena. John Wiley and Sons, Inc.,
New York, NY, 1960.

[10] W. Bub and P. Lugner. Systematik der Modellbildung. Teil 1: Konzeptionelle Modellbildung.
In I. Troch, editor, Modellbildung für Regelung und Simulation: Methoden-Werkzeuge-Fallstudien,
GMA Aussprachetag, Langen, BRD, 25-26 März 1992, VDI Ber. Nr. 925, pages 1–18. VDI Verlag,
Düsseldorf, 1992.

[11] P. Cicarelli and R. Aris. Continuous reactions in a non-isothermal CSTR- I. Multiplicity of steady
states. Chemical Engineering Science, 49:621–631, 1994.

[12] M. M. Denn. Modeling for process control. In C. T. Leondes, editor, Control and Dynamic
Systems, volume 15, pages 147–194. Academic Press, New York, NY, 1979.

[13] M. M. Denn. Process Modeling. Longman Scientific and Technical, Harlow, Essex, UK, 1986.

[14] G. Eigenberger. Modelling and simulation in industrial chemical reaction engineering. In

K. H. Ebert and P. Deuflhard, editors, Modelling of Chemical Reaction Systems. Proc. Int. Workshop,
Heidelberg, FRG, September 1-5, 1980, pages 284–304. Springer Verlag, Berlin, 1981.

[15] R. G. E. Franks. Mathematical Modeling in Chemical Engineering. John Wiley and Sons, Inc.,
New York, NY, 1967.

[16] R. G. E. Franks. Modeling and Simulation in Chemical Engineering. John Wiley and Sons, Inc.,
New York, NY, 1972.
42 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

[17] J. C. Friedly. Dynamic Behavior of Processes. Prentice Hall, Inc., Englewood Cliffs, NJ, 1972.

[18] C. Georgakis. On the use of extensive variables in process dynamics and control. Chemical
Engineering Science, 41:1471–1484, 1986.

[19] S. W. Golomb. Mathematical models: uses and limitations. IEEE Transactions on Reliability,
20:130–131, 1971.

[20] J. L. Guy. New CE Refresher: Process Dynamics. McGraw-Hill Book Co., New York, NY, 1983.

[21] K. Hartmann and K. Kaplick. Analysis and Synthesis of Chemical Process Systems. Elsevier
Publ., Amsterdam, The Netherlands, 1990.

[22] G. F. Hewitt, G. L. Shires, and T. R. Bott. Process Heat Transfer. CRC Press, Boca Raton, FL,
USA, 1994.

[23] D. M. Himmelblau and K. B. Bischoff. Process Analysis and Simulation. Deterministic Systems.
John Wiley and Sons, Inc., New York, NY, 1968.

[24] V. Hlavacek, M. Kubicek, and K. Visnak. Modelling of chemical reactors XXVI. Multiplicity
and stability analysis of a continuous stirred tank reactor with exothermic consecutive reactions
A → B → C. Chemical Engineering Science, 27:719–742, 1972.

[25] C. J. Hoogendoorn and T. H. Van der Meer. Fysische Transportverschijnselen II. Delftse Uitgevers
Maatschappij, Delft, Nederland, 1991.

[26] J. Ingham, I. J. Dunn, and E. Heinzle. Chemical Engineering Dynamics. Modelling with PC
Simulation. VCH Verlagsgesellschaft, Weinheim, W.Germany, 1994.

[27] R. E. Kalman. Comments on the scientific aspects of modeling. In M. Marois, editor, Towards a
Plan of Actions for Mankind, pages 493–505. North-Holland Publishing Company, Amsterdam, 1974.

[28] V. Kecman. State-Space Models of Lumped and Distributed Systems, volume 112 of Lecture
Notes in Control and Information Sciences. Springer Verlag, Berlin, 1988.

[29] H. G. Kwatny and V. E. Mablekos. The modeling of dynamical processes. In Proceedings of the
1975 IEEE Conference on Decision and Control, Houston, Texas, USA, 10-12 December 1975, pages
271–281. IEEE, New York, NY, 1975.

[30] H. G. Kwatny and V. E. Mablekos. On the modeling of dynamical processes with applications in
power systems engineering. In L. H. Fink and K. Carlsen, editors, Systems Engineering for Power:
Status and Prospects, Engineering Foundation Conference, Henniker, NH, USA, August 17-22, 1975,
pages 368–392. NTIS, Springfield, 1975.

[31] L. Ljung and T. Glad. Modeling of Dynamic Systems. Prentice Hall, Inc., Englewood Cliffs, NJ,
1994.

[32] B. Lohmann and W. Marquardt. On the systematization of the process of model development.
Computers and Chemical Engineering, 20:S213–S218, 1996.

[33] W. L. Luyben and L. A. Wenzel. Chemical Process Analysis: Mass and Energy Balances. Prentice
Hall, Inc., Englewood Cliffs, NJ, 1988.
1.8. LITERATURE FOR THIS CHAPTER 43

[34] W. Marquardt. Dynamic process simulation: Recent progress and future challenges. In Y. Arkun
and W. H. Ray, editors, Chemical Process Control - CPC IV, Proceedings of the 4th Engineering
Foundations Conf., South Padre Island, TX, USA, feb. 17-22, 1991, pages 131–180. American Institute
of Chemical Engineers, New York, 1991.

[35] W. Marquardt. Towards a process modeling methodology. In R. Berber, editor, Methods of Model
Based Process Control, Proceedings of the NATO Advanced Study Institute, Antalya, Turkey, august
7-17, 1994, pages 3–40. Kluwer Publ., Dordrecht, the Netherlands, 1995.

[36] W. Marquardt and M. Zeitz. Rechnergestützte Modellbildung in der Verfahrenstechnik. In

I. Troch, editor, Modellbildung für Regelung und Simulation: Methoden-Werkzeuge-Fallstudien, GMA
Aussprachetag, Langen, BRD, 25-26 März 1992, VDI Ber. Nr. 925, pages 307–341. VDI Verlag,
Düsseldorf, 1992.

[37] B. A. Ogunnaike and W. H. Ray. Process Dynamics, Modeling, and Control. Oxford University
Press, Inc., New York, NY, USA, 1994.

[38] H. M. Paynter. Thermofluid process dynamics in boiler modeling. In D. A. Berkowitz, editor,

Proceedings of the Seminar on Boiler Modeling, Bedford, MA, USA, nov. 6-7, 1974, pages 187–206.
The MITRE Corporation, Bedford, MA, USA, 1975.

[39] J. W. Polderman and J. C. Willems. Introduction to Mathematical Systems Theory. A Behavioural

Approach. Springer Verlag, New York, NY, 1998.

[40] J. W. Ponton and P. J. Gawthrop. Systematic construction of dynamic models for phase equilibrium
processes. Computers and Chemical Engineering, 15:803–808, 1991.

[41] H. A. Preisig, M. Kimmich, and D. W. T. Rippin. A study of dynamic system modelling.

Computers and Chemical Engineering, 12:455–460, 1988.

[42] P. Profos. Some considerations on theoretical process modeling. Trans. ASME / Journal of
Dynamic Systems, Measurement, and Control, 94:282–284, 1972.

[43] P. Profos. Modellbildung und ihre Bedeutung in der Regelungstechnik. In R. Isermann, editor,
Prozessmodelle 1977. Modellbildung und Identifikation technischer Prozesse, Tagung Wiesbaden, 25
und 26 April 1977, VDI-Berichte Nr. 276, pages 5–12. VDI Verlag, Düsseldorf, 1977.

[44] B. Roffel and J. E. Rijnsdorp. Process Dynamics, Control, and Protection. Ann Arbor Science
Publishers, Ann Arbor, MI, USA, 1982.

[45] R. W. H. Sargent. Advances in modelling and analysis of chemical process systems. Computers
and Chemical Engineering, 7:219–237, 1983/84.

[46] H. E. A. Van den Akker and R. F. Mudde. Fysische Transportverschijnselen I. Delftse Universitaire
Pers, Delft, Nederland, 1996.

[47] A. Varma and N. R. Amundson. Some observations on uniqueness and multiplicity of steady
states in non-adiabatic chemically reacting systems. Canadian Journal of Chemical Engineering,
51:206–226, 1973.

[48] M. Weiss and H. A. Preisig. Structural analysis in the dynamical modelling of chemical engineer-
ing systems. Mathematical and Computer Modelling of Dynamical Systems, 6:325–364, 2000.
44 CHAPTER 1. BASIC MODELLING OF DYNAMIC ENGINEERING SYSTEMS

[49] J. R. Welty, C. E. Wicks, and R. E. Wilson. Fundamentals of Momentum, Heat, and Mass Transfer.
3rd Edition. John Wiley and Sons, Inc., New York, NY, 1984.

[50] M. R. Westerweele. Five steps for building consistent dynamic process models and their imple-
mentation in the computer tool Modeller. PhD thesis, Technische Universiteit Eindhoven, 2003.

[51] J. C. Willems. System theoretic models for the analysis of physical systems. Ricerche di
Automatica, 10:71–106, 1979.

[52] J. C. Willems. Paradigms and puzzles in the theory of dynamical systems. IEEE Transactions on
Automatic Control, 36:259–294, 1991.