Chapter 13

Dynamic Simulation
of Control Systems

Dynamic simulation is a most convenient method for analyzing the performance of

process control systems. This chapter presents the basics for leaming how to develop
and use dynamic simulations of processes and instruments-controllers, valves, sen-
sors, etc. Most of the simulations presented here will be based on the dynamic modeling
technique presented in Chapters 3 and 4, as well as the instrument models presented
in Chapter 5.
Two methods of dynarnic simulation are presented:

• Simulation of linear systems through transfer function representation.

• Simulation of control systems from basic mecbanistic models.

The first of these methods has the advantage of not requiring concem about the initial
conditions because, as indicated in Chapter 2, transfer functions assume the variables
are deviations from steady-srate initial conditions, that is, the initial conditions are
zero. However, it requires the development of the transfer functions and, in the case
of nonlinear systems, the linearization of the model equations.
The second method is more fundamental and is most suitable for the modeling of
nonlinear systems. The major difficulty is that the modeler must be very concerned
with using the correct initial conditions of ali the variables. When simulating control
systems for continuous processes it is also necessary to ensure that the initial conditions
represen! a steady state.
The bulk of the cbapter assumes that the processes are lumped, that is, the variables
are functions only of time. A brief discussion of the simulation of stiff systems is
included at the end.
The examples presented use the program MATLAB with Simulink, because it
is one of the most powerful block-oriented simulation prograrns that is commonly
available to students and engineers. Simulations using VisSim are also included in
the web page that accompanies this book.

13-1 USES AND TOOLS OF DYNAMIC SIMULATION

This section presents the uses of dynarnic simulation and sorne of the simulation
prograrns currently available in the market.

Principies and Practice o/ Automatic Process Cor11roVfhirrl Edition, by C. A. Smith and A. B. Corripio
ISBN 0-471-4319().7 Copyright e 2006 John Wiley & Sons, lnc.

452
 13· I Uses and Tools of Dynamic Simulation 453

13-1.1 Uses of Dynamic Simulation

Sorne common uses of process dynamic simulation are
• Study the dynamic behavior of processes.
• Design, tune, and evaluate process control systems.
• Train operators on the operation of processes.
• Design batch processes and their operation.
Notice that the design of continuous processes is left out of this Iist. The reason is
that such design is best carried out by steady-state simulators such as ASPEN PLUS,
HySys, Proll, Chemcad, and similar prograrns.
The preceding chapters of this book have presented pencil-and-paper methods
of studying process dynamic behavior, and of designing and tuning control systems.
However, the methods presented have required considerable algebraic manipulation
to obtain the results-linearization, Laplace transforms, combination of simultaneous
equations, etc. Also, little has been required in terms of analyzing the final
response of the process and the control system. In contrast, dynamic simulation
cannot only produce the same results with little or no algebraic manipulations, but
it can also produce the time responses of the control systems to various inputs,
Nevertheless, the analytical procedures presented so far are most useful for
understanding the concepts of stability, the root causes of instability, the need for
dynarnic compensators and how to design them, and the like.
As pointed out in the list above, dynamic simulation is an excellent tool for training
operators in the operation of highly integrated processes and their control systems. This
practice, modeled after the training of astronauts for space travel, has been practiced
in industry for at least four decades. Such training is sometimes carried out with a
simulator interfaced to a control console that looks and feels like the console the
operators use to control the real process.

13-1.2 Tools for Dynamic Simulation

Before the personal computer revolution, dynamic simulation required considerable
prograrnming skills from the modeler. Initially simulation was carried out in ana-
log computers. The differential equations were solved by wiring electrical compo-
nents that were analogous to the physical systems being simulated. Not only was the
required wiring of electrical components cumbersome, but the requirement of scaling
the problem variables to the set ranges of the simulator was a slcill difficult to master.
However, analog simulation had the advantage of graphically connecting components
in a way that matched the way that control system components are connected, as, for
example, in block diagrams.
When analog simulation was replaced in the 1970s and 1980s by digital computer
programs, the set of skills required from the modeler changed, but the degree of diffi-
culty was not much reduced. Procedural languages such as FORTRAN, BASIC, Pascal,
and C were required to put together and solve the differential equations representing
the dynamics of the system. Fortunately, a number of numerical integration routines
became available that eliminated the need for the modeler to develop his or her own.
Currently there are a good number of computer programs that make dynarnic
simulation extremely accessible to modelers with little or no programming skills. A
list of the most commonly used prograrns is given in Table 13-1.1. The ones that
are most useful have the following features:
454 Chapter 13 Dynamic Simulation of Control Systems

Table 13-1.1 Dynamic simulation programs

Program Cornpany Web page

MATLAB/ MathWorks [Link] Procedural la

Simulink block-orien
programmi
VisSim Visual [Link] Block-oriente
Solutions, programm
lnc. with sorne
MathCad
TK Solver UTS [Link] Equatio
handl
symbolic m
MathCad MathSoft [Link] Equation-forr
handles un
symbolic m
Mathematica Wolfrann [Link] Equation-form
Research, handles un
lnc. symbolic m
LabVIEW National [Link] Block-oriente
Instruments collection
program th
for simula
Maple Waterloo [Link] Comrnand la
Maple, !ne. includes in
plotting su
Aspen Plus AspenTech www .[Link] Provides unit
Dynamic modules a
therrnodyn
packages
HySys AspenTech www .[Link] Provides unit
modules and
thermodynamic property
packages

• Graphical interface for block-oriented prograrnming.

• Blocks to generate input signals (sources or signa] generators) such as steps,
ramps, and sine waves.
• Time plots of variables.
• Blocks to si mu late the most common control devices-transportation lag,
trans- fer functions, signa! lirniters, sample-and-hold, etc.
• Ability to develop custom simulation blocks using a procedural language.
The examples in this chapter use MATLAB with Simulink because it is one of the
most powerful dynamic simulation programs that is readily available to students and
engineers. Toe web page accompanying this book also has VisSim versions of most
simulations. VisSim is currently available as part of MathCad.
As pointed out in the introduction to this chapter, there are two basic meth-
ods to carry out the dynamic simulation of a system. The next section presents the
simulation of linear systems by simulating the transfer functions of their components.
 13-2 Simulation of Linear Transfer Functions 455

The rest of the sections of this chapter present the simulation of models developed
from basic principies.

13-2 SIMULATION OF LINEAR TRANSFER FUNCTIONS

ln Chapters 2. 3. 4, and 5 we learned to develop transfer functions of processes, con-
trollers, control valves, and sensor-transmitters; we did this by application of the
Laplace transform to linear dynamic models. Processes represented by nonlinear dif-
ferential equations required linearization, so that the transfer functions were valid only
in a narrow range of operating conditions over which the linear approximations were
valid. When such transfer function models are available, the transfer functions can be
simulated to study the response of the control system to various inputs. The following
example shows how to carry out such a simulation.

EXAMPLE 13-2.1 In Example 6-1. l we developed the transfer functions for the feedback control loop to control the
temperature of a continuous stirred tank heater. Simulate the loop using the transfer functions
Simu/ation of developed there. Then use the simulation to verify the ultimate gain and period of the loop
Continuous Stirred determined in Example 6-2.2. Check also the response of a proportional controller tuned for
Tank Heater quarter decay ratio response to step change in process How.

SOLUTION Figure 13-2.1 shows the block diagram for the simulation of this problem in MATLAB with
Simulink. The blocks used are as follows:
• Step change blocks are u sed for the inputs: the set point SP and the process flow F.
• Gain blocks are used for the proportional controller Kc and the set point gain Ksp·

�
u
F
r
l.079s+2.06
2.583s2+5.454s+0.617
Gf
;c..�������
- T

M, %CO 1.652 0.730

0.2s+l 2.583s2+5.454s+0.617
SP Ksp Kc Gv Gs

C, %TO
0.75s+l
H

Seo pe

.... ""'" c::::::J •""'4S

Figure 13-2.1 Simulink diagram for simulation of stirred tank heater using transfer functions.
456 Chapcer 13 Dynamic Simulation of Control Systems

• Summers are used for the error computation E and the sum of the process flow and
steam flow effects on the outlet temperature T.
• Transfer function blocks are used for the valve Gv, the sensor-transmiuer H, and thc
process transfer functions G F and G s.
• A scope block is used to plot the temperature response versus time.
The values of the transfer functions are obtained from the results of Example 6-1. l. Most of the
numerical values are from Table 6-1.1, except for the transfer functions G F and G s that
muse be computed into the fonn required by the Simulink block. These calculations, using
values from Table 6-1.1, are
2.06(0.524s + 1)
=-----'----'-----
{r s+ l)(r,s + 1) - K, (4.93s +l)(0.524s 1) -0.383
l .079s + 2.06
= 2.583s2 + 5.454s + 0.617
KwK, l.905 X 0.383
G,(s) = = ----'-------
(rs + l)(r,s + 1) - K, (4.93s + 1)(0.524s + 1) -0.383
0.730
- 2.583s2 + 5.454s + 0.617
To obtain the ultimate gain and period the controller gain is set and the simulation is run, checking
the response in the scope after each run. The gain is then increased until the response is a
sustained oscillation. When sustained oscillations are obtained, the controller gain is the ultimate
gain and the ultimate period is the period of the response. Figure 13-2.2 shows the response
of the transmitter and controller ourputs when the controller gain is set at 10.4 %CO/%TO.
From the figure we see that three complete oscillations occur in 13.9 min, for a period of

Scope

C%TO
2
l.8
I '\ / '\ /
l.6
l.4
I \ I \ I \ I \
l.2
I
I \
\ I
I 1
1
' I
I 1
1
l
I \ I 1 I 1
I'
0.8
0.6
I \ I 1 \ I \
1 I I \ \
0.4
\ / \ / \ \
0.2
I ' '-' '-'
o M%CO
20

l5

10
'\
I\ ""'\ ¡, /
5

o
\ I \ I \ / \ I
\\_ I \ / \
' \ /
-5

-10
/
'l/ V \.., J

o 2 4 6 8 10 12 14 16 18 20
Time offset O

Figure 13-2.2 Response of transmitter and controller outputs for the ultimate controller gain.
 4 Chapter Dynamic Simulation of Control 13 Process Simulation

R•Scope ¡;;J[¡¡Jég
11 • lll UPIJO fD H� lll lllUjllBl-1>1 1
C%TO
o
-0.05
\
\
--0.1
--0.15
--0.2
1 / <,
-0.25
\
'
j
--0.3 <;
-0.35
\ / /

--0.4
\ /
\ /
--0.45
--0.5 '
M%CO
2.5

2
I
r\
1.5
I
/ \ \
/ <,

,_/
1

)
0.5

o
o
/ 2 4 6 8 10 12 14 16 18 20
lime offset O

Figure 13-2.3 Response of transmitter and controller outputs to a step change in process flow.

13.9/3 = 4.6 min. These are the values obtained by direct substitution in Example 6-2.2. The
responses to a step change in process flo� of 1.0 lb/min are shown in Fig. 13-2.3. Toe
controller
gain is set to 10.4/2 = 5.2 %CO/%TO, as indicated by the quarter decay ratio tuning formulas
of Table 7-1.1 for a proportional controller. Such responses, easy to obtain from the simulation,
require significant aJgebraic manipulations to obtain by analytical meaos.

Comparison of the Simulink diagram of Fig. 13-2.1 to the block diagram of Fig. 6-1.7
shows a one-to-one correspondence of the blocks. This is one of the major advantages
of block-oriented simulation for studying process control systems.
The initial conditions of ali the variables in the diagram are zero because, being
transfer functions, the variables are deviations from initial steady-state conditions (see
Section 2-3.1 and Exarnple 6-1.1).
Justas the simulation of Exarnple 13-2.1 was used to tune the controller by the
ultimate gain method, it could have been used to determine the unit step response of
the process. To do this, the controller is simply replaced by a step input block.
Having demonstrated how to simulate linear transfer function models of processes,
the next section presents a more fundamental method of simulating differential equation
models of the process.

