- Email: [email protected]
Distillation calculations using a locally paramaterized continuation method
Computers and Chemtcal Printed in Great Bntam
Engineerrng Vol. 9,
No. 6, pp. 593-599,
1985 0
0098-1354i85 $3.00 + .oO 1985 Pergamon Press Lfd
DISTILLATION CALCULATIONS USING A LOCALLY PARAMATERIZED CONTINUATION METHOD+ GEORGE Computing
and Information
D. BYRNE and
Support
LYNN A. BAIRD
Division, Exxon Research Park, NJ 07932, U.S.A.
and Engineering
Company,
Florham
(Received 21 November 1983; revision received 22 December 1984; received for publication 12 February 1985) Abstract-Distillation models are large scale nonlinear systems of equations which represent heat and material balances and phase equilibria for fractionation systems. These models are often solved by Newton’s method or one of its variants. Unfortunately, the method may not converge unless the starting point or initial guess for the solution is quite close to the solution which is sought. To resolve this difficulty, at least partially, we have turned to a convergence technique called the locally parameterized continuation method. Continuation is a mechanism for tracing a path along which lie the solutions of a family of problems. In the case of distillation, the starting point on the path corresponds to the solution of an easy, but somewhat unrealistic, distillation model, and the stopping point is the solution of the difficult more realistic model, which was of principal concern. The family of problems, and hence the path, is essentially generated by the variation of an artificial or physical continuation parameter. The continuation process then corresponds to the projection of a tangent along the path, a Newton correction to find the next point on the path and so on, until the terminus is reached. Turns are accomodated through the interchange of the continuation parameter and the appropriate member of the state vector. Scope-The steady state simulation of distillation towers is used for their design and operation. Exxon Research and Engineering Company uses several different proprietary simulation models for various distillation systems. A program called CHEMDIST[lS] has been used since 1967 on the most difficult of these systems. More particularly, it is used for such applications as multifeed demethanizers, butadiene extractive distillation towers, NMP-lube oil towers and sour water strippers. Each of these applications involves highly nonideal systems. The commercial and environmental impact of this class of problems is substantial. The difficulty with CHEMDIST is that is had been observed to fail to converge for some of these problems. These failures can result in the loss of both design and production time, as well as the time of engineering staff members. Consequently, we have developed a convergence strategy which we believe will reduce the lost time. The strategy, called continuation, has been described elsewhere[9,12] and popularized in recent advertisements[5]. In principle, it is not new to chemical engineering[l3,14]. However, the use of locally parameterized continuation (to be described later), in the setting of a distillation tower, is evidently new. The main body of this paper includes the equations describing the steady state operation of a distillation tower, background information on locally parameterized continuation, the linear algebraic method used and some numerical results. Conclusions and Significance-We believe that we have demonstrated that locally parameterized continuation can be used with a commercial distillation code. We do not believe that we have solved all convergence problems for such a code. In the course of the project we have learned several things that may or may not be well known to our readers. We do not believe we can sufficiently stress the potentially dominant role of a thermodynamic properties package. For distillation of highly nonideal compositions, the time spent in execution within a thermodynamics package, the function evaluation and Jacobian evaluation can consume more than 80% of the computing time. The number of significant digits given and the smoothness of the approximations can evidently have considerable effect on distillation codes. We believe it might be more effective to design a distillation code with continuation from scratch than it would be to retrofit continuation to an existing code. Generally speaking, calculations on an IBM machine should all be in double precision. The choice of a linear algebra package can indeed influence computer storage requirements and code speed. However, the potentially predominant role of the thermodynamic properties package limits the gains markedly. Finally, it seems clear to us that a synergistic effort is required for a very high quality code. This effort should involve specialists in numerical mathematics, chemical engineering and computer science.
1. THE MATHEMATICAL
MODEL
The mathematical model for the steady state operation of a distillation column can be most easily described
‘Prepared for presentation puter Design of Distillation
at the Symposium Columns, Spring
on ComNational
AIChE Meeting, Houston, TX, March 27.-31, 1983.
by referring to an annotated schematic diagram for a distillation column, Fig. 1. Figure 2 then corresponds to an input-output diagram for a single tray or stage for that column. Thus, for each of the NTRAY trays or stages involving NCOMP components, we can obtain the following system of nonlinear equations developed
earlier
by Niedzwiecki
[6], where similar equations
bered from the bottom 593
and Wolfe[
151 (see also
are given). Trays are numof the column.
594
GEORGE
D. BYRNE and
Condenser
LYNN A. BAIRD
1 OVERALL L ,+, -
MATERIAL
BALANCE:
(L, + v,, + v,m, + Ff + F,Y,
i =:NTRAY-1 Reflux
1 VAPOR TION:
W; -
PHASE
W,Y, -
MOLE
R, = 0
FRACTION
(1.3) SUMMA-
NCOMP
c
K,,,x,,, -
1.0 = 0
(1.4)
,=,
1
1
Sidestream
j=3
j=2
Bottoms Product
Fig. 1. Distillation
N COMPONENT Lj+
lxj+
+ Fti -
x i=,
+
V,KjiXji)
+
Rj,i + FjY_i,i -
wjy,,i=
1 HEAT
MATERIAL
BALANCES:
I.i
(LjXJi
-
tower schematic.
0
y-,Kj_,,iXj-1 ’i
(1.1)
Wf-i,j
.
BALANCE:
Lj+Ihj+I,iXj+I.r NCOMP
NCOMP
-
+
Lj hj, Xj,i +
1 i=,
1 i=l
(1.2)
q H,i K,i Xji
Vj_1Hj-,,iKj-1.iXj-l.i
+ Qjj
+ @)_I + Qj + QR,j - Q!$ -
Q”w.j-1= 0 .
2. BACKGROUND
Lj+lxi+l, i
ViYj,i
For simplicity,
t
_
z =
*
Fi,iL*
9 1
Rj,i
QR,i
I
t
iYj-1.
Fig. 2. Input/output
i
ON THE CONTINUATION METHOD
4
vj-l,
The unknowns for each stage are NCOMP components: x,,, , temperature: T, liquid flow rate: L,, vapor flow rate: V,. The enthalpies H,,, and h,,, and equilibrium coefficients K,,, are computed from the other unknowns in an Exxon Research and Engineering Co. proprietary library called EDL (Exxon Data Library)[4]. Thus we are concerned with a system of NEQN = NTRAY x (NCOMP + 3) nonlinear equations in as many unknowns. Of course, this system could be developed from a standard reference such as [8]. It has been noted elsewhere[6,8,12] that Newton’s method can be applied to such systems. It has also been noted and observed that Newton’s method can and does fail for highly nonideal systems. These failures led to our work with continuation in the setting of a distillation code. We simply note the following about the above. For a total condenser, the heat duty is fixed. So the heat balance is taken to be invariant (its derivatives vanish). Similarly, the derivatives of 2: Kx are taken to be zero. For a partial condenser, either the temperature or the vapor distillates are fixed and the derivatives of the heat balance are taken to be zero. These remarks apply to the top tray, the site of the condenser, only. Consequently, no Newton correction is made for either the temperature or the vapor flow rate in the case of a total condenser. For a partial condenser the temperature is similarly held constant.
Lixj, i
diagram for stage j (steady state).
let
[x,.,,x,,z ,... ,x I.NCOMPI~,,v,~-L...~
T NTRAY~
XNTRAY,I,...‘XNTRAY.NCOMP,
v NTRAY,
LNTRAYlr
9
where the superscript T denotes transpose. The trays or stages (see Fig. 1) are numbered from the bottom to the top. Next let F(Z) denote the system of equations developed for the entire column. That is, the system (l.l)-(1.4) is written forj = 1, thenj = 2,..., then j = NTRAY. It has been noted earlier that it is
Distillation
difficult
sometimes
calculations
using a locally
to find Z = Z* to satisfy
F(z*)=o
(2.11
the basic idea is to imbed this hard problem in a family of problems as follows. Construct a homotopy H&g) so that it is relatively easy to find Z,, = Z that satisfies H(Z,,O) = 0
parameterized
continuation
method
pletely satisfactory, as positive Ag does not accommodate turns (or switchbacks) in the solution curve. Finally, it is worth mentioning that it is generally neither practical nor possible to pick Z, = Z* or Z, = Z*. To test a code such choices might be made, but we did not make them in this work. Fairly standard shortcut calculations are used to obtain the starting profile, Z,.
(2.21 3. LOCALLY PARAMETERIZED
or so that (2.2) is automatically construct H so that H (Z*,l)
= F(Z*)
satisfied.
Further,
(2.3)
= 0
That is, H(Z,l) = F(Z). This can be done in a variety of ways. The actual homotopy used by us and by Wolfe and Niedzwiecki[lS] is proprietary. Wolfe and Niedzwiecki[l5], and later Salgovic et a[.[1 31, used the following ad hoc device for continuation. * Estimate a state vector Z,. * Compute Z, to satisfy (2.2) by, say, Newton’s method. - Set g, = 0 + Ag. * Use (Z,,,Ag) as a starting point to compute Z, to satisfy H(Z,,g,) = 0. * Continue this marching process until g = 1 and H(Z*,l) = 0. This process is depicted in Fig. 3. The difficulties with this procedure are as follows. A fixed step length Ag will sometimes result in convergence failure of the underlying nonlinear system solver. This is attributable to the notion that the starting point (Z,,g,+ ,) = [Z,,(j+ l)Ag] may be too far from the solution point (Z,+ ,,g,+ ,) for convergence to occur. General use of a very small Ag can lead to excessive computations. Furthermore, automatic step size reduction is not com-
2
A
Legend Solution
Curve-
Solution
Point
595
.
CONTINUATION
One technique which avoids some of the problems cited is locally parameterized continuation. The basic idea is to again use the homotopy or one parameter family of problems H(Z,g) as before. We assume that a starting point (ZO,gO1 is prescribed so that H(Z,,g, ) = 0. If this condition is not satisfied, then Newton’s method can be used to find (Z,,g, 1, because g, is prescribed. A continuation step from (Z, ,g,) to (Z,, ,,g,, ,) can be described as follows[l2]. Step I. Compute the tangent vector r to the curve H(Z,g) = 0 at the point (Z,,g,). Step 2. Pick the local continuation parameter 0, = max(lT:(:j + 1,2 ,..., NVARI. Step 3. Compute the step length h,. Step 4. Compute the predicted solution point
(Z;+,>g;+,,
= (Z,,g,) + h,T,.
Step 5.
Apply Newton’s method to find (Z,+i,g,+,) to satisfy H(Z, +,,g,+,> = 0. If Newton’s method converges, then i+i+l, else h,-qh,(O
w 0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.6
0.9
1.0
Fig. 3. A schematic diagram of the ad hoc, marching continuation method, with Ag = 0.1.
In both step 1 and step 5 of the algorithm sketch of Section 2, a linear algebraic system is solved. Here we give that solution procedure, which is due to Keller[ 161.
GWRGE D.BYRNE~~~LYNNA.BAIRD
596
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.6 0.9 1.0
Fig. 4. Schematic diagram of locally parameterized continuation.
In either step, the general problem has the following form: Solve for X dX=B
,
with G?’ an NVAR x NVAR matrix and with X and B both column vectors with NVAR components. In fact the system can be partitioned and is
where el,, is a row vector with NVAR componentsall of which are zero except the LCTN-th, and it is unity. Both x and b are column vectors with NEQN + NVAR1 components. Moreover,
J
=
JH(Z,g) JH tZ,g> ~~ aZ Jg I
.
of .G/.
need not be filled with nonzero elements. The procedure can now be given in quasi-Fortran as follows: Step I. Decompose aHl JZ into block lower and block upper triangular factors. Solve (JH/JZ)a = b for a by a block Step 2. solver. Solve (JHIJZ)P = (JHiJg) for P by a Step 3. block solver. If LCTN = NVAR, Step 4. then x (NVAR) = b (NVAR), else x (NVAR) = [a(LCTN)b (NVAR)] IP(LCTN). For Z = l,NEQN, x(1) = a(l) Step 5. x (NvAR)*PU). Here, a and /3 are column vectors with NEQN components. Moreover, as in Fortran, an argument cf an array denotes the component of the array. The block tridiagonal decomposition and solve routines we use are revised versions of earlier codes by Hindmarsh[ 171. The newer versions use LINPACK subroutines and Basic Linear Algebra Subroutine (BLAS)[3].
5.NUMERICALRESULTS
The element in the Zth row and mth column of the NEQN x NEQN matrix JH/JZ is JH’/JZ”. The superscripts denote the component of a vector. Note that aHlag is a column vector with NEQN components. The overall structure of the matrix & is depicted in Fig. 6. Note that since the balance equations for a given tray are coupled only to those equations for trays immediately above and immediately below, aHl JZ is block tridiagonal. Further, the column vector aHlag
H(Z,g) = 0 Fig. 5. A schematic diagram showing one locally terized continuation step.
Fig. 6. The structure
parame-
The locally parameterized continuation method was incorporated in a developmental version of the CHEMDIST code. This version of CHEMDIST was then tested using three cases supplied by our Engineering Technology Department. These three test cases are simulations of highly nonideal distillation systems which, due to the lack of a good initial guess for the nonideal system, had required using the ad hoc stepping procedure described in Section 1. We expected that the new continuation method would provide some improvement in convergence reliability due to its capability for varying the step size and turning corners on the solution curve. Each test run is discussed separately. Easy sour water stripper The purpose of a sour water stripping column is to separate the sour water feed, containing impurities such as carbon dioxide, ammonia and sulfur compounds, into two product streams: a “clean” water
Distillation calculations using a locally parameterized continuation method bottoms product which can safely be discharged to the environment and a vapor distillate stream containing the impurities. The latter stream can then be sent to a sulfur recovery process or to treatment. Sour water strippers are complex towers, since the impurities in the water are ionizable species with highly nonideal behavior. We label this problem “easy” only because the previously used convergence method of CHEMDIST results in a successful solution. Thus, as a minimum requirement, our development CHEMDIST must be able to solve this problem. The tower consists of eight stages, including a partial condenser at the top and a reboiler at the bottom. Sour water, containing large amounts of ammonia and hydrogen sulfide, is fed to the column at stage 7 (first stage below the condenser). The purpose of the tower is to reduce the amount of ammonia and hydrogen sulfide to acceptable levels. The simulation of the tower using CHEMDISTIPECON, the developmental version of CHEMDIST, was as follows. The continuation process started with the imbedded continuation parameter g equal to 0.1 for this test run and converged after nine steps and a total of 36 Newton steps. Initial increments in g were relatively large, advancing g from 0.1 to 0.8 in two steps. Subsequent g increments were smaller, illustrating the need for a variable step length. Comparison of this solution with the traditional CHEMDIST version showed the solutions to be identical, within machine accuracy. Several subsequent runs in which the continuation input parameters (initial step size and continuation parameter, error tolerances) were varied also converged in approximately the same number of steps. All computations reported here were on an IBM 3081. The thermal properties package, the Exxon Data Library, uses single precision arithmetic as does the basic CHEMDIST code. PECON, on the other hand, is in double precision. We used an interface routine to lengthen and shorten precision of data. Hard sour water stripper This sour water stripping column is labeled “hard” because the previously used method of convergence in CHEMDIST failed. The program marched forward from g equal to 0.1-0.9 with successful intermediate solutions; however, when taking the final step to 1.0, the program could not converge. This tower is much like the previous sour water stripper, except that instead of a reboiler, a steam feed to the bottom state of the tower provides the vapor boilup. Initial runs of this test case resulted in g being increased from 0.1 to 0.999 in 35 steps, but the predicted step size then became smaller than the minimum step size, which was set to approximately the machine accuracy, and calculations were stopped. Since g was very nearly 1.0, we decided to modify the code so that in such a case, the program would print a warning that the case was not converged, but then would increment g to 1.O (leaving all other variables unchanged) and calculate the residual values from the heat and material balance equations. Experienced users can then examine this information to see if the problem has
597
been solved sufficiently for engineering purposes, even though a mathematical convergence has not been reached. Examination of the hard-sour-water-stripper residual values indicate that, in fact, the solution is good enough for engineering use. Thus we have demonstrated an improvement of the locally parameterized continuation method in CHEMDIST. Secondary butyl alcohol recovery tower The secondary butyl alcohol (SBOH) recovery tower problem is another simulation that was successfully solved by the traditional CHEMDIST convergence method. Our main reason for choosing this test case was the different dimensions of the problem. Though still a relatively small problem, this distillation is of four components over 12 stages. The purpose of the tower is to recover the secondary butyl alcohol present in the feed stream, mostly water, which enters the tower near the bottom. A steam feed below stage 1 provides the boilup for the tower, and a total condenser at the top sends a portion of the top product back to the tower as reflux. Interestingly, the product stream of primary importance is a liquid sidestream from the middle of the tower which has a higher concentration of SBOH than either the top product or the bottoms product. This condition exists due to the nonideality of the system (azeotrope formation); if the components acted ideally, the highest concentration of SBOH would be at the top of the tower.+ Our development CHEMDIST solved this problem easily, needing only three continuation steps and a total of 12 Newton iterations to reach convergence. Comparison of this solution to that by standard CHEMDIST showed no differences. Ill-conditioning experiments Examination of the development CHEMDIST test case results show the nonlinear systems to be rather highly ill-conditioned; that is, the elements in the state vector vary widely in magnitude. For example, mole fractions are often lOE-8 and smaller, while temperatures are the order of lOE2 and loading (vapor and liquid flow rates) are lOE3. We suspected that the illconditioned nature of these systems was at least partly responsible for computational difficulties and conducted additional numerical experiments with a small, but similar system. We developed a single stage, two-component (water, propanol) flash model which we could run interactively. The model used the same set of material and heat balance equations as CHEMDIST and similar correlations for component properties. In the model, a feed stream of specified rate, composition and temperature is flashed. Increasing the feed rate (and thus the vapor and liquid product rates) results in an increase in the difference in magnitude of the elements of the state vector, which we estimate with a quantity called the condition number. In this way we were able to study the computational difficulties encountered as ‘The boiling point of SBOH and H,O are so similar that separation would be virtually impossible for the ideal system.
GEORGE D. BYRNE and LYNN A. BAIRD
598
the system became more ill-conditioned, i.e. as the condition number became large. Cases run at low flow rates converged quite easily and had condition factors ranging as high as lOE14. Since increasing the total feed rate simply results in a solution where the vapor and liquid product rates are likewise increased (temperature, product compositions and vapor-liquid ratio remain constant), no computational problems might be expected. But, in fact, increasing the flow rate in the flash model eventually resulted in nonconvergence, due to a failure in the linear algebra routines or because the predicted step sizes became smaller than the minimum size. For the nonconverged cases, the condition number was on the order of lOE16, indicating a highly ill-conditioned system. Though not able to solve the ill-conditioning problems we encountered, the tests with the small system convinced us that the problems did not lie specifically with the CHEMDIST program, but rather with the nature of the system of equations being solved [cl. l)(1.4)] The
and
x Y
h H J F
Q R W F
NCOMP NSTAGE Subscripts f i i R L V ;
test
results
may
be summarized
1. D. F. Davidenko, On a new method of numerical solution of systems of nonlinear equations. Dokl. Akod. Nook USSR 88, 925-947 (1953).
as
follows:
CHEMDIST / PECON
CHEMDIST Easy sour water stripper Hard sour water stripper SBOH recovery SBOH recovery 33: TCHEMDIST did not converge dividual problems. fA 33 tray SBOH tower.
33 188 12 24
16 25t 19 23 with coaxing.
CHEMDIST/PECON
Note added in prooj The paper by J. D. Seader, R. Chavez C., and T. L. Wayburn, “Multiple Solutions to Systems of Interlinked Distillation Columns by Differential Homotopy Continuation,” presented at the Annual AIChE Meeting, San Francisco, Nov. 1984 may be of interest to the reader. Acknowledgements-The authors are indebted to Dr. L. A. Barnstone for his encouragement, enthusiasm and support. J. L. Niedzwiecki and Dr. R. G. Wolfe answered many questions and provided us with a good insight toward distillation columns and distillation calculations. Dr. R. S. Stepleman believed continuation would help, and Dr. M. J. Eisner saw the opportunities. We also acknowledge early collaboration with Professor W. C. Rheinboldt of the University of Pittsburgh. Finally, we acknowledge the very helpful discussions with Dr. R. G. Segers, Communications and Computer Sciences Department, Exxon Corporation. These discussions were invaluable in the development of the bordered, block tridiagonal linear algebra package. NOMENCLATURE L V T
Liquid flow rate Vapor flow rate Temperature
Total continuation
Total Newton iterations
Problem
Feed Stage Component Reaction Liquid phase Vapor phase ith iteration Feed
REFERENCES
its scaling.
numerical
Liquid mole fraction Vapor mole fraction Molal liquid enthalpy Molal vapor enthalpy Equilibrium constant Feed rate Heat duty Moles reacted Sidestream Molal feed rate or function No. of components No. of stages
CHEMDIST 9 9 9 9
steps CHEMDIST /PECON 9 35 3 6
was not tuned for these in-
On steplength al2. C. Den Heijer & W. C. Rheinboldt, gorithms for a class of continuation methods. SIAM J. Numer. Anal. 18, 925-947 (1981). C. B. Moler, J. R. Bunch & G. W. 3. J. J. Dongarra, Stewart, LZNPACK User’s Guide. Society for Industrial and Applied Mathematics, Philadelphia (1979). and Exxon Research and Engineer4. Exxon Corporation ing Co., Exxon Data Library Monuol. Proprietary Information (1978). method. Scienti/ic 5. General Motors, The continuation American 23-25 (June 1982) (Advertisement). 6. R. P. Goldstein & R. B. Stanfield, Flexible method for the solution of distillation design problems using the Newton-Raphson technique. Znd. Engng Chem. Process Des Develop. 9, 78-84 (1970). 7. A. C. Hindmarsh, Solution of Block-Tridiagonal Systems of Linear Algebraic Equations. Report No. UCID-30150, Lawrence Livermore National Laboratory, Livermore, CA (Feb. 1977). 8. C. D. Holland, Fundamentals of Multicomponent Distillotion. McGraw-Hill, New York (1981). 9. H. B. Keller, Global homotopies and Newton methods. In C. deBoor & G. H. Golub (eds.), Recent Advances in Numerical Analysis, pp. 73-94. Academic Press, New York (1978). Iterative Solution of 10. J. M. Ortega & W. C. Rheinboldt, Nonlinear Equations in Several Variables. Academic Press, New York (1970).
Distillation
calculations
using a locally parameterized
11. R. H. Perry & C. H. Chilton (eds.), Chemical Engineers’ Handbook, Fifth Edition, McGraw-Hill, New York (1973). 12. W. C. Rheinboldt, Solution fields of nonlinear equations and continuation methods. SIAM .I Numer. Anal. 17, 221-237 (1980). 13. A. Salgovic, V. Hlavacek, &J. Ilavsky, Global simulation of countercurrent separation processes via one-parameter imbedding techniques. Chem. Engng Sci. 36, 1599-1604 (1981). 14. W. E. Schiesser, Application of Davidenko’s method in chemical engineering. In Proceedings of the Second World Congress of Chemical Engineering, Montreal, Vol. 5, pp. 3 15- 3 19, American Institute of Chemical Engineers, New York (1981).
continuation
method
599
CHEMDIST User’s Munual, Program 2525 (Revised), Report No. EE.53E.74. Exxon Engineering Technology Dept., Exxon Research and Engineering Co., Florham Park, NJ, Proprietary Information (19741. Near Singular 16 H. B. Keller, The Bordering Algorithm Points of Higher Nullity. SIAM .I. Sci. and Slat. Compufing (in press). Solution of Block-Tridiagonal Systems 17 A. C. Hindmarsh, of LinearAlgebraic Equation, Report UCID-30150. Lawrence Livermore National Laboratory, Livermore, CA (February 1977).
15 R. G. Wolfe & J. L. Niedzwiecki,
Engineerrng Vol. 9,
No. 6, pp. 593-599,
1985 0
0098-1354i85 $3.00 + .oO 1985 Pergamon Press Lfd
DISTILLATION CALCULATIONS USING A LOCALLY PARAMATERIZED CONTINUATION METHOD+ GEORGE Computing
and Information
D. BYRNE and
Support
LYNN A. BAIRD
Division, Exxon Research Park, NJ 07932, U.S.A.
and Engineering
Company,
Florham
(Received 21 November 1983; revision received 22 December 1984; received for publication 12 February 1985) Abstract-Distillation models are large scale nonlinear systems of equations which represent heat and material balances and phase equilibria for fractionation systems. These models are often solved by Newton’s method or one of its variants. Unfortunately, the method may not converge unless the starting point or initial guess for the solution is quite close to the solution which is sought. To resolve this difficulty, at least partially, we have turned to a convergence technique called the locally parameterized continuation method. Continuation is a mechanism for tracing a path along which lie the solutions of a family of problems. In the case of distillation, the starting point on the path corresponds to the solution of an easy, but somewhat unrealistic, distillation model, and the stopping point is the solution of the difficult more realistic model, which was of principal concern. The family of problems, and hence the path, is essentially generated by the variation of an artificial or physical continuation parameter. The continuation process then corresponds to the projection of a tangent along the path, a Newton correction to find the next point on the path and so on, until the terminus is reached. Turns are accomodated through the interchange of the continuation parameter and the appropriate member of the state vector. Scope-The steady state simulation of distillation towers is used for their design and operation. Exxon Research and Engineering Company uses several different proprietary simulation models for various distillation systems. A program called CHEMDIST[lS] has been used since 1967 on the most difficult of these systems. More particularly, it is used for such applications as multifeed demethanizers, butadiene extractive distillation towers, NMP-lube oil towers and sour water strippers. Each of these applications involves highly nonideal systems. The commercial and environmental impact of this class of problems is substantial. The difficulty with CHEMDIST is that is had been observed to fail to converge for some of these problems. These failures can result in the loss of both design and production time, as well as the time of engineering staff members. Consequently, we have developed a convergence strategy which we believe will reduce the lost time. The strategy, called continuation, has been described elsewhere[9,12] and popularized in recent advertisements[5]. In principle, it is not new to chemical engineering[l3,14]. However, the use of locally parameterized continuation (to be described later), in the setting of a distillation tower, is evidently new. The main body of this paper includes the equations describing the steady state operation of a distillation tower, background information on locally parameterized continuation, the linear algebraic method used and some numerical results. Conclusions and Significance-We believe that we have demonstrated that locally parameterized continuation can be used with a commercial distillation code. We do not believe that we have solved all convergence problems for such a code. In the course of the project we have learned several things that may or may not be well known to our readers. We do not believe we can sufficiently stress the potentially dominant role of a thermodynamic properties package. For distillation of highly nonideal compositions, the time spent in execution within a thermodynamics package, the function evaluation and Jacobian evaluation can consume more than 80% of the computing time. The number of significant digits given and the smoothness of the approximations can evidently have considerable effect on distillation codes. We believe it might be more effective to design a distillation code with continuation from scratch than it would be to retrofit continuation to an existing code. Generally speaking, calculations on an IBM machine should all be in double precision. The choice of a linear algebra package can indeed influence computer storage requirements and code speed. However, the potentially predominant role of the thermodynamic properties package limits the gains markedly. Finally, it seems clear to us that a synergistic effort is required for a very high quality code. This effort should involve specialists in numerical mathematics, chemical engineering and computer science.
1. THE MATHEMATICAL
MODEL
The mathematical model for the steady state operation of a distillation column can be most easily described
‘Prepared for presentation puter Design of Distillation
at the Symposium Columns, Spring
on ComNational
AIChE Meeting, Houston, TX, March 27.-31, 1983.
by referring to an annotated schematic diagram for a distillation column, Fig. 1. Figure 2 then corresponds to an input-output diagram for a single tray or stage for that column. Thus, for each of the NTRAY trays or stages involving NCOMP components, we can obtain the following system of nonlinear equations developed
earlier
by Niedzwiecki
[6], where similar equations
bered from the bottom 593
and Wolfe[
151 (see also
are given). Trays are numof the column.
594
GEORGE
D. BYRNE and
Condenser
LYNN A. BAIRD
1 OVERALL L ,+, -
MATERIAL
BALANCE:
(L, + v,, + v,m, + Ff + F,Y,
i =:NTRAY-1 Reflux
1 VAPOR TION:
W; -
PHASE
W,Y, -
MOLE
R, = 0
FRACTION
(1.3) SUMMA-
NCOMP
c
K,,,x,,, -
1.0 = 0
(1.4)
,=,
1
1
Sidestream
j=3
j=2
Bottoms Product
Fig. 1. Distillation
N COMPONENT Lj+
lxj+
+ Fti -
x i=,
+
V,KjiXji)
+
Rj,i + FjY_i,i -
wjy,,i=
1 HEAT
MATERIAL
BALANCES:
I.i
(LjXJi
-
tower schematic.
0
y-,Kj_,,iXj-1 ’i
(1.1)
Wf-i,j
.
BALANCE:
Lj+Ihj+I,iXj+I.r NCOMP
NCOMP
-
+
Lj hj, Xj,i +
1 i=,
1 i=l
(1.2)
q H,i K,i Xji
Vj_1Hj-,,iKj-1.iXj-l.i
+ Qjj
+ @)_I + Qj + QR,j - Q!$ -
Q”w.j-1= 0 .
2. BACKGROUND
Lj+lxi+l, i
ViYj,i
For simplicity,
t
_
z =
*
Fi,iL*
9 1
Rj,i
QR,i
I
t
iYj-1.
Fig. 2. Input/output
i
ON THE CONTINUATION METHOD
4
vj-l,
The unknowns for each stage are NCOMP components: x,,, , temperature: T, liquid flow rate: L,, vapor flow rate: V,. The enthalpies H,,, and h,,, and equilibrium coefficients K,,, are computed from the other unknowns in an Exxon Research and Engineering Co. proprietary library called EDL (Exxon Data Library)[4]. Thus we are concerned with a system of NEQN = NTRAY x (NCOMP + 3) nonlinear equations in as many unknowns. Of course, this system could be developed from a standard reference such as [8]. It has been noted elsewhere[6,8,12] that Newton’s method can be applied to such systems. It has also been noted and observed that Newton’s method can and does fail for highly nonideal systems. These failures led to our work with continuation in the setting of a distillation code. We simply note the following about the above. For a total condenser, the heat duty is fixed. So the heat balance is taken to be invariant (its derivatives vanish). Similarly, the derivatives of 2: Kx are taken to be zero. For a partial condenser, either the temperature or the vapor distillates are fixed and the derivatives of the heat balance are taken to be zero. These remarks apply to the top tray, the site of the condenser, only. Consequently, no Newton correction is made for either the temperature or the vapor flow rate in the case of a total condenser. For a partial condenser the temperature is similarly held constant.
Lixj, i
diagram for stage j (steady state).
let
[x,.,,x,,z ,... ,x I.NCOMPI~,,v,~-L...~
T NTRAY~
XNTRAY,I,...‘XNTRAY.NCOMP,
v NTRAY,
LNTRAYlr
9
where the superscript T denotes transpose. The trays or stages (see Fig. 1) are numbered from the bottom to the top. Next let F(Z) denote the system of equations developed for the entire column. That is, the system (l.l)-(1.4) is written forj = 1, thenj = 2,..., then j = NTRAY. It has been noted earlier that it is
Distillation
difficult
sometimes
calculations
using a locally
to find Z = Z* to satisfy
F(z*)=o
(2.11
the basic idea is to imbed this hard problem in a family of problems as follows. Construct a homotopy H&g) so that it is relatively easy to find Z,, = Z that satisfies H(Z,,O) = 0
parameterized
continuation
method
pletely satisfactory, as positive Ag does not accommodate turns (or switchbacks) in the solution curve. Finally, it is worth mentioning that it is generally neither practical nor possible to pick Z, = Z* or Z, = Z*. To test a code such choices might be made, but we did not make them in this work. Fairly standard shortcut calculations are used to obtain the starting profile, Z,.
(2.21 3. LOCALLY PARAMETERIZED
or so that (2.2) is automatically construct H so that H (Z*,l)
= F(Z*)
satisfied.
Further,
(2.3)
= 0
That is, H(Z,l) = F(Z). This can be done in a variety of ways. The actual homotopy used by us and by Wolfe and Niedzwiecki[lS] is proprietary. Wolfe and Niedzwiecki[l5], and later Salgovic et a[.[1 31, used the following ad hoc device for continuation. * Estimate a state vector Z,. * Compute Z, to satisfy (2.2) by, say, Newton’s method. - Set g, = 0 + Ag. * Use (Z,,,Ag) as a starting point to compute Z, to satisfy H(Z,,g,) = 0. * Continue this marching process until g = 1 and H(Z*,l) = 0. This process is depicted in Fig. 3. The difficulties with this procedure are as follows. A fixed step length Ag will sometimes result in convergence failure of the underlying nonlinear system solver. This is attributable to the notion that the starting point (Z,,g,+ ,) = [Z,,(j+ l)Ag] may be too far from the solution point (Z,+ ,,g,+ ,) for convergence to occur. General use of a very small Ag can lead to excessive computations. Furthermore, automatic step size reduction is not com-
2
A
Legend Solution
Curve-
Solution
Point
595
.
CONTINUATION
One technique which avoids some of the problems cited is locally parameterized continuation. The basic idea is to again use the homotopy or one parameter family of problems H(Z,g) as before. We assume that a starting point (ZO,gO1 is prescribed so that H(Z,,g, ) = 0. If this condition is not satisfied, then Newton’s method can be used to find (Z,,g, 1, because g, is prescribed. A continuation step from (Z, ,g,) to (Z,, ,,g,, ,) can be described as follows[l2]. Step I. Compute the tangent vector r to the curve H(Z,g) = 0 at the point (Z,,g,). Step 2. Pick the local continuation parameter 0, = max(lT:(:j + 1,2 ,..., NVARI. Step 3. Compute the step length h,. Step 4. Compute the predicted solution point
(Z;+,>g;+,,
= (Z,,g,) + h,T,.
Step 5.
Apply Newton’s method to find (Z,+i,g,+,) to satisfy H(Z, +,,g,+,> = 0. If Newton’s method converges, then i+i+l, else h,-qh,(O
w 0.0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.6
0.9
1.0
Fig. 3. A schematic diagram of the ad hoc, marching continuation method, with Ag = 0.1.
In both step 1 and step 5 of the algorithm sketch of Section 2, a linear algebraic system is solved. Here we give that solution procedure, which is due to Keller[ 161.
GWRGE D.BYRNE~~~LYNNA.BAIRD
596
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.6 0.9 1.0
Fig. 4. Schematic diagram of locally parameterized continuation.
In either step, the general problem has the following form: Solve for X dX=B
,
with G?’ an NVAR x NVAR matrix and with X and B both column vectors with NVAR components. In fact the system can be partitioned and is
where el,, is a row vector with NVAR componentsall of which are zero except the LCTN-th, and it is unity. Both x and b are column vectors with NEQN + NVAR1 components. Moreover,
J
=
JH(Z,g) JH tZ,g> ~~ aZ Jg I
.
of .G/.
need not be filled with nonzero elements. The procedure can now be given in quasi-Fortran as follows: Step I. Decompose aHl JZ into block lower and block upper triangular factors. Solve (JH/JZ)a = b for a by a block Step 2. solver. Solve (JHIJZ)P = (JHiJg) for P by a Step 3. block solver. If LCTN = NVAR, Step 4. then x (NVAR) = b (NVAR), else x (NVAR) = [a(LCTN)b (NVAR)] IP(LCTN). For Z = l,NEQN, x(1) = a(l) Step 5. x (NvAR)*PU). Here, a and /3 are column vectors with NEQN components. Moreover, as in Fortran, an argument cf an array denotes the component of the array. The block tridiagonal decomposition and solve routines we use are revised versions of earlier codes by Hindmarsh[ 171. The newer versions use LINPACK subroutines and Basic Linear Algebra Subroutine (BLAS)[3].
5.NUMERICALRESULTS
The element in the Zth row and mth column of the NEQN x NEQN matrix JH/JZ is JH’/JZ”. The superscripts denote the component of a vector. Note that aHlag is a column vector with NEQN components. The overall structure of the matrix & is depicted in Fig. 6. Note that since the balance equations for a given tray are coupled only to those equations for trays immediately above and immediately below, aHl JZ is block tridiagonal. Further, the column vector aHlag
H(Z,g) = 0 Fig. 5. A schematic diagram showing one locally terized continuation step.
Fig. 6. The structure
parame-
The locally parameterized continuation method was incorporated in a developmental version of the CHEMDIST code. This version of CHEMDIST was then tested using three cases supplied by our Engineering Technology Department. These three test cases are simulations of highly nonideal distillation systems which, due to the lack of a good initial guess for the nonideal system, had required using the ad hoc stepping procedure described in Section 1. We expected that the new continuation method would provide some improvement in convergence reliability due to its capability for varying the step size and turning corners on the solution curve. Each test run is discussed separately. Easy sour water stripper The purpose of a sour water stripping column is to separate the sour water feed, containing impurities such as carbon dioxide, ammonia and sulfur compounds, into two product streams: a “clean” water
Distillation calculations using a locally parameterized continuation method bottoms product which can safely be discharged to the environment and a vapor distillate stream containing the impurities. The latter stream can then be sent to a sulfur recovery process or to treatment. Sour water strippers are complex towers, since the impurities in the water are ionizable species with highly nonideal behavior. We label this problem “easy” only because the previously used convergence method of CHEMDIST results in a successful solution. Thus, as a minimum requirement, our development CHEMDIST must be able to solve this problem. The tower consists of eight stages, including a partial condenser at the top and a reboiler at the bottom. Sour water, containing large amounts of ammonia and hydrogen sulfide, is fed to the column at stage 7 (first stage below the condenser). The purpose of the tower is to reduce the amount of ammonia and hydrogen sulfide to acceptable levels. The simulation of the tower using CHEMDISTIPECON, the developmental version of CHEMDIST, was as follows. The continuation process started with the imbedded continuation parameter g equal to 0.1 for this test run and converged after nine steps and a total of 36 Newton steps. Initial increments in g were relatively large, advancing g from 0.1 to 0.8 in two steps. Subsequent g increments were smaller, illustrating the need for a variable step length. Comparison of this solution with the traditional CHEMDIST version showed the solutions to be identical, within machine accuracy. Several subsequent runs in which the continuation input parameters (initial step size and continuation parameter, error tolerances) were varied also converged in approximately the same number of steps. All computations reported here were on an IBM 3081. The thermal properties package, the Exxon Data Library, uses single precision arithmetic as does the basic CHEMDIST code. PECON, on the other hand, is in double precision. We used an interface routine to lengthen and shorten precision of data. Hard sour water stripper This sour water stripping column is labeled “hard” because the previously used method of convergence in CHEMDIST failed. The program marched forward from g equal to 0.1-0.9 with successful intermediate solutions; however, when taking the final step to 1.0, the program could not converge. This tower is much like the previous sour water stripper, except that instead of a reboiler, a steam feed to the bottom state of the tower provides the vapor boilup. Initial runs of this test case resulted in g being increased from 0.1 to 0.999 in 35 steps, but the predicted step size then became smaller than the minimum step size, which was set to approximately the machine accuracy, and calculations were stopped. Since g was very nearly 1.0, we decided to modify the code so that in such a case, the program would print a warning that the case was not converged, but then would increment g to 1.O (leaving all other variables unchanged) and calculate the residual values from the heat and material balance equations. Experienced users can then examine this information to see if the problem has
597
been solved sufficiently for engineering purposes, even though a mathematical convergence has not been reached. Examination of the hard-sour-water-stripper residual values indicate that, in fact, the solution is good enough for engineering use. Thus we have demonstrated an improvement of the locally parameterized continuation method in CHEMDIST. Secondary butyl alcohol recovery tower The secondary butyl alcohol (SBOH) recovery tower problem is another simulation that was successfully solved by the traditional CHEMDIST convergence method. Our main reason for choosing this test case was the different dimensions of the problem. Though still a relatively small problem, this distillation is of four components over 12 stages. The purpose of the tower is to recover the secondary butyl alcohol present in the feed stream, mostly water, which enters the tower near the bottom. A steam feed below stage 1 provides the boilup for the tower, and a total condenser at the top sends a portion of the top product back to the tower as reflux. Interestingly, the product stream of primary importance is a liquid sidestream from the middle of the tower which has a higher concentration of SBOH than either the top product or the bottoms product. This condition exists due to the nonideality of the system (azeotrope formation); if the components acted ideally, the highest concentration of SBOH would be at the top of the tower.+ Our development CHEMDIST solved this problem easily, needing only three continuation steps and a total of 12 Newton iterations to reach convergence. Comparison of this solution to that by standard CHEMDIST showed no differences. Ill-conditioning experiments Examination of the development CHEMDIST test case results show the nonlinear systems to be rather highly ill-conditioned; that is, the elements in the state vector vary widely in magnitude. For example, mole fractions are often lOE-8 and smaller, while temperatures are the order of lOE2 and loading (vapor and liquid flow rates) are lOE3. We suspected that the illconditioned nature of these systems was at least partly responsible for computational difficulties and conducted additional numerical experiments with a small, but similar system. We developed a single stage, two-component (water, propanol) flash model which we could run interactively. The model used the same set of material and heat balance equations as CHEMDIST and similar correlations for component properties. In the model, a feed stream of specified rate, composition and temperature is flashed. Increasing the feed rate (and thus the vapor and liquid product rates) results in an increase in the difference in magnitude of the elements of the state vector, which we estimate with a quantity called the condition number. In this way we were able to study the computational difficulties encountered as ‘The boiling point of SBOH and H,O are so similar that separation would be virtually impossible for the ideal system.
GEORGE D. BYRNE and LYNN A. BAIRD
598
the system became more ill-conditioned, i.e. as the condition number became large. Cases run at low flow rates converged quite easily and had condition factors ranging as high as lOE14. Since increasing the total feed rate simply results in a solution where the vapor and liquid product rates are likewise increased (temperature, product compositions and vapor-liquid ratio remain constant), no computational problems might be expected. But, in fact, increasing the flow rate in the flash model eventually resulted in nonconvergence, due to a failure in the linear algebra routines or because the predicted step sizes became smaller than the minimum size. For the nonconverged cases, the condition number was on the order of lOE16, indicating a highly ill-conditioned system. Though not able to solve the ill-conditioning problems we encountered, the tests with the small system convinced us that the problems did not lie specifically with the CHEMDIST program, but rather with the nature of the system of equations being solved [cl. l)(1.4)] The
and
x Y
h H J F
Q R W F
NCOMP NSTAGE Subscripts f i i R L V ;
test
results
may
be summarized
1. D. F. Davidenko, On a new method of numerical solution of systems of nonlinear equations. Dokl. Akod. Nook USSR 88, 925-947 (1953).
as
follows:
CHEMDIST / PECON
CHEMDIST Easy sour water stripper Hard sour water stripper SBOH recovery SBOH recovery 33: TCHEMDIST did not converge dividual problems. fA 33 tray SBOH tower.
33 188 12 24
16 25t 19 23 with coaxing.
CHEMDIST/PECON
Note added in prooj The paper by J. D. Seader, R. Chavez C., and T. L. Wayburn, “Multiple Solutions to Systems of Interlinked Distillation Columns by Differential Homotopy Continuation,” presented at the Annual AIChE Meeting, San Francisco, Nov. 1984 may be of interest to the reader. Acknowledgements-The authors are indebted to Dr. L. A. Barnstone for his encouragement, enthusiasm and support. J. L. Niedzwiecki and Dr. R. G. Wolfe answered many questions and provided us with a good insight toward distillation columns and distillation calculations. Dr. R. S. Stepleman believed continuation would help, and Dr. M. J. Eisner saw the opportunities. We also acknowledge early collaboration with Professor W. C. Rheinboldt of the University of Pittsburgh. Finally, we acknowledge the very helpful discussions with Dr. R. G. Segers, Communications and Computer Sciences Department, Exxon Corporation. These discussions were invaluable in the development of the bordered, block tridiagonal linear algebra package. NOMENCLATURE L V T
Liquid flow rate Vapor flow rate Temperature
Total continuation
Total Newton iterations
Problem
Feed Stage Component Reaction Liquid phase Vapor phase ith iteration Feed
REFERENCES
its scaling.
numerical
Liquid mole fraction Vapor mole fraction Molal liquid enthalpy Molal vapor enthalpy Equilibrium constant Feed rate Heat duty Moles reacted Sidestream Molal feed rate or function No. of components No. of stages
CHEMDIST 9 9 9 9
steps CHEMDIST /PECON 9 35 3 6
was not tuned for these in-
On steplength al2. C. Den Heijer & W. C. Rheinboldt, gorithms for a class of continuation methods. SIAM J. Numer. Anal. 18, 925-947 (1981). C. B. Moler, J. R. Bunch & G. W. 3. J. J. Dongarra, Stewart, LZNPACK User’s Guide. Society for Industrial and Applied Mathematics, Philadelphia (1979). and Exxon Research and Engineer4. Exxon Corporation ing Co., Exxon Data Library Monuol. Proprietary Information (1978). method. Scienti/ic 5. General Motors, The continuation American 23-25 (June 1982) (Advertisement). 6. R. P. Goldstein & R. B. Stanfield, Flexible method for the solution of distillation design problems using the Newton-Raphson technique. Znd. Engng Chem. Process Des Develop. 9, 78-84 (1970). 7. A. C. Hindmarsh, Solution of Block-Tridiagonal Systems of Linear Algebraic Equations. Report No. UCID-30150, Lawrence Livermore National Laboratory, Livermore, CA (Feb. 1977). 8. C. D. Holland, Fundamentals of Multicomponent Distillotion. McGraw-Hill, New York (1981). 9. H. B. Keller, Global homotopies and Newton methods. In C. deBoor & G. H. Golub (eds.), Recent Advances in Numerical Analysis, pp. 73-94. Academic Press, New York (1978). Iterative Solution of 10. J. M. Ortega & W. C. Rheinboldt, Nonlinear Equations in Several Variables. Academic Press, New York (1970).
Distillation
calculations
using a locally parameterized
11. R. H. Perry & C. H. Chilton (eds.), Chemical Engineers’ Handbook, Fifth Edition, McGraw-Hill, New York (1973). 12. W. C. Rheinboldt, Solution fields of nonlinear equations and continuation methods. SIAM .I Numer. Anal. 17, 221-237 (1980). 13. A. Salgovic, V. Hlavacek, &J. Ilavsky, Global simulation of countercurrent separation processes via one-parameter imbedding techniques. Chem. Engng Sci. 36, 1599-1604 (1981). 14. W. E. Schiesser, Application of Davidenko’s method in chemical engineering. In Proceedings of the Second World Congress of Chemical Engineering, Montreal, Vol. 5, pp. 3 15- 3 19, American Institute of Chemical Engineers, New York (1981).
continuation
method
599
CHEMDIST User’s Munual, Program 2525 (Revised), Report No. EE.53E.74. Exxon Engineering Technology Dept., Exxon Research and Engineering Co., Florham Park, NJ, Proprietary Information (19741. Near Singular 16 H. B. Keller, The Bordering Algorithm Points of Higher Nullity. SIAM .I. Sci. and Slat. Compufing (in press). Solution of Block-Tridiagonal Systems 17 A. C. Hindmarsh, of LinearAlgebraic Equation, Report UCID-30150. Lawrence Livermore National Laboratory, Livermore, CA (February 1977).
15 R. G. Wolfe & J. L. Niedzwiecki,