An Interactive Simulation Package For the Dynamic Behaviour of Binary Distillation Columns

An Interactive Simulation Package For the Dynamic Behaviour of Binary Distillation Columns

Copni)(hl © IF .-\C COlllrol o f Di slill alioll COIUIlIIl S and Chemical Re anors. I\ollrn e molllh . l · K I ~ 1 ~ l i AN INTERACTIVE SIMULATION PA...

1MB Sizes 1 Downloads 81 Views

Copni)(hl © IF .-\C COlllrol o f Di slill alioll COIUIlIIl S and Chemical Re anors. I\ollrn e molllh . l · K I ~ 1 ~ l i

AN INTERACTIVE SIMULATION PACKAGE FOR THE DYNAMIC BEHAVIOUR OF BINARY DISTILLATION COLUMNS G. Guardabassi and R. Scattolini Dipartimfllto di Elettroniw, Politeenico di Milano , Piazw Leonardo da Vinei 32 , 20 I 33 M ilano, Italy

Abstract This paper presents a simulation package for the dynamics of plate-type distillation columns. The model underneath basically consists of material balance equations . The relationship between the liquid and vapour compositions in the flows leaving any single plate is described in terms of the Murphree vapour efficiency and relative volatility. Dependence of the latter on operating pressure and composition is taken into ac c ount. A coordinate "closed loop" control of both the size of the time step and the value of the averaging parameter used in the implicit integration algorithm improves the performance and t he numerical robustness of the code. The package, designed to run on a Digital VAX 11/750 under VMS operating system, is constituted by a set of FORTRAN programs coordinated by a unique command file. Keywords. Chemical industry; numerical anal ysis; computer-aided design.

1 . INTRODUCTION The transient behav ior of plate-type distillation columns has long been studied, both theoretically and experimentally ([1)-[5) .

the second stripping section, the feed section, the first and (possibly) the second enriching (or rectifying) section, the condenser (Fig.l). The external (exogenous) variables are denoted as follows.

For binary (or almost binary) systems it appears t hat ma ss balance equati ons provide a description o f the plant that is a c curate enough for many contro l design purposes .

F xF

In this paper, we present an interactive simulati o n package by which a variety of "open loop" experiments c an e as i ly be simulated, visualized, and re corded. Th i s is extremely useful not only to hel p tuning the model on the behaviour of a real plant , but also to produce data best suited for further c ontrol- o riented processing; e . g.: c omputer-aided linearization, order reduction, "n i cely nonlinear" modeling, etc . An enlarged vers i on of the package, including an ample class of control system structures and algorithms is under development .

w D xD R B xB PT PB

:= feed flow rate : = concentration of the more volatile nent in the feed stream := reboiler manipulated variable : = distillate flow rate : = concentration of the more volatile nent in the distillate stream : = reflux flow rate := bottom flow rate : = concentration of the more volatile nent in the bottom stream : = top pressure : = bottom pressure

compo-

compo-

compo-

The plates of the column (condenser and reboiler included) are numbered bottom-up from 1 to n . The variables relevant to the i-th plate are denoted as follows .

The paper is organized as follows . The basic model is described in Section 2 . Section 3 presents the numer i cal integration alg o rithm. The "expert mechanism" for adjusting both the step size and the value o f the averag i ng parameter used in the impl i cit integration algorithm is described in Section 4. In Section 5, the structure and the o rganization of the pac kage are presented. Finally, a few e xamples given in Section 6, serve primarily the purpose of illustrating the code, more than the accuracy by which a real system can in fact be des c ribed by the model of Section 2 .

Li Vi xi Yi Qi zi

2 . THE DYNAMIC DISTILLATION MODEL

Hi Ei Pi

I n o rder to fit with the geometr y of several real plants, without augmenting too much the number of parameters needed t o describe the column, this is assumed to consist of no more than 7 uniform sections: the rebo iler, the first and (possibly)

liquid flow rate to plate i-I vapor flow rate to plate i+l concentration of the more volatile component in the liquid hold-up concentration of the more volatile component in the vapor hold-up flow rate of the liquid side-stream entering the plate concentration of the more volatile component in the liquid side-stream entering the plate liquid hold-up vapor Murphree efficiency pressure

Assumptions: 1) Ei and Hi are constant, each section

III

and uniform

within

G. Guardabassi and R. Scattolini

112 2)