13-3 PROCESS SIMULATION

An altemative way to simulate the process to the one presented in the preceding
section is to numerically solve the differential equations that represent the model. In
4 Chapter Dynamic Simulation of Control 13 Process Simulation

this manner we not only save having to develop the transfer function of the process,
but also having to linearize the process model if it happens to be nonlinear. In fact
both linear and nonlinear process models can be simulated with equal ease.
Today's simulation software offers a variety of very powerful methods for the
numerical solution of differential equations. Further, it is completely transparent to the
user of the software which method is used. Ali that is needed is a model in which each
differential equation is in first-order form, such as:
d x( t)
- -- ; ¡¡ - - = f[x(t), u(t), t] (13-
3.1)

where
x(t) = state variable
u(t) = input variable
t = time
and f is any linear or nonlinear function of the variables. Fortunately, process
models developed from basic principies are exactly in the form of Eq. 13-3.1, as you
can easily verify by examining ali of the models developed in Chapters 3 and 4.
Block-diagrarn oriented simulation software provides an integrator block designed
to numerically integrate equations in the form of Eq. 13-3.1. In addition to the input
derivative, the integrator block requires the initial condition of the state variable x.
The symbol lis is used in Simulink and VisSim to indicate an integrator block. This
does not mean that the block uses Laplace transforms to carry out the integration. lt
is only a convenient label for the block.
The following examples show how to simulate processes using MATLAB with
Si mu link.

EXAMPLE 13-3.1 This is the gas process modeled in Section 3-6 to obtain a response of the pressure of the gas
in the tank to a step input in the signal to the fan pushing the gas into the tank. The model
Simulation of Gas developed in Section 3-6 is Eq. 3-6.1 O:
Process
V dp(t)
lif,(t)-pf.(t)=--
RT dt
1
°
f;(t) = 0.16m,(t)
/0 (1) = 0.00506m0(t)j p(t)[p(t) - p¡ (1)]

where
,n;(t) = signal to the input fan, %
ni0(t) = signal to outlet valve, %
p(t) = pressure in the tank, psia
PI (t) = pressure downstream of valve, psia
p = gas density, 0.00263 lbmole/scf
V = tank volume, 20 ft3
R = ideal gas law constant, 10.73 psia-ft3/lbmole-ºR
T = gas temperature, 520°R
and the constants are in units of scf/min-%. The tank model consists of one differential equation
and two algebraic equations. The initial condition of the pressure, p(O), is needed to define
the state of the tank, so the pressure is said to be the state variable. There are three inputs or
forcing functions, the two signals m;(t) and m0{l), and the downstream pressure Pi(t).
4 Chapter Dynamic Simulation of Control 13 Process Simulation

For a simulation to properly represent the response of the process to the input signals it
is necessary that the initial conditions be at steady state. This steady state is usually obtained
from the process design conditions. However, we must realize that the steady-state requirement
imposes a restriction on the values the design conditions can take. Of the design conditions
given in Section 3-6, the minimum required specifications are
m;(O) = m; = 50 % m0(0) = m = 50 %
0 P1(0) = p1 = 14.7 psia
Using the three model equations, with dp(O)/dt = O at steady state, the other three initial
conditions can be determined,
/;(O)= 0.16m;(O) = (0.16)(50) = 8.0 scf/min
/0(0) =/;(O)= 8.0 scf/min
-------
/.(0) = 0.00506m0(0)J p(O)[p(O) - p¡ (O)]
.. p(O) = 39.8 psia
Having the initial conditions, the differential equation is solved for the derivative as follows:
p( ) pRT
d t =
v [f;(t) -
fo(/)]
= (0.00263);i·73)(520) [f;(t) _ fo(I)]

= 0.734[/;(1) - /0(1)]
Figure 13-3.l shows the Simulink block diagram used to simulate the gas process. Step
input blocks are used to generate the three input signals, and multiplication, addition, and a
square-root function block are used to carry out the calculation of the outlet flow. Gain
blocks are used to introduce the constant coefficients in the model, and an integrator block
integrates

file �d1t yiew �imulation F2'fTlill !ools �lp

, O. i6S::::,-,,.--,,-�-+f
mi,% fi, scf/min dp/dt
mi lntegrator

p(t)

+
Function pl, psia

pi

Scope

'""'
Figure 13-3.1 Simulink diagram for the simulation of the gas process.
4 Chapter Dynamic Simulation of Control

1v Scope ¡;Jlil]IEJ
11• !i IIP\JO >"A!.::1121 llll�la\'l'I
1
p(t)
43
42.5 _,,.....
42
V
41.S
41 /
40.S
/
40
/
39.S

39
55 m(t)

54.5
54
53.S
53
52.S

Sl " .
S
51
50.5
so
o 5 10 15 20 25 30 35 40
Time offset O

Figure 13-3.2 Response of gas pressure to step change in the signal to the fan (inlet How).

the differential equation to obtain the pressure in the tank. The initial condition on the integrator
is set to 39.8 psia.
Figure J 3-3.2 shows the response of the pressure to a step change of 5 % in the signal to
the inlet fan. The step change takes place at t = 5 min, to show that the pressure is indeed at
steady state before the inlet How changes. The response is complete in about 25 min, which is
consistent with the time constant of 5.24 min determined in Section 3-6.
You might ask what determines that the units of time are minutes. Toe fact that the units of
d p(t)/dt are psia/min, is what determines that the time variable is in minutes, as the
integration with respect to time gives the tank pressure in psia.

The following example shows how to simulate a process represented by multiple

differential equations.

EXAMPLE 13-3.2 This is the stearn heater of Example 6-1.1, the same heater that was simulated from the lin-
ear transfer functions in Example 13-2.1. Here we will simulate it from the basic
Simulation of differential equations derived in Example 6-1.1.
Stirred-Tank Heater
dT(t) 1 UA
-- = - f(t)[T;(t) - T(t)]+ --[T,(t) - T(t)]
dt V Vpc,
dT,( t)
- - =
dt
1
-(J.w(t)
CM
- U A[T,(t) - T(t)]}

where we have solved for the derivatives. The explanation of the variables and the values of
the system parameters and of the design conditions are given in Example 6-1.1. There it was
determined that the initial steady-state conditions are T(O) = 150°F and T,(O) = 230°F. Also
the initial design conditions are f(O) = 15 ft3/ntin, T;(O) = IOOºF, and w(O) = 42.2·lb/min.
 13-3 Process Sirnulation 461

CM

'""'
Figure 13-3.3 Simulink diagrarn for stirred tank heater.

Figure 13-3.3 shows the Simulink diagram for the simulation of the heater. Step input
signa) generators are used for the three inputs: the process and steam flows, and the inJet tem-
perature. Two integrators (1/s) are used to integrate the two differential equations and arithmetic
blocks-summers, multipliers, gains, and constant blocks are used to calculate the derivative
functions from the equations given above. A scope is used to plot the response of the input
signal and the two temperatures.
Figure 13-3.4 shows the responses of the temperatures to a step change in process flow.
The process gain can be determined from its definition (136.9 - 150)°F/(16 - 15)(ft3 /min) =
-3. J°F/(ft3 /min). This compares with the analytical result obtained in Example 6-1.1, K F(
(1 - K,) = [-2.06ºF/(ft3 /min)]/(J - 0.383) = -3.3°F/(ft /min). The slight difference is due
3

to the linear approximation.

EXAMPLE 13-3.3 Many important specialty chemical products are produced in bioreactors by processes such as
fermentation. Most of these processes are carried out in batch mode by filling a tank with a
Simulation of a substrate solution and inoculating it with a smalJ amount of biomass. The biomass, feeding
Batch Bioreactor on the substrate, reproduces to produce the desired product, until the substrate is consumed. This
example is presented here to show sorne of the special characteristics of biochemical processes.
A dynamic model of the growth of the biomass concentration x(I) and of the consumplion
of the substrate concemration, s(t), is given on a per unit volume basis as follows:
x( )
d t = µ,
(t)x(t)

ds(t) 1
-- = --µ,(t)x(t)
dt y
462 Chapter 13 Dynamic Simulation oí Control Systems

jFloatinS scopej rm, ºF

100.---- �,---.-======'...._.------ � :'.._---.------,------,
149.5f-------'"-s;,:--f--------tf--------t------t------t------l
1491----....::,.....c+------+------t-----+------+------I
148.5f-------f',....=---+-------,-----+------+------i
1,af-------t--:::,,,-..-1::::------l-----+----+------t
14174.57�r:=====t=====t::::'.=:'.:::=�===�;;i::====t====
::J
146.5�----�----�----c-c' -c-----'-------'------'
TslO, ºF
230.5,=.::...._:----,-----,---- --,-----.-----,------,
230f-----""-sc:--t----+-------l-----+----+------t
29.51----'- "'=f------t------t------+------+------I
2229f------�"-_c ,....._,.. ._:---t------+-------+-----I------�
228.5f-------f----.. : :,,,_<:::l:------t------t------t------l
22 7 .5
22 2 8t========�==========���====���::::�::::::������
��
1113,��---�----�---- - --------�----�
l(t ),

/mi n
- -
15.81-1----+-----+-------,-----+-----+-----;
15.6t-t----+------+-------<-----+------+-------l
15.41-t----+------+-------<-----+------+-------l
15.21-1----+-----+-------,-----+-----+-----;
15o��---�5 ----�10----�15-----2�0 ----�25
----�30
Time offset O

Figure 13-3.4 Response of heater outlet temperature and steam chest temperature to a step
change in process flow.

where y is the yield in biomass per unit mass of substrate and µ(t) is the biomass growth rate
function, usually in h-1• This growth rate function is analogous to the kinetic models used to
model chentical reactors. 1t is designed to match experimental reactor data. Here we will use
the Monod model with adaptability which has the following form:
µ.( ) [ s(t) ]
d t = e, µ.m k + s(t) - µ.(t)

where a is the adaptability parameter, and k and [Link] are the parameters oí the model. These
three differential equations are solved in the Simulink diagram of Fig. 13-3.5 with a= 15 h-1,
k = 0.5 g/liler, s(O) = 2.5 g/liter, µ.(O)= µ.m = 1.2 h-1 and x(O) = 0.001 g/liter,
Figure 13-3.6 shows a plot of the biomass and substrate concentrations with time in
hours. It takes about 3 h for the growth rate of the biomass to become noticeable, but it then
takes off exponentially while the substrato concentration plot is a mirror image of the biomass
plot. Eventually, about 10 h from the start, the substrate has been converted into biomass, the
growth stops, and the batch is complete.
Biomass reactors have been the subject of numerous control studies that have been reported
in the literature.
 The preceding examples showed how to simulare dynamic models consisting of
multiple differential equations. Each derivative function is calculated from the state
and input variables and fed to a separate integrator block.
A convenient feature of Simulink and other block-oriented simulation programs is
the ability to combine the blocks for the simulation of a system into a subsystem. In
 13-3 Process Simulation 463

file �dil �ew §.1mulation FQfmat Iools !::Jelp

!I Eill3

Batch Monod Model with Adaptability

X
x, glliter
Mu*x X

s, glliter

Seo pe

s/(k+s)

�-------< 0.5
k

Figure 13-3.5 Simulink diagram for batch bioreactor.

2r---r---r---r---r---r---r---r-----:P-----r-----1
1.51---+---+----+---f---+---+---+/-//---t---+-----I

[Link]-----+-----+--+---+----IL-----±,,-" '-----+--t-----
/
+-----t
ol_�J_�...L�....L�dc::::::::::::::I:...__JL_�L �J_�...L�-1
s, glliter
·•1-T---i--r�"T=::::::¡=--T-1--r-T-1
2

,1------1----t---+--+---+- - - ::,,,.k: -- --+---+---

- - -- < ,
+-----I
[Link]---+----+----+--->---+----+---->,'sr+----t---<-----1
r-,

0.5f---+--+---f---+--+---f---+--'<"--+---t-----l
-,
oe----l----l-----l----1-----1-----e----+---l'�--t-----l
�·�!:---_¡_---,'----�--_J,'-----�--_J,,-----:----�---,:--------!.o
Time offset O

Figure 13-3.6 Plot of biomass and substrato concentrations versus time in hours.
4 Chapter 13 Dynamic Simulation of Control

file t_dit �iew �imulation FQrmat !ools ljelp

!ll'!ilEI

f(t), ft3/min
f(t)

f(t)
T(t)
Ti(t), ºF D
Ti(t) Ts(t)
Seo pe
CST Heater

w(t), lb/min

w(t)

....,
Figure 13·3.7 Subsystem block for the continuous stirred tank heater.

this manner tbe simulation of the heater shown in Fig. 13-3.3 can be represented as a
single block with three inputs and two outputs, as in Fig. 13-3.7. lf we now replace
the block of the heater for the transfer function blocks of Fig. 13-2.1, we can study the
temperature control system with the more fundamental nonlinear model of the heater.
However, this would present a problem: the transfer function models of the valve and
the transmitter in Fig. 13-2.1 have zero initial conditions, but we need nonzero initial
conditions of the steam flow for the fundamental model of the heater. To resolve this
problem we need more fundamental models of the instrumentation. This is the subject
of the next section.

13-4 SIMULATION OF CONTROL INSTRUMENTATION

When simulating a control system using linear transfer functions, as in Section 13-2,
the instrumentation can also be simulated using linear transfer functions, Chapters 5
and 11 presented the linear transfer functions for sensor-transmitters, control valves,
feedback controllers, and lead-lag units. However, if we would limit ourselves to these
types of simulations, we would not be able to study by simulation the effect of process
nonlinearities. Such study requires simulating the process from the original differential
equations, as in Section 13-3. As we saw there, an important characteristic of these
simulations is that the initial conditions are not zero and must be carefully taken into
consideration. As simple transfer function models do not allow for nonzero initial
conditions, we rnust develop special simulation models of the various instruments.
This is the subject of this section.
 13-4 Simulation of Control lnstrumentation 465

Most of the instrument models to be developed here will require the simulation of
a first-order lag with a nonzero initial condition. The transfer function of a first-
order lag with unity gain is
Y(s)
= (13-4.1)
X(s) rs +1
Where Y(s) is the output signal. X(s) is the input signa], and r is the time constant.
To model the response of the output, the transfer function is rearranged thus:
(rs + l)Y(s) = X(s)

