- Email: [email protected]
An “efficient” continuation method for the solution of difficult equilibrium stage separation process problems
0098-I 354/88
Compuf. them. Engng, Vol. 12, No. 1, pp. 9%103, 1988 Printed in Great Britain. All rights reserved
$3.00 + 0.00
Copyright Q 1988 Pergamon Journals Ltd
AN “EFFICIENT” CONTINUATION METHOD FOR THE SOLUTION OF DIFFICULT EQUILIBRIUM STAGE SEPARATION PROCESS PROBLEMS D. J. VICKERY, J. J. FERRARI and R. TAYL~R~ Department
of Chemical Engineering, Clarkson University,
Potsdam, NY 13676, U.S.A.
(Received 20 May 1986; final revision received 8 June 1987; received for publication
16 June 1987)
Abstract-A continuation method is proposed for solving equilibrium stage separation process problems that are extremely difficult to solve using Newton’s method. The Murphree stage efficiency is used as the continuation parameter. The procedure used to generate initial estimates of the independent variables is completely automatic. Numerical results show the method to be very effective at converging to the solution of the MESH equations. The method can be implemented very quickly in an existing code that solves the equilibrium stage model equations using Newton’s method.
equations exists. How we would deal with other specifications is discussed at the end of the manuscript. An extremely large number of methods have been devised for solving these equations ([l-5] for reviews and/or descriptions). The single most versatile and powerful approach is to solve all of the MESH equations for all stages simultaneously using Newton’s method [6-91. In Newton’s method, we repeatedly solve a linear approximation to equations (l-5)
INTRODUCFION
Multicomponent separation modelled using the MESH librium stage model: The total material balance:
processes
equations
usually
are
of the equi-
MfE(l+ry)I:+(l+rf)Lj - %,I -Lj_$-((Fy+F:)=O. The component
material
(1)
balances:
VI(Xk+L- Xd = -(ww)
M; = (1 + ry) Vjyd + (1 + r~)L,x,, -~+,~~,j+I--j-,xi,,-I-_(f~+Sf;)=o.
The equilibrium
or, if the Murphree
(3)
vapor efficiency is used,
Eg E E,MvK,j~, - y,j + (1 - E,?‘)yi,j+,
The summation
for (X,,,). (X) denotes the vector of independent variables (Lj, vj, T,, yv, xij for each stage) and k indicates the current iteration; thus (X,) is the current (X) vector and (F) is the vector of discrepancy functions (all the MESH equations for all stages). [JI] is the Jacobian matrix with elements that are partial derivatives of (F) with respect to the independent variables (X). In order to obtain convergence, Newton’s method requires reasonable initial estimates of the independent variables. Unfortunately, there is no guarantee that Newton’s method will converge from any particular initial estimate. Techniques that have been used to improve the reliability of Newton’s method include: damping the Newton step, line searching and so on; but none of these techniques has proven completely satisfactory. The methods most recently proposed for assisting convergence of Newton’s method are homotopy continuation methods [l&13]. In homotopy continuation methods, a sequence of related problems is solved beginning with a problem that is easy to solve and ending with the difficult problem whose solution is actually desired. In a homotopy method the transformation of the easy problem into the difficult problem is done by introducing into the model equations an additional parameter that is adjusted from some initial value (usually zero) giving the easy problem to a final value
t2)
relations: E,7 E Kiixii - yii = 0,
= 0.
(3a)
equation: s, 3 ,$, (Yij - x0) = 0
(4)
and the heat balance equations: H,=(l
+r~)~Hjv+(l
+ri)L,H,!-
- Lj_ ,HF_, - (F;H;’
V,+,Hy+,
+ F;HkF) + Qj = 0. (5)
This paper addresses the problem of solving the equations with the “standard” specifications [I]: number and type of all stages; flow rate, composition, temperature, pressure and stage location of all feeds; pressure on all stages, heat load on all stages. For distillation columns the specification of the heat load on the condenser and reboiler is replaced by the specification of the reflux ratio and bottoms or distillate flow rate. This set of specifications ensures that a solution of the MESH MESH
tTo whom all correspondence
(6)
should be addressed. 99
100
D. J. VICKERY et al.
(usually unity) which gives the difficult problem. Each problem in the sequence is solved using some iterative technique-Newton’s method normally being the method of choice-with the initial guess for that problem obtained using the solution to the previous problem. So, at least for small (and sometimes quite large) changes in the parameter, Newton’s method should have no difficulty converging. Homotopy continuation methods are quite closely related to parametric continuation methods, a class of methods in which a parameter that occurs naturally in the model equations is varied. Parametric continuation methods are most often employed to investigate the sensitivity of a model to a particular quantity: the reflux ratio or bottoms flow rate are parameters that have been used in parametric solutions of the MESH equations 1141. In this respect parametric continuation methods differ from homotopy continuation methods which are employed “merely” to solve difficult problems. A parameter occurring naturally in the MESH equations which we thought would be a good candidate for a continuation parameter is the stage efficiency, E w . The maximum separation possible on a given stage is obtained with an efficiency of unity, a true equilibrium stage. On the other hand, for a vanishingly small stage efficiency the stage performs no separation worth mentioning and the streams leaving the stage have essentially the same flow rates, composition and temperature as the combined feeds. This fact can be exploited in a simple yet efficient (pun intended) continuation method for solving difficult equilibrium stage separation process problems. Although the method we describe belongs to the class of parametric continuation methods, we advocate the method for solving difficult problems, not for parametric sensitivity studies. In this sense it is more closely related to homotopy continuation methods. When we submitted this paper, we were not aware of any other similar work. Since then we have learnt that Miiller [15] has used the stage efficiency as a continuation parameter. He used a stage-to-stage method for solving the MESH equations at each parameter value. Sereno [ 161 has also used the efficiency as a continuation parameter for solving liquid-liquid extraction problems. THE “EFFICIENT”
CONTINUATION
METHOD
steps of the efficient continuation method are quickly stated:
ence tolerance for intermediate problems is 1000 times less demanding than that used to converge the final problem in the sequence. Step 4. If EM”= Efi”., stop. Otherwise, continue with Step 5. Step 5. Increase EM” by specified amount. Step 6. Generate initial estimate for next problem in sequence and return to Step 3. Steps 1, 4 and 5 are discussed in more detail later in the paper. However, Steps 2 and 6 deserve further comment before proceeding further. Step 2. Initialization The initialization procedure is completely automatic; nothing, not even the end stage temperatures, need be guessed. The procedure makes use of the method for solving the MESH equations described in the Appendix. The procedure for using this method to initialize the efficient continuation method is as follows: (1) calculate interstage flowrates assuming constant molar flows. Take account of sidestreams and the thermal condition of each feed stream, (2) set E,!” = 0.0 for all stages j, (3) for each species i, compute xii and y, by solving the quad-diagonal system (A2). Note that with E$“”= 0.0, K-values do not have to be estimated. The diagonal element B, must be unity as though the column were equipped with a total reboilersee the Appendix. (4) For each stage j, normalize x0 and y, and compute temperatures Tj from bubble point calculation. Store last computed set of K-values: &=K~(T,, xkj*Ykj, P). (5) Set Ey” = EzLti, # 0. (6) For each species i, compute xij and yij by solving the quad-diagonal system (A2). Repeat (4)-(6) as many times as desired or until the temperature profile has converged. Step 6 Initial estimates of the independent variables for the problems other than the first in the sequence can be obtained from the solution to the previous problem in one of two ways: (a) use the solution at the previous parameter value unchanged or, (b) solve the linear System [J]
The
Step 1. Specify initial and final values of EM” (0.1 and 1.0 say). Specify intermediate values of EM” at which calculations are to be performed. step 2. Generate initial estimates of flow rates, temperatures and mol fractions of each stream. step 3. Solve MESH equations at current value of EM” using Newton’s method. The converg-
(AX) = - AEM”@‘),
(7)
where: (AX) is the (calculated) change in the (X) vector from one problem to the next, AEM” is the prescribed change in EM” from one problem to the next, (F’) is the derivative of (F) with respect to EM”. The only nonzero elements of (F’) are obtained quite straightforwardly by differentiating equation (3a)
Stage separation process problems with respect to E MV.[J] in equation (7) is the Jacobian of(F) and was calculated when we solved the MESH equations at the present value of EM”. Equation (7) has exactly the same form as equation (6) from which we calculate successive iterates of (X) at any given value of EM”. Thus, to estimate (X) for the next problem in the sequence, we need do nothing more than solve a single set of linear equations with two different right hand side vectors. NUMERICAL
STUDIES
Test problems Specifications for seven problems that we used to test the efficient continuation method are listed in Table 1. All of these problems, which are taken or adapted from the literature (see Vickery and Taylor, [ 17,181 for sources and for a more complete commentary on the problems), are hard to converge using Newton’s method. Only problem 1 can be solved by Newton’s method from the initial estimate recommended by Henley and Seader [ 1] (unless we use some tricks to try to get the method to converge). Problem 7 (and a number of others that we have solved [19]) involves an interlinked system of columns. The linking between the columns gives rise to off band elements in the Jacobian matrix, [J], as well as in the quad-diagonal matrix in equation (A2). These off-band elements are handled using Kubicek’s [20] method. The initial value of EM” was 0. I in all problems. We took steps in EMV of 0.3. That is, following the solution of the initial problem at EM’ = 0.1, we try again at EM” = 0.4, 0.7 and finally, 1.0. RESULTS AND DISCUSSION The
efficient continuation method solved all of these problems (as well as others [19]) without any difficulty. The number of Newton iterations required at each of three values of EM” is summarized in Table 2. The first number in the column headed EM” = 0.1 is the number of Newton iterations required after converging the quad-diagonal algorithm for temperatures and compositions. The second number is the number of iterations required after just one pass through the quad-diagonal algorithm. In the other two columns in Table 2 (headed EM” = 0.4 and 1.O) the first number is the number of Newton iterations needed to converge from an initial (X) vector obtained from the solution to the previous problem by way of equation (7) and the second is the number of iterations needed when the initial (X) vector is the solution to the previous problem. The first number in the last column in Table 2 is the sum of the first numbers in the three preceding columns; the second is the total number of Newton iterations required if equation (7) is not used. None of these totals is very high; indeed, they are comparable to the number of iterations that Newton’s method would be expected to take on very difficult problems that it can solve.
101
Table 1. Specifications of test problems (all flows in kmol h-‘, oreasures in barl 1. 5 = 42 c = 4: methanol, acetone, water, ethanol P = 1, total condenser, R = 3, partial reboiler, B = 124 liquid fad 1 at 321.9 K to Stage 6,1; = 0, 0, 50, 0 liquidfeedzat 310.8K toStage21,f;=65, 25, 5, 5 2. s = I5 c = 3: water, acetone, furfural P = 1, partial condenser, R = 2.5, partial reboiler, B = 50 liquid feed I at 334.7 K to Stage 8, X = 45, 50, 5 3. s = 22 c = 3: methanol, ethanol, water P = 1, total condenser, R - 3, partial reboiler, B = 50 liquid feed 1 at 333.3 K toStage12, 1; = 10, 20, 70 4. s = 25 c = 3: methanol, acetone, chloroform P = 1, partial condenser, R = IO,partial reboiler, B - 75 liquid feed I at 330.7 K to Stage 15,1; = 23, 30, 47 5. s = 29 c = 3: n-heptane, toluene, methanol P = 1.013, total condenser, R = 3.64, partial rcboiler, B = 161.8 feed 1 at 333.8 K to Stage 18, f; = 53.7, 46.4, 356.8 6. s = 34 c = P = 1.013, liquid feed liquid feed
3: n-heptane, toluene, methyl ethyl ketone total condenser, R = 1.0, partial &oiler, 1 at 352.5 K to Stage 15,1; = 100, 15, 85 2 at 352.0 K to Stage 27, X = 170,0, 0
B = 14.7
7. Column 1: s = 32 c = 3: benzene, phenol, cyclohexane P = 1.013, total condenser, R = 5.0, partial reboiler, B = 749.9 liquid feed 1 at 298.1 K to Stage 5.1; = 0, 700, 0 liquid feed 2 at 298. I K to Stage 24, J; = 50, 0, 50 Column 2: s = 12 c = 3: benzene, phenol, cyclohcxanc P = 1.013, total condenser, R = 3.0, partial t&oiler, B = 699.9 feed I to Stage 10, bottoms of Column I
Convergence of the initial problem is always very rapid (two to five iterations is quite typical) at values of EM” less than or equal to 0.1. Converging the quad-diagonal (QD) system leads to a small reduction in the number of Newton iterations required to converge the initial problem. However, this method sometimes requires a great many QD iterations. With hindsight, it would appear better to attempt only two or three QD iterations. Since the initialization procedure described above and the Wang and Henke method used by Waybum and Seader [lo] are equivalent at EM” = 1.0, it is pertinent to ask why the initial estimates appear to be so much better at very low values of the efficiency. There are two related reasons for this: (a) at extremely low values of the efficiency the composition and temperature profiles are virtually flat (with small step changes at feed stages, reboilers and total condensers). Sharp gradients are avoided; thus, there is no serious tendency for Newton’s method to make ridiculous corrections to the set of independent variables;
Table 2. Summary of results Iteration count efficiency value
Problem 1 2 3 4 5 6 7
0.1 214 313 313 213 315 3/5 5/6
0.4 3/4 3/4 617 313 515 414 518
1.0 615 415 617 415 516 315 7/Fail
Total iterations ll/ll IO/l2 15117 9110 13/14 IO/l2 17/-
102
D. J. VICKERYet ol.
(b) the nonlinearity of the MESH equations is reduced quite considerably by multiplying the K-values by a small positive number, EMV. It is better to use equation (7) to generate initial estimates of (X) in subsequent problems. The intermediate problem at EMV= 0.7 was omitted in the calculations reported here, thereby demonstrating that the solution at EMV= 0.4 and equation (7) provides a sufficiently good estimate of the solution at EMV= 1 (at least in these problems). The simpler method of setting the initial (X) vector equal to the solution of the previous problem would have been successful in problem 7 if the solution at EMV= 0.7 had been obtained. On rare occasions, we have experienced difficulty in going from the problem at a parameter value of 0.1 to a parameter value of 0.4 in one step. Examples of problems in this category involved strongly nonideal systems in interlinked systems of columns [19]. In these cases, we have obtained convergence quite easily to talking more (smaller) steps in the parameter, EM”. Another case in point is the dehydration of ethanol using benzene as posed by Magrmssen et al. [I 91. The value of EM” at which the calculations were halted was 1.0 in all cases. However, if the stage efficiency is known to be 60% (for example), we could stop the calculations when we reach that point. If the stage efficiency differs in different parts of the column or even from tray to tray and component to component, that too can be taken into account simply by preventing further increases in the efficiency when the prescribed maximum efficiency for that stage/component has been reached. Other
notes
we have described is not capable of solving problems in which a product purity is specified. To deal with problems of this kind, we would have to convert the problem to one with standard specifications (using, for example, a shortcut method), obtain the solution to that problem and then use reflux ratio and/or bottoms flow rate as parameters in a parametric continuation method in order to find the values that met the product purity specifications.
whether by trying to solve equilibrium stage models of real “far-from-equilibrium” operations, we have been making multicomponent separation process problems harder to solve than they really needed to be. Acknowledgements-DJV was the recipient of an NSF fellowship while this work was undertaken. We would like to thank the reviewer who drew our attention to the work of Sereno [15] and Professor D. W. T. Rippin who kindly supplied us with a copy of Miiller’s thesis. NOMENCLATURE
EM”= Murphrec vapor efficiency E = Equilibrium or efficiency equation F = Feed flow rate (F) = Column matrix of mass balance and efficiency equations (F) = Column matrix of MESH equations H = Enthalpy H = Heat balance equation [J] = Jacobian matrix K = K-value L = Liquid flow rate M = Component material balance Q = Heat load on stage r = Ratio of side stream flow to intcrstage flow s = Number of stages S = Summation equation V = Vapor flow rate (X) = Column matrix of mol fractions (X) = Column matrix of independent variables x = Liquid phase mol fractions y = Vapor phase mol fractions .z = Feed mol fractions Subscripts i = Component index j = Stage index k = Iteration counter
The method
CONCLUDING
REMARKS
The “efficient continuation method” is an effective method for solving equilibrium stage separation process problems that are diflicult to solve in other ways. The method is very easily implemented in an existing code that solves the MESH equations using Newton’s method. The initialization procedure is completely automatic; no initial estimates need be supplied by the user of the method. The fact that “inefficient equilibrium stage” separation calculations are much easier to converge than true equilibrium stage calculations makes us wonder
REFERENCES
1. E. J. Henley and J. D. Seader, Equilibrium Sruge Separafion Operations in Chemical Engineering. Wiley, New York (1981). 2. C. J. Ring, Separation Processes. McGraw-Hill, New York (1980). 3. C. D. Holland, Fundamentals of Multico~ponent Distillation. McGraw-Hill, New York (1981). 4. Y. L. Wang and J. C. Wang, In Proceedings of First Conference on Foundations of Computer-Aided Process Design, (W. D. Seider and R. S. H. Mah, Eds). AIChE
(1980). 5. J. D.‘Seader, C/rem. Eng. Ed. 19, 88 (1985). 6. P. A. Whitehouse, PhD Thesis in Chemical Engineering, University of Manchester Institute of Science &d Technoloav (1964). 7. L. Gy Naphtali, Paper presented at AIChE S6th Nat1 Mtg, San Francisco (1965). 8. L. M. Naphtali and D. P. Sandholm, AIChE .Jll7,148 (1971). 9. Goldstein and Stanfield, Ind. Engng Chem. Process Des. Dev. 9, 78 (1970). 10. T. L. Waybum and J. D. Scader, In Proceedings of Second Conference on Foundations of Computer-Aided Process Design. (A. W. Westerberg and H. H. Chien, Eds). CACHE (1983).
Stage separation 11. T. L. Waybum and J. D. Seader, Comput. them. Engng 11, 7 (1987). 12. J. D. Seader, Computer MoaWing of Chemical Processes. AIChE Monograph (1985). 13. G. Byrne and L. Baird, Comput. them. Engng 9, 593 (1985). 14. J. Jehnek, V. Hlavacek and M. Kubicck, Chem. Engng Sci. 28, 1555 (1973). 15. F. R. Mtiller. PhD Thesis, ETH Zurich (1979). 16. A. M. Sereno, Doctoral Thesis, University of Porto, Portugal (1985). 17. D. J. Vickery and R. Taylor, AIChE JI 32, 547 (1986). 18. D. J. Vickery and R. Taylor, Paper presented at New Orleans AIChE Nat1 Mtg (1986). 19. D. J. Vickery, T. L. Waybum and R. Taylor, Proc. IChemE Symp. DistiIlation and Absorption 2, B305 (1987). 20. M. Kubicek, V. Hlavacek and F. Prochaska, Chem. Engng Sci. 31, 277 (1976). 21. T. Magnussen, M. Michelson and A. Fmdenslund, Proc. IChemE Symp. Distillation, London (1979). 22. W. F. Huber, Hydrocarbon Process. August 121 (1977).
process
103
problems
we may write: (F) = [AlBCO](X) -(R)
B,
C,
D,
A,
Bz
Cr
Dr
A,
B,
Cr
J’,+,&j+,xi,j+,
--L,~,~c~,~~,=+~.
(Al)
If the flow rates and K-values are known we may solve the resulting linear equations for the xi . What makes these methods so elegant is that the toe $. crent matrix of each linear system is tridiagonal; linear systems with a tridiagonal coefficient matrix can be solved very easily using the well known Thomas algorithm [l, 21. Murphree efficiencies are not normally included in the Wang and Henke method (or other related methods) because if equation (4) is used to eliminate the y,,, the upper triangular portion of the coefficient matrix is filled out [22]. We could use equation (4) to eliminate the xy from the component mass balances and solve the resulting tridiagonal system for the. vapor phase compositions instead but this is not recommended ([2], p. 472). Rather than combine equations (2) and (4) into a single set of equations, we define a column matrix of discrepancy functions (QT = (MT. I, El.I
. MY,j_IvEn,]-t1M>vEu, x M;j+,,E,,+,...M;,,~,,),
where M$ is the component material balance equation (2) and E,. represents equation (3a). Each pair of equations depends on only four mol fractions: .x~,~_,, y,, xi/ and y,,, , . Thus, if we define a column matrix of mol fractions (X) by: (XY = (Yr I, xf,l~~.Y~.j-l~xi.j-l, x Yij9 +Yi,j+
D,
A,_,’
B,_;
C,-;
A,-,
BN--I
D,_, G-1
A,
BN
where N (=2s) is the order of the matrix. or! the odd numbered rows of [ABCD] the elements are given by: = --&.I
B2j_1=(1+r~)~
A Quad-diagonal Algorithm for Solving Equilibrium Stage Separation Process Problems Incorporating Stage Eficiencies Many of the methods that have been developed for solving the MESH equations (the Wang and Henke method for distillation, the Bumingham and Otto algorithm for absorption and others-see Henley and Seader [l] Chap. 15; or King [2] Chap. 10) make use of the fact that when the Equilibrium relations (3) are used to eliminate the y,, from the component mass balances (2) we have: (Vjx;I+Lj)xij-
(A2)
[ABCD]=
A,_, APPENDIX
= (0).
With the equations and variables ordered in this way, the coefficient matrix [ABCD] has four adjacent diagonals with the additional diagonal in the upper triangular portion of the matrix:
I3 4,j+ I ’ ’ ’ Yi,.* x,3),
C,_,
=(l
+r:)Lj
D~1-r = -J’j+l
(A3)
(these elements being derived from the component materialbalances) whereas on the even numbered rows we have: Arj= -1 Bz, = E?‘“K. c,
= (;‘-
&)
D,=O
(A4)
which belong to the equilibrium relations. For a total condenser set EM” = 0. For a total n&oiler, set EMV= 0 and B,= 1. The right hand side matrix (R) has elements: R, _ , =
4.2,
R2, = 0.
(A5)
The indexj corresponds to the stage number and stages are numbered down from the top. This linear system of equations can be solved for the x,~ and y,, very easily using Gaussian elimination. We have used this method of solving the component material balances and efficiency equations in place of the corresponding steps in the Wang and Henke method for distillation simulations. While not as fast as the conventional tridiagonal matrix algorithm (the quad-diagonal matrix has twice the number of rows as the tridiagonal matrix in the Wang and Henke method), it permits the straightforward solution of problems in which the stage efficiency is not unity. We estimate that the quad-diagonal algorithm will be faster than Huber’s [22] method of dealing with stage efficiencies (with an upper Hessenberg matrix of order s) if the number of stages exceeds nine. An interactive computer program implementing this algorithm for distillation calculations is available from the authors.
Compuf. them. Engng, Vol. 12, No. 1, pp. 9%103, 1988 Printed in Great Britain. All rights reserved
$3.00 + 0.00
Copyright Q 1988 Pergamon Journals Ltd
AN “EFFICIENT” CONTINUATION METHOD FOR THE SOLUTION OF DIFFICULT EQUILIBRIUM STAGE SEPARATION PROCESS PROBLEMS D. J. VICKERY, J. J. FERRARI and R. TAYL~R~ Department
of Chemical Engineering, Clarkson University,
Potsdam, NY 13676, U.S.A.
(Received 20 May 1986; final revision received 8 June 1987; received for publication
16 June 1987)
Abstract-A continuation method is proposed for solving equilibrium stage separation process problems that are extremely difficult to solve using Newton’s method. The Murphree stage efficiency is used as the continuation parameter. The procedure used to generate initial estimates of the independent variables is completely automatic. Numerical results show the method to be very effective at converging to the solution of the MESH equations. The method can be implemented very quickly in an existing code that solves the equilibrium stage model equations using Newton’s method.
equations exists. How we would deal with other specifications is discussed at the end of the manuscript. An extremely large number of methods have been devised for solving these equations ([l-5] for reviews and/or descriptions). The single most versatile and powerful approach is to solve all of the MESH equations for all stages simultaneously using Newton’s method [6-91. In Newton’s method, we repeatedly solve a linear approximation to equations (l-5)
INTRODUCFION
Multicomponent separation modelled using the MESH librium stage model: The total material balance:
processes
equations
usually
are
of the equi-
MfE(l+ry)I:+(l+rf)Lj - %,I -Lj_$-((Fy+F:)=O. The component
material
(1)
balances:
VI(Xk+L- Xd = -(ww)
M; = (1 + ry) Vjyd + (1 + r~)L,x,, -~+,~~,j+I--j-,xi,,-I-_(f~+Sf;)=o.
The equilibrium
or, if the Murphree
(3)
vapor efficiency is used,
Eg E E,MvK,j~, - y,j + (1 - E,?‘)yi,j+,
The summation
for (X,,,). (X) denotes the vector of independent variables (Lj, vj, T,, yv, xij for each stage) and k indicates the current iteration; thus (X,) is the current (X) vector and (F) is the vector of discrepancy functions (all the MESH equations for all stages). [JI] is the Jacobian matrix with elements that are partial derivatives of (F) with respect to the independent variables (X). In order to obtain convergence, Newton’s method requires reasonable initial estimates of the independent variables. Unfortunately, there is no guarantee that Newton’s method will converge from any particular initial estimate. Techniques that have been used to improve the reliability of Newton’s method include: damping the Newton step, line searching and so on; but none of these techniques has proven completely satisfactory. The methods most recently proposed for assisting convergence of Newton’s method are homotopy continuation methods [l&13]. In homotopy continuation methods, a sequence of related problems is solved beginning with a problem that is easy to solve and ending with the difficult problem whose solution is actually desired. In a homotopy method the transformation of the easy problem into the difficult problem is done by introducing into the model equations an additional parameter that is adjusted from some initial value (usually zero) giving the easy problem to a final value
t2)
relations: E,7 E Kiixii - yii = 0,
= 0.
(3a)
equation: s, 3 ,$, (Yij - x0) = 0
(4)
and the heat balance equations: H,=(l
+r~)~Hjv+(l
+ri)L,H,!-
- Lj_ ,HF_, - (F;H;’
V,+,Hy+,
+ F;HkF) + Qj = 0. (5)
This paper addresses the problem of solving the equations with the “standard” specifications [I]: number and type of all stages; flow rate, composition, temperature, pressure and stage location of all feeds; pressure on all stages, heat load on all stages. For distillation columns the specification of the heat load on the condenser and reboiler is replaced by the specification of the reflux ratio and bottoms or distillate flow rate. This set of specifications ensures that a solution of the MESH MESH
tTo whom all correspondence
(6)
should be addressed. 99
100
D. J. VICKERY et al.
(usually unity) which gives the difficult problem. Each problem in the sequence is solved using some iterative technique-Newton’s method normally being the method of choice-with the initial guess for that problem obtained using the solution to the previous problem. So, at least for small (and sometimes quite large) changes in the parameter, Newton’s method should have no difficulty converging. Homotopy continuation methods are quite closely related to parametric continuation methods, a class of methods in which a parameter that occurs naturally in the model equations is varied. Parametric continuation methods are most often employed to investigate the sensitivity of a model to a particular quantity: the reflux ratio or bottoms flow rate are parameters that have been used in parametric solutions of the MESH equations 1141. In this respect parametric continuation methods differ from homotopy continuation methods which are employed “merely” to solve difficult problems. A parameter occurring naturally in the MESH equations which we thought would be a good candidate for a continuation parameter is the stage efficiency, E w . The maximum separation possible on a given stage is obtained with an efficiency of unity, a true equilibrium stage. On the other hand, for a vanishingly small stage efficiency the stage performs no separation worth mentioning and the streams leaving the stage have essentially the same flow rates, composition and temperature as the combined feeds. This fact can be exploited in a simple yet efficient (pun intended) continuation method for solving difficult equilibrium stage separation process problems. Although the method we describe belongs to the class of parametric continuation methods, we advocate the method for solving difficult problems, not for parametric sensitivity studies. In this sense it is more closely related to homotopy continuation methods. When we submitted this paper, we were not aware of any other similar work. Since then we have learnt that Miiller [15] has used the stage efficiency as a continuation parameter. He used a stage-to-stage method for solving the MESH equations at each parameter value. Sereno [ 161 has also used the efficiency as a continuation parameter for solving liquid-liquid extraction problems. THE “EFFICIENT”
CONTINUATION
METHOD
steps of the efficient continuation method are quickly stated:
ence tolerance for intermediate problems is 1000 times less demanding than that used to converge the final problem in the sequence. Step 4. If EM”= Efi”., stop. Otherwise, continue with Step 5. Step 5. Increase EM” by specified amount. Step 6. Generate initial estimate for next problem in sequence and return to Step 3. Steps 1, 4 and 5 are discussed in more detail later in the paper. However, Steps 2 and 6 deserve further comment before proceeding further. Step 2. Initialization The initialization procedure is completely automatic; nothing, not even the end stage temperatures, need be guessed. The procedure makes use of the method for solving the MESH equations described in the Appendix. The procedure for using this method to initialize the efficient continuation method is as follows: (1) calculate interstage flowrates assuming constant molar flows. Take account of sidestreams and the thermal condition of each feed stream, (2) set E,!” = 0.0 for all stages j, (3) for each species i, compute xii and y, by solving the quad-diagonal system (A2). Note that with E$“”= 0.0, K-values do not have to be estimated. The diagonal element B, must be unity as though the column were equipped with a total reboilersee the Appendix. (4) For each stage j, normalize x0 and y, and compute temperatures Tj from bubble point calculation. Store last computed set of K-values: &=K~(T,, xkj*Ykj, P). (5) Set Ey” = EzLti, # 0. (6) For each species i, compute xij and yij by solving the quad-diagonal system (A2). Repeat (4)-(6) as many times as desired or until the temperature profile has converged. Step 6 Initial estimates of the independent variables for the problems other than the first in the sequence can be obtained from the solution to the previous problem in one of two ways: (a) use the solution at the previous parameter value unchanged or, (b) solve the linear System [J]
The
Step 1. Specify initial and final values of EM” (0.1 and 1.0 say). Specify intermediate values of EM” at which calculations are to be performed. step 2. Generate initial estimates of flow rates, temperatures and mol fractions of each stream. step 3. Solve MESH equations at current value of EM” using Newton’s method. The converg-
(AX) = - AEM”@‘),
(7)
where: (AX) is the (calculated) change in the (X) vector from one problem to the next, AEM” is the prescribed change in EM” from one problem to the next, (F’) is the derivative of (F) with respect to EM”. The only nonzero elements of (F’) are obtained quite straightforwardly by differentiating equation (3a)
Stage separation process problems with respect to E MV.[J] in equation (7) is the Jacobian of(F) and was calculated when we solved the MESH equations at the present value of EM”. Equation (7) has exactly the same form as equation (6) from which we calculate successive iterates of (X) at any given value of EM”. Thus, to estimate (X) for the next problem in the sequence, we need do nothing more than solve a single set of linear equations with two different right hand side vectors. NUMERICAL
STUDIES
Test problems Specifications for seven problems that we used to test the efficient continuation method are listed in Table 1. All of these problems, which are taken or adapted from the literature (see Vickery and Taylor, [ 17,181 for sources and for a more complete commentary on the problems), are hard to converge using Newton’s method. Only problem 1 can be solved by Newton’s method from the initial estimate recommended by Henley and Seader [ 1] (unless we use some tricks to try to get the method to converge). Problem 7 (and a number of others that we have solved [19]) involves an interlinked system of columns. The linking between the columns gives rise to off band elements in the Jacobian matrix, [J], as well as in the quad-diagonal matrix in equation (A2). These off-band elements are handled using Kubicek’s [20] method. The initial value of EM” was 0. I in all problems. We took steps in EMV of 0.3. That is, following the solution of the initial problem at EM’ = 0.1, we try again at EM” = 0.4, 0.7 and finally, 1.0. RESULTS AND DISCUSSION The
efficient continuation method solved all of these problems (as well as others [19]) without any difficulty. The number of Newton iterations required at each of three values of EM” is summarized in Table 2. The first number in the column headed EM” = 0.1 is the number of Newton iterations required after converging the quad-diagonal algorithm for temperatures and compositions. The second number is the number of iterations required after just one pass through the quad-diagonal algorithm. In the other two columns in Table 2 (headed EM” = 0.4 and 1.O) the first number is the number of Newton iterations needed to converge from an initial (X) vector obtained from the solution to the previous problem by way of equation (7) and the second is the number of iterations needed when the initial (X) vector is the solution to the previous problem. The first number in the last column in Table 2 is the sum of the first numbers in the three preceding columns; the second is the total number of Newton iterations required if equation (7) is not used. None of these totals is very high; indeed, they are comparable to the number of iterations that Newton’s method would be expected to take on very difficult problems that it can solve.
101
Table 1. Specifications of test problems (all flows in kmol h-‘, oreasures in barl 1. 5 = 42 c = 4: methanol, acetone, water, ethanol P = 1, total condenser, R = 3, partial reboiler, B = 124 liquid fad 1 at 321.9 K to Stage 6,1; = 0, 0, 50, 0 liquidfeedzat 310.8K toStage21,f;=65, 25, 5, 5 2. s = I5 c = 3: water, acetone, furfural P = 1, partial condenser, R = 2.5, partial reboiler, B = 50 liquid feed I at 334.7 K to Stage 8, X = 45, 50, 5 3. s = 22 c = 3: methanol, ethanol, water P = 1, total condenser, R - 3, partial reboiler, B = 50 liquid feed 1 at 333.3 K toStage12, 1; = 10, 20, 70 4. s = 25 c = 3: methanol, acetone, chloroform P = 1, partial condenser, R = IO,partial reboiler, B - 75 liquid feed I at 330.7 K to Stage 15,1; = 23, 30, 47 5. s = 29 c = 3: n-heptane, toluene, methanol P = 1.013, total condenser, R = 3.64, partial rcboiler, B = 161.8 feed 1 at 333.8 K to Stage 18, f; = 53.7, 46.4, 356.8 6. s = 34 c = P = 1.013, liquid feed liquid feed
3: n-heptane, toluene, methyl ethyl ketone total condenser, R = 1.0, partial &oiler, 1 at 352.5 K to Stage 15,1; = 100, 15, 85 2 at 352.0 K to Stage 27, X = 170,0, 0
B = 14.7
7. Column 1: s = 32 c = 3: benzene, phenol, cyclohexane P = 1.013, total condenser, R = 5.0, partial reboiler, B = 749.9 liquid feed 1 at 298.1 K to Stage 5.1; = 0, 700, 0 liquid feed 2 at 298. I K to Stage 24, J; = 50, 0, 50 Column 2: s = 12 c = 3: benzene, phenol, cyclohcxanc P = 1.013, total condenser, R = 3.0, partial t&oiler, B = 699.9 feed I to Stage 10, bottoms of Column I
Convergence of the initial problem is always very rapid (two to five iterations is quite typical) at values of EM” less than or equal to 0.1. Converging the quad-diagonal (QD) system leads to a small reduction in the number of Newton iterations required to converge the initial problem. However, this method sometimes requires a great many QD iterations. With hindsight, it would appear better to attempt only two or three QD iterations. Since the initialization procedure described above and the Wang and Henke method used by Waybum and Seader [lo] are equivalent at EM” = 1.0, it is pertinent to ask why the initial estimates appear to be so much better at very low values of the efficiency. There are two related reasons for this: (a) at extremely low values of the efficiency the composition and temperature profiles are virtually flat (with small step changes at feed stages, reboilers and total condensers). Sharp gradients are avoided; thus, there is no serious tendency for Newton’s method to make ridiculous corrections to the set of independent variables;
Table 2. Summary of results Iteration count efficiency value
Problem 1 2 3 4 5 6 7
0.1 214 313 313 213 315 3/5 5/6
0.4 3/4 3/4 617 313 515 414 518
1.0 615 415 617 415 516 315 7/Fail
Total iterations ll/ll IO/l2 15117 9110 13/14 IO/l2 17/-
102
D. J. VICKERYet ol.
(b) the nonlinearity of the MESH equations is reduced quite considerably by multiplying the K-values by a small positive number, EMV. It is better to use equation (7) to generate initial estimates of (X) in subsequent problems. The intermediate problem at EMV= 0.7 was omitted in the calculations reported here, thereby demonstrating that the solution at EMV= 0.4 and equation (7) provides a sufficiently good estimate of the solution at EMV= 1 (at least in these problems). The simpler method of setting the initial (X) vector equal to the solution of the previous problem would have been successful in problem 7 if the solution at EMV= 0.7 had been obtained. On rare occasions, we have experienced difficulty in going from the problem at a parameter value of 0.1 to a parameter value of 0.4 in one step. Examples of problems in this category involved strongly nonideal systems in interlinked systems of columns [19]. In these cases, we have obtained convergence quite easily to talking more (smaller) steps in the parameter, EM”. Another case in point is the dehydration of ethanol using benzene as posed by Magrmssen et al. [I 91. The value of EM” at which the calculations were halted was 1.0 in all cases. However, if the stage efficiency is known to be 60% (for example), we could stop the calculations when we reach that point. If the stage efficiency differs in different parts of the column or even from tray to tray and component to component, that too can be taken into account simply by preventing further increases in the efficiency when the prescribed maximum efficiency for that stage/component has been reached. Other
notes
we have described is not capable of solving problems in which a product purity is specified. To deal with problems of this kind, we would have to convert the problem to one with standard specifications (using, for example, a shortcut method), obtain the solution to that problem and then use reflux ratio and/or bottoms flow rate as parameters in a parametric continuation method in order to find the values that met the product purity specifications.
whether by trying to solve equilibrium stage models of real “far-from-equilibrium” operations, we have been making multicomponent separation process problems harder to solve than they really needed to be. Acknowledgements-DJV was the recipient of an NSF fellowship while this work was undertaken. We would like to thank the reviewer who drew our attention to the work of Sereno [15] and Professor D. W. T. Rippin who kindly supplied us with a copy of Miiller’s thesis. NOMENCLATURE
EM”= Murphrec vapor efficiency E = Equilibrium or efficiency equation F = Feed flow rate (F) = Column matrix of mass balance and efficiency equations (F) = Column matrix of MESH equations H = Enthalpy H = Heat balance equation [J] = Jacobian matrix K = K-value L = Liquid flow rate M = Component material balance Q = Heat load on stage r = Ratio of side stream flow to intcrstage flow s = Number of stages S = Summation equation V = Vapor flow rate (X) = Column matrix of mol fractions (X) = Column matrix of independent variables x = Liquid phase mol fractions y = Vapor phase mol fractions .z = Feed mol fractions Subscripts i = Component index j = Stage index k = Iteration counter
The method
CONCLUDING
REMARKS
The “efficient continuation method” is an effective method for solving equilibrium stage separation process problems that are diflicult to solve in other ways. The method is very easily implemented in an existing code that solves the MESH equations using Newton’s method. The initialization procedure is completely automatic; no initial estimates need be supplied by the user of the method. The fact that “inefficient equilibrium stage” separation calculations are much easier to converge than true equilibrium stage calculations makes us wonder
REFERENCES
1. E. J. Henley and J. D. Seader, Equilibrium Sruge Separafion Operations in Chemical Engineering. Wiley, New York (1981). 2. C. J. Ring, Separation Processes. McGraw-Hill, New York (1980). 3. C. D. Holland, Fundamentals of Multico~ponent Distillation. McGraw-Hill, New York (1981). 4. Y. L. Wang and J. C. Wang, In Proceedings of First Conference on Foundations of Computer-Aided Process Design, (W. D. Seider and R. S. H. Mah, Eds). AIChE
(1980). 5. J. D.‘Seader, C/rem. Eng. Ed. 19, 88 (1985). 6. P. A. Whitehouse, PhD Thesis in Chemical Engineering, University of Manchester Institute of Science &d Technoloav (1964). 7. L. Gy Naphtali, Paper presented at AIChE S6th Nat1 Mtg, San Francisco (1965). 8. L. M. Naphtali and D. P. Sandholm, AIChE .Jll7,148 (1971). 9. Goldstein and Stanfield, Ind. Engng Chem. Process Des. Dev. 9, 78 (1970). 10. T. L. Waybum and J. D. Scader, In Proceedings of Second Conference on Foundations of Computer-Aided Process Design. (A. W. Westerberg and H. H. Chien, Eds). CACHE (1983).
Stage separation 11. T. L. Waybum and J. D. Seader, Comput. them. Engng 11, 7 (1987). 12. J. D. Seader, Computer MoaWing of Chemical Processes. AIChE Monograph (1985). 13. G. Byrne and L. Baird, Comput. them. Engng 9, 593 (1985). 14. J. Jehnek, V. Hlavacek and M. Kubicck, Chem. Engng Sci. 28, 1555 (1973). 15. F. R. Mtiller. PhD Thesis, ETH Zurich (1979). 16. A. M. Sereno, Doctoral Thesis, University of Porto, Portugal (1985). 17. D. J. Vickery and R. Taylor, AIChE JI 32, 547 (1986). 18. D. J. Vickery and R. Taylor, Paper presented at New Orleans AIChE Nat1 Mtg (1986). 19. D. J. Vickery, T. L. Waybum and R. Taylor, Proc. IChemE Symp. DistiIlation and Absorption 2, B305 (1987). 20. M. Kubicek, V. Hlavacek and F. Prochaska, Chem. Engng Sci. 31, 277 (1976). 21. T. Magnussen, M. Michelson and A. Fmdenslund, Proc. IChemE Symp. Distillation, London (1979). 22. W. F. Huber, Hydrocarbon Process. August 121 (1977).
process
103
problems
we may write: (F) = [AlBCO](X) -(R)
B,
C,
D,
A,
Bz
Cr
Dr
A,
B,
Cr
J’,+,&j+,xi,j+,
--L,~,~c~,~~,=+~.
(Al)
If the flow rates and K-values are known we may solve the resulting linear equations for the xi . What makes these methods so elegant is that the toe $. crent matrix of each linear system is tridiagonal; linear systems with a tridiagonal coefficient matrix can be solved very easily using the well known Thomas algorithm [l, 21. Murphree efficiencies are not normally included in the Wang and Henke method (or other related methods) because if equation (4) is used to eliminate the y,,, the upper triangular portion of the coefficient matrix is filled out [22]. We could use equation (4) to eliminate the xy from the component mass balances and solve the resulting tridiagonal system for the. vapor phase compositions instead but this is not recommended ([2], p. 472). Rather than combine equations (2) and (4) into a single set of equations, we define a column matrix of discrepancy functions (QT = (MT. I, El.I
. MY,j_IvEn,]-t1M>vEu, x M;j+,,E,,+,...M;,,~,,),
where M$ is the component material balance equation (2) and E,. represents equation (3a). Each pair of equations depends on only four mol fractions: .x~,~_,, y,, xi/ and y,,, , . Thus, if we define a column matrix of mol fractions (X) by: (XY = (Yr I, xf,l~~.Y~.j-l~xi.j-l, x Yij9 +Yi,j+
D,
A,_,’
B,_;
C,-;
A,-,
BN--I
D,_, G-1
A,
BN
where N (=2s) is the order of the matrix. or! the odd numbered rows of [ABCD] the elements are given by: = --&.I
B2j_1=(1+r~)~
A Quad-diagonal Algorithm for Solving Equilibrium Stage Separation Process Problems Incorporating Stage Eficiencies Many of the methods that have been developed for solving the MESH equations (the Wang and Henke method for distillation, the Bumingham and Otto algorithm for absorption and others-see Henley and Seader [l] Chap. 15; or King [2] Chap. 10) make use of the fact that when the Equilibrium relations (3) are used to eliminate the y,, from the component mass balances (2) we have: (Vjx;I+Lj)xij-
(A2)
[ABCD]=
A,_, APPENDIX
= (0).
With the equations and variables ordered in this way, the coefficient matrix [ABCD] has four adjacent diagonals with the additional diagonal in the upper triangular portion of the matrix:
I3 4,j+ I ’ ’ ’ Yi,.* x,3),
C,_,
=(l
+r:)Lj
D~1-r = -J’j+l
(A3)
(these elements being derived from the component materialbalances) whereas on the even numbered rows we have: Arj= -1 Bz, = E?‘“K. c,
= (;‘-
&)
D,=O
(A4)
which belong to the equilibrium relations. For a total condenser set EM” = 0. For a total n&oiler, set EMV= 0 and B,= 1. The right hand side matrix (R) has elements: R, _ , =
4.2,
R2, = 0.
(A5)
The indexj corresponds to the stage number and stages are numbered down from the top. This linear system of equations can be solved for the x,~ and y,, very easily using Gaussian elimination. We have used this method of solving the component material balances and efficiency equations in place of the corresponding steps in the Wang and Henke method for distillation simulations. While not as fast as the conventional tridiagonal matrix algorithm (the quad-diagonal matrix has twice the number of rows as the tridiagonal matrix in the Wang and Henke method), it permits the straightforward solution of problems in which the stage efficiency is not unity. We estimate that the quad-diagonal algorithm will be faster than Huber’s [22] method of dealing with stage efficiencies (with an upper Hessenberg matrix of order s) if the number of stages exceeds nine. An interactive computer program implementing this algorithm for distillation calculations is available from the authors.