- Email: [email protected]
The analytical modelling and dynamic behaviour of tray-type binary distillation columns
Mathl. Comput. Modelling Vol. 16, No. 2, pp. 105-129, Printed in Great Britain. All rights Reserved
THE
0895-7177192 $5.00 + 0.00 Copyright@1992 Pergamon Press plc
1992
ANALYTICAL MODELLING AND DYNAMIC BEHAVIOUR OF TRAY-TYPE BINARY DISTILLATION COLUMNS MOHAMMAD Department
H.
N.
TABRIZI
of Math/Computer
Science
East Carolina University Greenville, NC 27858, U.S.A. JOHN Department
B.
EDWARDS
of Mining, McGill University
3480 University St., Montreal Quebec, Canada H3A 2A7 and Canadian Centre for Automation (Received
and Robotics in Mining
September
1991)
Abstract A completely original, parametric transfer-function matrix (TFM) model for a binary distillation column, is derived. From the ordinary differential equations describing the individual trays of a multitray column, a partial differential equation (PDE) re p resentation is derived together with the necessary boundary conditions. Large-signal steady-state solutions are derived to provide parametric data for a small-signal PDE model obtained by linearization. The small-signal PDEs are solved analytically for sinusoidal inputs to produce the parametric TFM model for the system. The plant parameters and operating conditions are chosen to produce a physically symmetrical system which, in turn, yields a completely diagonal TFM choosing input and output vectors as the “sum and difference” of the circulating vapour and liquid flows within the column and the “difference and sum” of the top and bottom compositions, respectively. This deduction accords well with the findings of previous authors McMoran [l], Shinskey [2], and Stainthorp [3] based on part impirical, approximate analytical by Edwards and Jassim [4], or purely simulation methods [5]. Near firstorder-lag behaviour is also a feature common to the findings of these earlier investigations and the precise derivation herein reported.
1. INTRODUCTION
In earlier research, a number of important behavioural features have been noted in the simulation of distillation columns by the authors and by others including Armstrong and Wood [6], Radmaker and Rijndorp [7], th e most important being the near symmetrical responses obtained in the stripping and rectifying section to changes in both vapour rate V and distillate rate (D = V - R) where R = the reflux rate. Furthermore, in the approximate analysis carried out in previous studies [$I, it was realized that progress to a theoretical solution was only practicable if approximations were made in an attempt to produce symmetry of the system equations. Furthermore, the solutions seemed to be more readily obtained if the deviations E, and cS of the linearized equilibrium curve from the 45 degree line (Figure 1) were (a) small and (b) nearly equal. Despite the attempt to set up a perfectly symmetrical system in terms of operating conditions and geometry, the results obtained were still not completely symmetrical but nearly so. The simulation results did not agree precisely with prediction, but did so approximately. It is therefore encouraging to search for a completely symmetrical system which will hopefully then generate a completely symmetrical response and allow accurate analysis to be accomplished fairly easily. Symmetry would seem to be the key to initial success in making predictions of column behaviour. The problem is threefold:
105
106
M.H.N.
TABRIZI, J.B. EDWARDS
_.___-_----
--
true
Y h5”
equilibrium
x stripping section
I
\
I
‘\
I
/ ‘\I
I z
1.0
%I
Figure 1. Typical a-curve and its straight line approximation.
(i) In precisely (ii) What (iii) What
what sense should
are the precise dependent are the input
variables
the column variables
be symmetrical? between
which should
which symmetry
be set quiescently
should
be sought?
in a symmetrical
way?
These problems are answered by finding the “axes” of symmetry, since the column inputs from the end vessels (Figure 2) are L, and V,, the column should be operated with these nominally equal. Similarly, the outputs from the column V, and L, likewise should be operated at nominally equal values. This obviously implies that a mixed feed of liquid Fl and vapour F, should enter the column as shown in Figure 3, where Fj = F, nominally. With regard to the dependent variables, symmetry should be sought, between rectifier vapour and stripper liquid compositions Y(h) and X/(-h), respectively. More precisely, since X’(-h) is the deviation from its extreme value, zero, and 1 -Y(h) is th e d eviation of Y(h) from its extreme values, 1 .O, symmetry should be expected between X’( -h) and 1 - Y(h), (here X , X’ are liquid compositions in rectifier and stripper and Y, Y’ are vapour composition in rectifier and stripper). Returning to fluid flows, it is clear that for symmetry, if the top product is removed as liquid, D, the bottom product should be removed as vapour, W, and nominally D and W should be set equal. Finally, the column geometry should be symmetrical, i.e., an equal number of trays in rectifier and stripping section and the end vessels should be physically symmetrical. In addition to the column and its operating variables being symmetrical, piecewise linear approximation of the equilibrium curve in Figure 1 should be symmetrical. (Note from Figure 1 that the true, ideal equilibrium curve is symmetrical about the 45O line.) If the straight line slope is tan-’ o for the stripper, then for the rectifier it should be (1 - tan-’ o) as indicated. This produces a knee on the true axis of symmetry, i.e., on the 45’ line. This has implications for the feed mivtnre which &nnld he in m-tllilihrillm and in some sense svmmetrical. A little thought
Parametric model for distillation columns
107
reveals that system symmetry should be produced by setting the feed mixture compositions at the knee. Thus, if z and Z denotes the vapour and liquid feed compositions, respectively, then 2=1-z
and
z=CYz.
hccumul.ntor
Feed (vapour l&id)
h’ = 0
and
Stripping section
Bottom product l., - vs
4
4
“S
UUU
vapour stream
stream
Figure 2. Illustrating the complete system. WI 16:2-H
M.H.N.
108
vapour
TABRIZI,
J.B. EDWARDS
rectifier
tray 1
rectifier
tray 0
feed, Fv,s
Y'(O)
liquid feed, Fe,2
stripper
stripper
tray 0
tray -1
\XI(-1) Figure 3. Showing composition and flow variablesassociated with the feed section.
Additionally, if one input is to be a perturbation in D, i.e., now in V, - L, for symmetry of response, the other input should be V, + L, (also = V, + Ld). This, however, emerges from inspection of the model equations derived by the authors and now to be developed. Despite the difficulties of derivation, the need for parametric TFM remain for the following main reasons: (a) To test detailed
column
(b) As a base for control (c) To check coarse TFM
simulations
system
(which are fraught
with numerical
pitfalls).
design.
models obtained
theoretically
(d) To allow control design to influence plant ships between plant and TFM parameters.
or by experiment.
design by making
transparent
the relation-
In Section 2 of this paper, we develop the large signal partial differential equations (PDEs) for a long binary tray column from the discrete tray equations for mass and (implied) heat balance based on the usual idealizing assumptions, i.e., (a) Equal
latent
(b) Adiabatic
heats/mole/component,
conditions,
(c) Tray equilibrium, (d) Fixed molar holdup
(liquid),
(e) Zero vapour
holdup,
and
(f) A piecewise
linear equilibrium
curve.
109
Parametric model for distillation columns
The large scale model is solved for symmetrical operating conditions to yield an analytical steady-state composition profile which is used to generate the plant/model parameter relationships for the small perturbation model. The latter is derived in Section 3 by implicit differentiation of the large-signal PDE model and its derived boundary conditions for the feed- and end-vessels. The parametric transfer-function matrix (TFM) model between input changes to circulating flow (vapour and liquid) and terminal composition (top and bottom) is then found for a symmetrical plant by precise solution of the double Laplace transformed equations in the space-frequency domain. In Section 4 the general forms of the inverse Nyquist loci (generated by computation of the TFM elements) are examined and checked against analytically derived high and low frequency limits and against simulated steps responses. The TFM model was presented by the authors Edwards and Tabrizi [9] at the IFAC Symposium DYCORD 86 in Bournmouth, UK, in December 1986, but the limited length of a symposium paper precluded full and proper development and appraisal of such a model with its associated complexity of derivation. The paper was very well received, encouraging that the full paper is presented herein. 2. LARGE
2.1.
Vapour/Liquid
SIGNAL
MODEL
Equilibrium
The vapour above and liquid on each tray are assumed to be in equilibrium with one another respectively, at all times, so that if Y(n) and X(n) d enote the vapour and liquid mole-fraction, of the more volatile component on tray n of the rectifier, then for an ideal mixture: 1 - X(n) Y(n)X(n)
[l -Y(n)]
(1)
= @,
where p is the constant relative volatility of the mixture. Similarly, at the n th tray in the stripper: 1 -X’(n) Y’(n)X’(n)
[l - Y’(n)]
(2)
= p,
this relationship while the primes distinguishing variables in the strinning section. Linearizing rectifying section, and Y’ = 1 -X’, for retaining symmetry about the line Y = 1 - X,-for-the the stripping section, one obtains:
Y(n) Y’(n) where CYis a constant
= y
+ e,
= crX’(n),
where
Such linearization has often been used in earlier and Wilkinson and Armstrong [lo]. General
(4)
and is given by: o=l+c,
2.2.
(3)
c > 0.
dynamic
analysis
(5) by Armstrong
and Wood
[6]
Tray Equation
In common with others including Armstrong and Wood [6], Edwards and Jassim [4], Jawson of continuous and Smith [ll], Stainthorp and Searson [3], and Voetter [12], the assumption equilibrium effectively implies zero vapour capacitance (or infinite mass transfer rates between liquid and vapour phases). Thus, from Figure 4, the material balance for tray n of the rectifier may be written as HI
6/Q
-dX(n) = L, {X(n + 1) A
X(n)}
+ v, {Y(n - 1) -Y(n)},
M.H.N.
110
X(n+l)
TABRIZI, J.B. EDWARDS
-
tray
n+l
Y(n) -X(n) tray
---.
---_-
n Y (n-l)
(I
&
Lr --
x(n-1) -
-----
tray
I
I
I
I I Figure
4. Showing
composition
and flow variables
n-l
T.-.
-r
associated
’
I
with a general rectifier tray.
where Hl is rectifier liquid capacitance per unit length of column and 6h’ is total length of cell n, while L, and V, denote molar liquid and vapour flow rate within the rectifier’. By eliminating Y(n) in terms of X(n) using Equation (3), one can achieve: Hl
6h’
-dX(n) = L, dt
If X was a continuous
(X(72 + 1) - X(n)}
function
+ (%)
of height h’, by Taylor’s
X(n + 1) - X(n) = %SV+
$$$
{X(n theorem
(6)
- 1) - X(n)}. one can write:
(Sh’)2 + higher powers of bh’,
(7)
and X(n - I) - X(n) If n is sufficiently closely approximated by:
= -T
(Sh’)2 + higher powers of 6h’.
6h’ + $$$&
large, the truly discrete function X(n) may be, from Equations (6)-(8), by the continuous functions X(h’, t) so that the rectifier is now described
~ If the column
(8)
is operated
under
a particular
working
(6h’)2 + higher powers of bh’. condition,
where
v, = Lpa, then the odd powers of 6h’ disappear. second may be ignored, then in steady
(9)
PO)
Also, if h’ > 6h’ such that state
all powers of 6h’ above
the
2
(;/,;2 =‘I ‘Liquid and vapour flow rata, will be spatially components have equal latent heats per mole.
independent
(11) if operating
conditions
are adiabatic
and the two
Parametric model for distillation columne
111
thereby yielding a constant composition gradient through the rectifier (i.e., all trays perform the same duty: a good design criteria). Similarly, for the stripper, if HI, L, and V, denote liquid capacitance per unit length, liquid, and vapour molar flow rates (Figure 2), then msss balance considerations, together with Equation (5), give:
v =
H’6h’ 1
L, {X’(n
Again by using a continuous
H,! iih’ F
+ 1) -X’(n)}
approximation,
= (L, - Qr)
F
+ v,cr {X’(n + 1) - X’(n)}. X’(h’,t)
to X’(n,t),
we obtain:
6h’ (c5h’)2 + higher powers of 6h’,
and if the column
is operated
such that L, = V,a,
then in steady
(12)
(13)
state 82X’ (14
(ah/)2 = "
To solve general PDEs (9) and (12) or their special cases (11) and (14), it is necessary consider the four boundaries of the system, i.e., the terminating vessels and the feed trays. 2.3’.
Feed Boundary
to
Conditions
The feed section of the column is illustrated in Figure 3. Here F denotes the feed rate of the vapour and also of liquid, their composition being z and 2, respectively. The tray above the feed point is considered first. The mass balance is: H, 6h’ -dX(o) dt If for symmetry
the column
= Fz + V,Y’(O) - KY(O) is operated
+ L, {X(l)
- X(O)}.
such that
v, = L,
and
(16)
K = L, then using the relationships
(16), (17), (3), (lo),
(37) and (13), one obtains:
“(‘) = Fr f K {X’(O) - Y(O)} + V,Bh’
HI bh’ ax(o) at
(15)
ij/,,
(
(18)
where T
6h’ = X(1) - X(0).
Now F = V, - V, . Considering
the special
case, by substituting
the relations
(19) (10) and (17), one can obtain:
F=V,-L,=V, so that F=V-;.
(20)
M.H.N.
112
Furthermore,
J.B. EDWARDS
TABRIZI.
if the feed mixture is in equilibrium so that .z = oz,
(21)
and the composition such that the point Z, 2 lies on the -45O line of the vapour/liquid diagram for symmetry, i.e., r=l-2, (22) then it follows that ff
and
%=l$-cr Z=-. Substituting yields:
(23)
1
(24)
l+o
the relations (20) and (23) into (18) an d e 1iminating X(0) in favour of Y(0) therefore o HI 8Y(O) 6/$ = ---& K at
+ X’(0) -t (1 - Y(O)} + F
6/Z’.
(25)
Expressing normalized distance as h=-$
and normalized time as
then the Equation
(26)
tv, r=cuH16h’1
(27)
(25) reduces to the simplified form: 8Y(O) = &
--&
+X’(0) + (1 - Y(O)} + F.
Now, taking the tray immediately below the feed point and taking a mass balance, we obtain H’I Sh’ q
= F 2 +
and, again under symmetrical
H; ax’(o) v, dt
v, {Y’(-1)
- Y’(O)} + L,X(O)
- -&X’(O),
(29)
operation, Equation (29) reduces to: &,
-- 2
=
o+l
If, therefore, the tray capacitances
v
6h’ L (1 - Y(O)} - X’(0).
(30)
are such that H&Y&
(31)
then, in terms of normalized time and distance, Equation (30) can be reduced to a simpler form as
ax’(o)
-
Equations 2.4.
67
= -
2
(Y-i-1
- qp
- (1 -
Y(O)}-X’(O).
(32)
(28) and (32) clearly share a marked degree of symmetry.
Terminal
Boundary
Conditions
The variables associated with the accumulator and reboiler ends of the process are shown in Figure 5, the top tray being the Nth and N + 1 denoting accumulator quantities. If H, is the molar capacitance of the accumulator, then on taking a mass balance on this vessel, we obtain:
H
a
ax(;t+ ‘)=v,
{Y(N) - X(N
+ 1)).
113
Parametric model for distillation columns
Accumulator
(top
product) +
“r-L,
X(N)
Accumulator-end variables.
X' (-d,+ 6h’ )
(bottom product)
Reboiler
Heater
I
I
\
‘b
Reboiler-end variables.
Figure 5.
Regarding Y ss a spatially continuous variable, Y(N) using a truncated Taylor expansion to give: H 0 ax(;t+
l) = v,
{
Y(N+
1) -
can be eliminated in terms of Y(N -t 1)
OY(N + 1) 6h’ - X(N + 1) ._#
, >
so that, in terms of Y (IV + 1) only, we have H oBY(N+ 0 at
1)
Y(N+l)(l-o)+o--l-
(33)
And if Y(N + 1) is now replaced by Y(Ll), where L1 is the normalized length of the rectifier, then in terms of r and h we finally obtain the general result: T a-dY(Ll) A7
= E{l _ Y(&))
- 8y(Lr) --E---’
(34)
M.H.N.
114
TABRIZI,
J.B. EDWARDS
where T, is the normalized time-constant of the accumulator and is given by HCZ HI 6h’
Ta=-.
(35)
A similar treatment of the reboiler mass balance equation gives H*
dX’(-J52)
L, {X/(-L2
=
dt
+ 6h’) - Y’(-L2)},
where Hb is the molar capacitance of the reboiler and Ls, the normalized length of the strippingsection, produces the general results: T* 8X’(-L2)
8X’(-L2) dh
= -cX’(-L*)+
a7
)
(36)
where Ta is the normalized reboiler time constant given by
Ha
Tb=-. 2.5.
(37)
HI 6h’
Steady-State Solutions
General PDEs (9) and (12), under the particular operating conditions (10) and (13), yield the steady-state differential Equations (11) and (14). And, under the further operating conditions (16), (17), and (23), these must be solved subject to special feed boundary conditions (28) and (32), together with the general terminal conditions (34) and (36). Therefore in our symmetrical steady-state situation the system for solution is as follows: if then
Ls = L1 = L,
(36)
d2X’ &=
(39)
(dh>2 =01
X’(0) + 1 - Y(0) + 7
dY(0)
2 = o+l
(40)
dX’(0) X’(0) + 1 - Y(0) + --&-
2 = o+l’
(41)
and
(42)
7=&-Y(L)} dX’(-L)
= cX’(-L). dh As shown in the Appendix this system has the following solutions: X’(-h)
(43)
= 1 -Y(h),
(44)
Y(h) = Y(0) + G h,
(45)
where profile gradient G is constant and given by G=
(46)
(o+l)(2t22+(1+1)’
and initial values X(0) and Y(0) are given by 2(cL+l) (CY+ 1)(26 L + cr + 1)’
X’(0) = 1- Y(0) =
(47)
From (44) to (47) it is also readily shown that the terminal values X’(-L) by: X’(-L)
= 1 - Y(L)
=
Figure 6 illustrates the steady-state composition even tray loading: a desirable design condition.
(a + 1)(2A
+ (Y+ 1) *
profiles graphically.
and Y(L)
are given (48)
Their linearity indicates
Parametric model for distillation columns
115
: ’ (= 1 - Y)
____
T
--
__-------
I f I
2cL GL = (a+1) (zcL+a+l)
;
T t 2(1+ EL) I , (a+l) (2cL +a+l)
1 t GL=(a+l)
I i
$cL+a+l)
I
I I
1
:.
I
L-h
Figure 6. Steady-state composition profile.
3. SMALL 3.1.
Partial
Diferential
PERTURBATION
MODEL
Equations
The small perturbation PDEs for the system may be derived by implicit differentiation of the general large signal PDEs (9) and (12). The lower case symbols denote small changes in the variables. Substituting the steady state solutions for the upper-case symbols, one can obtain: H,iih’$-:G+
which in terms of normalized
time and distance BY -_= ar
Also from the special
symmetrical
$
cl
H;6h’F
r and y, rather I
v,
8% .
(49)
(W2
(50)
case (17), PDE (50) reduces
and output
than 2, becomes
=
va)G
67 if input
(Sh’)‘,
(lo - v)G
a;c’ = (I-
In particular,
&
+ a2x’ dhz’
v,
vectors are defined, ‘=
[
;;;:
to the normalized
form: (51)
thus
1
(52)
and (53)
Then the small signal PDEs (49) and (51) may be grouped
to form the simple matrix equation
as follows: (54)
M.H.N.
116
Laplace
(2 -p)&$)
-Z(O)
‘W denotes where superscript and “.“, the spatial derivative. Terminal
3.2.
Boundary
Differentiating reboiler we obtain
transforms
(33) implicitly
to obtain
the coefficient
treatment
1 ] 21,
to h and 7, “N” with respect
applies
boundary
-cy(~1)
6h’
dh’
of v in (56) is zero, so that (56) reduces
1
aY(Ll) ar
EY(L1)
a-=-
to the reboiler
+
F}
to r only
- p}
equation
for the
.
to the normalized
(56) form:
.
and yields the similar &‘(-L2) dh
= --m’(-L~)+
aT
Inverted
a;
the small-perturbation
BY(h)
T* al’(---La)
3.3.
[ -@0-l)
to r yields
Conditions
T
Similar
= s-1
with respect
Y(LI)c+c---
In steady-state
to h and in p with respect
(54) in s with respect
transforming
J.B. EDWARDS
TABRIZI,
(57)
result:
.
(58)
U-Tube Model
In the symmetrical case L1 = L2 = L, solutions of the small signal equations are considerably simplified by redefining h, h’ from the origin at the ends of the column rather than the central feed point, by first bending the column conceptually into an inverted U-tube as shown in Figure 7, so that co-ordinate h is replaced by L - h in the rectifier and h by L + h in the stripper. Oddpowered spatial derivatives in the rectifier are consequently reversed in sign, while even-powered derivatives in the rectifier and all derivatives in the stripping section are unaffected in sign. The terminal boundary conditions (57) and (58) thus become T
-ado) = -
a
T* a%‘(O) =
ar
+ y
(59)
and
ax'(o) --Ed(O)+ 7.
(60)
If T, = Tb = T,
(61)
then one can obtain TpC(O) = -q(o) -
+2(o).
(62)
The PDE (54) is unaffected by a change of distance base, since this involves only even-powered spatial derivatives so that its transformed version (55) is also unaffected provided z(O) and l(O) are now interpreted as terminal conditions rather than feed-point conditions. may therefore be eliminated from (55) using boundary conditions (62) to give
(s2-
p);- (s +
E+
Tp)i(O)
= s-l
[ -(ao-
‘)
(y ;
The unknown
a
1 ] 4,
after inversion now involving only the one unknown vector n(O), which may be eliminated substitution of the feed conditions not yet involved in this small signal analysis.
by
Parametric
model for distillation columns
117
Feed
h =L A
Stripper
Rectifier
T h
Reboiler
end
kcumulator
end
Figure 7. Inverted U-tube model.
3.4. Feed Boundary Condition By implicit differentiation of the large signal Equations (15) and (29) one can achieve the small signal equations for the feed boundary condition. The upper case symbols are replaced by the special case steady-state solutions. Retaining the original distance base for the moment, we thus obtain for the rectifier:
H,
6h’
-‘do) = v {Y’(O) at
Normalizing
-
(64) and defining aY(o> - 67 = +
By using the steady-state
Y'(O) Y(O)}+ v, yin y(0) rather
r
~(0)
1a’
ax(o) + 3
oh’
o!
6h’
‘do)
dh”
(64)
than r(O), one can achieve
Gv) - Y(O)}+ 1a g solutions
+
+ fgp
+ c’(0) - y(0).
(65)
P
of Section 2.5, Y’(0) - Y (0) can be eliminated,
thus giving
Y’(0) - Y(0) = &X’(O) + 1 - Y(0) - 1, so
Y’(0) - Y(0)
(a + 1)v 2
= -G
+crl
F +T+z’(O)-y(0).
>
(66)
118
M.H.N.TABRIZI, J.B.EDWARDS
An identical result as
treatment
applied
to the first tray of the stripping
ax’(o)_ 5 Transforming
(a + 111 -
_-av +
v, {
ar
2
(66) and (67) to the new distance aY(L) -=67
G
(a + l)v
1 -
v,
1
produced
the similar
+ Y(0)-z'(O).
base for the inverted
+ol
2
ax’(o)
r
section
(67)
U-tube
therefore
yields:
(68)
-Y.z’(L)-y(L), >
and ar’(L)=G a7 Adding yields:
and subtracting
&+2 0
A+ or in terms of Laplace
ax’(L) _crv + (o + I)! - T + Y(L)- x'(L). 2 v, { > the (68) and (69) and grouping the results in matrix
form therefore
0 (70)
$+7!i
transforms:
P+& 1@)=[Oir-0&+1)1k. 0
[ p+$+20 3.5.
(69)
(71)
Solutions
From the transformed
Inverting g(h)
Equation
from the s, p domain = { &cosh
&h
(63), one can now achieve
to the h, p domain
gives
C(O) ~rii
+ (6 + Tp) sinh ,/$h}
-(ao-
>[
‘)
’ a+1
1 f.
(72)
In order to substitute obtained by differentiating as’(L) = --$ {psinh &L ah
the boundary condition (71), one must obtain w (72) with respect to h and setting h = L, thus giving + ~;i;(c + up) cash ,@L} g(O) + sinh$L
‘)
O a+1
may be
] h, (73)
where G(O) @>
The elements of (74) are now calculated For ql(O) we have (,+2+ and, on substituting
[ -(aoe
which
for &(L)
and w,
(P + 2) {,/S-h {psinh &L
=
1
(74)
. [ Gz(O) individually from (71), (72) and (73).
-3
K(L)
= O.%iG,
we get &L
+ (6 + T’P) sinh JFLI
+ fi(c
{ 1 - coshJi?L}
+ TP) co& fiL} c Gil -
(5)
G(0) G(O)
sinh,/j?Lr
Gi = 0.5c zi,
Parametric model for distillation columns
119
[m
Im
0.4 lPo.91
w=o.
6
_.0.02
0.2
w=o.72
a
I I 0.4
I 0.2
54
b Re
! E - 0.1 L - 10 T = 1.0
I
I\!
’L,---0
2
I
4
6
8
Figure 8. Inverse Nyquist loci of gr’(O,
yielding
10
jw) for the tray column.
the final solution:
w=
6 {(P +
Ql
((1-t
‘J)(coshJiiL- 1)/p+ (sinhJiTL)/Jii+ conditions,
( > ,+;
so that on substituting
G(L)=
for &i(L) and w
+(a + l)(cosh&L
(75)
we obtain -0.5 (3a + 1) &,
we get
+ ,/Xc + TP) sinhJjjL
{p coshfiL
0.5)
{(p + 2)(c + Tp) + p} sinh@L/&’
TIP i- 2 + 6) coshfil+
For &(O), from the feed boundary
yielding
b Re
I
I
I
+ fisinh&L
+ (c + Z’p)coshJiiL}
- 1) ii2 + (CY+ l)(sinh&L)
$
@2(O)
= -0.5 (3a + 1) i?z,
the final solution:
m=_ z2
(o + l)(coshfiL
- 1) + (o + l)sinh@/,/F
{p( I+ T) + e} coshfiL
+ 0.5(3a + 1) + fi( 1+ c + Q)sinh&L *
(76)
Therefore, setting p = ju, Equations (75) and (76) may be used to compute the two Nyquist or Inverse Nyquist loci for this diagonal system. Despite the complexity of the hyperbolic functions involved, the general form of the loci may be deduced with some precision analytically by comiderine
the high md law frenuencv
behaviour
of these functions.
M.H.N.
120
0.4
0.5
TABRIZI, J.B. EDWARDS
0.3
I
0.2
II
0.1
1 I
I
I I
I
-
Re
I 0.18 --
0.28
0.1
--
0.38--
0.2
0.48 _.
0.58--
0.3
0.68--
0.78-m
0.88
--
0.99-. Figure 9. Inverse Nyquist loci of g,‘(O,
4.
The system
GENERAL
FORM
may be described
$&P) where transfer
function
matrix
=
NYQUIST
(77)
takes the diagonal
glc;pp) g2(; 1 p)
[
LOCI
form
= G(h,P)Z(h>P),
G(h,p)
(TFM)
G(hP)
function
0.5
jw) for the tray c&mm.
OF INVERSE
in transfer
0.4
form
’
’
or alternatively G(h,P) where the inverse
TFM G*(h,p)
= G’(h,P)pbP)l
is also diagonal, G’(h,p)
=
taking
(78) the form:
Cgi(;yp)gl(; p)
’
1
1
(79)
where d(hP)
= L7,‘hP)
and
&(hP)
= g,‘@,P).
(80)
Attention here will be confined to the system outputs and therefore to G*(O, p), the elements of which may be calculated directly from (75) and (76) by substituting the known plant parameters CY,L, E and T. Therefore setting p = jw, Equations (75) and (76) may b e used to compute the two Nyquist or Inverse Nyquist loci for this diagonal system. Typical locus shapes are given in Figures 8 and 9 for L = 10. E = 0.1. of= 1 + c) = 1.1 and T = 1.0.
Parametric model for distillation columns
121
4.0
of
model
3.5
3.0 Invalid
2.C
0.4 1.0.
0.2
I
c
ot1
I
I I
0.3
012
0
I
I I
.4
o!s
+E
Figure 10. Parameter space showing static perfo rmance loci and region of model validity.
4.1.
Zero-Frequency
By considering deduced that
Behaviour
the limits
as p +
0 of the right-hand
sides of (75) and
2CL + cx + 1
!Z(O,O) =
and
EL2 + CL + 0.5c
g2+(oYO) = -(o
+ l)L +6.5(3o
+ 1) *
(76), it is quickly
031) (82)
(These relations may be otained by alternatively solving the small signal PDEs (54), subject to the boundary conditions, with operator $$ set to zero.) L = actual length of rectifier (or stripper)/bh = number of trays per section which, of course, must be >> 1 for the spatial approximation to hold good. An important feature to note is that the static gain gr(0,0) for the long tray column can only be positive unlike packed columns [13,14], where similar analysis reveals the possibilitv of a negative value, depending on plant and operating parameters.
M.H.N.
122
TABRIZI, J.B. EDWARDS
P,(’ 3,( 1)
(E =0.2)
(lz=0.15) 3.0
t (c=O.l)
2.0 (E =0.05)
6
4
2
10
8
Figure 11. Comparing values of the static gain gl(O, 0) with those calculated from the simulated step responses.
4.2.
High-Frequency
Behaviour
As w is increased from zero, it becomes sions (75) and (76). C onsidering frequencies
necessary in band
to consider
the functions
of p in expres-
(83) which is a wide band,
recalling
that L > 1.0,
it is clear that hyperbolic
functions
coshfil
(84)
and sinh&iL
may be replaced
by 0.5 exp &L,
Jir = (o.5w)“.5(1 -l-j),
since (85)
and, therefore, has a positive real part allowing the terms in exp(-@L) to be neglected. Because of the upper bound on p set by (83), only the lowest powers of p outside the hyperbolic functions need be considered, so that from (75) we obtain: &(%P)
@lP
w*
2++2c/@
and C2(%P) 82(P)
the exponential terms clearly the lowest powers of p, since
-(a
~
+ I/&
+ I)(1 + l/l/T)
(87)
p(l+T)+c+&+c+Tp)’
cancelling.
However, it is insufficient
merely
to now ignore
c < 1.0, if (84) is to be satisfied.
as insnection
of the static
all but (88)
performance
curves
of Figure
10 reveals.
Parametric model for distillation cohmms
123
Rl(O,““)
0.2 --
0
0
0
0
0
0
A
A
A
A
A
0
0
0
0
0
.x
x
x
x
0
0 v
0 0.1 --
A
A
A
0
O
A A
0 0 0
x
x.X
x
X
Y I
I
I
I
4
2
I
t
6
8
wL
10
Figure 12. Comparing high frequency gains gl(O, oo) (analytical vdues) with those calculated from the simulated step responses.
Terms reduce
having to:
c as the coefficient Cl
[
are therefore
non-dominant
241 + (YIP
(0,P)lG (PI
T2iO,PV~2(P)
and (86) and (87) consequently
1 =
L-l <
[
-(I
+ a&p
(P".51< 1.0,
1 ’
(89)
(91) L-l <
Ip”.gl < 1.0.
(92) Thus the system approaches a purely integrating process above w = L-l. Above w = 1.0 the dynamic response changes, but of course all real processes exhibit highly complex behaviour if the frequency of excitation is raised sufficiently. However, the band L-’ < w < 1.0 is a wide band since the inverse gain moduli are of order: I!?i(O,0011= :
and
&(O, =))I = ;9
(93)
where L-l < w”.5 cg: 1.0, compared
with their approximate
zero-frequency
l~;(O,O)l = ;
and
(94)
values from (81)-(82): 19x0,
ON= &.
In approximating the high-frequency behaviour of the tray-column become apparent2, which is not surprising since this model is confined
(95) no more phenomena have to long columns and would
‘No wave phenomena at the ends of the cohunns ((LBcan occur in packed color [13] for instance.) have been revealed in this andvsis. A consideration of h # 0 might, however, produce such effects.
124
M.H.N.
TABRIZI,
EDWARDS
T
IR
A 10
J.B.
0.95
0.038
a
0.76
0.019
0.57
6
t------+
0'.2
0.4
Re
0.38
4
2
r-
I
I I
I
L_-_,
I
1
1 200
I
400
I
I
1
I 1000
I
a00
600
Figure 13. Inverse Nyquist loci of jT’(O,
b
Rl?
jw).
be expected to drastically attenuate any reflected waves well before their arrival at the top and bottom of the column. Furthermore, we note that the signs of the high and low frequency gains are identical in both cases; this would seem to validate the use of a multivariable first order lag approximation to tray column dynamics. Such a model is readily derived from the system PDEs without the necessity of tedious dynamic solution as is now demonstrated.
4.3.
Effect of Terminal
Capacitance
As it was shown in Section 4.2, T drops out of the high-frequency analysis and it follows that only the final portion of the step-response is affected by changes in T. This is confirmed by the simulation results of Figure 17.
5. VALIDATION
OF TFM
MODEL
The validity of the analytically derived TFM model is clearly confirmed by Figure 11, which from Equation (81) with simulated step compares value of the static gain, gi(O,O), calculated responses. Figure 12 makes a similar comparison of high frequency gain, the analytical values having been obtained from the first equation of set (89). It is clear that agreement is good, even for L down to as little as 3 units of normalized length, i.e., much better than we might have expected from the restriction L >> 1.0 assumed in the derivations. A typical locus obtained from to note the analytical model for gi(0,jw) (Equation (75)) is shown in Figure 8. It is interesting the nearly straight line nature of this locus in passing between the zero to high frequency region. Such loci closely resemble those of first order lags. In particular, assuming a first order lag fit:
Parametric
model for distillation columns
125
dRe
0.016
__
0.032-e
-- (I.01
-- (1.02
t3.03
0.049-.
t 0.066.-
0.04
0.05
Figure 14. Inverse Nyquist loci of gl’(O,jw)
for the tray column.
91(0,0) ih(O,P)
=
l+.$$$$P to the high and low frequency transfer function obtained from Equations (89) and (81), the straight locus is also presented in Figure 13. The similarity of the frequency calibration of Figure 15 compares the open-loop step responses of the sl’(O,G) and Kr(O,j w ) is reassuring. accurately simulated system with the first-order lag model derived from Equations (89) and (81), again confirming the theoretical analysis. First-order lag approximation of Owens [15] may in fact be obtained much more readily directly from the PDEs by setting g = 0 and 00, and form an excellent basis for controller design for tray-columns as has been demonstrated by Tabrizi [8]. It is interesting to note that, for the small terminal capacitance, the frequency calibrations on the inverse Nyquist locus (Figures 8 and 14) change nearly proportionally to frequency. Multiplying w by a high-frequency gain = 10 (from 91) on Figure 8 yields the length of the imaginary component. Similar multiplication by the high frequency gain of 0.5 (again from 91) on Figure 14 does likewise. This confirms that the model with small terminal hold-ups is approximately of the multivariable first-order lag type having the high frequency gains predicted. The theoretical static gains are calculated from (81) and (82) to be:
From time response
of Figures
gi(0) = 2.69,
(96)
g2(0) = -231.5.
(97)
15 and 16 the steady-state
gains are:
gl(0) = 2.7,
(98)
g2(0) = -230.
(99)
126
M.H.N.
TABRIZI, J.B.
Full process simulation First order lag model
0000
E =
EDWARDS
0.1
L = 10 T = 1.0
1
Figure 15. Unit step response of g~(O,p) for the tray column.
From frequency
responses
Figures
8 and 14, the static
gains are:
gr(0) = 2.63,
(100)
gs(0) = -230.
(101)
It is clear that the static gains (98) and (99) compare favourably with (100) and (101) and also there is good agreement with theoretical values (96) and (97). 45’ lines intersect the loci of 13 and 14 at w = 0.035 and 0.007, respectively, yielding predicted time constants of 25 and 143. These values compare well with figures of 27 and 150 obtained from the simulated step responses of Figures 15 and 16, again confirming model validity. 6. CONCLUSIONS
A parametric TFM model has been derived completely analytically for long, symmetricallyoperated, tray-type distillation columns separating binary mixtures. The system has been shown to enjoy a completely diagonal structure3 provided the selected outputs are y(h, 0) - z’(h, 0) and y(h, 0) +z’(h, 0) and selected inputs are v + 1 and v - 1, where y denotes the change in top vapour composition, 2’ that in bottom liquid composition and v and 1 are perturbations in the flow rates of vapour and reflux, respectively. Increase of v + 1 drives y - I’ positively while v - 1 drives y + x’ negatively and with much greater gain. This accords with common experience. The gains of both elements of the TFM are shown to have the same sign at high and low frequency. Unlike the packed columns [13,14,16], non-minimum phase effects are not anticipated in tray-column behaviour, thus permitting confidence in the general application of a multivariable first-order 3The deduction accords well with the findings of previous authors based on part imperical McMoran [l], Shinskey [2], and Stainthorp [3] approximate analytical by Edwards and Jassim [4], or purely simulation methods [5].
Parametric model for distillation col-a
127
-100 __,
g*(O,P)
E - 0.1 L = 10 T = 1.0
--
-x0--
-200--
T Figure 16. Unit step response of 92(0, p) for the tray column.
E = 0.1 L = 10
,
I
I
t
,
b 2000 Figure 17. Unit step response of the long tray column and the effect of terminal capacitance. 200
600
1000
I’ 1400
1800
128
M.H.N.
TABRIZI, J.B.
EDWARDS
lag approximation of Owens [15] to long tray-type columns. The analysis has indicated the absence of significant effects from the reflection of travelling waves in the tray-type column on the composition changes at the ends of the column. However, such effects might nevertheless occur closer to the feed tray, and this possibility should be investigated since control measurements are frequently taken from points well away from the end of distillation columns. Since it has relied on the approximation of discretely-changing spatial functions by approximate continuous functions, the analysis here reported is, of course, restricted strictly to columns of many trays separating mixtures of low relative volatility. Predicted and measured response parameters seem to accord closely however down to normalized column length = 3. The substantially vertical shape calibration and position of the inverse Nyquist loci have been confirmed by low-and (fairly) highfrequency analysis and by simulated step response tests, provided terminal capacitances are small. Large terminal capacitances cause a slow tail on the step responses and the frequency calibrations of the loci are distorted accordingly. One extension of the research reported in this paper was carried out by Tabrizi and Edwards [17] w h ere it involved the reduced order modelling of large terminal capacitance columns. This revealed that a lumped, second-order lead/lag structure could closely approximate the dynamic response predicted (a) by the detailed TFM derived in the present paper and (b) by simulation. The same structure also carries over to packed columns where it also allows the possibility of non-minimum phase behaviour. At very high frequency, the loci depart from the vertical, the step responses being strictly of the Bessel function rather than exponential type. Further work is needed in this area to predict more accurately the extreme initial and final behaviour of long tray columns by simple formulae. The precise and original TFM model derived in this paper will provide a sound test for future simple models (Bessel function ielated or otherwise) that might emerge in future research.
REFERENCES
1. P.D. McMoran, Application of the inverse Nyquist method lJ.h’.A.C. Control Convention, page 122, (1971). 2. F.G. Shinskey, Process 3. F.P. Stainthorp
Control System, McGraw-Hill,
and H.M. Searson,
Tran~. Inst.
4. J.B. Edwards and H.J. Jassim, An analytical Inst. Chem. Eng., 55 (l), 17-29, (1977). 5. H.H. Rosenbrock and C. Storey, Computational (1966).
of 4fh
to a distillation column model, In Proc.
New York, (1967).
Chem.
Eng., 51 (l), 42-50,
(1973).
study of the dynamics of binary distillation Techniques
for Chemical
columns,
Trans.
Engineers, Pergamon Press, London,
6. W.D. Armstrong and R.M. Wood, An introduction to the theoretical evaluation of the frequency response of a distillation column to a change in reflux flow rate, Trans. Inst. Chem. Engrs., 39 (2), 80-90 (1961). 7. 0. R&maker, J.E. Rijindorp and A. Mearleved, Elsevier, Amsterdam, (1975). a. M.H.N. Tab&i, U.K., (1982).
Dynamics
and Conirol
of Coniinuous
Dislillation
Units,
The dynamics and control of binary distillation columns, Ph.D. Thesis, Sheffield University,
9. J.B. Edwards and M.H.N. Tabrizi, Parametric transfer function for binary tray-distillation columns, In Proc. of IFAC Symposium on Dynamics and Control of Chemical Reactors and Distillation Columns, (DYCORD-86), pp.lSS-194, B ournmouth, U.K., (Dec. 1986). 10.
W.L. Wilkinson and W.D. Armstrong, An approximate method fractionating column, Chem. Eng. Sci. 7 (l/2), l-7 (1957).
11.
M.A. Jawson and W. Smith, Counter transfer processes in non-steady A225, page 226, London, (1965).
12. H. Voetter, 13.
Plant
and Process
Dynamics
Characteristics,
of predicting
Butterworth,
composition
state, In PTOC. Royal
response Sot.,
of a
(London),
London, (1957).
M.H.N. Tabrizi, Towards an insight into the dynamic behaviour of tray distillation columns through analysis and simulation of dynamically asymmetric packed columns, Research Report No. 113, Sheffield University, U.K., (1979).
14. J.B. Edwards and M. Guilandoust, Parametric transfer function for packed distillation columns, In Proc. of IFAC Symposium on Dynamica and Conlrol of Chemical Reactors and Diatilllalion Columns, (DYCORD-86), pp. 195-202, Bo urnmouth, U.K., (1986). 15. D.H. Owens, First- and second-order-like IEE, 122, (S), 935-941, (1975).
structures
in linear multivariable
control system design, In Proc.
129
Parametric model for distillation columns
16. M.H.N. Tabrizi, An investigation into the dynamics of the distillation column, search for a wave effects, In Proc. IEEE Southeastcon 90, pp. 599-602 (1990). 17. M.H.N. Tabriti and J.B. Edwards, Comparison of lumped and distributed parameter model performance, pp. 305-308; Modelling for controller design in tray distillation development of a parametric, lumped parameter, transfer function model, pp. 301-304, In Proc. IEEE Soufheasfcon-88, (1988). APPENDIX Of Steady-State
Calculation
Composition Profile.
The large-signal steady-state system is described by the differential equations: d2Y
d2X’
(A.11
(dh)2=(dh)2=Ov which are subject to the boundary conditions
:
X’(0) + 1 -Y(O) x’(o)+1
+ y
-Y(o)+
y dY(L)
-
dh dX’( L) dh From the system symmetry
= -&
(A-2)
= --&,
(A-3)
= e{l-Y(L)},
(A.4)
= cX’(-L).
(A.5)
it is clear that
X’(-h)
= 1 -Y(h),
(A.61
and from (A.l) z
= G =
dh
dx’
= constant.
G’
-=
constant,
dh
(A-7)
From (A.6) it is obvious that G=G’.
(A.6)
Hence,
Y(h)
= Y(0) + G h, = X’(0)
X’(h)
+ G h.
(A-9)
Now from (A.2) or (A.3) and (A.6) we get 2X’(O) + G = 2 a+1
(A.lO)
and from (A.4) or (A.5) and (A.9) G = cX’(-L)
= c {X’(O)
- GL} ,
therefore, X’(0) Hence, eliminating X’(0)
between (A.lO)
G(cL+l)
=
and (A.ll),
2G(cL+l)
we let +G=
2 -I a+1
c G= stated as Equation
2c
(A.12)
a+1(2eLtatl)’
(46) in the main text. Hence from (A.ll)
X’(0) =
(A.ll)
2(eL t 1) a+l(IcLtatl)
and (A.12)
This is stated as (48) in the main text.
= 1 -Y(L)
=
(A.13)
= 1 -Y(O).
This is stated as Equation (47) in the main text. X’( -L) and Y(L) may now be calculated from general Equations (A.9)
X’(-L)
one can get:
to give:
2 a + 1 (2eL + (Y t 1)’
(A.14)
THE
0895-7177192 $5.00 + 0.00 Copyright@1992 Pergamon Press plc
1992
ANALYTICAL MODELLING AND DYNAMIC BEHAVIOUR OF TRAY-TYPE BINARY DISTILLATION COLUMNS MOHAMMAD Department
H.
N.
TABRIZI
of Math/Computer
Science
East Carolina University Greenville, NC 27858, U.S.A. JOHN Department
B.
EDWARDS
of Mining, McGill University
3480 University St., Montreal Quebec, Canada H3A 2A7 and Canadian Centre for Automation (Received
and Robotics in Mining
September
1991)
Abstract A completely original, parametric transfer-function matrix (TFM) model for a binary distillation column, is derived. From the ordinary differential equations describing the individual trays of a multitray column, a partial differential equation (PDE) re p resentation is derived together with the necessary boundary conditions. Large-signal steady-state solutions are derived to provide parametric data for a small-signal PDE model obtained by linearization. The small-signal PDEs are solved analytically for sinusoidal inputs to produce the parametric TFM model for the system. The plant parameters and operating conditions are chosen to produce a physically symmetrical system which, in turn, yields a completely diagonal TFM choosing input and output vectors as the “sum and difference” of the circulating vapour and liquid flows within the column and the “difference and sum” of the top and bottom compositions, respectively. This deduction accords well with the findings of previous authors McMoran [l], Shinskey [2], and Stainthorp [3] based on part impirical, approximate analytical by Edwards and Jassim [4], or purely simulation methods [5]. Near firstorder-lag behaviour is also a feature common to the findings of these earlier investigations and the precise derivation herein reported.
1. INTRODUCTION
In earlier research, a number of important behavioural features have been noted in the simulation of distillation columns by the authors and by others including Armstrong and Wood [6], Radmaker and Rijndorp [7], th e most important being the near symmetrical responses obtained in the stripping and rectifying section to changes in both vapour rate V and distillate rate (D = V - R) where R = the reflux rate. Furthermore, in the approximate analysis carried out in previous studies [$I, it was realized that progress to a theoretical solution was only practicable if approximations were made in an attempt to produce symmetry of the system equations. Furthermore, the solutions seemed to be more readily obtained if the deviations E, and cS of the linearized equilibrium curve from the 45 degree line (Figure 1) were (a) small and (b) nearly equal. Despite the attempt to set up a perfectly symmetrical system in terms of operating conditions and geometry, the results obtained were still not completely symmetrical but nearly so. The simulation results did not agree precisely with prediction, but did so approximately. It is therefore encouraging to search for a completely symmetrical system which will hopefully then generate a completely symmetrical response and allow accurate analysis to be accomplished fairly easily. Symmetry would seem to be the key to initial success in making predictions of column behaviour. The problem is threefold:
105
106
M.H.N.
TABRIZI, J.B. EDWARDS
_.___-_----
--
true
Y h5”
equilibrium
x stripping section
I
\
I
‘\
I
/ ‘\I
I z
1.0
%I
Figure 1. Typical a-curve and its straight line approximation.
(i) In precisely (ii) What (iii) What
what sense should
are the precise dependent are the input
variables
the column variables
be symmetrical? between
which should
which symmetry
be set quiescently
should
be sought?
in a symmetrical
way?
These problems are answered by finding the “axes” of symmetry, since the column inputs from the end vessels (Figure 2) are L, and V,, the column should be operated with these nominally equal. Similarly, the outputs from the column V, and L, likewise should be operated at nominally equal values. This obviously implies that a mixed feed of liquid Fl and vapour F, should enter the column as shown in Figure 3, where Fj = F, nominally. With regard to the dependent variables, symmetry should be sought, between rectifier vapour and stripper liquid compositions Y(h) and X/(-h), respectively. More precisely, since X’(-h) is the deviation from its extreme value, zero, and 1 -Y(h) is th e d eviation of Y(h) from its extreme values, 1 .O, symmetry should be expected between X’( -h) and 1 - Y(h), (here X , X’ are liquid compositions in rectifier and stripper and Y, Y’ are vapour composition in rectifier and stripper). Returning to fluid flows, it is clear that for symmetry, if the top product is removed as liquid, D, the bottom product should be removed as vapour, W, and nominally D and W should be set equal. Finally, the column geometry should be symmetrical, i.e., an equal number of trays in rectifier and stripping section and the end vessels should be physically symmetrical. In addition to the column and its operating variables being symmetrical, piecewise linear approximation of the equilibrium curve in Figure 1 should be symmetrical. (Note from Figure 1 that the true, ideal equilibrium curve is symmetrical about the 45O line.) If the straight line slope is tan-’ o for the stripper, then for the rectifier it should be (1 - tan-’ o) as indicated. This produces a knee on the true axis of symmetry, i.e., on the 45’ line. This has implications for the feed mivtnre which &nnld he in m-tllilihrillm and in some sense svmmetrical. A little thought
Parametric model for distillation columns
107
reveals that system symmetry should be produced by setting the feed mixture compositions at the knee. Thus, if z and Z denotes the vapour and liquid feed compositions, respectively, then 2=1-z
and
z=CYz.
hccumul.ntor
Feed (vapour l&id)
h’ = 0
and
Stripping section
Bottom product l., - vs
4
4
“S
UUU
vapour stream
stream
Figure 2. Illustrating the complete system. WI 16:2-H
M.H.N.
108
vapour
TABRIZI,
J.B. EDWARDS
rectifier
tray 1
rectifier
tray 0
feed, Fv,s
Y'(O)
liquid feed, Fe,2
stripper
stripper
tray 0
tray -1
\XI(-1) Figure 3. Showing composition and flow variablesassociated with the feed section.
Additionally, if one input is to be a perturbation in D, i.e., now in V, - L, for symmetry of response, the other input should be V, + L, (also = V, + Ld). This, however, emerges from inspection of the model equations derived by the authors and now to be developed. Despite the difficulties of derivation, the need for parametric TFM remain for the following main reasons: (a) To test detailed
column
(b) As a base for control (c) To check coarse TFM
simulations
system
(which are fraught
with numerical
pitfalls).
design.
models obtained
theoretically
(d) To allow control design to influence plant ships between plant and TFM parameters.
or by experiment.
design by making
transparent
the relation-
In Section 2 of this paper, we develop the large signal partial differential equations (PDEs) for a long binary tray column from the discrete tray equations for mass and (implied) heat balance based on the usual idealizing assumptions, i.e., (a) Equal
latent
(b) Adiabatic
heats/mole/component,
conditions,
(c) Tray equilibrium, (d) Fixed molar holdup
(liquid),
(e) Zero vapour
holdup,
and
(f) A piecewise
linear equilibrium
curve.
109
Parametric model for distillation columns
The large scale model is solved for symmetrical operating conditions to yield an analytical steady-state composition profile which is used to generate the plant/model parameter relationships for the small perturbation model. The latter is derived in Section 3 by implicit differentiation of the large-signal PDE model and its derived boundary conditions for the feed- and end-vessels. The parametric transfer-function matrix (TFM) model between input changes to circulating flow (vapour and liquid) and terminal composition (top and bottom) is then found for a symmetrical plant by precise solution of the double Laplace transformed equations in the space-frequency domain. In Section 4 the general forms of the inverse Nyquist loci (generated by computation of the TFM elements) are examined and checked against analytically derived high and low frequency limits and against simulated steps responses. The TFM model was presented by the authors Edwards and Tabrizi [9] at the IFAC Symposium DYCORD 86 in Bournmouth, UK, in December 1986, but the limited length of a symposium paper precluded full and proper development and appraisal of such a model with its associated complexity of derivation. The paper was very well received, encouraging that the full paper is presented herein. 2. LARGE
2.1.
Vapour/Liquid
SIGNAL
MODEL
Equilibrium
The vapour above and liquid on each tray are assumed to be in equilibrium with one another respectively, at all times, so that if Y(n) and X(n) d enote the vapour and liquid mole-fraction, of the more volatile component on tray n of the rectifier, then for an ideal mixture: 1 - X(n) Y(n)X(n)
[l -Y(n)]
(1)
= @,
where p is the constant relative volatility of the mixture. Similarly, at the n th tray in the stripper: 1 -X’(n) Y’(n)X’(n)
[l - Y’(n)]
(2)
= p,
this relationship while the primes distinguishing variables in the strinning section. Linearizing rectifying section, and Y’ = 1 -X’, for retaining symmetry about the line Y = 1 - X,-for-the the stripping section, one obtains:
Y(n) Y’(n) where CYis a constant
= y
+ e,
= crX’(n),
where
Such linearization has often been used in earlier and Wilkinson and Armstrong [lo]. General
(4)
and is given by: o=l+c,
2.2.
(3)
c > 0.
dynamic
analysis
(5) by Armstrong
and Wood
[6]
Tray Equation
In common with others including Armstrong and Wood [6], Edwards and Jassim [4], Jawson of continuous and Smith [ll], Stainthorp and Searson [3], and Voetter [12], the assumption equilibrium effectively implies zero vapour capacitance (or infinite mass transfer rates between liquid and vapour phases). Thus, from Figure 4, the material balance for tray n of the rectifier may be written as HI
6/Q
-dX(n) = L, {X(n + 1) A
X(n)}
+ v, {Y(n - 1) -Y(n)},
M.H.N.
110
X(n+l)
TABRIZI, J.B. EDWARDS
-
tray
n+l
Y(n) -X(n) tray
---.
---_-
n Y (n-l)
(I
&
Lr --
x(n-1) -
-----
tray
I
I
I
I I Figure
4. Showing
composition
and flow variables
n-l
T.-.
-r
associated
’
I
with a general rectifier tray.
where Hl is rectifier liquid capacitance per unit length of column and 6h’ is total length of cell n, while L, and V, denote molar liquid and vapour flow rate within the rectifier’. By eliminating Y(n) in terms of X(n) using Equation (3), one can achieve: Hl
6h’
-dX(n) = L, dt
If X was a continuous
(X(72 + 1) - X(n)}
function
+ (%)
of height h’, by Taylor’s
X(n + 1) - X(n) = %SV+
$$$
{X(n theorem
(6)
- 1) - X(n)}. one can write:
(Sh’)2 + higher powers of bh’,
(7)
and X(n - I) - X(n) If n is sufficiently closely approximated by:
= -T
(Sh’)2 + higher powers of 6h’.
6h’ + $$$&
large, the truly discrete function X(n) may be, from Equations (6)-(8), by the continuous functions X(h’, t) so that the rectifier is now described
~ If the column
(8)
is operated
under
a particular
working
(6h’)2 + higher powers of bh’. condition,
where
v, = Lpa, then the odd powers of 6h’ disappear. second may be ignored, then in steady
(9)
PO)
Also, if h’ > 6h’ such that state
all powers of 6h’ above
the
2
(;/,;2 =‘I ‘Liquid and vapour flow rata, will be spatially components have equal latent heats per mole.
independent
(11) if operating
conditions
are adiabatic
and the two
Parametric model for distillation columne
111
thereby yielding a constant composition gradient through the rectifier (i.e., all trays perform the same duty: a good design criteria). Similarly, for the stripper, if HI, L, and V, denote liquid capacitance per unit length, liquid, and vapour molar flow rates (Figure 2), then msss balance considerations, together with Equation (5), give:
v =
H’6h’ 1
L, {X’(n
Again by using a continuous
H,! iih’ F
+ 1) -X’(n)}
approximation,
= (L, - Qr)
F
+ v,cr {X’(n + 1) - X’(n)}. X’(h’,t)
to X’(n,t),
we obtain:
6h’ (c5h’)2 + higher powers of 6h’,
and if the column
is operated
such that L, = V,a,
then in steady
(12)
(13)
state 82X’ (14
(ah/)2 = "
To solve general PDEs (9) and (12) or their special cases (11) and (14), it is necessary consider the four boundaries of the system, i.e., the terminating vessels and the feed trays. 2.3’.
Feed Boundary
to
Conditions
The feed section of the column is illustrated in Figure 3. Here F denotes the feed rate of the vapour and also of liquid, their composition being z and 2, respectively. The tray above the feed point is considered first. The mass balance is: H, 6h’ -dX(o) dt If for symmetry
the column
= Fz + V,Y’(O) - KY(O) is operated
+ L, {X(l)
- X(O)}.
such that
v, = L,
and
(16)
K = L, then using the relationships
(16), (17), (3), (lo),
(37) and (13), one obtains:
“(‘) = Fr f K {X’(O) - Y(O)} + V,Bh’
HI bh’ ax(o) at
(15)
ij/,,
(
(18)
where T
6h’ = X(1) - X(0).
Now F = V, - V, . Considering
the special
case, by substituting
the relations
(19) (10) and (17), one can obtain:
F=V,-L,=V, so that F=V-;.
(20)
M.H.N.
112
Furthermore,
J.B. EDWARDS
TABRIZI.
if the feed mixture is in equilibrium so that .z = oz,
(21)
and the composition such that the point Z, 2 lies on the -45O line of the vapour/liquid diagram for symmetry, i.e., r=l-2, (22) then it follows that ff
and
%=l$-cr Z=-. Substituting yields:
(23)
1
(24)
l+o
the relations (20) and (23) into (18) an d e 1iminating X(0) in favour of Y(0) therefore o HI 8Y(O) 6/$ = ---& K at
+ X’(0) -t (1 - Y(O)} + F
6/Z’.
(25)
Expressing normalized distance as h=-$
and normalized time as
then the Equation
(26)
tv, r=cuH16h’1
(27)
(25) reduces to the simplified form: 8Y(O) = &
--&
+X’(0) + (1 - Y(O)} + F.
Now, taking the tray immediately below the feed point and taking a mass balance, we obtain H’I Sh’ q
= F 2 +
and, again under symmetrical
H; ax’(o) v, dt
v, {Y’(-1)
- Y’(O)} + L,X(O)
- -&X’(O),
(29)
operation, Equation (29) reduces to: &,
-- 2
=
o+l
If, therefore, the tray capacitances
v
6h’ L (1 - Y(O)} - X’(0).
(30)
are such that H&Y&
(31)
then, in terms of normalized time and distance, Equation (30) can be reduced to a simpler form as
ax’(o)
-
Equations 2.4.
67
= -
2
(Y-i-1
- qp
- (1 -
Y(O)}-X’(O).
(32)
(28) and (32) clearly share a marked degree of symmetry.
Terminal
Boundary
Conditions
The variables associated with the accumulator and reboiler ends of the process are shown in Figure 5, the top tray being the Nth and N + 1 denoting accumulator quantities. If H, is the molar capacitance of the accumulator, then on taking a mass balance on this vessel, we obtain:
H
a
ax(;t+ ‘)=v,
{Y(N) - X(N
+ 1)).
113
Parametric model for distillation columns
Accumulator
(top
product) +
“r-L,
X(N)
Accumulator-end variables.
X' (-d,+ 6h’ )
(bottom product)
Reboiler
Heater
I
I
\
‘b
Reboiler-end variables.
Figure 5.
Regarding Y ss a spatially continuous variable, Y(N) using a truncated Taylor expansion to give: H 0 ax(;t+
l) = v,
{
Y(N+
1) -
can be eliminated in terms of Y(N -t 1)
OY(N + 1) 6h’ - X(N + 1) ._#
, >
so that, in terms of Y (IV + 1) only, we have H oBY(N+ 0 at
1)
Y(N+l)(l-o)+o--l-
(33)
And if Y(N + 1) is now replaced by Y(Ll), where L1 is the normalized length of the rectifier, then in terms of r and h we finally obtain the general result: T a-dY(Ll) A7
= E{l _ Y(&))
- 8y(Lr) --E---’
(34)
M.H.N.
114
TABRIZI,
J.B. EDWARDS
where T, is the normalized time-constant of the accumulator and is given by HCZ HI 6h’
Ta=-.
(35)
A similar treatment of the reboiler mass balance equation gives H*
dX’(-J52)
L, {X/(-L2
=
dt
+ 6h’) - Y’(-L2)},
where Hb is the molar capacitance of the reboiler and Ls, the normalized length of the strippingsection, produces the general results: T* 8X’(-L2)
8X’(-L2) dh
= -cX’(-L*)+
a7
)
(36)
where Ta is the normalized reboiler time constant given by
Ha
Tb=-. 2.5.
(37)
HI 6h’
Steady-State Solutions
General PDEs (9) and (12), under the particular operating conditions (10) and (13), yield the steady-state differential Equations (11) and (14). And, under the further operating conditions (16), (17), and (23), these must be solved subject to special feed boundary conditions (28) and (32), together with the general terminal conditions (34) and (36). Therefore in our symmetrical steady-state situation the system for solution is as follows: if then
Ls = L1 = L,
(36)
d2X’ &=
(39)
(dh>2 =01
X’(0) + 1 - Y(0) + 7
dY(0)
2 = o+l
(40)
dX’(0) X’(0) + 1 - Y(0) + --&-
2 = o+l’
(41)
and
(42)
7=&-Y(L)} dX’(-L)
= cX’(-L). dh As shown in the Appendix this system has the following solutions: X’(-h)
(43)
= 1 -Y(h),
(44)
Y(h) = Y(0) + G h,
(45)
where profile gradient G is constant and given by G=
(46)
(o+l)(2t22+(1+1)’
and initial values X(0) and Y(0) are given by 2(cL+l) (CY+ 1)(26 L + cr + 1)’
X’(0) = 1- Y(0) =
(47)
From (44) to (47) it is also readily shown that the terminal values X’(-L) by: X’(-L)
= 1 - Y(L)
=
Figure 6 illustrates the steady-state composition even tray loading: a desirable design condition.
(a + 1)(2A
+ (Y+ 1) *
profiles graphically.
and Y(L)
are given (48)
Their linearity indicates
Parametric model for distillation columns
115
: ’ (= 1 - Y)
____
T
--
__-------
I f I
2cL GL = (a+1) (zcL+a+l)
;
T t 2(1+ EL) I , (a+l) (2cL +a+l)
1 t GL=(a+l)
I i
$cL+a+l)
I
I I
1
:.
I
L-h
Figure 6. Steady-state composition profile.
3. SMALL 3.1.
Partial
Diferential
PERTURBATION
MODEL
Equations
The small perturbation PDEs for the system may be derived by implicit differentiation of the general large signal PDEs (9) and (12). The lower case symbols denote small changes in the variables. Substituting the steady state solutions for the upper-case symbols, one can obtain: H,iih’$-:G+
which in terms of normalized
time and distance BY -_= ar
Also from the special
symmetrical
$
cl
H;6h’F
r and y, rather I
v,
8% .
(49)
(W2
(50)
case (17), PDE (50) reduces
and output
than 2, becomes
=
va)G
67 if input
(Sh’)‘,
(lo - v)G
a;c’ = (I-
In particular,
&
+ a2x’ dhz’
v,
vectors are defined, ‘=
[
;;;:
to the normalized
form: (51)
thus
1
(52)
and (53)
Then the small signal PDEs (49) and (51) may be grouped
to form the simple matrix equation
as follows: (54)
M.H.N.
116
Laplace
(2 -p)&$)
-Z(O)
‘W denotes where superscript and “.“, the spatial derivative. Terminal
3.2.
Boundary
Differentiating reboiler we obtain
transforms
(33) implicitly
to obtain
the coefficient
treatment
1 ] 21,
to h and 7, “N” with respect
applies
boundary
-cy(~1)
6h’
dh’
of v in (56) is zero, so that (56) reduces
1
aY(Ll) ar
EY(L1)
a-=-
to the reboiler
+
F}
to r only
- p}
equation
for the
.
to the normalized
(56) form:
.
and yields the similar &‘(-L2) dh
= --m’(-L~)+
aT
Inverted
a;
the small-perturbation
BY(h)
T* al’(---La)
3.3.
[ -@0-l)
to r yields
Conditions
T
Similar
= s-1
with respect
Y(LI)c+c---
In steady-state
to h and in p with respect
(54) in s with respect
transforming
J.B. EDWARDS
TABRIZI,
(57)
result:
.
(58)
U-Tube Model
In the symmetrical case L1 = L2 = L, solutions of the small signal equations are considerably simplified by redefining h, h’ from the origin at the ends of the column rather than the central feed point, by first bending the column conceptually into an inverted U-tube as shown in Figure 7, so that co-ordinate h is replaced by L - h in the rectifier and h by L + h in the stripper. Oddpowered spatial derivatives in the rectifier are consequently reversed in sign, while even-powered derivatives in the rectifier and all derivatives in the stripping section are unaffected in sign. The terminal boundary conditions (57) and (58) thus become T
-ado) = -
a
T* a%‘(O) =
ar
+ y
(59)
and
ax'(o) --Ed(O)+ 7.
(60)
If T, = Tb = T,
(61)
then one can obtain TpC(O) = -q(o) -
+2(o).
(62)
The PDE (54) is unaffected by a change of distance base, since this involves only even-powered spatial derivatives so that its transformed version (55) is also unaffected provided z(O) and l(O) are now interpreted as terminal conditions rather than feed-point conditions. may therefore be eliminated from (55) using boundary conditions (62) to give
(s2-
p);- (s +
E+
Tp)i(O)
= s-l
[ -(ao-
‘)
(y ;
The unknown
a
1 ] 4,
after inversion now involving only the one unknown vector n(O), which may be eliminated substitution of the feed conditions not yet involved in this small signal analysis.
by
Parametric
model for distillation columns
117
Feed
h =L A
Stripper
Rectifier
T h
Reboiler
end
kcumulator
end
Figure 7. Inverted U-tube model.
3.4. Feed Boundary Condition By implicit differentiation of the large signal Equations (15) and (29) one can achieve the small signal equations for the feed boundary condition. The upper case symbols are replaced by the special case steady-state solutions. Retaining the original distance base for the moment, we thus obtain for the rectifier:
H,
6h’
-‘do) = v {Y’(O) at
Normalizing
-
(64) and defining aY(o> - 67 = +
By using the steady-state
Y'(O) Y(O)}+ v, yin y(0) rather
r
~(0)
1a’
ax(o) + 3
oh’
o!
6h’
‘do)
dh”
(64)
than r(O), one can achieve
Gv) - Y(O)}+ 1a g solutions
+
+ fgp
+ c’(0) - y(0).
(65)
P
of Section 2.5, Y’(0) - Y (0) can be eliminated,
thus giving
Y’(0) - Y(0) = &X’(O) + 1 - Y(0) - 1, so
Y’(0) - Y(0)
(a + 1)v 2
= -G
+crl
F +T+z’(O)-y(0).
>
(66)
118
M.H.N.TABRIZI, J.B.EDWARDS
An identical result as
treatment
applied
to the first tray of the stripping
ax’(o)_ 5 Transforming
(a + 111 -
_-av +
v, {
ar
2
(66) and (67) to the new distance aY(L) -=67
G
(a + l)v
1 -
v,
1
produced
the similar
+ Y(0)-z'(O).
base for the inverted
+ol
2
ax’(o)
r
section
(67)
U-tube
therefore
yields:
(68)
-Y.z’(L)-y(L), >
and ar’(L)=G a7 Adding yields:
and subtracting
&+2 0
A+ or in terms of Laplace
ax’(L) _crv + (o + I)! - T + Y(L)- x'(L). 2 v, { > the (68) and (69) and grouping the results in matrix
form therefore
0 (70)
$+7!i
transforms:
P+& 1@)=[Oir-0&+1)1k. 0
[ p+$+20 3.5.
(69)
(71)
Solutions
From the transformed
Inverting g(h)
Equation
from the s, p domain = { &cosh
&h
(63), one can now achieve
to the h, p domain
gives
C(O) ~rii
+ (6 + Tp) sinh ,/$h}
-(ao-
>[
‘)
’ a+1
1 f.
(72)
In order to substitute obtained by differentiating as’(L) = --$ {psinh &L ah
the boundary condition (71), one must obtain w (72) with respect to h and setting h = L, thus giving + ~;i;(c + up) cash ,@L} g(O) + sinh$L
‘)
O a+1
may be
] h, (73)
where G(O) @>
The elements of (74) are now calculated For ql(O) we have (,+2+ and, on substituting
[ -(aoe
which
for &(L)
and w,
(P + 2) {,/S-h {psinh &L
=
1
(74)
. [ Gz(O) individually from (71), (72) and (73).
-3
K(L)
= O.%iG,
we get &L
+ (6 + T’P) sinh JFLI
+ fi(c
{ 1 - coshJi?L}
+ TP) co& fiL} c Gil -
(5)
G(0) G(O)
sinh,/j?Lr
Gi = 0.5c zi,
Parametric model for distillation columns
119
[m
Im
0.4 lPo.91
w=o.
6
_.0.02
0.2
w=o.72
a
I I 0.4
I 0.2
54
b Re
! E - 0.1 L - 10 T = 1.0
I
I\!
’L,---0
2
I
4
6
8
Figure 8. Inverse Nyquist loci of gr’(O,
yielding
10
jw) for the tray column.
the final solution:
w=
6 {(P +
Ql
((1-t
‘J)(coshJiiL- 1)/p+ (sinhJiTL)/Jii+ conditions,
( > ,+;
so that on substituting
G(L)=
for &i(L) and w
+(a + l)(cosh&L
(75)
we obtain -0.5 (3a + 1) &,
we get
+ ,/Xc + TP) sinhJjjL
{p coshfiL
0.5)
{(p + 2)(c + Tp) + p} sinh@L/&’
TIP i- 2 + 6) coshfil+
For &(O), from the feed boundary
yielding
b Re
I
I
I
+ fisinh&L
+ (c + Z’p)coshJiiL}
- 1) ii2 + (CY+ l)(sinh&L)
$
@2(O)
= -0.5 (3a + 1) i?z,
the final solution:
m=_ z2
(o + l)(coshfiL
- 1) + (o + l)sinh@/,/F
{p( I+ T) + e} coshfiL
+ 0.5(3a + 1) + fi( 1+ c + Q)sinh&L *
(76)
Therefore, setting p = ju, Equations (75) and (76) may be used to compute the two Nyquist or Inverse Nyquist loci for this diagonal system. Despite the complexity of the hyperbolic functions involved, the general form of the loci may be deduced with some precision analytically by comiderine
the high md law frenuencv
behaviour
of these functions.
M.H.N.
120
0.4
0.5
TABRIZI, J.B. EDWARDS
0.3
I
0.2
II
0.1
1 I
I
I I
I
-
Re
I 0.18 --
0.28
0.1
--
0.38--
0.2
0.48 _.
0.58--
0.3
0.68--
0.78-m
0.88
--
0.99-. Figure 9. Inverse Nyquist loci of g,‘(O,
4.
The system
GENERAL
FORM
may be described
$&P) where transfer
function
matrix
=
NYQUIST
(77)
takes the diagonal
glc;pp) g2(; 1 p)
[
LOCI
form
= G(h,P)Z(h>P),
G(h,p)
(TFM)
G(hP)
function
0.5
jw) for the tray c&mm.
OF INVERSE
in transfer
0.4
form
’
’
or alternatively G(h,P) where the inverse
TFM G*(h,p)
= G’(h,P)pbP)l
is also diagonal, G’(h,p)
=
taking
(78) the form:
Cgi(;yp)gl(; p)
’
1
1
(79)
where d(hP)
= L7,‘hP)
and
&(hP)
= g,‘@,P).
(80)
Attention here will be confined to the system outputs and therefore to G*(O, p), the elements of which may be calculated directly from (75) and (76) by substituting the known plant parameters CY,L, E and T. Therefore setting p = jw, Equations (75) and (76) may b e used to compute the two Nyquist or Inverse Nyquist loci for this diagonal system. Typical locus shapes are given in Figures 8 and 9 for L = 10. E = 0.1. of= 1 + c) = 1.1 and T = 1.0.
Parametric model for distillation columns
121
4.0
of
model
3.5
3.0 Invalid
2.C
0.4 1.0.
0.2
I
c
ot1
I
I I
0.3
012
0
I
I I
.4
o!s
+E
Figure 10. Parameter space showing static perfo rmance loci and region of model validity.
4.1.
Zero-Frequency
By considering deduced that
Behaviour
the limits
as p +
0 of the right-hand
sides of (75) and
2CL + cx + 1
!Z(O,O) =
and
EL2 + CL + 0.5c
g2+(oYO) = -(o
+ l)L +6.5(3o
+ 1) *
(76), it is quickly
031) (82)
(These relations may be otained by alternatively solving the small signal PDEs (54), subject to the boundary conditions, with operator $$ set to zero.) L = actual length of rectifier (or stripper)/bh = number of trays per section which, of course, must be >> 1 for the spatial approximation to hold good. An important feature to note is that the static gain gr(0,0) for the long tray column can only be positive unlike packed columns [13,14], where similar analysis reveals the possibilitv of a negative value, depending on plant and operating parameters.
M.H.N.
122
TABRIZI, J.B. EDWARDS
P,(’ 3,( 1)
(E =0.2)
(lz=0.15) 3.0
t (c=O.l)
2.0 (E =0.05)
6
4
2
10
8
Figure 11. Comparing values of the static gain gl(O, 0) with those calculated from the simulated step responses.
4.2.
High-Frequency
Behaviour
As w is increased from zero, it becomes sions (75) and (76). C onsidering frequencies
necessary in band
to consider
the functions
of p in expres-
(83) which is a wide band,
recalling
that L > 1.0,
it is clear that hyperbolic
functions
coshfil
(84)
and sinh&iL
may be replaced
by 0.5 exp &L,
Jir = (o.5w)“.5(1 -l-j),
since (85)
and, therefore, has a positive real part allowing the terms in exp(-@L) to be neglected. Because of the upper bound on p set by (83), only the lowest powers of p outside the hyperbolic functions need be considered, so that from (75) we obtain: &(%P)
@lP
w*
2++2c/@
and C2(%P) 82(P)
the exponential terms clearly the lowest powers of p, since
-(a
~
+ I/&
+ I)(1 + l/l/T)
(87)
p(l+T)+c+&+c+Tp)’
cancelling.
However, it is insufficient
merely
to now ignore
c < 1.0, if (84) is to be satisfied.
as insnection
of the static
all but (88)
performance
curves
of Figure
10 reveals.
Parametric model for distillation cohmms
123
Rl(O,““)
0.2 --
0
0
0
0
0
0
A
A
A
A
A
0
0
0
0
0
.x
x
x
x
0
0 v
0 0.1 --
A
A
A
0
O
A A
0 0 0
x
x.X
x
X
Y I
I
I
I
4
2
I
t
6
8
wL
10
Figure 12. Comparing high frequency gains gl(O, oo) (analytical vdues) with those calculated from the simulated step responses.
Terms reduce
having to:
c as the coefficient Cl
[
are therefore
non-dominant
241 + (YIP
(0,P)lG (PI
T2iO,PV~2(P)
and (86) and (87) consequently
1 =
L-l <
[
-(I
+ a&p
(P".51< 1.0,
1 ’
(89)
(91) L-l <
Ip”.gl < 1.0.
(92) Thus the system approaches a purely integrating process above w = L-l. Above w = 1.0 the dynamic response changes, but of course all real processes exhibit highly complex behaviour if the frequency of excitation is raised sufficiently. However, the band L-’ < w < 1.0 is a wide band since the inverse gain moduli are of order: I!?i(O,0011= :
and
&(O, =))I = ;9
(93)
where L-l < w”.5 cg: 1.0, compared
with their approximate
zero-frequency
l~;(O,O)l = ;
and
(94)
values from (81)-(82): 19x0,
ON= &.
In approximating the high-frequency behaviour of the tray-column become apparent2, which is not surprising since this model is confined
(95) no more phenomena have to long columns and would
‘No wave phenomena at the ends of the cohunns ((LBcan occur in packed color [13] for instance.) have been revealed in this andvsis. A consideration of h # 0 might, however, produce such effects.
124
M.H.N.
TABRIZI,
EDWARDS
T
IR
A 10
J.B.
0.95
0.038
a
0.76
0.019
0.57
6
t------+
0'.2
0.4
Re
0.38
4
2
r-
I
I I
I
L_-_,
I
1
1 200
I
400
I
I
1
I 1000
I
a00
600
Figure 13. Inverse Nyquist loci of jT’(O,
b
Rl?
jw).
be expected to drastically attenuate any reflected waves well before their arrival at the top and bottom of the column. Furthermore, we note that the signs of the high and low frequency gains are identical in both cases; this would seem to validate the use of a multivariable first order lag approximation to tray column dynamics. Such a model is readily derived from the system PDEs without the necessity of tedious dynamic solution as is now demonstrated.
4.3.
Effect of Terminal
Capacitance
As it was shown in Section 4.2, T drops out of the high-frequency analysis and it follows that only the final portion of the step-response is affected by changes in T. This is confirmed by the simulation results of Figure 17.
5. VALIDATION
OF TFM
MODEL
The validity of the analytically derived TFM model is clearly confirmed by Figure 11, which from Equation (81) with simulated step compares value of the static gain, gi(O,O), calculated responses. Figure 12 makes a similar comparison of high frequency gain, the analytical values having been obtained from the first equation of set (89). It is clear that agreement is good, even for L down to as little as 3 units of normalized length, i.e., much better than we might have expected from the restriction L >> 1.0 assumed in the derivations. A typical locus obtained from to note the analytical model for gi(0,jw) (Equation (75)) is shown in Figure 8. It is interesting the nearly straight line nature of this locus in passing between the zero to high frequency region. Such loci closely resemble those of first order lags. In particular, assuming a first order lag fit:
Parametric
model for distillation columns
125
dRe
0.016
__
0.032-e
-- (I.01
-- (1.02
t3.03
0.049-.
t 0.066.-
0.04
0.05
Figure 14. Inverse Nyquist loci of gl’(O,jw)
for the tray column.
91(0,0) ih(O,P)
=
l+.$$$$P to the high and low frequency transfer function obtained from Equations (89) and (81), the straight locus is also presented in Figure 13. The similarity of the frequency calibration of Figure 15 compares the open-loop step responses of the sl’(O,G) and Kr(O,j w ) is reassuring. accurately simulated system with the first-order lag model derived from Equations (89) and (81), again confirming the theoretical analysis. First-order lag approximation of Owens [15] may in fact be obtained much more readily directly from the PDEs by setting g = 0 and 00, and form an excellent basis for controller design for tray-columns as has been demonstrated by Tabrizi [8]. It is interesting to note that, for the small terminal capacitance, the frequency calibrations on the inverse Nyquist locus (Figures 8 and 14) change nearly proportionally to frequency. Multiplying w by a high-frequency gain = 10 (from 91) on Figure 8 yields the length of the imaginary component. Similar multiplication by the high frequency gain of 0.5 (again from 91) on Figure 14 does likewise. This confirms that the model with small terminal hold-ups is approximately of the multivariable first-order lag type having the high frequency gains predicted. The theoretical static gains are calculated from (81) and (82) to be:
From time response
of Figures
gi(0) = 2.69,
(96)
g2(0) = -231.5.
(97)
15 and 16 the steady-state
gains are:
gl(0) = 2.7,
(98)
g2(0) = -230.
(99)
126
M.H.N.
TABRIZI, J.B.
Full process simulation First order lag model
0000
E =
EDWARDS
0.1
L = 10 T = 1.0
1
Figure 15. Unit step response of g~(O,p) for the tray column.
From frequency
responses
Figures
8 and 14, the static
gains are:
gr(0) = 2.63,
(100)
gs(0) = -230.
(101)
It is clear that the static gains (98) and (99) compare favourably with (100) and (101) and also there is good agreement with theoretical values (96) and (97). 45’ lines intersect the loci of 13 and 14 at w = 0.035 and 0.007, respectively, yielding predicted time constants of 25 and 143. These values compare well with figures of 27 and 150 obtained from the simulated step responses of Figures 15 and 16, again confirming model validity. 6. CONCLUSIONS
A parametric TFM model has been derived completely analytically for long, symmetricallyoperated, tray-type distillation columns separating binary mixtures. The system has been shown to enjoy a completely diagonal structure3 provided the selected outputs are y(h, 0) - z’(h, 0) and y(h, 0) +z’(h, 0) and selected inputs are v + 1 and v - 1, where y denotes the change in top vapour composition, 2’ that in bottom liquid composition and v and 1 are perturbations in the flow rates of vapour and reflux, respectively. Increase of v + 1 drives y - I’ positively while v - 1 drives y + x’ negatively and with much greater gain. This accords with common experience. The gains of both elements of the TFM are shown to have the same sign at high and low frequency. Unlike the packed columns [13,14,16], non-minimum phase effects are not anticipated in tray-column behaviour, thus permitting confidence in the general application of a multivariable first-order 3The deduction accords well with the findings of previous authors based on part imperical McMoran [l], Shinskey [2], and Stainthorp [3] approximate analytical by Edwards and Jassim [4], or purely simulation methods [5].
Parametric model for distillation col-a
127
-100 __,
g*(O,P)
E - 0.1 L = 10 T = 1.0
--
-x0--
-200--
T Figure 16. Unit step response of 92(0, p) for the tray column.
E = 0.1 L = 10
,
I
I
t
,
b 2000 Figure 17. Unit step response of the long tray column and the effect of terminal capacitance. 200
600
1000
I’ 1400
1800
128
M.H.N.
TABRIZI, J.B.
EDWARDS
lag approximation of Owens [15] to long tray-type columns. The analysis has indicated the absence of significant effects from the reflection of travelling waves in the tray-type column on the composition changes at the ends of the column. However, such effects might nevertheless occur closer to the feed tray, and this possibility should be investigated since control measurements are frequently taken from points well away from the end of distillation columns. Since it has relied on the approximation of discretely-changing spatial functions by approximate continuous functions, the analysis here reported is, of course, restricted strictly to columns of many trays separating mixtures of low relative volatility. Predicted and measured response parameters seem to accord closely however down to normalized column length = 3. The substantially vertical shape calibration and position of the inverse Nyquist loci have been confirmed by low-and (fairly) highfrequency analysis and by simulated step response tests, provided terminal capacitances are small. Large terminal capacitances cause a slow tail on the step responses and the frequency calibrations of the loci are distorted accordingly. One extension of the research reported in this paper was carried out by Tabrizi and Edwards [17] w h ere it involved the reduced order modelling of large terminal capacitance columns. This revealed that a lumped, second-order lead/lag structure could closely approximate the dynamic response predicted (a) by the detailed TFM derived in the present paper and (b) by simulation. The same structure also carries over to packed columns where it also allows the possibility of non-minimum phase behaviour. At very high frequency, the loci depart from the vertical, the step responses being strictly of the Bessel function rather than exponential type. Further work is needed in this area to predict more accurately the extreme initial and final behaviour of long tray columns by simple formulae. The precise and original TFM model derived in this paper will provide a sound test for future simple models (Bessel function ielated or otherwise) that might emerge in future research.
REFERENCES
1. P.D. McMoran, Application of the inverse Nyquist method lJ.h’.A.C. Control Convention, page 122, (1971). 2. F.G. Shinskey, Process 3. F.P. Stainthorp
Control System, McGraw-Hill,
and H.M. Searson,
Tran~. Inst.
4. J.B. Edwards and H.J. Jassim, An analytical Inst. Chem. Eng., 55 (l), 17-29, (1977). 5. H.H. Rosenbrock and C. Storey, Computational (1966).
of 4fh
to a distillation column model, In Proc.
New York, (1967).
Chem.
Eng., 51 (l), 42-50,
(1973).
study of the dynamics of binary distillation Techniques
for Chemical
columns,
Trans.
Engineers, Pergamon Press, London,
6. W.D. Armstrong and R.M. Wood, An introduction to the theoretical evaluation of the frequency response of a distillation column to a change in reflux flow rate, Trans. Inst. Chem. Engrs., 39 (2), 80-90 (1961). 7. 0. R&maker, J.E. Rijindorp and A. Mearleved, Elsevier, Amsterdam, (1975). a. M.H.N. Tab&i, U.K., (1982).
Dynamics
and Conirol
of Coniinuous
Dislillation
Units,
The dynamics and control of binary distillation columns, Ph.D. Thesis, Sheffield University,
9. J.B. Edwards and M.H.N. Tabrizi, Parametric transfer function for binary tray-distillation columns, In Proc. of IFAC Symposium on Dynamics and Control of Chemical Reactors and Distillation Columns, (DYCORD-86), pp.lSS-194, B ournmouth, U.K., (Dec. 1986). 10.
W.L. Wilkinson and W.D. Armstrong, An approximate method fractionating column, Chem. Eng. Sci. 7 (l/2), l-7 (1957).
11.
M.A. Jawson and W. Smith, Counter transfer processes in non-steady A225, page 226, London, (1965).
12. H. Voetter, 13.
Plant
and Process
Dynamics
Characteristics,
of predicting
Butterworth,
composition
state, In PTOC. Royal
response Sot.,
of a
(London),
London, (1957).
M.H.N. Tabrizi, Towards an insight into the dynamic behaviour of tray distillation columns through analysis and simulation of dynamically asymmetric packed columns, Research Report No. 113, Sheffield University, U.K., (1979).
14. J.B. Edwards and M. Guilandoust, Parametric transfer function for packed distillation columns, In Proc. of IFAC Symposium on Dynamica and Conlrol of Chemical Reactors and Diatilllalion Columns, (DYCORD-86), pp. 195-202, Bo urnmouth, U.K., (1986). 15. D.H. Owens, First- and second-order-like IEE, 122, (S), 935-941, (1975).
structures
in linear multivariable
control system design, In Proc.
129
Parametric model for distillation columns
16. M.H.N. Tabrizi, An investigation into the dynamics of the distillation column, search for a wave effects, In Proc. IEEE Southeastcon 90, pp. 599-602 (1990). 17. M.H.N. Tabriti and J.B. Edwards, Comparison of lumped and distributed parameter model performance, pp. 305-308; Modelling for controller design in tray distillation development of a parametric, lumped parameter, transfer function model, pp. 301-304, In Proc. IEEE Soufheasfcon-88, (1988). APPENDIX Of Steady-State
Calculation
Composition Profile.
The large-signal steady-state system is described by the differential equations: d2Y
d2X’
(A.11
(dh)2=(dh)2=Ov which are subject to the boundary conditions
:
X’(0) + 1 -Y(O) x’(o)+1
+ y
-Y(o)+
y dY(L)
-
dh dX’( L) dh From the system symmetry
= -&
(A-2)
= --&,
(A-3)
= e{l-Y(L)},
(A.4)
= cX’(-L).
(A.5)
it is clear that
X’(-h)
= 1 -Y(h),
(A.61
and from (A.l) z
= G =
dh
dx’
= constant.
G’
-=
constant,
dh
(A-7)
From (A.6) it is obvious that G=G’.
(A.6)
Hence,
Y(h)
= Y(0) + G h, = X’(0)
X’(h)
+ G h.
(A-9)
Now from (A.2) or (A.3) and (A.6) we get 2X’(O) + G = 2 a+1
(A.lO)
and from (A.4) or (A.5) and (A.9) G = cX’(-L)
= c {X’(O)
- GL} ,
therefore, X’(0) Hence, eliminating X’(0)
between (A.lO)
G(cL+l)
=
and (A.ll),
2G(cL+l)
we let +G=
2 -I a+1
c G= stated as Equation
2c
(A.12)
a+1(2eLtatl)’
(46) in the main text. Hence from (A.ll)
X’(0) =
(A.ll)
2(eL t 1) a+l(IcLtatl)
and (A.12)
This is stated as (48) in the main text.
= 1 -Y(L)
=
(A.13)
= 1 -Y(O).
This is stated as Equation (47) in the main text. X’( -L) and Y(L) may now be calculated from general Equations (A.9)
X’(-L)
one can get:
to give:
2 a + 1 (2eL + (Y t 1)’
(A.14)