Mikles J., Fikar M. Process Modelling, Identification, and Control
Подождите немного. Документ загружается.
138 4 Dynamical Behaviour of Processes
0 500 1000 1500 2000
−1.5
−1
−0.5
0
t [s]
x
1
[d
m
]
nonlinear model
linearised model
Fig. 4.12. Response of the tank to step change of q
0
where
t
- time variable,
c
A
- molar concentration of A (mole/volume) in the outlet stream,
c
Av
- molar concentration of A (mole/volume) in the inlet stream,
V - reactor volume,
q - volumetric flow rate,
ϑ - temperature of reaction mixture,
ϑ
v
- temperature in the inlet stream,
F - heat transfer area,
α - overall heat transfer coefficient,
ρ - liquid density,
c
p
- liquid specific heat capacity,
H
r
- heat of reaction,
ϑ
c
- cooling temperature,
ϑ
cv
- cooling temperature in the inlet cooling stream,
q
c
- volumetric flow rate of coolant,
V
c
- coolant volume in the jacket,
ρ
c
- coolant density,
c
pc
- coolant specific heat capacity,
k
0
- frequency factor,
g - activation energy divided by the gas constant.
Compared to the CSTR model from page 31 we included the equa-
tion (4.69) describing the behaviour of cooling temperature ϑ
c
. The state-
space variables are c
A
,ϑ,ϑ
c
, the input variable is q
c
, and the output variable
is ϑ. We assume that the variables c
Av
,ϑ
v
,ϑ
cv
are held constant and the other
parameters of the reactor are also constant.
The reactor model can be rewritten as
4.2 Computer Simulations 139
dx
1
dt
=1− x
1
− f
1
(x
1
,x
2
) (4.71)
dx
2
dt
= b
1
− x
2
− f
2
(x
2
)+f
3
(x
3
)+f
1
(x
1
,x
2
) (4.72)
dx
3
dt
= f
4
(x
2
) −f
5
(x
3
) −f
6
(u
,x
3
)+f
7
(u
) (4.73)
k = k
0
exp
−
gρc
p
c
s
Av
H
r
x
2
(4.74)
where
x
1
=
c
A
c
s
Av
,x
2
=
ρc
p
ϑ
c
s
Av
H
r
x
3
=
ρ
c
c
pc
ϑ
c
c
s
Av
H
r
,t=
q
V
t
u
=
q
c
q
s
c
,b
1
=
ρc
p
ϑ
s
v
c
s
Av
H
r
f
1
(x
1
,x
2
)=
V
q
x
1
k
0
exp(−
gρc
p
c
s
Av
H
r
x
2
) ,f
2
(x
2
)=
αF
Vρc
p
x
2
f
3
(x
3
)=
αF
Vρ
c
c
pc
x
3
,f
4
(x
2
)=
αF V
qV
c
ρc
p
x
2
f
5
(x
3
)=
αF V
qV
c
ρ
c
c
p
x
3
,f
6
(u
,x
3
)=
Vq
s
c
qV
c
u
x
3
f
7
(u
)=
ρ
c
c
pc
Vq
s
c
ϑ
s
cv
c
s
Av
H
r
qV
c
and variables with superscript
s
denote steady-state values.
In the steady-state for state variables x
1
,x
2
,x
3
holds
0=1− x
s
1
− f
1
(x
s
1
,x
s
2
) (4.75)
0=b
1
− x
s
2
− f
2
(x
s
2
)+f
3
(x
s
3
)+f
1
(x
s
1
,x
s
2
) (4.76)
0=f
4
(x
s
2
) −f
5
(x
s
3
) −f
6
(u
s
,x
s
3
)+f
7
(u
s
) (4.77)
We define deviation variables and functions
x
1
= x
1
− x
s
1
(4.78)
x
2
= x
2
− x
s
2
(4.79)
x
3
= x
3
− x
s
3
(4.80)
u = u
− u
s
(4.81)
f
1
(x
1
,x
2
)=f
1
(x
1
+ x
s
1
,x
2
+ x
s
2
) −f
1
(x
s
1
,x
s
2
) (4.82)
f
2
(x
2
)=f
2
(x
2
+ x
s
2
) −f
2
(x
s
2
) (4.83)
f
3
(x
3
)=f
3
(x
3
+ x
s
3
) −f
3
(x
s
3
) (4.84)
f
4
(x
2
)=f
4
(x
2
+ x
s
2
) −f
4
(x
s
2
) (4.85)
f
5
(x
3
)=f
5
(x
3
+ x
s
3
) −f
5
(x
s
3
) (4.86)
f
6
(u, x
3
)=f
6
(u + u
s
1
,x
3
+ x
s
3
) −f
6
(u
s
,x
s
3
) (4.87)
f
7
(u)=f
7
(u + u
s
) −f
7
(u
s
) (4.88)
140 4 Dynamical Behaviour of Processes
The original mathematical model of the reactor given by (4.67) - (4.69)
can finally be rewritten as
dx
1
dt
= −x
1
− f
1
(x
1
,x
2
) (4.89)
dx
2
dt
= −x
2
− f
2
(x
2
)+f
3
(x
3
)+f
1
(x
1
,x
2
) (4.90)
dx
3
dt
= f
4
(x
2
) −f
5
(x
3
) −f
6
(u, x
3
)+f
7
(u) (4.91)
Figure 4.13 shows the Simulink block-scheme that uses the program rea7m1.m
www
(Program 4.5) as its S-function for the solution of the differential equa-
tions (4.89) - (4.91).
Program 4.5 (MATLAB file rea7m1.m)
function [sys,x0,str,ts] = rea7m1(t,x,u,flag)
% Non-linear deviation model of the CSTR;
% 3 equations, Reaction of the type A ----> B;
% k0 = a1 [1/min];
% g = a2 [K];
% ro = a3 [kg/m3];
% cp = a4 [kJ/kg.K];
% V = a5 [m3];
% Hr = a6 [kJ/kmol];
% F = a7 [m2];
% al = a8 [kJ/m2.min.K];
% roc = a9 [kg/m3];
% cpc = a10 [kJ/kg.K];
% Vc = a11 [m3];
% qs = a12 [m3/min];
% Thvs = a13 [K];
% cavs = a14 [kmol/m3];
% Thcvs = a15 [K];
% qcs = a16 [m3/min];
% cas = a17 [kmol/m3];
% Ths = a18 [K];
% Thcs = a19 [K];
% x(1) - concentration of A, dimensionless, deviation;
% x(2) - temperature of the mixture, dimensionless, deviation;
% x(3) - jacket temperature, dimensionless, deviation;
% u(1) - cooling flowrate, dimensionless, deviation;
% u=qc/qcs ;
% u(1)=<-0.2, 0.2> = cca<-0.001m3/min,0.001m3/min>;
4.2 Computer Simulations 141
switch flag,
case 0
[sys,x0,str,ts] = mdlInitializeSizes;
case 1,
sys = mdlDerivatives(t,x,u);
case 3,
sys = mdlOutputs(t,x,u);
case 9
sys = [];
end
function [sys,x0,str,ts] = mdlInitializeSizes()
sizes = simsizes;
sizes.NumContStates = 3;
sizes.NumDiscStates = 0;
sizes.NumOutputs = 6;
sizes.NumInputs = 1;
sizes.DirFeedthrough = 0;
sizes.NumSampleTimes = 1;
sys = simsizes(sizes);
x0 = [0; 0; 0];
str = [];
ts = [0 0];
function sys = mdlDerivatives(t,x,u)
% Calculation of f1c, f2c,...,f7c in the steady-state;
a1=10000000000000;
a2=11078;
a3=970;
a4=4.05;
a5=1.8;
a6=149280;
a7=5.04;
a8=130;
a9=998;
a10=4.18;
a11=0.7;
a12=0.25;
a13=370;
a14=0.88;
a15=285;
a16=0.05;
a17=0.0345;
142 4 Dynamical Behaviour of Processes
a18=385.877;
a19=361.51753;
x1cs = a17/a14;
x2cs = (a3*a4*a18)/(a14*a6);
x3cs = (a9*a10*a19)/(a14*a6);
u1cs = a16/a16;
f1cs = ((a5/a12)*x1cs*a1)*(exp((-a2*a3*a4)/(a14*a6*x2cs)));
f2cs = ((a7*a8)/(a12*a3*a4))*x2cs;
f3cs = ((a7*a8)/(a12*a9*a10))*x3cs;
f4cs = ((a7*a8*a5)/(a12*a11*a3*a4))*x2cs;
f5cs = ((a5*a7*a8)/(a12*a11*a9*a10))*x3cs;
f6cs = ((a5*a16)/(a12*a11))*u1cs*x3cs;
f7cs = ((a9*a10*a5*a16*a15)/(a14*a6*a12*a11))*u1cs;
x1c = x(1) + x1cs;
x2c = x(2) + x2cs;
x3c = x(3) + x3cs;
u1c = u(1) + u1cs;
f1c = ((a5/a12)*x1c*a1)*(exp((-a2*a3*a4)/(a14*a6*x2c)));
f2c = ((a7*a8)/(a12*a3*a4))*x2c;
f3c = ((a7*a8)/(a12*a9*a10))*x3c;
f4c = ((a7*a8*a5)/(a12*a11*a3*a4))*x2c;
f5c = ((a5*a7*a8)/(a12*a11*a9*a10))*x3c;
f6c = ((a5*a16)/(a12*a11))*u1c*x3c;
f7c = ((a9*a10*a5*a16*a15)/(a14*a6*a12*a11))*u1c;
f1 = f1c - f1cs;
f2 = f2c - f2cs;
f3 = f3c - f3cs;
f4 = f4c - f4cs;
f5 = f5c - f5cs;
f6 = f6c - f6cs;
f7 = f7c - f7cs;
% Right hand sides of ODEs;
sys(1)=-x(1) - f1;
sys(2)=-x(2) - f2 + f3 + f1;
sys(3)=f4 - f5 - f6 + f7;
function sys = mdlOutputs(t,x,u)
a3=970;
a4=4.05;
a6=149280;
a9=998;
4.2 Computer Simulations 143
a10=4.18;
a14=0.88;
sys(1)=x(1);
sys(4)=a14*x(1);
sys(2)=x(2);
sys(5)=(a14*a6*x(2))/(a3*a4);
sys(3)=x(3);
sys(6)=(a14*a6*x(3))/(a9*a10);
Step Input
Scope2
Scope1
Scope
S−Function
rea7m1
Fig. 4.13. Simulink block scheme for the nonlinear CSTR model
Linearised mathematical model of the CSTR with steady-state x
1
= x
2
=
x
3
= 0 is of the form
dx
1
dt
= −x
1
−
∂f
1
(0, 0)
∂x
1
x
1
−
∂f
1
(0, 0)
∂x
2
x
2
(4.92)
dx
2
dt
= −x
2
−
∂f
2
(0)
∂x
2
x
2
+
∂f
3
(0)
∂x
3
x
3
+
∂f
1
(0, 0)
∂x
1
x
1
+
∂f
1
(0, 0)
∂x
2
x
2
(4.93)
dx
3
dt
=
∂f
4
(0)
∂x
2
x
2
−
∂f
5
(0)
∂x
3
x
3
−
∂f
6
(0, 0)
∂x
3
x
3
−
∂f
6
(0, 0)
∂u
u +
∂f
7
(0)
∂u
u (4.94)
or
dx
dt
= Ax + Bu (4.95)
where
x =(x
1
,x
2
,x
3
)
T
A =
⎛
⎝
a
s11
a
s12
a
s13
a
s21
a
s22
a
s23
a
s31
a
s32
a
s33
⎞
⎠
, B =
⎛
⎝
0
0
b
s31
⎞
⎠
144 4 Dynamical Behaviour of Processes
a
s11
= −1 −
∂f
1
(0, 0)
∂x
1
,a
s12
= −
∂f
1
(0, 0)
∂x
2
a
s13
=0,a
s21
=
∂f
1
(0, 0)
∂x
1
a
s22
= −1 −
∂f
2
(0)
∂x
2
+
∂f
1
(0, 0)
∂x
2
,a
s23
=
∂f
3
(0)
∂x
3
a
s31
=0,a
s32
=
∂f
4
(0)
∂x
2
a
s33
= −
∂f
5
(0)
∂x
3
+
∂f
6
(0, 0)
∂x
3
,b
s31
=
∂f
7
(0)
∂u
−
∂f
6
(0, 0)
∂u
Figures 4.14, 4.15, 4.16 show responses of the CSTR to step change of q
c
(Δu = 10). Numerical values of all variables are given in Program 4.6.
File kolire8.m (Program 4.6) calculates for the given steady-state of the
CSTR parameters of the corresponding linearised mathematical model.
Program 4.6 (MATLAB file kolire8.m)
% Linearised model of the CSTR, 3 equations;
% state-space equations, Reaction A ----> B;
% k0 = a1 [1/min];
% g = a2 [K];
% ro = a3 [kg/m3];
% cp = a4 [kJ/kg.K];
% V = a5 [m3];
% Hr = a6 [kJ/kmol];
% F = a7 [m2];
% al = a8 [kJ/m2.min.K];
% roc = a9 [kg/m3];
% cpc = a10 [kJ/kg.K];
% Vc = a11 [m3];
% qs = a12 [m3/min];
% Thvs = a13 [K];
% cavs = a14 [kmol/m3];
% Thcvs = a15 [K];
% qcs = a16 [m3/min];
% cas = a17 [kmol/m3];
% Ths = a18 [K];
% Thcs = a19 [K];
% x(1) - concentration of A, dimensionless, deviation;
% x(2) - temperature of the mixture, dimensionless, deviation;
% x(3) - jacket temperature, dimensionless, deviation;
% u(1) - cooling flowrate, dimensionless, deviation;
% u=qc/qcs ;
4.3 Frequency Analysis 145
% u(1)=<-0.2, 0.2> = cca<-0.001m3/min,0.001m3/min>;
a1=10000000000000; a2=11078;
a3=970; a4=4.05;
a5=1.8; a6=149280;
a7=5.04; a8=130;
a9=998; a10=4.18;
a11=0.7; a12=0.25;
a13=370; a14=0.88;
a15=285; a16=0.05;
a17=0.0345; a18=385.877;
a19=361.51753;
% Calculation of the state-space coefficients;
as11=-1-((a5*a1)/a12)*(exp(-a2/a18));
as12=-((a5*a1*a6*a2*a17)/(a12*a3*a4*((a18)^2)))*(exp(-a2/a18));
as13=0;
as21=((a5*a1)/a12)*(exp(-a2/a18));
as22=-1-((a7*a8)/(a12*a3*a4));
as221=((a5*a1*a6*a2*a17)/(a12*a3*a4*((a18)^2)))*(exp(-a2/a18));
as22=as22 + as221;
as23=(a7*a8)/(a12*a9*a10);
as31=0;
as32=(a7*a8*a5)/(a12*a11*a3*a4);
as33=-(a7*a8*a5)/(a12*a11*a9*a10)-(a5*a16)/(a12*a11);
bs11=0;
bs21=0;
bs31=-(a5*a16*a9*a10*a19)/(a12*a11*a14*a6);
bs31=bs31+(a9*a10*a5*a16*a15)/(a14*a6*a12*a11);
4.3 Frequency Analysis
4.3.1 Response of the Heat Exchanger to Sinusoidal Input Signal
Consider a jacketed heat exchanger described by the differential equation (see
Example 3.13)
T
1
dy
dt
+ y = Z
1
u (4.96)
where y = ϑ − ϑ
s
is the deviation of the liquid temperature from its steady
state value, u = ϑ
p
− ϑ
s
p
is the deviation of the jacket temperature from its
146 4 Dynamical Behaviour of Processes
0 2 4 6 8 10
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
t
x
1
nonlinear model
linearised model
Fig. 4.14. Responses of dimensionless deviation output concentration x
1
to step
change of q
c
0 2 4 6 8 10
−3.5
−3
−2.5
−2
−1.5
−1
−0.5
0
t
x
2
nonlinear model
linearised model
Fig. 4.15. Responses of dimensionless deviation output temperature x
2
to step
change of q
c
steady state value, T
1
is the process time constant, and Z
1
is the process
steady-state gain.
The process transfer function of (4.96) is given as
G(s)=
Y (s)
U(s)
=
Z
1
T
1
s +1
(4.97)
Let the following sinusoidal signal with amplitude A
1
and frequency ω
influence the heat exchanger:
u(t)=A
1
sin ωt 1(t) (4.98)
4.3 Frequency Analysis 147
0 2 4 6 8 10
−9
−8
−7
−6
−5
−4
−3
−2
−1
0
t
x
3
nonlinear model
linearised model
Fig. 4.16. Responses of dimensionless deviation cooling temperature x
3
to step
change of q
c
Then the Laplace transform of the process output can be written as
Y (s)=
Z
1
T
1
s +1
A
1
ω
s
2
+ ω
2
(4.99)
Y (s)=
Z
1
A
1
ω
T
1
s
3
+ s
2
+ T
1
ω
2
s + ω
2
(4.100)
Y (s)=
Z
1
A
1
ω
T
1
s
3
+
1
T
1
s
2
+ ω
2
s +
ω
2
T
1
(4.101)
The denominator roots are:
s
1
= −jω, s
2
=jω, s
3
= −1/T
1
After taking the inverse Laplace transform the output of the process is
described by
y(t)=
Z
1
A
1
T
1
ω(−jω)
−jω[3(−jω)
2
+
2
T
1
(−jω)+ω
2
]
e
−jωt
+
+
Z
1
A
1
T
1
ω(jω)
jω[3(jω)
2
+
2
T
1
(jω)+ω
2
]
e
jωt
+
+
Z
1
A
1
T
1
ω(−
1
T
1
)
−
1
T
1
[2(−
1
T
1
)
2
+
2
T
1
(−
1
T
1
)+ω
2
]
e
−
t
T
1
(4.102)
y(t)=
Z
1
A
1
−2ωT
1
− 2j
e
−jωt
+
Z
1
A
1
−2ωT
1
+2j
e
jωt
+ Ke
−
t
T
1
(4.103)
where K =(2ωT
1
A
1
)/(2T
1
− 2+ω
2
T
2
1
).