¡
rsY(s) = X(s) - Y(s)

Y(s) = � �[X(s) - (13-4-2)

Y(s)]}

Equation 13-4.2 is simulated by feeding an integrator block (1/s) the difference between
the input and the output signals divided by the time constant. To start at steady state,
the initial condition of the integrator must be set to the initial value of the input.
This model of a first-order lag will appear in every instrument model presented in
this section.

13-4.1 Control Valve Simulation

Figure 13-4.l shows the Simulink block diagram for the simulation of a control valve
with constant pressure drop. This diagram assumes that the valve model is a subsystem

file �d1t '[iew iimulatiofl F�mat !ools ljelp

m(O

fmax/100
Valve Characteristic function: General model of a control valve.
Linear: f(u) = u(l)
Eq.%, f(u) = 100•50•10.01•u11H) Constant pressure drop.
Q.O., f(u) = 10o•sqrt(O.O! •u(I)) lntegrator initial condltion is the
initiat o/o vatve position

''""
Figure 13-4.1 Simulink block diagram for a general control valve with constant pressure drop.
4 Chapter Dynamic Simulation of Control 13 Simulation of Control 4

with an input m(t}, the controller output in% CO, andan output f(t), the flow through
the valve in the most convenient units.
The first item in the valve model is a switch to select the valve as air-to-open
(fail closed}, in which case the signal 111(1) is fed through, or air-to-close (fail opened),
in which case the signa] is changed to (100 - m(t)]% CO. Notice that the controller
output signaJ m(t) varíes in the range O to 100 % CO and is not a deviation variable.
Therefore, changing the sign of the valve gain is not sufficient to change its action,
because that would only work on deviation variables.
The next part of the val ve model is the first-order lag with the val ve time constant
rv- This is simply solving Eq. 13-4.2. The output of the first-order lag is the valve
position vp(t), in % VP. To start at steady state the initial condition of the integrator
must be set to the initial value of the input. Notice that this could be m(O) or 100 -
111(0), depending on the action of the valve. If the val ve time constant is negligible, the
entire simulation of the first-order lag rnust be deleted or bypassed, because the time
constant cannot be set to zero.
The valve dynarnic term is followed by a general purpose "Math Function" block
to enter the valve characteristics function. Such a block has, in this case, a single input
u(l), and allows any calculation using this input. To simulate a linear valve, enter
simply u(l), and for an equal percentage valve the statement 100 * a" (0.01 * u(]) - 1)
is required, where a is the rangeability parameter of the valve, for example, 50 (see
Eq. 5-2.7). A quick-opening valve could be simulated by the square root function.
Finally, the capacity of the valve is entered as a constant multiplier, fma,/100, and
the 100 because the valve position is in % VP. When the pressure drop across the
valve is not affected by the simulation variables, the maximum flow through the valve
is a constant. This is usually the case in the simulation of temperature or composition
control loops. The maximum flow through the valve is the flow when the valve is fully
opened. For a liquid valve, from Eq. 5-2. l,

Ímax = Cv,maxy
¡¡:;;:
G¡ (13-4.3)

where fmax is in gallons per minute. However, for the purposes of the simulation of
the valve, fmax can be entered in the diagram of Fig. 13-4. l in any convenient units
(lb/hr, ft3/min, etc.). The units of the flow out of the model, f(t), will be the sarne as
the units of fmax· For a gas valve, fmax are calculated from either Eq. 5-2.3 or 5-2.4,
as required. ln these equations, Cv,max is used for Cv lo obtain the maximum flow.
For liquid val ves, when the pressure drop across the val ve is affected by the process
variables, as in the simulation of level or pressure control loops, the block diagram
of Fig. 13-4.2 should be used. The only differences with the diagram of Fig. 13-4.1
are that the flow through the valve is calculated using Eq. 5-2.1, andan additional
input is required for the pressure drop across the valve t:,.Pv· A similar diagram
could be developed for simulating a gas or steam valve.

Simulation of Flow Control Loops

The diagram of Fig. 13-4. l can be used to simulate a flow control loop when the
pressure drop across the valve is not affected by the process variables. When this
is done,

• the input m(t} is the flow set point in % TO of the flow transmitter,
• the switch is in the air-to-open position,
 4 Chapter Dynamic Simulation of Control 13 Simulation of Control 4
