COMPARISON OF “STEPPING” AND GENERAL COMPLEX MATRIX INVERSION TECHNIQUES
IN CALCULATING THE FREQUENCY RESPONSE OF BINARY DISTILLATION COLUMNS The dynamics of a binary distillation column are calculated in the frequency domain by the “stepping” technique and a general complex matrix inversion method. The computation times for each method are compared as a function of column size.
YNAMIC behavior of distillation columns is one of the most studied areas in chemical process control, because of its importance in the chemical and petroleum industries. Mathematical modeling of complex chemical processes usually requires simplifying assumptions, so that the equations in the simulation can be solved without exceeding economically justifiable computer costs. illso, methods of computation have been sought that make more efficient use of computer time. This paper compares methods of determining the dynamics of binary distillation columns in the frequency domain in order to assess their relative efficiencies quantitatively. Distillation column dynamics are often computed in the frequency domain because the large number of equations makes solution for transfer functions intractable in the Laplace domain. The “stepping” technique for determining dynamics in the frequency domain was developed in the classic work by Lamb and Rippin (1960). This method is very efficient because it involves no matrix inversion and no iteration. Other workers (Wood, 1967) have approached the solution via a general complex matrix inversion technique. A special-purpose complex matrix inversion technique that takes advantage of the structure of the system equations has been discussed by Bollinger (1968). The program has, however, not been released for public use. No comparisons of the computation efficiencies of these techniques have appeared in the literature. This paper presents a direct comparison of the stepping technique and general complex matrix inversion technique. Using the same assumptions as Lamb and Rippin, the linearized equations for the nth tray are: dxn/dt = AnlLn+l - An2V
+
Ansxn+l
-
the reboiler equations in the same form as Equation 1, the equations are transformed into the frequency domain. A value of frequency is specified. The next step is to set one of the variables a t the base or a t the feed tray (q, F , V, B , X B ) to the unit impulse, 1 Oi, and the remainder to zero. The compositions and flow rates are calculated tray by tray until the top of the column is reached. The procedure is repeated for the other base or feed tray variables. The final simultaneous equations can be solved for XD and XB as functions of F , xf, V , and R a t a particular frequency, yielding a single point on each of the eight open-loop transfer function Bode plots ( x D / F , X D / Z ~ ,. . . ). Then a new value of frequency is chosen and the procedure is repeated. The program can be altered to calculate the response of any composition or flow in the column, if these transfer functions are desired. The stepping technique is applicable to any staged system where one can solve explicitly for the variables from stage to stage in a stepwise manner.
+
Complex Matrix Inversion
The state variable equation derived from Equation 3 is: d x / d t = Ax(t)
By(t)
(41
The sans serif denotes a column vector and the boldface denotes a general matrix. By transforming Equation 4 into the frequency domain and gathering like terms, we have:
(iwI- A ) x ( i w ) = B y ( i w ) The frequency response of by Equation 6. ~ ( i w= )
The A’s are constants made up of the steady-state holdups, flows, and compositions. The equations for all the trays can be expressed more compactly in matrix notation.
+
x(iw)
(5)
to an input y ( i w ) can be solved
(iw1 - A ) - ’ B y ( i w ) = P ( i w ) y ( i w )
(6 1
P(iw)is the plant transfer matrix which contains the t’ransfer functions relating all the inputs y ( i w ) and all the outputs x(iw). The transfer vector, Ph, between the kth input y k and the out’puts x can be determined by setting y k equal t o the unit impulse. Define a column vector Y k whose element’s y j are: y3 = 0
+ oi
j=l,m
(7) L # k
yj=l+Oi
j=k
(81
Substituting Yk into Equation 6,
(3 1
Stepping Technique
As t’hename implies, this method involves stepwise calculations from the base of the column to the top. Starting with 838
l&EC
FUNDAMENTALS
Pk(iw)
=
(iwI - A ) - ’ B Y k ( i , )
(9 1
The complex matrix inversion subrouting for solving Equation 9 was obtained from Temple University computing center. The program employs the technique of Wilkinson and workers at Stanford University (McKeeman, 1962), which is based on forming associated real matrices (Lanczos, 1964). Although other inversion methods might also have been studied, only the general technique was used, because this is
the type of program that would be found in most computer libraries. 5 .o
Results and Discussion
Four column sizes were investigated: 6, 10,20, and 30 trays. The calculations were done on a Control Data CDC 6400 digital computer and respective C.P.U. times for the two methods were compared. Table I shows the computation times for the matrix inversion, Ti, and stepping technique, Ts, for 65 values of frequency. The number of additions and multiplications for each technique can be estimated from the execution times per operation given in Table I. Figure 1 shows a linear relationship between stepping time and column size. The intercept on the ordinate is the time required for peripheral calculations such as log modulus, etc., and is nearly constant for either procedure. The equation for the curve is: Ts = 0.11 ( X T ) 2 (10)
4.5
4.0
i u-
3.5
5 t-
," 3.0
a a
Y v) 2.5 t-
+
2.0
Figure 2 shows that the computation time for the general complex matrix inversion method is a nonlinear function of column size. The equation is:
I.5
T i = 1.45(NT)2.37 Equation 11 can be used for estimating computation t,imes for matrix inversion methods of the general type, provided of course, the user knows the relat'ive speed of his computer as compared to the CDC 6400. In addition to excessive computer times, the general matrix inversion method has a second disadvant'age of requiring a large amount of computer memory. The 30-tray case was the largest column size that could be run on the CDC 6400 using only core storage. Over 145K (octal) words of st'orage were required. The data also indicate that doubling t'he column size increases the matrix inversion time by a factor of 5 to 6. Lanczos (1964) st,ated that doubling the column size increases the computation time by a factor of 8. It is believed that t,he shortcutt,ing techniques in the program when operating on a sparse A matrix account for this difference. Both techniques gave five significant figure accuracy down to -80 db. on the 60-bit-word CDC 6400. Bollinger (1968) reported that. the dynamics for a 75-tray column mere calculated for 8 frequencies in less than 5 minut,es on an I B X 7094 computer using a streamlined matrix inversion technique. Basing this example on 65 frequencies, as the cases in this paper, and estimat'ing t'he CDC 6400 to be 3 times as fast as the IBRl 7094, gives a rough time estimate of 800 seconds for a 75-tray column. Comparable stepping and general complex matrix inversion times are about 10 and 40,000 seconds, respectively.
Table 1.
Stepping and Matrix Inversion Computation Timesa K O . Of
Stepping Time,
Trays
Set.*
Matrix Inversion Time, Sec.
6
2.64 3.16 4.25 5.37
100.4 336.2 1891.3 5654
10
20 30
' Cases for 65 values of frequency on CDC 6400 computer. The CDC 6400 has a cycle time of 0.1 microsecond and typical execution times: Addition = 0.3 microsecond. Multiplication = 1.0 microsecond. Division = 2.9 microseconds. Execution time only.
*
I
I .o
I
I
5
I
10 15 20 NUMBER O F T R A Y S
I
I
25
30
Figure 1. "Stepping" technique time as a function of column size
IOK
F
E K
t-
a 5
10
5 10 N U M B E R O F TRAYS
I
Figure 2. size
50
Matrix inversion time as a function of column
Conclusions
The stepping method has been quantitatively shown to be much faster than the general matrix inversion method for calculating column dynamics in the frequency domain. This difference is more marked as column size increases, because the time required for the mat,rix inversion method is a higherVOL.
8
NO.
4
NOVEMBER
1969
839
order function of column size, while the stepping method is a linear one.
PA
= transfer function vector relating x and specific
x
disturbance y k
Nomenclature
Yk
= vector of n state variables (compositions, flows) = vector of forcing functions in which one is equal
F = feed flow rate, moles/unit time L = liquid flow rate, moles/unit time N T = number of trays in column R = reflux flow rate, moles/unit time
y
= vector of m forcing functions-cg., feed rate, feed
t
Ti Ts V w
x XB
XD xj
to unit impulse, and remainder are zero composition SUBSCRIPTS
= time = matrix inversion computation time, seconds
stepping method computation time, seconds vapor flow rate, moles/unit time frequency, radians/unit time mole fraction more volatile component in liquid mole fraction more volatile component in bottoms stream = mole fraction more volatile component in distillate stream = mole fraction more volatile component in feed stream
= = = = =
GREEKLETTERS
i, k = indexes referring to elements of y n = general tray n T = top tray literature Cited
Bollinger, R. E., “Predicting Fractionator Dynamics Using a Frequency Domain Solution Technique,” A.1.Ch.E. meeting, St. Louis, hlo., February 1968. Lamb, D. E., Rippin, D. W. T., A.1.Ch.E. meeting, Washington, December 1!60. Lanczos. C.. Auulied Analvsis.” Prentice-Hall., Ennlenood “ , Cliffs,‘ N. J., 19’64. McKeeman, W. hl., C. A . C. M . 6, 553 (1962). Wood, R. X., Trans. Insl. Chem. Engrs. (London) 46, 190 (1967).
= hydraulic constant, time
@
J. P. SHUNTB W. L. LUYBEN
M ATRTc E s
Lehigh Cniversity Bethlehem, Pa. 18015
A, B = coefficient matrices I = identity matrix P = plant transfer function matrix
RECEIVED for review September 4, 1968 ACCEFTEDMay 27, 1969
DISPERSION IN SMOOTH PIPES WITH ADSORBING WALLS Two methods of deriving the dispersion coefficient for a solute flowing in a smooth pipe with adsorbing walls are presented: direct derivation using Aris’ method of moments and interpretation of Aris’ work on diffusion in a concentric flow system. The possible usefulness of the results to represent an adsorbing packed bed is discussed.
RECEKT paper (Dayan and Levenspiel, 1968) considered
A the axial dispersion of
flowing fluid in packed beds of porous adsorbing solids. One special case of the general solution in this paper (Case 4) suggested that the dispersion coefficient in smooth pipes with adsorbing walls is
D =D
+ K (UZa2/D)
(1
where K accounts for the adsorption as well as the radial diffusion and axial convection effects. As a n example of the generality and flexibility of Ark’ method of moments (1956) we show here that this result can be obtained in two other ways: by direct derivation of the dispersion coefficient for this particular system, and by reinterpreting Aris’ findings (1959a) on the distribution of solute by diffusion, convection, arid exchange between the flowing phases in concentric laminar flow. Direct Derivation of Equations
Consider a carrier fluid containing adsorbing solute flowing through a section of infinite pipe as shown in Figure 1. A material balance for the flowing solute gives
840
l&EC
FUNDAMENTALS
+
where ( r ) measures the variation from the mean velocity of flow as given by
C ( r ) = T+(r)
(3)
An over-all material balance for all solute adsorbed or flowing in a section of pipe of length dx gives
aC 2aw a26 =D -- C+(r) ( d C / d r ) - a at ax2 at
_,+.here
-
2
1
C=2
”
rCdr
(4)
( 51
and (6 1
Introducing the coordinate z = x - U&, which moves a t the of the solute, the material balances of mean velocitJ7 LTeff Equations 2 and 4 become, in terms of C = C ( Z ,r, t ) and = W ( 2 ,t ) ,
w