in each plate, the vapor hold-up is negligeable (as compared with the liquid hold-up)

for all Finally, given by

i = 1,2, ... ,n-2 the vapour flow rates are supposed to be

Vi(t) = V1 (t) + ~ (Vi - V1 ) V1 (t) = V- + k1 w(t) where Vi is the value of Vi at the nominal equilibrium . This may be obtained, of course, by running an accurate steady-state simulation program . Since no energy balance equation is (explicitly) included in the dynamic model, the assumption underneath the equation above is that the profile of the vapour flow rates is constant, and independent of the operating condition (if the nominal profile of the vapour flow rates is not known, set ~=O ) . Fig. 1 :

The distillation column

Under the above assumptions, total and partial material balances at each plate lead to writing:

The digital simulation of the model given in Section 2 calls for the numerical solution of a set of stiff nonlinear differential equations. In fact, the stiffness ratio, namely the ratio between the maximum and the minimum absolute eigenvalue of the Jacobian matrix, usually ranges from a few tens to a few thousands with values of some hundreds common.

As for the single-plate sections, we set f'

reboiler (i=1) :

,', feed plate

*

condenser (i=n):

Vo

L1 = 0

zl

xl

Qi

F

Vn

Ln+1 = 0

zn

xn , Ln = R

-B

Q1 El = 1 zi

xF Qn

-D

, En

The relationships between liquid and vapor concentrations, across each plate, are described by: Yi - Yi-1 Y1

3. THE NUMERICAL INTEGRATION ALGORITHM

Y1 *

Yi *

standard implicit method of numerical A integration has been used to compute the state increments on each time-step. In this sense, our integration algorithm is a fairly straightforward variation of Rosenbrock's "second programme" [1]. A distinctive feature of our routine, as compared with [1], is the closed-loop control of the main integration parameters. Any discussion of this important issue is however postponed to the next section. Here, we consider a single time-step. First, the essentials of the implicit numerical integration method used in the sequel will be recalled. Then, the specific formulas relevant to the considered class of distillation columns (Section 2) will be given in detail. 3.1 The numerical integration method

a(x,p) x ",(x,p)

Consider a dynamic system described by

1 + (a(x,p)-l) x

x(t) in what follows, the relative volatility is supposed to be given by: a(x,p) = AO(p) (l-x) + A1(P) x where in turn, for each

k=O,l

A little thought reveals that the 4 parameters needed to specify a( . ,.) are uniquely determined by the values of a in 4 points of the (x,p) plane. These 4 values of a are in fact the only primary information needed to describe the vapourliquid equilibrium of the considered mixture. The program, then, automatically computes the above parameters on the basis of this "simple" information. As for the pressures, it is supposed that Pn-1(t) = Pn(t) = PT PT + (PB - PT) (n-i-1)/(n-2)