!::: c,al,o02JG,(sl' IJ�0
Eile �dit �iew ,S.1mulation F21mat !ools !jelp

2
DeltaPv psi

º·ªr�·-.
Gf
X
Product2 Math
Function

m(t)

Constant

Tauv Cv,max/100

Valve Characteristic function:

Linear: f(u) = u(l) General model of a control valve.
Eq.%, f(ul = 100•5o•co.01 •u(l)-]I Variable pressure drop. Constant [Link].
lntegrator initial condition is the
a.o., f(u) = 10o•sqrt(0.01 •u(l)) initial % valve position

,..,,
Figure 13�4.2 Simulink block diagram for a control valve with variable pressure drop and
constant specific gravity.

• the time constan! <v is the closed-loop time constan! of the flow loop,
, the valve characteristic function is usually linear, and
• fmax is the upper limit of the flow transmitter in the loop.
This way it is not necessary to simulate and tune the flow controller. lf, on the other
hand, the pressure drop across the valve is a function of the process variables, the
diagram of Fig. 13-4.2 must be used to simulate the valve in a loop that will include
a flow controller and a flow transmitter.

13-4.2 Simulation of Feedback Controllers

This subsection develops a model of a proportional-integral-derivative (PID) controller
that provides for nonzero initial conditions, direct limits of O and 100 % in the con-
troller output, and a dynamic gain limit on the derivative mode. This model is highly
recommended when there is a possibility that the controller output may saturate. The
model will be presented in two parts, the proportional-integral (PI) and the proportional-
derivative (PD). Then the two parts will be combined to simulate a series PID controller.
A model of a parallel PID controller is also presented.
Proportional-Integral (PI) Controller
The transfer function of a PI controller was given by Eq. 5-3.14,

G,(s) = M((s) = K, (1 + _I_) = K, (r¡s + (5-3.14)

1) s) x¡s t is
468 Chapter 13 Dynamic Sirnulation of Control Systems

where K¿ is the proportional gain and r,

is the integral time. The problem with the
direct simulation of this transfer function is that a limiter on the controller output <loes
not prevent the integrator output from winding up, while it is not easy to directly limit
the integrator. A more convenient model is developed as follows:
r¡sM(s) = KcE(s)(r¡s + 1)
r¡sM(s) + M(s) = KcE(s)(r,s + 1) + M(s)
M(s) = KcE(s) + x¡s 1+ 1 M(s) (13-4.4)

Equation 13-4.4 is programmed in the Simulink block diagram of Fig. 13-4.3. The
diagram shows that the controller output m(t) is the sum of the terms KcE and the
output lagged by a first-order time constan! r1. A limiter with limits O and 100 % CO
is inserted before the ourput is fed back to the calculation. Thus there is no way for
any of the variables in the loop to go outside those limits. To start at steady state, the
integrator initial condition must be the desired initial value of the controller output,
m(O), counting on the initial error to be zero. The only undesirable feature of the PI
model given here is that it is not possible to combine it with a derivative terrn to
simulate a parallel PID controller.

Proportional-Derivative (PD) Unit

We learned in Section 5-3 that the derivative mode requires a filter in practice:
U(s) rvs
= (13-4.5)
X(s) at os +1
1!11 l:il 13

X +
+ m(t), %CO
OP
Product Saturation
0.8
%TO

The integrator initial condition is the

initial value of the controller output. lntegrator

'""
..., . "
'"
Figure 13-4.3 Simulink block diagram for proportional-integral (PI) controller with
saturation limits.
4 Chapter Dynamic Simulation of Control 13 Simulation of Control 4
where U(s) is the output, X(s) is the input, and ro is the derivative time. The filter
consists of a first-order lag with a time constant ar0. The value of a is selected
small enough-0.05 to 0.2-so that the filter does not affect the performance of the
controller. The filter limits the output of the derivative mode on a sudden input, such
as a step change, the limit being equal to l/a times the magnitude of the step. This is
why the filter is also called the dynamic gain limit. The model of the filter derivative
of Eq. 13-4.5 is developed thus:
arosU(s) + U(s) = rosX(s)

U(s) = -1 -(rosX(s) - U(s)) = -1 ( X(s) -

U(s) ) (13-4.6)
-1
aros a ros
This equation gives a filtered derivative mode. To obtain the PD unit we simply add
the input variable.
Y(s) = X(s) + U(s) (13-4.7)
Where Y(s) is the output. Figure 13-4.4 shows the Simulink block diagram for the
PD unit. The constan! multiplier 1/a is set egua! to 8, which is then the dynamic gain
limit. To prevent a bump at the start the integrator initial condition must be the initial
value of the input x(O).
Combining Eqs. 13-4.5 and 13-4.7, the transfer function for the PD unit is
Y(s) ros (a+ !)ros+ l
--= l+ =------ (13-4.8)
X(s) aros+ 1 aros+ l

Proportional Plus Derivative Unit

+
x(t), %TO + y(tl, % TO
In Out

+ 8
u(t)

1/Alpha

U(s)fTauD*s s +

lntegrator Product
Tau O

The integrator initiar condition must be

the initial value of the input.

'""' .
Figure 13·4.4 Simulink block diagram for proportional-plus-derivative (PD) unit.""This also
serves to introduce a net lead of TauD with an atrenuation factor of Alpha.
4 Chapter Dynamic Simulation of Control 13 Simulation of Control 4
We see that the PO unit is a lead-lag unit with a net lead of ro and a ratio of the lead
to the lag of (a+ 1)/a. This means that the diagrarn ofFig. 13-4.4 can also be used to
insert a net lead with a negligible lag into other control schemes such as feedforward
controllers and decouplers.

Simulation of a Series Proportlonal-Integral-Derrvatíve (PID) Controller

The diagrams of Figs. 13-4.3 and 13-4.4 are Simulink subsystem blocks that can now
be put together into a series PID controller as in Fig. 13-4.5. Notice that the PO unit
is connected on the measurement input instead of on the error signa!. This is the
recommended arrangement. If derivative action on the error signa! is desired, the PO
block should be moved into the PI block of Fig. 13-4.3 and inserted just after the error
calculation. The diagram of Fig. 13-4.5 is itself a Simulink subsystem that can be used
to simulate the PID controllers in the simulation of a control system.
The diagram of Fig. 13-4.5 includes an auto-manual switch to introduce a step
change into the controller output. This is where the step change must be applied to
carry out the process step test procedure described in Section 7-2.1. To obtain the step
response the switch is put in the bottom position and the step function M (t) is set
up to introduce a step change. Notice that the initial value of the step function must
match the initial value of the controller output, and its final value is the initial value
plus the magnitude of the step change.

!!I 1§113

1
� -� %TO
SP

�-2� %TO
In
PV

P+D

The initial value of the step Step

block must match the initial
value of the controller cutput.

....,
Figure 13·4.5 Block diagram for series P[O controller with derivative on the measurement.
The P + l and P + D blocks are subsystems containing the diagrarns of Figs. 13-4.2 and
13·4.4, respectively. The auto-manual switch allows a step change in controller output.
 13-4 Simulation of Control Instrumentation 471

Simulation of a Parallel Proportional plus Integral plus Derivative (P +I+

D) Controller
Figure 13-4.6 shows a Simulink block diagram for a parallel proportional-integral-
derivative (P + 1 + D) controller. This model was developed by adding the integral
terrn to the PD block of Fig. 13-4.4, and it has the following features:
• The derivative term is filtered with a first-order lag with a time constant equal
tO TD/8.
• The derivative term acts on the process variable and not on the error. This is
the recommended arrangement.
• lf the initial conditions are not zero, the diagrarn explains how to set the initial
conditions for the two integrators.
The transfer function of the controller of Fig. 13-4.6 was presented in Chapter 5 as
Eq. 5-3.17.

13-4.3 Simulation of Sensors-Transmitters

The Simulink block diagram of a sensor-transmitter model shown in Fig. 13-4.7 has
two parts, a dynamic terrn consisting of first-order lag with time constant rr, and the
change of scale from engineering units to %TO. Also shown in the model is the change
of scale of the set point from engineering units to %TO, which is identical to that
of

file gd1t �iew �!mulation fQOMI Iools !:!elp

!I 1::1113

Proportional Plus Integral Plus Derivative Controller

Product
OP
PV %TO

Y(sl
Taul

lnitial condition on 'Integral' is the X

initial condition of the output Y(O).
lnitial condition on 'lntPD' is Kc*PV(Ol. +
Product2
Taud

Ready

Figure 13-4.6 Simulation of paraJlel proportional-integral-derivative (P + l + D) controller

with derivative on the measurement.
4 Chapter Dynamic Simulation of Control

!ll':ilEI

b(t), % TO
PV
Productl Gain

Measurement scale change to % TO

TauT
The integrator initiat Tlow
condition is the initial
value of the input Set point scale change to % TO
variable, c{O)
THi

Tlow and TH i are the lower and

%TO
SP
I l--------------'
cset(tl E.U.

. '
Figure 13-4.7 Simulink block diagram to simulate a sensor transmitter and change the set
point scale.

the process variable. lf the transmitter is not connected to a feedback controller, the
set point scale change portion of the diagram can be deleted.
To start at steady state, the integrator initial condition must be the same as the
initial value of the measured variable, c(O). The scale change calculation shown in
Fig. 13-4.7 is for a linear transmitter. lt is derived by knowing that the transmitter
output signa] b(t) is zero when the process variable c(I) = Ttow- and 100 %TO
when the process variable c(t) = TM. So the calculation is
c(I) - T¡
b(t) = OW 100 %TO (13-4.9)
Thi - Ttow
Where T1ow and TM are the lower and upper limits of the transmitter range, in the same
units as c(I). The same calculation is performed on the set point input to produce the
controller reference signa) r(t) in %TO.
lf the transmitter is nonlinear, a Simulink "Math Function" block must be used to
program the nonlinear function. This is the block type used in Figs. 13-4.1 and 13-4.2
for the valve characteristics.

13-4.4 Simulation of Lead-Lag Dynamic Compensation

Section 11.3 showed that lead-lag units are commonly used to introduce dynamic
compensation in feedforward control schemes. In many situations, the tuning of the
lead-lag device can be simplified considerably by using either the required net lead
 13-4 Simulation of Control lnstrumentation 473

without a lag or the required net lag without a lead. When this is the case, no
special device is required to simulate the lead-lag unit. For example,
• A net lead with negligible lag can be simulated using the PD unit of Fig. 13-
4.4 with rv = T/ead - T/ag· If a nonnegligible lag is desired, then set a=
T/aglrv.
• A net lag without a lead can be simulated using a first-order lag with the time
constant set to r1ag - Ttead·

The only situation not covered by the above two cases is when a net lag with a
nonnegligible lead is desired.
Figure 13-4.8 shows the Simulink block diagram for a lead-lag unit. It is developed
as follows:
Y(s) T/eadS + 1
= (13-4.10)
X(s) T/ags+l
T/agsY(s) + Y(s) = T/eadsX(s) + X(s)
Tfead I
Y(s) = -X(s) +-[X(s) - Y(s)]
Ttag r¡08s

To prevent an initial sudden change in output, the integrator initial condition must be
set to ( T/ag - <1ead)lr1ag times the initial value of the input, x(O).

Simulation of Dead Time

Simulink, VisSim, and most other block-oriented simulation programs provide a dead-
time block with the ability to set nonzero initial values.

[ile �d1t Yiew �1mulat1on [Link] Jools �elp

X
X
+ (Lead/Lag)•xcs)
Product

Y(s)
Out
X(s)
In [X(s)-Y(s))/(Lag's)
lntegrator

2 lnitial condition on the integrator must be

(lag-lead)/lag times the initial value of
Lag the input.

''"'
Figure 13-4.8 Simulink block diagram for a lead-lag unit.
4 Chapter Dynamic Simulation of Control 13 Simulation of Control 4
The following examples wiU demonstrate how to use the models of the instruments
presented here for simulating of feedback and feedforward control of the temperature
in the stirred tank heater of Example 13-3. l.

EXAMPLE 13-4.1 Simulate the feedback temperature control of the stirred tank heater of Example 6-1. l. Use a
series PID controller with the tuning parameters determined in Example 7-1.1 and show the
Feedback Temperature responses to step changes in process flow.
Control of Stirred Tank
Heater
SOLUTION Figure 13-4.9 shows the Simulink block diagram of the control loop. Step input signa) generators
are used for the three inputs, the temperature set point rSP(t), the process flow /(r), and the
process in Jet temperature T; (t). Subsystem blocks are used for the temperature sensor-transmitter,
K,p and H(s), the series PID controller, the control valve Gv(s), and the process (STHtr).
The first three blocks are those detailed, respectively, in Figs. 13-4.7, 13-4.5, and 13-4. 1,
while the process block is detailed in Fig. 13-3.3. The parameters for the instrument blocks
are those given in Examples 6-1.1 and 7-1.1, namely,
Sensor-transmitter: rr = O. 75 min T10w = 1 OOºF
PID controller: K, = 6. 1 %CO/%T O t¡ = 2.3 min
Control valve: 'fv = 0.20 min /max = 84.4lb/min
Here [Link] maximum flow [Link] [Link] valve is obtained by doubling the steady-state flow of
42.2 lb/min determined in Example 6-1.1, as the valve is sized for 100 % overcapacity .

./
� csthfb' �@]0
�ile �dit Yiew �imolation FQr!Tlal Ioo1s !::!elp

ft3/min
f(t)

TI(TJ
f(O T(t) T(t), 'F

Ts(t) Ts(t), 'F

m(t), % CO '---�
PV STHtr
Ksp & H(s)
Series PID Gv(sl

T(t), 'F

Ready
. '
Figure 13-4.9 Simulink block diagram of the temperature control loop for [Link] continuous
stirred tank heater.
 4 Chapter Dynamic Simulation of Control 13 Simulation of Control 4

-., Stope r;JIC]IEJ

ll • i!iHP!.I' .a>HC:. 111 llll�HB!irl 1
155
no ºF

150
r>: � -

V
1---"

145
Ts(t) °F
300
I '\ ,...--....,
250

200
r----... /

150

100
m(t),%CO
100
95
I .....
'
I
90
,
85
I
I
'
758

657
'
'
' , -
60
o 5 10 15 20 25 30
Time offset O

Figure 134.10 Responses of stirred tank heater to step changes in process flow.

The initial condition of the valve position, which is also that of the controller output, is
calculated so that the initial flow is 42.2 lb/min, for an equal percentage val ve with a rangeability
parameter of 50,
UJ = Wmax • 5QVfJ- I
_ ln(w/w,ruu) ln(42.2/84.4)
vp = 1 + = 1 + = 0 . 823
ln50 ln50
or 82.3 % VP. The other initial conditions are the design conditions from Example 6-1.1, that
is T(O) = 150°F, T,(O) = 230°F, b(O) = 50 % TO.
Figure 13-4. 1 O shows the responses of the outlet and steam temperatures, and of the con-
troller output, to a step change in process flow from 15 to 10 ft3hnin at t = 1 min, followed
by another step change from 10 to 20 ft3/min at t = 15 min. Notice that the response to the
increase in flow is much more oscillatory than the step in the decrease of the flow. This is proba-
bly due to the increase in the gain of the equal percentage valve which is not really needed here,
given tbe assumptions of the problem. Notice also that the controller output saturates briefly at
100 % CO for the increase in process ñow. These effects cannot be detected when the linear
transfer function simulation of Fig. 13-2.1 is used. There is another effect that could require
a more detailed simulation of the control valve. Notice how much the steam temperature rises
for the increase in process flow. As the steam temperature increases, the pressure in the stearn
chest increases, because it is the vapor pressure of water at the condensing temperature. This
could affect the pressure drop across the valve and reduce the capacity of the valve, depending
on the supply pressure of the steam. This difficulty is addressed in one of the problems at tbe
end of this chapter.

EXAMPLE 13-4.2 Design a feedforward controller for the stirred tank heater of Example 13-4.1. Toe controller
is to compensate for changes in process flow and inlet temperature. Simulate the feedforward
Feedforward Control control scheme with and without dynamic compensation. Assume a flow controller is instalJed
o/ Stirred Tank Heater to control the steam flow and the feedforward controller is to determine the set point of the
476 Chapter J 3 Dynamic Simulation of Control Systems

flow controller. The process flow transmitter has a range of O to 25 ft3/min anda negligible lag,
and the temperature transmitter has a range of 50 to 150ºF and a lag of 0.75 min. The
steam flow transmitter has a range of O to 75 lb/min, and the time constant of the flow
control loop is 0.2 min.

SOLUTION To design the feedforward controller we follow the procedure outlined in Section 11-2. The
control objective is to maintain the process temperature at the set point, T(t) = rsct(t). The
disturbances are the process flow f(t) and the inlet temperature T;(t), and the manipulated vari-
able is the steam flow set point, wse1(1). A steady-state enthalpy balance on the tank establishes
a simple relationship between these variables.
w(t)!. = f(t)pcp(T(t) - T;
(t))
PCp
:. w'"Ctl = T J(1)(T''"(1) -
T;Ctll
where p is the fluid density in lb/ft3, Cp is the specific heat in Btu/Jb-ºF, and A is the steam
latent heat of condensation in Btu/lb.
Figure 13-4.11 shows the feedforward controller based on the equation derived above. It
consists of a summer to compute the change in temperature and a multiplier to obtain the
product of the flow, the temperature change, and the constant parameter. For simplicity, the
scale conversion to the fraction of the transmitter outputs has been omitted. ln practice the
process flow and inlet temperature transmitter ranges must be taken into consideration in the
calculation. The inlet temperature transmitter H, (s) is simulated as a first-order lag with
unity gain, and the steam flow controller as a linear valve.

l!I r;¡¡ 13
file �dil �iew §.1mulatioo F21mat Iools !jelp

f(t), ft3/min
f(t)
f(t), ft3/min f(t) T(t)
�---------l-�----------1.!r,(t)
w(t) Ts(t)
X
fsp(t) f(t) w(t), lb/min CSTHtr

Ht(s) Product Steam FC

Tsp(t) �------�
68'0.8'100/(966'75)
R ho *Cp* l 00/( lambda* f
max)

The conversión of the inlet signals and of the

set point to % TO has been omitted for simplicity.

"""'
200% ode23s . "
Figure 13-4.11 Simulink block diagram of feedforward control of stirred tank heater without
dynamic compensation.
 13 Simulation of Control 4

!11':1113

f(t),ft3/min
f(t)
CSTHtr

LeadF
f(t) T(t)
Ti(t),ºF
Ti(t)
Ts(t)
In Out X
w(t),lb/min
Ht(s) LeadT Product
Steam FC
Tsp(t)
68'0.8' 100/(966'75)
Rho*Cp* 100/{Lambda*fmax)

The conversion of the inlet signals and of the

set point to% TO has been omitted for simplicity.

jcxie23s
. "
Figure 13ª4.12 Block diagram of feedforward controller with dynamic compensation.

Figure 13-4.12 shows the feedforward controller with dynamic compensation. The lead
blocks have been simulated as the PD unit of Fig. 13-4.4. After sorne tuning, the net lead for
the flow disturbance was set at 1.0 min with a= 1/1.5-this is equivalent to a Jead of 1.67 min
with a lag of 0.67 min. A net Jead of 1.7 min with a = 1 /4 was used on the inlet temperature
disturbance, equivalent to a lead of 2.12 min with a lag of 0.42 min. The requirement for a
longer lead on the temperature is to compensate for the transmitter lag.
Figure 13-4.13 shows the responses of the heater variables to a step change in process
ftow from 15 to 10 ft3/min at t = 1 min, followed by a step change in inJet temperature
from
100 to I IOºF at t = 30 min. The responses shown are for feedback control and feedforward
control with and without dynamic compensation. Compared with the feedback response, the static
feedforward controller reduces the initial deviation in temperature by about half and dynamic
compensation further halves the initial deviation. This is because the static feedforward controller
takes action as soon as the process ftow changes, while the feedback controller must wait for the
deviation in temperature to take action. The dynamic lead units give additional initial correction,
as Fig. 13-4. 13 shows.
Although slower than for the feedback controller, the feedforward controllers retum the
process temperature to its set point of I 50ºF after the initial transient. This is because in
this ideal simulation situation, the values of the process parameters-density, specific heat,
and latent heat-are known exactly. In practice a feedback controller must be installed in
combination with the feedforward controller to compensate for errors in the feedforward model.
When this is done the output of the feedback controlJer can be connected to the set point of the
feedforward controller.

The preceding examples show how to put together the models presented in
Sections 13-3 and 13-4 to simulate complete process control systems. One
advantage
478 Chapter 13 Oynamic Simulation of Control Systems

l!JIJ Scope ¡;J[g]l8]

11 a lll HPI ,o Jl> Hc:i lll ll'll@IIBl'ir I 1

T,°F
153
152.5
J \
152 \
151.5
151 \ ,. '
150.5 ,
150 ' "
149.5 ' '
149
To,ºF
240
230
1\
220
\\
210
200 v\/ \
190
180
60 m%CO

50
1\
 H / �
30
{{
20
10
o
o 10 20 30 40 50 60
lime offset O

Figure 13-4.13 Comparison of responses of the stirred tank heater for feedback control and
feedforward control with and without dynamic compensation.

of setting the models as subsystems is that they can be easily copied and pasted
from one simulation to another.

13-5 OTHER SIMULATION ASPECTS

So far we have discussed just the basic aspects of dynamic simulation as it applies
to continuous systems. Most process simulators provide a number of features and
specialized blocks that are extremely useful to the engineer in modeling processes.
Sorne of these will be briefly discussed in this section.

Stiffness
When numerically solving a system of differential equations, sorne equations may
respond much faster than others. This characteristic is known as stiffness, and it requires
special numerical integration methods. For exarnple, in the simulations above, the PD
and lead units introduce the time constant aro which is much smaller than the
other time constants in the simulation. When a standard integration method is used to
solve a stiff system, such as MATLAB' s ODE45 or ODE23, the response of sorne
variables will be noisy due to integration errors. When this happens, select a method
designed to handle stiffness-such as ODE23s-the "s" at the end of the method's
name in MATLAB designates the method as appropriate for stiff systems.

Complex Models
Sorne dynamic models are too complex to simulate using computation blocks. Examples
are models of tray distillation columns and distributed parameter systerns, that is,
 13-5 Other Simulation Aspects 479

systems in which the variables are functions of both space and time. To program the
many equations, procedure-oriented programs such as C++, Pascal, Basic, or MAT-
LAB m-files, are much more convenient than programming with blocks. However,
block programming is still desirable to simulate the control and instrumentation com-
ponents and to put together the control system. It is therefore convenient to be able
to create customized blocks that solve the equations in a procedural language while
accepting inputs from and generating outputs to the control components. Many block-
oriented programs provide such a feature. In MATLAB the S-function block provides
that feature.
The S-function block accepts a vector of input variables and generales a vector
of output variables. Procedural statements contained in a function stored in an m-file,
which is associated with the block, perform the calculations.

Discrete Time Operations

In sorne control systems the process variables are not continuously available to the
controllers. For example, when a gas chromatograph is used to measure the composition
of a process stream, the analyzer takes a sample and analyzes it overa period of time,
the cycle time. The controller does not see the composition until the end of the cycle.
Also, the controller cannot tell if the composition changes during the cycle.
To simulate discrete time operations, simulation programs provide a set of blocks,
known as discrete blocks. Two of the most useful are the sample-and-hold and the delay
blocks. The sample-and-hold block samples its input variable and maintains it constant
for one sample time, just Jike the analyzer. However, to fully simulate the analyzer
a delay block of one sample must be connected after the sample-and-hold block to
simulate the fact that the composition is not available until the end of the cycle.

Other Features
Modero simulation programs provide many other features, including the storage of
the response data in files to use with other programs, the manipulation of vectors and
matrices, blocks for the simulation of special effects such as hysteresis and dead bands,
and so on. The user of such programs should carefully review the on-line help they
provide for learning how to use these features when they are needed.

EXAMPLE 13-5.1 Although distillalion columns are sorne of the most complex operations modeled, a simple
model is presented here to demonstrate the use of the S-function in Simulink. Figure 13-5.1
Simple Model of a is a schematic of the column.
Distillation Column A saturated liquid feed, consisting of 55 mole% benzene (xF) and the balance
toluene, flows at the rate F of 3500 Jbmole/h into a distillation column. For the purpose of this
model it may be assumed that the relative volatility a is constant and equal to 2.46, and that the
latent heat of the mixture A is constant and equal to 13,860 Btu/lbmole. The column can be
represented by
6 ideal plates with the feed entering between plates 3 and 4. Each plate has a constara hold-up
of 450 lbmole. A condenser may be assumed to condense ali the vapors leaving the column to
rnaintain the column pressure constant, and the condensate leaves the condenser at its bubble
poinl. In the condenser accumulator drum, the hold-up varies from 50 to 550 lbmole as the
level transmitter signa( varies from O LO 100 % TO. It may be assumed that the rate Vs of the
vapors leaving the reboiler is proportional to the steam ttow Ws. Saturated steam at 15 psig
is used. At this pressure the Iatent heat of condensation As is 952 Btu/lb. The partial reboiler
behaves asan additional equilibrium stage. The column bottom has a hold-up of 100 lbrnole
when the
480 Chapter 13 Dynamic Simulation of Control Systems

Fx,
Q

Figure 13-5.1 Sketch of the distillation column for

Example 13-5.1.

level is at the bottom of the rransmitter range and 1100 lbmole when the leve) is at the top
of the range.
At the design conditions the distillate composition xo is 92 mole % benzene and the
bottoms composition XB is 10 mole%. Toe reftux ratio L,/ D is approximately 2.5.
Model the column assuming equimolaJ overflow, that is, no change in liquid and vapor
rates from stage to stage, and perfectly mixed trays, accumulator, and bottom.
SOLUTION In modeling processes it is a good idea to first define the input variables to the process. For the
columns these are the feed rate F, benzene mole fraction XF, and enthalpy condition q, that is,
the fraction that is liquid at the conditions of the feed tray. For saturated liquid feed q = 1. The
variables to be manipulated by the control systern are the reftux Lr, the steam ftow to reboiler
Ws, the distillate ftow D, and the bottoms product flow B.
By the assumption of equimolal overflow, the liquid leaving each tray is equal to L, in the
rectifying section-above the feed tray-and Ls in the stripping section-below the feed tray.
Also the vapor rates are V, and Vs in the rectifying and stripping sections, respectively.
From an enthalpy balance on the reboiler, neglecting sensible heat effects, we obtain the
proportionality between the steam rate and the stripping vapor rate:

V,(1) = -x, W,(1) = 952 Btu/lb

W,(1) = 0.0687W,(1)
x 13,860 Btu/lbmole
From the definition of the feed enthalpy condition, total mass balances on the feed tray give
l,(1) = L,(1) + q F(1)
V,(1) = V,(1) + (! - q)F(t)
A benzene balance on each stage, assuming constant hold-up and perfect mixing, gives the
following:
dx, (1)
M;-- = l,(1)[xo(1) - x1(1)]
di + V,(1)[)'2(1) - y,(1)]

dx2(1)
M; - = l,(1)[x,(1) - x2(1)] + V,(l)[y,(1)- )'2(1)]
di
dx3(1)
M;-- = l,(1)[x2(I) - x3(1)] + V,(1)(y4(1) - )'3(1)]
di
dx4(1)
M;-- = l,(1)x3(1) V,(1)y5(1) - V,(l)y4(1) - L,(l)x4(1)
di +FxF(I)
 13-5 Other Simulation Aspects 481

dxs(t)
M;---¡¡--
=L,(1)[x4(t) - x5(t)] V,(t)[y6(t) - ys(t)]

dx6(1)
M¡-- = L,(1)[x5(t) - X6(1)]
dt + V,(t)[ys(t) - Y6(1)]
6 eq., 14 unk. (6x;, 6 y;, xo, YB)
The equilibrium relationships, based on the relative volatility. give us the vapor compositions:
ax;(t)
y¡ ()t =------ i = 1, ... ,6
1 (a - l)x;(t)
+
12 eq., 14 unk.
The equilibrium relationship also applies to the reboiler, as it is assumed to be an equilibrium
stage.

ys(t) = axs(t)
1 + (a- l)xs(t)
13 eq., 15 unk. (x8)
Total and benzene mole balances are then canied out on the condenser accumulator.
dMv(t) = V,(t)- L,(t)- D(t)
dt
d[M v(t)xv(t)]
--d--- = V,(t)y, (t) - [L,(r)
+D(t)]xv(t)
15 eq., 16 unk. (Mo)
Similarly, total and benzene mole balances are written around the column bottoms and reboiler,
assumed perfectly mixed:
dMs(t)
-d-- = L,(r) - B(t) - V,(t)

r z : : ' : = L,(t)x6(t) - B(t)xs(t) - V,(t)ys(t)

[M s (t)x s ( t) ]

17 eq., 17 unk. (Solved!)

This completes the simple model of the column. It must be kept in mind that this model
requires two level controllers to maintain the material balances around the accumulator and the
column bottoms.
Because the hold-ups are not constant, the benzene balances on the accumulator and bonoms
do not give the derivative of the compositions directly. This problem may be resolved with the
aid of a little algebra. Expanding the derivative of the accumulator, and substituting the total
mole balance, we obtain
d[M o(t)xv(t)] dxv(t) d M o(t)
---d--- = Mo(t)-d-- +xv(t)-d-,-
dxo(t)
= Mo(t) - +xo(t)[V,(t)- D(t)-L,(t)]
dt
= V,(t)y1 (t) - [D(r) + L,(t)]xv(t)
dxv(t)
Mo(t)-- = V,(t)[y,(t)-xv(t)]
dt
Similarly, the benzene balance around the column bottoms is converted to:
dxs(t)
Ms(t) - = L,(t)[x6(t)-xs(t)]- V,(t)[ys(t)-xs(t)]
dt
The moclel of the column consists of I O differential equations and the 7 equilibrium relationships,
plus the 2 total mole balances around the feed stage. As in most control system simulations, it is
desirable to start from the initial conditions at steady state. With modem computation tools such
4 Chapter Dynamic Simulation of Control

Table 13-5.1 Distillation column conditions

Design conditions Initial conditions

F = 3500
lbmole/h XF = 0.55 xo =
D = 1921 lbmole/h q = 1.00 x, =
B = 1579 lbmole/h x2 =
L, = 4802 lbmolefu V, = 6723 lbmole/h X3 =
L, = 8302 lbmole/h V, = 6723 lbmole/h X4 =
w, = 97873 lb/h X5 =
xo =
XB =

file [Link] Y'iew �imulation F2"flll !ools t!elp

Feed

xF

q
Reflux
Scopel
D Block Para111eters S functmn
,
:
:::::�-;;;:-é;ks ���c:-;-��!
' � ,:J ��
AOla'd l'llllli conb"' IOS.f� SUlldatds. V,u IAd 11'8:

- -- - - - -
r �= s:.=·- . � : l d !
)o@ñí3
�::, t i
1

I
OP BLCSP 1
f(t) m(t)
PVI+----' ·- ··- - - -- . . . .J
Bvalve BLC

[Link] �IIIIIM. ¡
Di:J�Dii]c:¡;;;¡¡::J

Figure 13-5.2 Sirnulink schematic for the distilJation column using an S-function.

as MATLAB and MathCAD, the equations can be solved simultaneously for the steady-state
conditions by setting the derivatives with respect to time equal to zero. We must also set the
feed rate, composition, and enthalpy condition to their design values, as well as the bouoms
composition and the reflux ratio L,/D. The solution of the 19 equations with the remaining
19 unknowns produce Lhe steady-state conditions listed in Table 13-5.1. Notice Lhat the
initial values of the hold-ups in the condenser accumulator and column bottoms cannot be
calculated from the steady-state column equations, so we assume them to correspond to 50 % of
the level transmitter outputs. Toe initial value of the steam rate is obtained from the enthalpy
balance on the column bottoms presented above.
 13-5 Other Simulation Aspects 483

The simulation of the distillation column model using Simulink computing blocks would
be very complex and difficult to manage. A more ccnvenient method is the use of the Simulink
S-function block, which is like building a custom block to simulate the column. Figure 13-5.2
shows the Simulink block diagram for the column showing the S-function block. The block
refers to an m-file called Distillation.m, which contains the equations for the simulation of the
column. Figure 13-5.3 shows the contents of the ,n-file for the distillation model. The following
are the features of this program:

function [sys, xoJ = Distillation (t, x, flag)

% t = time (independent) variable; x = state variable vector

% u = input variable vector
% flag (See below. set by Simulink according to the stage of the problem)
%
% Binary distillation column with 6 equilibrium stages, feed on stage 4
o/o Equimolal overflow, constant relative volatility, constant pressure

xl = zeros (1, 6);

y =zeros(l,6);

% Set problem constants here

Alpha = 2.46; % Relative volatility

LambdaV = 13860; % latent heat of vapors, BTU/lbmole
Lambdas= 952; % latent heat of steam, BTU/lb (sat. 15 psig steam)
Mi =450; % Tray hold-up, lbmole
MOmin = 50; % minimum accumulator hold-up, lbmole
MOmax = 550; % maximum accumulator hold-up, lbmole
MBmin = 100; % mínimum bottoms hold-up, lbmole
MBmax = 1100; % maximum bottoms hold-up, lbmole

if flag== 1 % Calculation of the derivatives of the state variables

% Input variables

Feed = u(l)/60; % Feed rate, lbmole/min

xF = u(2); % Feed mole fraction
q = u(3); % Feed fraction liquid
Lr = u(4)/60; % Reflux flow, lbmole/min
D = u(5)/60; % Distillate rate, lbmole/min
Ws = u(6)/60; % Steam flow to reboiler, lb/min
8 = u(7)/60; % Bottoms flow, lbmo1e/min

% lnitial calculations

ls = Lr + q*Feed; % stripping liquid flow, lbmole/min

Vs = Ws*LambdaS/LambdaV; % stripping vapor rate, lbmole/min
Vr = Vs + (1-q) "Feed; % rectifying vapor rate, lbmole/min
% State variables

for i=l:6
xl(i) = x(i);
y(i) = Alpha*xl(l)/(l+(Alpha-1) -un».
end
xD = x{7);
xB = x(8);
MD = x{9);
M8 = x(lO);
yB = Alpha•xB/(l+(Alpha - 1) •xB);

·- .
Figure 13-5.3 m-file for distillation model S-function.
4 Chapter Dynamic Simulation of Control Proble 4

% Calculate the derivanves sys(i) = dx(i)/dt of the state variables

sys (1) = (Lr' CxD - xL(l)) + Vr' Cy(2) - y( 1)))/Mi; % Stage

l sys (2) = (Lr' (xl(l) - xLC2) + Vr'(y(3) -y(2)))/Mi; % Stage
2 sys (3) = (Lr' (xl(2) - xl(3) + Vr'(y(4) -y(3)))/Mi; % Stage
3 sys (4) = (Lr*xL(3) + Feed*xF - Ls*xl(4) + Vs*y(5) - Vr*y(4))/Mi;%
Feed sys (5) = (Ls' (xl(4) - xl(5) + Vs'(y(6) -y(5)))/mi; %
Stage 5 sys (6) = (Ls' (xLC5) - xLC6) + Vs'(yB - y(6)))/mi %
Stage 6

sys (7) = Vr' (y(!) - xD)/MD; % Distillate mole fraction

sys (8) = (Ls'(xL(6) - x8) - Vs'(yB - xB))/MB; % Bottoms mole fraction
sys (9) = Vr-Lr-D; % Accumulator hold-up
sys (10) = Ls-Vs-8; % Bottoms hold-up

elseif flag= = 3 % Calculate outputs

sys (1) = x(7); % Distil1ate mole fraction
sys (2) = x(8); % Bottoms mole fraction
Dlevel = (x(9) - MDmin)/(MDmax- MDmin) *100; % To accmulator level
Dlevel = min ( [Dlevel 100) );
sys (3) = max ( (Dlevel O] );
Blevel = (x(lO) - MBmin)/(MBmax- MBmin) % To bottoms level
*100; Blevel = min ( [Blevel 100] );
sys (4) = max ( [Blevel 01 );

elseif flag== O % lnitialize: [# cont. states, # disc. states = O,

sys = (10, O, 4, 7, O, O];% # outputs, # inputs]

MDo = MDmin + 0.5* (MDmax - MDmin); % lnitial ho1d-up at 50%TO

MBo = MBmin + 0.5* (MDmax - MBmin); % lntitial compostns. and hold-ups
xo = (0.8274, 0.7037, 0.570, 0.452, 0.3188, 0.1928, 0.9206, 0.0993, MDo,MBo];
else
sys = [ ];
end

Figure 13-5.3 (Co111inued)

• Starts by defining the parameters of the problem.

• Stores the seven input variables to the S-function iruo variable names that are easier to
recognize as model variables. Notice that the flows are assumed to be in lbmole/h and
lb/h, but they are divided by 60 to convert them to lbmole/min and lb/min. The unit of
the independent variable in the simulation is then minute.
• Tt moves the state variables, contained in vector x, into variable names that are easier to
identify as model variables, as was done with the input variables.
• lt calculates the derivatives of the state variables.
• Jt outputs the four output variables: xo, xs, and the levels in the accumulator and the
column bottom.
• It sets the initial conditions far the I O state variables.

The rest of the statements are requirements of the S-function that are described in the Simulink
help files.
The diagrarn of Fig. 13-5.2 shows the leve! controllers for the accurnulator drurn (DLC)
and the column bottoms (BLC). lt is important to realize that without these level controllers the
model would run away from the steady state, just as a real distillation column would.
Figure 13-5.4 shows the responses of the distillate and bottoms compositions to a step
change of 2000 lb/h in steam flow folJowed by a step in 140 lbmole/h in reflux flow. Notice
how the increase in reflux flow almost exactly cancels the effect of the stearn flow on both
 4 Chapter Dynamic Simulation of Control Proble 4

0.93,---.----.----.-----.----.----.------,-----,
0.09.9252l-;: :'\::==t====t====r====r====_.l:..:-:::=1::;:;�=J�=
=�
0. 9 15 1 -' '-c- - t - - - - t - - + -= -
0 .9 If ---- - -' <c + -- - - + --I - -+----:, r" ' ---- -
'\
c:::: :¡f-- - - -- t- - /

+ - - - + -
0.905f------f"'S::,---t-----t--
--+-
---- -J
�" --t------ct-----t----

0.91---+-"'"S--....c-+---+-/+-+---f----j---+---l
0.8951------f---"'t-<:::::-.. ---f-/-/----t------,t-----t----+----I
0.891-----l----+--==-1"---+----l-----+---+----l
[Link]�---�---�---�---�---�---�---�---�
0.1 \
0.095 \ ./ ....
0.09f-�---t---lc----+--+7"'-'--t----+--+---I
,/
0.085t----"'\rl-----t-----t-----,/ � ---+-----+-----+-----l
[Link]---t''c"-...---+---t-/--r--t---+---t----+-----j

0.070 100 200

r--
0.075!-----l--=-d-----hf----+------il-----+---+---I

300 400 500 600 700 800

Time offset o

Figure 13-5.4 Responses of the distillate and bottoms composition to a step change of
2000 lb/h in steam Oow followed by a step of 140 lbmole/h in reflux
flow.

compositions. If the steam and reflux flows were to be manipulated by a control system to control
both product compositions, the interaction between them would cause the controllers to cancel
each other's actions. Toe concept of loop interaction is studied in Chapter 12. Figure 13-5.4
also shows that the response of the column is very slow. It takes about 400 min, over 6 h, for
the column to reach a new steady state after a disturbance.

13-6 SUMMARY
Simulation is a very irnportant too) for analyzing, designing, can really save time and effort when analyzing specific sys-
and tuning control systems. We have seen a few examples on tems. lf you are planning to make a career in control systems
how to simulate simple processes and their control systems engineering, you should resolve to gain expertise in the use
using a block-oriented program. Although simulation can- of sorne dynamic simuJation program. We can tell you from
not completely replace the understanding you have acquired experience that you wiJJ find a satisfying and exciting future
through the mathematical analysis of dynamic responses, it in this field.

PROBLEMS to the center of the mass-L = 0.8 m. For small swing

angles and negligible air resistance, show that the dynamic
13-1. Simulate the response of the bird mobile of
Problem 2-9 and plot the response. Assume an initial con-
dition y(O) = 0.1 rn, and no externa! force, f(t) = O.
13-2. Model and simulare the response of a pendulum with
mass M = 0.5 kg, and a length-distance from the fulcrum
 4 Chapter Dynamic Simulation of Control Proble 4

response is represented
by

Md2x(t) = - Mg x(t)
dt2 L

where g = 9.8 m/s2 is the acceleration of gravity and x(t)

is the horizontal deviation of the pendulum from the verti-
cal position. Assume an initiaJ condition x(O) = 0.1 m. Note:
 4 Chapter Dynamic Simulation of Control Proble 4

Construct a pendulum with a string and a weight, or just use change of 0.10 mole fraction in the inlet composition z(t).
a yo-yo, and time the oscillations with a watch. See how it Assume that the flow rates are constant.
matches the response of the simulation when you match the
length of the string.
13-3. Simulate the punctured tank of Problem 2-23 and plot
the response of the pressure in the tank.
13-4. Simulate the response of the temperature of the turkey
of Problem 2-24. Assume the oven temperature is constant
at 800ºR, the turkey weighs 12 lb, is initially at 535°R,
and has an area of 3.5 ft2• Use a specific heat of 0.95
Btu/lb-ºR,
and an emissivity of 0.6.
13-5. Simulate the chemical reactor of Section 4-2.3 and
plot the responses of the reactor temperature to a step change
of 0.25 ft3/min in process flow, and of 0.1 ft3/min in
coolant flow. The reactor is initially at the design
conditions. Com- pare your responses with the responses of
the linear transfer functions given in Secuon 4-2.3 for those
responses by sim- ulating the transfer functions. Note: The
differences in the responses will be due to the nonlinear
nature of the reactor. You may decrease the size of the step
changes to observe how the responses of the nonlinear reactor
approach the responses of the transfer functions.
13-6. Simulate the mixing process of Problem 3-1 and plot
the response of the outlet concentration to a step change
of 5 gpm in flow /1. At the [Link] steady-state conditions
the flow from the tank is 100 gpm, and its concentra-
tion is 0.025 moles/cm3. The tank volume is 200 gallons,
and the feed compositions are [Link] and 0.05 moles/cm3.
Assume a tight- level controller keeps the volume in the tank
constant.
13-7. Add a feedback control loop for the concentration in
the tank of Problem 13-6. Use a transmitter with a time con-
stant of 1.0 min and a range of O to 0.1 moles/cm3, and a
linear control valve with a time constant of 0.1 min and sized
for a maximum flow of 100 gpm. The valve is installed in
the more concentrated inlet stream. Tune a PI controller and
obtain the plot of the outlet concentration to a 5-gpm step
change in the flow of the other inlet stream.
13-8. Simulate the isothermal rector of Problem 3-2 and
plot the response of the concentration of reactants at point 3
for a unit step change of 5 ft3/min in reactants flow. At
the
initial steady-state conditions the reactant flow is 50 ft3/min,
the feed concentration is 2 lbmoles/ft3, and the outlet con-
centration is 0.5 lbmoles/ft3. The tank volume is 150 ft3
(constant), the pipe length between the tank exit and point 3
is 400 ft, and the inside diameter of the pipe is 5.500 in. Note:
You may use a delay block with variable delay to simulate
the transportation lag in the pipe.
13-9. Simulate the flash drum of Problem 3-11 and plot
the response of the liquid and vapor compositions for a step
 4 Chapter Dynamic Simulation of Control Proble 4
13-10. Simulate the distillation tray of Problem 3-12 and
obtain the response of the leve! h(I) and the flow from the
tray /0(1) for step changes of 10 and 20 ft3/min in the flow
into the tray /¡(!). Determine the effect of nonlinearity by
comparing the responses for the two step changes.
13-11. Simulate the blending tank of Problem 3-18 using
the additional parameters given for this process in Problem
6-11. Use a PID controller tuned for quarter decay ratio
response and obtain the responses of the transmitter and con-
troller outputs to a step change of 0.1 m3 /min in the flow of
the concentrated solution.
13-12. Simulate the tanks of Fig. 4-1.1 to obtain the
responses of the levels in the tanks and of the flow out of the
second tank, /2(1), to a step change in the flow to the first
tank, /;(1), of0.2 m3/min. At the inicial steady state the inlet
flow is 5 m3/min, the pump flow /0(1) is 2 m3/min, and the
leve! in each tank is 2.5 m. Each tank has a cross-sectional
area of 9 m2. Compare the gains and time constants mea·
sured from the simulation responses with those calculated
from the results given in Section 4-1.1.
13-13. Do Problem 13-12 for the interacting tanks of
Fig. 4-2.1. In this case the initial levels are 5.0 m in the
first tank and 2.5 m in the second tank.
13-14. Simulate the thermal tanks in series of Fig. 4-1.5
to obtain the responses of the temperatures in the tanks for
a unit step change of 10 K in the inlet temperature to the
first tank, T1 (1). Assume the flow into the first tank !A is
constant and equal to 1.0 m3/min, and ignore the flow Is.
Each tank has a constant volume of 5.0 m3 and the liquid is
water. Compare the gains and time constants from the simu-
lation responses with those calculated from the results given
in Section 4-1.2.
13-15. Do Problem 13-14 for the interacting thermal tanks
of Fig. 4-2.4. Obtain the responses for the following values
of the recycle flow /R: O, !, 5, 20, and 100 m3/min.
Discuss
what happens to the temperature responses as the recycle
flow increases.
13-16. Simulate the reactors of Problem 4-9 to obtain the
responses of the reactant concentrations in each tank to a step
change of 0.5 lbmole/ft3 in the feed concentration. Use the
parameters given for these reactors in Problem 6-15. Obtain
the responses for the following values of the recycle flow
/R: O, 10.0, 50.0, and 500 ft3/ntin. Discuss what happens
to
the concentration responses as the recycle flow is increased.
13-17. Simulate the extraction process of Problem 4-6 to
e,
obtain the response of the exit raffinate concentration (t)
to a step change of 2 m3/min in solvent rate /2(t). At the
[Link] steady-state conditions the inlet flow /1 is 5 m3/[Link]
and the solute concentration is 0.4 kmole/m3. It is desired to
recover 90 % of the solute. The total volume of the extractor
is 25 m3, the equilibrium slope is m = 3.95, and the mass
transfer coefficient is Ka= 3.646 min-1 •
 4 Chapter Dynamic Simulation of Control Problems

13-18. Simulate the temperature control loop for the stirred the cascade control loop of part (c). Obtain the responses of
tank cooler of Problem 4-7. Use the parameters given for the transmitter ourputs and of the signal to the control valve
this process in Problem 6-19. Simulate the controller as a for a step change of I OºF in ambient air temperature, and to
PID tuned for quarter decay ratio response and obtain the a step change of 0.3 mass % in set point.
responses of the transmitter and controller outputs to a step
change in the flow of cooling water of 0.02 m3 /min, and 13-26. Simulate the loop of Problem 9-3 to compare the
to responses of the simple feedback loop and the cascade loop
a set-point step change of 2°C. proposed in the problem. Obtain the responses of the trans-
mitter and controller outputs to unit step changes in distur-
13-19. Simulate the tank of Problem 4-5 and plot the
bance and in set point. Use a PI controller for the master
responses of the outlet temperature to step changes of 5 gpm
loop tuned for minimum lAE response.
in the process ftow and of 5 % in the signa! to the steam
control valve. 13-27. Simulate the jacketed reactor of Problem 9-4 to
compare the responses of the simple feedback loop and the
13-20. Modify the simulation of the control valve given
cascade loop using the jacket temperature as the interme-
in Section 13-4. l to calculate the flow of a gas in which
diate variable. Obtain the responses of the transmitter and
both the inlet pressure PI and the pressure drop across the
controller outputs to a step change of 5°C in set point. Use a
valve ó.pv can vary. Using the data for the steam valve of
PID controller for the reactor temperature tuned for minimum
Example 5-2.2, obtain a plot of the flow versus pressure drop
IAE response.
as the pressure drop is ramped from 5 to 25 psi. Keep the
inlet pressure and the valve position constant. 13-28. Show that when the inlet temperature disturbance
is neglected in Example 13-4.2, the feedforward controller
13-21. Simulate a simple temperature control loop for the
becomes a simple ratio controller. Modify the simulation
nonisothermal reactor of Section 4-2.3 that was sirnulated in
given in that example and compare the responses of the outlet
Problem 13-5. Use the data given for the control loop in
temperature to those given in the example for step changes
Problem 6-12, but assume the temperature transmitter has a
in both a process flow and inlet temperature.
time constant of 1.0 min, and the control valve has a time
constant of 0.1 min. Tune a PID controller for quarter decay 13-29. In the simulation of the blending tank in
ratio response and obtain the responses of the transmitter Problem 13-11, add a ratio controller between the concen-
and controller outputs to a step change of -0.2 ft3 /min in trated stream and the dilute stream and compare the responses
reactants flow and to a step change of 1 ºR in set point. to a step change in concentrated stream. Simulate a flow con-
trol loop on the dilute stream, instead of a control valve, as
13-22. In Example 13-4.1 the steam valve was simulated as
a constant pressure drop valve. However, as the temperature shown in Section 13--4. J.
of the condensing steam Ts varies, so does the pressure in the 13-30. The nonisothermal chemical reactor of Section 4-2.3
steam coil, so the pressure drop across the val ve is Jower the was simulated in Problem 13-21 with a simple feedback con-
higher the steam temperature. To simulare this effect we must troller. Using the data from that problem, simulate and tune
calculate the pressure in the coi) as a function of the tempera- a cascade control system for the reactor with the jacket out-
ture of the steam using the Antoine equation for water. Then let temperature Tc(t) as the intermediare variable. The jacket
the model of the val ve must be modified to calculate the ftow temperature transmitter has a time constant of 1.0 min and a
of steam w(t) using Eqs. 5-2.3 and 5-2.4. Assume the steam range of 560 to 660°R. Use a PI slave controller tuned by
is supplied at 30 psia and 250°F (saturated). Size the valve the synthesis formulas of Table 7-4.1 for 5 % overshoot on a
for 100% overcapacity and use C¡ = 0.9. Simulate the heater set-point change, and a master PID controller tuned for quar-
with the new valve and obtain the responses of the temper- ter decay ratio response. Obtain the closed-Ioop responses of
atures and the controller output to a step change in process the reactor temperature and the signal to the valve far a step
flow of -5 ft3 /min followed by a step of +5 ft3 /min. change of -0.2 ft3/min in reactants, and for a step change
Com- of 1 ºR in set point.
pare your responses to those of Fig. 13-4.1 O. 13·31. The nonisothermal chemical reactor of Section 4-2.3
13-23. Simulate the composition control loop of Problem was simulated in Problem 13-21 with a simple feedback con-
6-17 and obtain the responses of the transmitter and con- troller. Using the data from that problem, simulate and tune
troller outputs to a step change of 5 gpm in inlet flow. Use a feedforward control system for the reactor with the pro-
a PID controller tuned for quarter decay ratio response. cess flow fe(!) as the major disturbance, The process flow
13-24. Simulate the level and temperature control loops of transmitter has a time constant of 0.1 min and a range of O
Problem 6-24 and obtain the responses of the transmitter and to 2.0 ft3/min. Instead of the control valve on the coolant,
simulate a flow control loop (see Section 13-4.1). Obtain the
controller outputs for the temperature loop to a step change of
responses of the reactor temperature and the signa! to the
5 gpm in outlet oil ftow. Use a light leve! controller anda PLD
valve for a srep change of -0.2 ft3/min in reactant flow,
temperature controJJer tuned for quarter decay ratio response.
with and without feedback trim.
13-25. Simulate tbe drier of Problem 9-1 to compare the
responses of the simple feedback control loop of part (a) and
4 Chapter Dynamic Simulation of Control Problems

13-32. Do Problem 13-31 using a ratio controller with the to connect the controllers to the control valves. In Problem
feedback controller adjusting the ratio. 12-4 a general pairing strategy is developed for the evapora-
13-33. Design a cascade control scheme for the reactors in tor. You may want to pair the controller both ways and see
series of Problem 13-23 using the concentration out of the which pairing results in better control. Tune the composition
first reactor as the intennediate variable. Assume a transmit- controller for each pairing and plot the responses to a step
change of 1000 lb/h in steam flow.
ter range of O to 4 lb/gal of reactant with a time constant of
1.0 min. Toe slave controller is a PI controller tuned by the 13-36. Distillation column control.
synthesis formulas of Table 7-4.1. Compare the responses Use the simulation of the distillation column of
to those of the single feedback control loop obtained in Example 13-5.1 to study composition control. Assume PI
Problem 13-23. controllers tuned for minimum IAE. The composition trans-
13-34. Simulate the caustic blending tank of Problem 12-2 mitters have ranges of O to 1.0 mole fraction benzene, and
with PI controllers for both the product flow and composi- time constants of 1.0 min, and the reflux and steam flow con-
tion. Assume the tank has a constan! hold-up of 10,000 lb, trol loops have time constants of 0.1 min and have ranges
the flow transmitter has a negligible time lag and a range of O to 200 % of the respective design ñows. Obtain the
of O to 60 klb/h, and the composition transmitter has a time responses of the transmitter and controller ourputs to a step
constan! of 1.0 min and a range of O to 50 mass % NaOH. change of 200 lbmole/h in feed How. Do the problem for the
Toe flow control loops have ranges of O to 60 klb/h and time following cases:
constants of 0.1 min. Tune the controllers for quarter decay (a) Control only the distillate composition.
ratio response and obtain the responses to a l O klb/h step (b) Control only the bottoms composition.
change in product flow set point. Do the problem for
(c) Control both compositions.
(a) the correct pairing of the controlled and manipu-
lated variables, 13-37. Gasoline blending.
(b) the other pairing, Simulate the gasoline blending system of Example 12-2.5
(c) using a decoupler as the one in Example 12-3.1. with PI controllers or the product flow, octane, and Reid
vapor pressure (RVP). Assume constant densities and a con-
13-35. Simulation of an stant mixer volume of 450 barreis, instead of an on-line
Evaporator. mixer. The flow transmitter has a time constant of 0.1 rnin
The evaporator sketched in Fig. 12-3.4 is to concentrate a anda range of O to 100 kbl/day, and the analyzer transmit-
sugar solution by vaporizing the water using steam. Model ters have time constants of 1.0 min and ranges of 60 to 11 O
the evaporator to find the response of the level and product octane and O to 20 RVP. Toe flow control loops on the inlet
composition to changes in feed rate, feed composition, steam streams have time constants of 0.1 min and ranges of O to
flow, and product flow. Assume the contents of the evapora- 150 % of the design flows. Obtain the responses of the trans-
tor are perfeclly mixed and that the economy is constant. Toe mitter and controller outputs to a step change of I O kbl/day
economy is defined as the mass of water vaporized per unit in product flow set point. Do the following cases,
mass of steam supplied. By assuming a constant economy (a) Use the correct pairing.
you avoid having to write energy balances on the evaporator. (b) Use an altemate pairing.
Toe leve! is a linear function of the mass of solution in the (e) Use the decoupler of Example 12-3.3.
evaporator. At design conditions the feed is 50,000 lb/h with
a composition of 12 weight % sugars. Toe desired product 13-38. Three pump and tank problem.
composition is 70 weight % sugars. The economy is 0.95 lb This problem was created by chemical engineers at duPont in
of water evaporated per lb of steam. The evaporator is ini- the early I 960s, the heyday of process sirnulation. The tank
tially at steady state with the level at 50 % of the transmitter. shown in Figure Pl3-l is fed by three identical pumps and
When the level in the evaporator is at O % of the transmitter discharges through a very long pipe. Because of Lhe inertia
range the mass of solution in the evaporator is 394 lb, and of the liquid in the pipe, the tank leve) response is under-
when the leve) is at 50 %, the mass is 747 lb. Control sys- damped. This causes a problem because, although the tank is
tem: Simulate a proportional Jevel controller with a gain of taJJ enough at steady state, the overshoot in leve! can cause
20 %CO/%TO-tight leve! control-and a PI composition the tank to overflow when ali three pumps are tumed on.
controller with the feed and product flows as the manipulated Your assignment is to simulate the tank and determine if
variables. Toe control valves are to be linear with constant there is a time sequence for tuming the pumps on that pre-
pressure drop and sized for twice the design flow. Assume vents overflowing the tank. Each pump can be simulated as
negligible dynamics of the level transmitter, and a 0.1 min a step function with the step time being the time at which
time constant on the control valve actuators. Toe composi- the pump is tumed on. The step time for the first pump can
tion transmitter has a time constant of 0.6 min and a range be assumed to be 1 min, and the problem is to find the step
of 40 to 90 weight % sugars. Notice that there are two ways
 4 Chapter Dynamic Simulation of Control Problems

d,(t)

�i---t===========:�==�•r,c,i
h(t)

P1 P2 P3 .....����-L�����----+I1
Figure Pl3-l Schematic for Problem 13-38.
times for the other two pumps so that the tank leve! does not radiation from the sun is neglected. The overall heat transfer
exceed the height of the tank. A force balance on the exit coefficient between the flowing river and the barge is esti-
pipe results in the following equation: mated as 1000 kjoule/hr-m2-ºC. Assume the liquid chlorine
2
has a constant latent heat of vaporization of 288 kjoule/kg,
df. (t) rrd ( gal)
pL-- = _ 7.48 - pg[h(t) - 2.22 X I0-3 f0 (t) anda constan! specific heat of 0.946 kjoule/kg-ºC, and that
dt 4 ft
the pressure inside the barge is the vapor pressure of chlo-
-0.5184 X I0-6 f.;(t)] rine at the temperature of the chlorine in the barge. The How
where f0(t) is the flow out of the tank in gal/min, h(t) is of chlorine gas through the nozzle opening is given by the
the leve) in the tank in ft, L is the length of the pipe in ft, following equation:
dp is the inside diameter of the pipe in ft, p is the fluid 2
den- rrD
sity in lb/gal, and g is the acceleration of gravity in lf subcritical flow: Wc = 4/2Mw [P(P - P0)]
ft/min2. The term in brackets contains a correlation for the
friction
losses in the pipe in ft of head. Each pump has a capacity Wc = D2 {M;; p
v"'i
7r
lf critical flow:
of 4
750 gal/min, the tank is a cylindrical vertical tank 9.0 ft high
by 3.41 ft inside diameter. The pipe is 3170 ft long by 1.19 ft where Wc is the flow of chlorine gas through the opening, D
inside diameter. Simulate the system, obtain responses of the is the diameter of the hole, Mw is the molecular weight of
leve! in the tank and the outlet flow, and determine the time chlorine, R is the ideal gas constant, T is the absolute tern-
sequence for starting the pumps that prevents the leve! in the perature, P is the pressure in the tank, and P0 is atmospheric
tank from exceeding the height of the tank. pressure. Toe pressure P in the tank is the vapor pressure
13-39. Ecological interaction of host-parasite populations. of
chlorine that can be calculated from the temperature T using
A problem that illustrates the interaction between species in
the Antoine equation. Simulate the tank by writing the mass
an ecological system, also from the heyday of dynamic simu-
and enthalpy balances on the barge and determine the inicial
lation, considers only two species in the system, a host (e.g.,
and steady-state Hows of chlorine out of the barge. How long
rats) and parasite (e.g., fleas). Assume that the net growth rate
will it take for the barge to empty? Note: A simple way to
of rats-birth rate minus death rate-is 5 rat/rat-year, in the
simulate the critical/subcritical How is to calculate the flow
absence of fleas. The fleas infect the rats and increase their
using both formulas and select the smaller of the two flows
death rate by 0.05 rat/flea-rat-year. The fleas need the rats
using a minmax block or a min function.
to build their nests, and their growth rate is 0.2 flea/rat-flea-
year, while the normal death rate of fleas is 20 fleas/flea-year. 13-41. Control of semibatch reactor.
Set up the differential equations and determine the popula-
A semibatch reactor is to produce butyl propionate from
tion cycle for rats and fleas as a function of time. Assume
2-bucanol and propionic anhydride using a sulfuric acid cat-
there are initially 100 rats and 20 ñeas.
alyst. The reaction is given by
13-40. Environrnental impact of chlorine barge accident.
(C2HsC0)20 + 2C4H90H ---> 2C4H900C3Hs + H20
ln the summer of 1965 a barge loaded with liquid chlorine
sank in the Mississippi River near downtown Baton Rouge, Feliu et al. (2003) repon the following kinetic model for this
Louisiana, When they were preparing to bring it up, one of reaction:
the authors was asked to estimate the rates of release of chlo-
rine in the event one of the 4-in. nozzles on the barge were to
re = (kt CA + k2Ccm1 + k3Cco12)Cs
be knocked off. The barge originally contained 136,000 kg where r e is the rate of reaction of the anhydride, kmol/min-
of liquid chlorine at river ambient temperature, 25ºC. The m3, CA, Cs, Ccarl, and Cesa are the concentrations of the
area of contact between the tank and the ri ver is 212 m2 butano!, the anhydride, and two forms of the catalyst, in
and, for worse-case scenario, it is assumed that the nozzle is kmoUm3. A reaction is proposed for the conversion of one
knocked off when the barge is still submerged in the water fonn of the catalyst to the other, as follows
but the nozzle is jusr above the surface of the water. Heat of
Cat1 -+ Ca12
4 Chapter Dynamic Simulation of Control Problems

Catalyst

Hot water

Butanol

Cold water

Anhydride
Figure PlJ-2 Semibatch reactor for Problem 13.41.
with kinetics fed to the jacket to be able to heat and cool the reactor as
necessary. Assurne the reactor is initially filled with 600 kg
rcarl = k4lO-H'Ccat2Cs
of essentially pure 2-butanol at 20°C and the flows of cold
H, = -(p1 Cca,t + p2Ccat2) (p, + �) and hot water are adjusted to produce an inlet jacket tem-
perature of 30°C, which is also the initial temperature of the
where T is the reactor temperature in K. Toe heat of reac- jacket. The catalyst, consisting of 70 weight% 2-butanol,
tion is -80,000 kJ/kmol and the reaction rate coefficients are 20 weight % sulfuric acid, and the rest water, at 20ºC, is fed
given by the standard Arrhenius expression at the rate of 1200 kg/h for a period of 5 min. Then the
k; = A;e-E¡/RT anby- dride, assurned pure and also at 20°C, is added at a
rate that must not exceed 1200 kg/hr. An important
The parameters of the [Link] model are operational con- straint is that the reactor temperature must
not exceed 60ºC. Assume the heat transfer coefficient
between the reactor and the jacket is constant at 20 kJ/min-
Subscript i P; A; E;, m2-ºC, and the heat trans-
kJ/kmol 2
fer area is 15 m . The total volume of the reactor is 3 and
3 ,
1 0.2002 J.93 X lQtl 80,480 m
2 0.3205 J.01 X 1014 79,160 the jacket holds 415 kg of water. Get the molecular weights,
3 -21.38 J.42 X 1014 69,970 specific gravities, and specific heats from Perry's Chemical
4 1271 5.05 X 1011 76,620 Engineers' Handbook.
Simulate the reactor and design a control system and
operating procedure that minimizes the duration of the
batch cycle. Reference: Feliu, J. A., l. Grau, M. A. Alós, and
J. J. Macías-Hemández, "Match Your Process Constraints
Figure Pl3-2 shows a schematic of the jacketed reac- Using Dynamic Simulation." CEP, New York: AIChE:
tor. Hot and cold water, at 90 and ISºC, respectively, are December 2003, 42-48.