= f(x(t),u(t»

a(x,p) and a time interval [",+ot) . We assume that the input is constant over the considered time interval (Fig.2): u(t) = u(,)

V t

E

[",+ot)

If, as is commonly the case, the input is given in sampled form, u(,) will generally take on the value of the last sample preceding , . Then, in case of periodic sampling,

where T is the sampling period and k, largest value of k such that kT < , . Then, x(,+ot) = x(,) + ox (1-9) f(x(,),u(,»

ox

is

the

(x)av ot

+ 9 f(x(,)+ox,u(,»

where 9 E [0,1] will henceforth be referred as averaging parameter. The state increment

to ox

11 3

Dynamic Behavio ur of Binar)' Distillatioll Columns is therefore a solution of linear) algebraic equation:

the

implicit

(nonBi (,)

Ei + 1 ~x(xi+1(')'Pi+1('» a Li (,) + Hi/6t

+

This could be solved by a Newton-Raphson routine:

Ei ~x(xi(')'Pi('»

o

6x(0) 6x(k+1) A z(k)

a L i (,) (l-E i ) := ----------

6x(k) + z(k) b(k)

+

+av i (,)

a Li+1 (,) := - ---------

Ei+1 ~x(xi+1(')'Pi+1('»

where A := I - a fx(x(,),u(,»

6t

b(k) := {(l-a) f(x(,),u(,»

Gi (,) := Li+1(') xi+1(') - Li (,) xi(') +

+

+ Vi - 1 (,) Yi-1(') - Vi (,) Yi(') +

+ a f(x(,)+6x(k),u(,»} 6t - 6x(k) .

+ Qi(.) zi(')

As usual, one more simply sets

Now, by forward elimination, one obtains: in the implicit equation for the linear equation:

6x,

thus obtaining

6Yi +

~i(.)

6Yi+1 = r i (.) , i = 1,2, .. . ,n-1 6Yn

=

r n (,)

where : It is easy to check that this amounts to stopping the Newton-Raphson routine just after the first step. 3.2 Computation of the distillation column

~

increment

for

the

When applied to the model of Section 2, the method above yields, for the increments, the following set of linear algebraic equations. {L i +1 (,) xi+1(') - Li (,) xi(') +

Yi(')

where

+

Vi - 1 (,) Yi-1(') - Vi (,) Yi(') +

+

Qi(') zi(') +

+

a [L i + 1 (,) 6xi+1 - Li (,) 6x i +

+

Vi - 1 (,) 6Yi-1 - Vi (,) 6Yi]} 6t Ei ~(xi('),Pi('»

=

1,2, ... ,n

i

a(x,p)

+ (l-E i ) Yi-1(')

while = -------

4. FEEDBACK CONTROL OF THE INTEGRATION PARAMETERS

(1 + (a(x,p)-l) x)2

The time step 6t is bound to an interval the extremes of which are denoted 6tmin and 6tmax, respectively . The minimum 6tmin is set equal to 0.1 Tmin ' where

AO(p) (l-x) + A1 (p) x

and

Tmin

k = 0,1

By solving for 6xi the 3rd equation, substituting back in the first, one gets

where A1 (,) := Cn (,) := 0 ; while, values of i = 1,2, ... ,n , := -

Thus, by backward substitution, 6y is readily computed. Finally , for all i = 1,2, ... ,n ,

In conclusion, the basic computation cycle is as f ollows . Given the state x(,) , the inputs F(,), xF('),PB('),PT(,),w(,),D(,), and the integration parameters 6t,a, one computes first y(,) and V(,) , then A(,),B(,) , C(,),G(,),~(,),r(,),6 y and finally 6x . My this routine pluS llnear lncerpolaclon, the top and bottom concentrations at the sampling i nstants tk = kT , where T is the same sampling period used f or the input functions, are computed and temporarily recorded for possible further processing. The length of the time step at the sampling instants is also recorded and displayed on request.

a(x,p) + x (l-x) (A1(p)-AO(p»

=

Bi (,) - Ai(') ~i-1(')

Ei ~x(xi('),Pi('»

+ (a(x,p)-l) x ~x(x,p)

Gi (,) - Ai (,) ri-1(')

6Yi - (l-E i ) 6Yi-1

a(x,p) x ~(x,p)

Ci (,) ~i(') := - - - - - - - - Bi (,) - Ai (,) ~i-1(')

(a Li (,) + Hi /6t)(1-E i ) Ei ~x(xi(')'Pi('»

min Hi/(Li+V i ) i

and

is the min i mum plate-time-constant. On the other side, there is no real interest in giving 6t values much larger than 0.1 the average settling time of the considered system . Thus, 6tmax is often set equal to

for all other

-

a Vi - 1(,)

Within the above range, 6t is varied according to a number of rules meant to protect the process against inconsistencies and numerical instability. The first rule is obvious.

(;. (;uardabassi and R. ScallOiini

IH

RULE 1/1. If, for some i, xi(-t) + 6xi ;. [0,1) , then, if 6t > 6tmin, set 6t = 6tmin and recompute 6x; else, send a message and wait for instructions (abort, reduce 6tmin). The second rule helps facing a slightly less catastrophic event. To express it, we need some preliminary definition. Let

The so obtained step size updating rule can be further refined by currently controlling the maximum relative increment 6x "(,) . As long as it ranges between 0.01 and 0.1, we leave the updating rule unchanged, whereas we increase the growing rate up to 5 times its "natural" value if 6x"(,) is less than 0.01 , while we decrease it down even to negative values if 6x"(,) > 0,1 .

be the relative increment of the i-th state variable at time, Then,

In keeping with the fact . that the risk of numerical instability is known to decrease as long as 9 gets close to I, we have: RULE 1/2. If 6x"(,) > 0.2 if 6t > 6tmin, set pute 6x else, send instructions (abort, 6tmin, continue).

+ 0.1 9 2 , then 6t = 6t/3 and recoma message and wait for restart with a reduced

As long as the inputs remain constant, all state and output transients of a distillation column tend to be monotone and quasi-exponential. Let's consider,

then,

for

a moment,

a

linear

time-

invariant first order system and see, on this "toy example", how the step size should best be varied in time. Precisely, let

Fig. 2

The function

~l(')

) i-----,.

hence,

In view of the requirement that the integration accuracy be roughly time-invariant, it should be apparent that making 6t(,) proportional to the inverse of the second derivative of x with respect to , is a quite conceivable choice; then,

~~•.~,,~.,~.,~,~-,~.~~--,7.w~'--7,...~~.c,.,~.~-c,.~~~_~,.~ ~~.~ ,.~ ~~.~ , .~ .. ~ l og ( d x" I

Fig. 3

6th)

The function

~2(')

x(O)-u This is an "open loop" time step programming rule. In order to give it a "closed loop" form, let's take the derivative and approximate it with the incremental ratio: d

6th)/T

All the preceding discusssion is summarized by the following RULE 1/3. If at time ,+6t(,) all input (exogenous) variables remain unchanged, as compared with the immediately preceding time interval, the size of the time step is updated by the formula:

6th)/T where hence,

On the basis of this analysis, we tentatively set

and the functions ~1(')'~2(') are shown in Fig.2 and in Fig.3 , respectively. Otherwise, 6t(,+6t(,» = 6tmin .

6t(0) = 6tmin

where the function ~l(') is given the shape of Fig.2 The deviations of ~l(') from the identity are justified, on one side, by the requirement that, in order to keep the total number of steps within reasonable bounds, the growing rate should not be less than 107., even for very small values of 6t(,)/T. On the other side by the constraint:

In connection with Rule 1/3, it is worth noticing that the sampling period of the output (and input) variables is usually much larger than 6tmin. Of course, the rule is in itself most effective in computing responses to input functions which are constant over relatively long time intervals (e.g. step functions). Due to high stiffness ratio, the values of 6t produced by the updating formula of Rule 1/3 may readily become quite large relative to Tmin' As a result, numerical instability may show up in some state variable, with undamped or poorly damped

D\'Ilamic Bcha"iour of Binan' Distillation Columns

115

oscillations of period roughly equal to 2 6t . An effective remedy is provided the next rule.

structure of the package,the memory occupation is 30KB at most.

Rule 114. Make a distinction between say "regular" and "auxiliary" time steps, the purpose of the latter being to discharge, so to say, the fast dynamics excited by each regular time step. Precisely, each regular step is followed by a sequence of N auxiliary time steps, the first of which has length 6tmin , while the subsequent ones grow (slowly) in size, according to Rule 113. At the end of the auxiliary step sequence, the length of the next regular time step is computed again by Rule 113, but setting now 6t(-c) , i n the right hand side of the step size updating formula, just equal to the length of the last regular time step.

A "simulation experiment" proceeds according to the following steps: Step 1 The program COLDIS provides the display of a sket ch of the bynary distillation column together with a list of the main symbols used during the simulation. Step 1 A list is provided of the data-files already recorded and containing data relevant to different columns and mixtures previously considered. If the data concerned with the problem at hand are already stored, the corresponding file is copied into a new temporary file. Otherwise the program COLLEG provides the input in interactive mode of the process and mixture data, namely: - a label of the simulation; - the units considered; - the number of plates of each section of the column and the corresponding efficiencies; - the liquid hold-up of each element; - the vector Vand the scalars V- , ~ and kl; the values of pressure, concentration and relative volatility needed to specify a(. ,.). A temporary file is then created for future manipulations. On request this data file can be assigned a name specified by the user and permanently recorded. Step l The program COLVER allows to check and possibly modify the data-file generated at the previous step. Step ~ A list is provided of the data-files already recorded and containing different equi l ibrium points previously computed. If the data concerned with the equilibrium point of interest are already recorded, the corresponding data file is copied into a temporary file. Otherwise the program LEGSTI provides the input of the nominal values of the following exogenous signals: top (PT) and bottom (PB) pressures, feed (F), stream (w) and distillate (D) flow-rates, feed composition (xF)' The corresponding steadystate composition in each plate is then computed. On re4uest the data and the results are displayed on the screen and permanently stored into a new file specified by the user. Step ~ The program VERSTI allows to check the nominal values of PT' PB' F, w, D, xF' In the case that any modification of these values is made, the corresponding new equilibrium point is computed and recorded on request. Step ~ The user can decide either to come to an end of the "simulation experiment" or to continue. In this case, the two possibilities illustrated in the following steps 6a and 6b are available. Step 6a The program EQUIL allows to compute a number of equilibrium points corresponding to some constant values of PT' PB' F, w, D, xF interactively specified by the user. On request, the results are recorded into a new data-file. Step 6b First the user is requested to specify by means of FORTRAN-type statements the form of the exogenous signals PT' PB' F, w, D, xF' This can easisy be done by making reference to a set of FORTRAN functions recorded into an existing program library and simulating a number of "canonical signals" (step, ramp, sine, square waves, ... ) (program IN). Then the transient of the concentrations on each plate corresponding to the initial state determined at steps 4 and 5, and to the given input functions is computed (program TRA). The final results can be displayed in semigraphical mode and permanently recorded. A further option allows to have the transient of the main simulation variables plotted on a CALCOMP M84 (program PLDIS). Step?. The "simulation experiment" is terminated. A new "simulation experiment" can be performed simply moving back to anyone of the previous steps.

N may conceivabl y range from 0 to 30 . In our code it is set, by default, equal to 15. Finally a word on the averaging parameter 8 As is well known, the smaller is 8 the higher is the risk of numerical instability. As for the accuracy, we may go back to our "toy example" and easily find that the exact value of 8 is then given by

Since our primary interest is in the slow dynamics and, in view of Rule 113, 6t/Tmax £ [0,1], the conclusion can be drawn that the most accurate value of 8 ranges from 0.5 (Crank-Nicholson's trapezoidal rule) to about 0.6 On the other hand, whenever the inputs hold constant for long periods of time, relative to the process dynamics, the risk of numerical instability increases as long as the transient settles and the time step approaches its maximum value. It is then conceivable to make e slowly increasing in time, and asymptotically equal to 1. The initial value 80 of e is specified by the operator, or set equal to 0.6 by default. The time- c onstant of the 8 updating mechanism is 2 Tmax ' Of course, the value of 8 is reset to 80 as soon as there is a change in the value of one of the input variables. 5 . STRUCTURE AND ORGANIZATION OF THE PACKAGE The interactive package SIMCOL is designed to run on a VAX-ll/750 using the VAX/VMS operating system. The graphical data-display is presently designed for VT100-type terminals with semigraphic options. The package is constituted by a set of FORTRAN77 programs coordinated by an unique command file. This structure provides flexibility to the management of the programs and allows easy maintenance and upgrading of the package. As a minor drawback, during a "simulation experiment", some recorded data-files have to be read more than once. However, being these data-files of small size, the input/output operations involved do not substantially affect the efficiency of the package. SIMCOL has been explicitely structured to meet with the following requirements : * provide a friendly interactive environment to the engineer unfamiliar with the use of computers; * guarantee an immediate comprehension of the results during a "simulation experiment"; * allow the readability of the recorded datafiles concerning both process data and simulation data and results; * provide a final software product reasonably portable on different computers equipped with a FORTRAN compiler. The total size of SIMCOL in its current form is about 203KB. However, thanks to the segmented

G. Guardabassi and R. Scattolini

116

Mah, Seider (Eds.), Foundations of ComputerAided Chemical Process Design, American Institute of Chemical Engineers, New York, 1981F. De Lorenzo et aI., "Theoretical and experimental researches on the dynamics of plate distillation columns. Part II: Experimental tests", Proc. X International Automation and Instrumen~n- Convention and Exhibiti~ Milano (Italy), Nov. 1968.

[5]

6. EXAMPLES The simulation package presented in the preceding sections has been tested on several cases and is currently used to help control-oriented modeling of industrial columns. In this section some aspects of the package performance are shown by examples . Any illustration of the man-machine interaction has instead been avoided due to page constraints. It has however to be stressed that improving the friendliness of this interaction has been a major engagement in every stage of the conception, design and implememtation of the package . Though discussing the accuracy of the basic model presented in Section 2 is out of the scope of this paper, in Fig.4 computed versus experimental responses to step changes in the distillate flow rate are shown. The plant is here a laboratory 20 plates column operating on a mixture benzenetoluene [6]. The sampling period used for this computation is T = 1. The CPU time needed to compute the two transients is 23.03 sec. Fig.5 shows the step length as a function of time, produced by the step length updating mechanism (Section 4) . Note that the minimum (initial) value of 6t is 0.0011

[6]

I~

X

In Fig . 6, the computed response is shown of the same column to a sinusoidal variation of the distillate flow rate: t (min) D(t) = {

0.071

t

~

0.071 + 0 . 01 sin(t 2n/30)

t

> 20

20

Fig. 4 Computed (- ) and expe rimenta l '0 ) response o f x to stepchanges in D. D

starting from the same initial equilibrium and using the same sampling period T = 1 . Due to the nonconstant shape of the input functions, the CPU time now becomes 8 min and 48.02 sec. Yet it may be noted that, by increasing the sampling period to T 3, the computation accuracy is not significantly affected, while the CPU time drops to 4 min and 41.69 sec. 7. CONCLUDING REMARKS A number of extensions of the simulation package SIMCOL are currently underway. In particular, an enlarged version of SIMCOL will embed some FORTRAN programs simulating an ample call of control system structures and algorithms. During the simulation, the transients of the main process and control variables will be displayed on a TEK 4109 color graphic terminal. At any time it will also be possible to modify both the control design (structure and parameters of the controller) and the values of set-points and independent exogenous signals.

Fig . 5 Tr ansien t of tl,e s Lepl e ngh t .

REFERENCES [1]

[2] [3]

[4]

H.H . Rosenbrock, "Calculation of the transient behaviour of distillation columns", British Chemical Engineering, Part I, pp . 364-367, Part II, pp.432-435, Part III, pp.491-494, 1958.R.G.E.Franks, Mathematical Modeling in Chemical Engineering, John Wiley & Sons, New York, 1967. W.L.Luyben, Process Modeling Simulation and Control for Chemical Engineers, McGraw-Hill, New York, 1973 . C.D.Holland, Fundamentals and Modeling of Separation Processes, Prentice-Hall, New York, 1975.

s1 .~1~~~~~~~~~~~~~~~~~ -o.Qo ' ~ .oo !ID:" ' 120 . 00

. 80.00 '

2 10. 011'

2' 0.00 ' 2 10 .0 0 '

c (min) Fig . 6 Compute d r e sponse o f x t o a s inu so i dal variation

oPD.