Ind. Eng. Chem. Process Des. Dev. 1904, 2 3 , 609-618
609
Sato, T. J. Inorg. Nucl. Chem. 1982, 2 4 , 1267. Sato, T. J. Inorg. Nucl. Chem. 1983. 2 5 , 441. Sharp, E. M.; Smutza, M. Ind. Eng. Chem. ProcessDes. D e v . 1965, 4 , 49. Wallace, D., Jr. Nod. Sc/. Eng. 1962, 14, 159.
McKay, H. A., Jr. Trans. Farahy Soc. 1958, 5 2 , 39. Robinson, C. Q.; Paynter, J. C. “Proceedh~gs,International Solvent Extraction Conference”; The Hague, 1971; paper no. 214. Rozen, A. M.; Khorkhorina, L. P. J. Inorg. Chem. USSR 1057, 2(8). 1956. Rozen, A. M.; Khorkhorina, L. P.; Karpacheva, S. M.; Agashkina, Q. D. Rad/&h/mlya 1962. 4(5). 591. Rozen, A. M.; Khorkhorina, L. P.; Karpacheva, S. M. J . Inorg. Chem. USSR 1957, 2 . 1441.
Received for review November 4, 1982 Accepted September 20, 1983
Local Models for Representing Phase Equilibria I n Multicomponent, Nonideal Vapor-Liquid and Liquid-Liquid Systems. 2. Application to Process Design Eldred H. Chlmowitz,+ Sandro Macchletto, Thomas F. Anderson,” and Leroy F. Stutzman Department of Chemical Englneerlng, Unlvershy of Connectlcut, Storrs, Connecticut 06268
Local thermodynamic and physical property models are used in computer-aided design of chemical processes. These models are especially sulted to numerical methods that require partial derivatives to obtain a solution. The algorithms developed here substantiilly reduce computational time while maintaining good convergence properties. Numerous phase-equlllbrium problems have been solved wlth these algorithms, which show them to be superior to conventional ones. Although attention has been restricted to small-scale problems in this work, local models can be readily used in very large design problems.
Introduction The simulation of industrial chemical processes requires that thermophysical-property equations be solved together with mass balances, energy balances, design equations, etc. Traditionally, these thermophysical properties are provided by specialized subroutines; numerous subroutines are made available to cover different ranges of temperature and pressure and different types of mixtures. A collection of these subroutines is known as a thermodynamic and physical property data base (TPPD). Since the relationships between the thermophysical properties and the process variables (temperature, pressure, compositions, etc.) are usually highly nonlinear and in many cases include implicit equations (e.g., density in an equation state), subroutines are written to provide “point values” of thermophysical properties at specified conditions. Usually, only these “point values” are incorporated in the process calculations. This conventional approach has two disadvantages. First, since only point values are used in the computational process, the dependence of these properties on temperature, pressure, and composition tends to be neglected. In some cases this may lead to a solution algorithm which converges slowly or not at all. The second disadvantage is that the rigorous thermophysical properties subroutines must be repeatedly accessed; the number of times these routines are accessed is proportional to the number of iterations required for the computational method to converge. Experience has shown that up to 80% of the total simulation cost may be due entirely to the evaluation of thermophysical properties.
A more recent concept for process simulation combines the rigorous thermophysical property equations with the process model equations. The solution is still iterative (e.g., a Newton-Raphson method may be employed), but all the equations are solved simultaneously. This approach also has several problems. The number of equations is dramatically increased. The enlarged set of equations poses a difficult problem to solve due to the highly nonlinear nature of most thermophysical models, and the analytic evaluation of the numerous partial derivatives can become very cumbersome. In addition, the partial derivatives must be redeveloped each time a new or different property relation is used. Westeberg et al. (1979) have estimated that evaluation of thermodynamic properties by this approach still accounts for nearly 80% of the computational cost. A third approach to process simulation has been suggested in the recent literature by Hutchison and Shewchuk (1974),Leesley and Heyen (1977),Boston and Britt (1978), Barrett and Walsh (1979), and Boston (1979). It involves the use of approximate models for representing thermophysical properties and the restructuring of the calculational procedure into two levels, resulting in a two-tiered approach. With this approach a sequence of problems is solved which has, in the limit, the same solution as the original one posed. In particular, the rigorous correlations for thermophysical properties are approximated by simple, local models. Parameters in these models are obtained from rigorous values provided by the thermodynamic and physical property data base. These parameters are either estimated or calculated initially, then updated, if necessary, at each solution of the current simulation problem. The two-tiered approach possesses several important advantages. The total number of rigorous thermophysical property evaluations can be substantially reduced. The local models can easily be incorporated into the process
‘Department of Chemical Engineering,University of Rochester, Rochester, NY 14627. 0196-4305/84/1123-0609$01.50/0
0
1984 American Chemical Society
610
Ind. Eng. Chem. Process Des. Dev., Vol. 23,No. 3, 1984
Table I. Various Local Models for Vapor-Liquid Equilibrium Ratios or Relative Volatilitiea in Nonideal Mixtures model approximation equations parameters no. model type A r e f , j j= 1, ..., 3 I relative volatility A 1., J j = 1,..., 3 I1 relative volatility Aref,jj = I,..., 3 ~~~~
I11
IV V
~
Ai, Ai,j Ai,j A I,J ..
equilibrium ratios equilibrium ratios equilibrium ratios (ideal model)
Table 11. Equations for Single-Stage Multicomponent Vapor-Liquid Equilibrium Flash Process Process Equations phase equilibrium yi = Kixi ( i = 1, ..., n ) (1) component x ~+ L y i V = z i F ( i = 1, ..., n ) (2) mass balance total mass balance L + V = F (3) enthalpy balance HLL t H V V = H F F + Q (4) constitutive zxi = z y i (5) czi= 1 (6) TP Property Equations K i = Ki(T,P,x,y) ( i = 1, ..., n ) ( 7 ) HL = H L ( T,P,x1 (8) H v = Hv(T,P,Y) (9) HF = HF(TFJ'FJ) (10)
model equations and their form is independent of the particular rigorous method used to obtain values for thermodynamic and physical properties. The key to using this approach lies in the formulation of accurate yet simple local models to represent the various thermophysical properties. In part 1 (Chimowitz et al., 1983) we presented various simple models for the equilibrium ratios (K values) in multicomponent vapor-liquid and liquid-liquid systems. These local models are presented in Table I. Here, we will illustrate how these simple models may be applied to basic equilibrium calculations and compare the results with established methods for solving these problems. In part 3 (Chimowitz et al., 1984) we extend the use of local models to larger systems and describe in detail efficient methods for initializing and updating parameters.
Phase-Equilibrium Calculations The set of equations which form the process model for a single-stage multicomponent equilibrium flash unit is shown in Table 11. For a feed mixture consisting of n components, there are a total of 3n + 7 equations-2n + 4 for the process model itself, given by eq 1-6 and n 3 equations for the thermophysical properties, eq 7-10. There are 4n 11variables and therefore, n + 4 degress of freedom exist for the most general problem. Usually, feed characteristics (composition, flow rate, temperature, and pressure) are specified, which take up n 2 degrees of freedom. Thus, six common types of flash calculations can be formulated depending on which two of the remaining variables (temperature, pressure, heat load, or vapor fraction) are specified. The remaining variables then become dependent quantities. The bubble and dew point calculations are merely special cases of the flash problem where the vapor fraction is set either to zero or to one. Furthermore, a liquid-liquid equilibrium stage can be described by similar equations. Additional phase-equilibrium problems may be formulated by considering the underdetermined flash calculation,
+
+
+
j j j j
=
1,..., 4
= 1,..., 3 = 1,..., 3 = 1 , ..., 2
i.e., the case where the given specifications do not fix all of the degrees of freedom in the system, while at the same time some overall ojective is to be met. For example, we may desire to determine the temperature which leads to the best purity and/or recovery of a particular component. This is a typical optimization problem and is not easily solved by conventional algorithms for phase-equilibrium calculations. A Newton-Raphson based method for solving the equations outlined above requires the partial derivatives with respect to all of the independent variables, e.g., d K i / d T , dKi/axi. These quantities can be found in one of two ways: either by perturbation or by analytic differentiation of the rigorous thermodynamic correlations. Perturbation is an expensive process, especially if all the composition derivatives are required at each stage of the Newton-Raphson iteration. In addition, a suitable perturbation step size must be chosen; this may be difficult for components present in small quantities. Developing the complex expressions for analytic derivative expressions also has a serious disadvantage; it is a time-consuming and error-prone task. Because of these difficulties,solution procedures for the single-stage equilibrium process have generally avoided applying the Newton-Raphson method simultaneously to all variables and equations. Instead, the original equations are transformed into a modified set of nonlinear equations given below.
J/ = CY&
+ (1- CY)&
- HF - q = 0
(12) Equation 11 is known as the Rachford-Rice equation (Rachford and Rice, 1952) and eq 12 is merely an enthalpy balance; CY is the vapor to feed ratio and q = Q/F. The solution of these equations is achieved by assuming that the functions given by eq 11 and 12 are dependent only on the variables CY and T . When the Newton-Raphson method is used to solve this considerably simpler set of equations, derivatives of both functions with respect to CY are obtained analytically and derivatives with respect to T are obtained by perturbation. The K values are considered independent of composition during each Newton-Raphson iteration and are recomputed after new values of CY and Tare obtained. Similar seta of equations are also used in bubble-point and dew-point calculations as well as liquid-liquid equilibrium calculations. The reader is referred to the recent work of Prausnitz et al. (1980) for a more detailed description of these methods. The principal disadvantage to the above method occurs when it is applied to highly nonideal systems. Then, the composition effects on the K values can be sufficiently strong so as to lead to slow convergence, and in the case of liquid-liquid equilibria such an algorithm may not converge at all. To improve on even these simple equi-
Ind. Eng. Chem. Process Des. Dev., Vol. 23, No. 3, 1984 611
R I O O R O U S MODEL
NON-LINEAR
Figure 1. Linear and nonlinear approximationsto a rigorous model.
librium Calculations, a simultaneous solution method based on the concept of local models is proposed.
Local Models in Phase Equilibrium Calculations As an alternative to the conventional algorithms for solving phase equilibrium problem, we propose that local models be used to represent the equilibrium relations. These local models are combined with the other process model equations to yield a single set of equations which is solved simultaneously. This results in a two-tiered approach. Rigorous thermophysical properties are calculated in the outside tier of the algorithm. These properties are used to estimate the parameters in the local models. The insider tier is used to find the solution of what we term the approximate process problem, which can be either a design, a simulation, or an optimization problem. A quantitative understanding of this procedure can be obtained by considering Figure 1. The continuous line represents a nonlinear function, f ( x ) , whose zero, x * , is sought. This represents the solution of the rigorous problem. If an accurate, approximate function can be developed called g ( x ) (given by the discontinuous line), then the root of g(n),xl,will be very close to the solution of the original problem, x*. Finding x1 from a starting point xo may require that an iterative procedure be applied to g(x). However, the important point is that the original function and ita derivatives are not required during these iterations. At x1 the local model g(x) is updated giving a more accurate approximation to f ( x ) and the process continued until convergence is achieved. In the multivariate situation, f ( x ) becomes a vectorvalued function F ( x )of equations describing the process (e.g., equilibrium relations, energy balances, etc.) and the single derivative is replaced by the Jacobian matrix of partial derivatives. It is computationally expensive and difficult to evaluate the partial-derivative matrix (Jacobian) of F ( x )when the rigorous thermophysical property relations are included in the system of equations. However, if local models are used for these properties, the evaluation of these derivatives is quite feasible. Thus, a solution to the rigorous problem is obtained by the use of local thermodynamic approximations in an analogous manner
to the solution of the single equation f ( x ) = 0. A similar approach, however, not using derivatives of the local models, was proposed by Boston and co-workers (1978,1979). For a detailed description of their method (called the inside-out algorithm), the reader is referred to the original series of papers. Here we note a number of relevant points about the inside-out method. The algorithm is two-tier in structure utilizing rigorous thermodynamic properties in the outside tier only. The traditional role of the variables is reversed; the parameters of the local models become the iteration variables in the outside tier; the usual (primitive) variables, such as composition and temperature are evaluated in the inside loop as dependent variables. The rationale for this reversal is that the model parameters are not very sensitive to the primitive variables, and hence, an accurate initialization is not necessary to assure convergence. Drawbacks of the algorithm appear to be that the use of Broyden’s method for updating the parameters (Broyden 1965) prevents the use of this type of update for any large problem (e.g., distillation),as noted by Boston and Britt (1978). Furthermore, in their original article these authors assume fixed process specifications. In particular, inlet variables must be specified; a different modification of the basic algorithm was needed for each of the six common types of the single-stageflash problem. We believe that a transformation of the basic flash equations as described by Boston and Britt is not really required and introduces unnecessary complexity into the problem. The method we propose considers the process model equations in the original form without introducing additional variables. The simultaneous solution method permits treating all types of flash calculations the same way and does not limit the possible combinations of problem specifications. The succeas of our approach requires more than accurate local equilibrium models. Methods are required for initializing and updating the parameters in these models as well as methods for converging the system of nonlinear equations. The implementation of these concepts into useful algorithms is described in the following sections, leading to what we have called the P-DELTA method. P-DELTA is an acronym for Process Design by Limiting Thermodynamic Approximations. Figure 2 shows a schematic diagram which outlines the essential features of the P-DELTA method. Application to Dew-Point Calculations The first computational problem that we have applied our method to is the dew-point problem. Many algorithms that have been developed for the solution of complex processes (e.g. multicomponent distillation) require repetitive bubble-point and/or dew-point calculations to be performed as subproblems. Hence, any improvement in these basic algorithms can result in significant savings on more complex process calculations. For the dew-point problem the local model used to represent the equilibrium ratio of each component in a mixture is given by In Ki = In
(:)
+ ailxi + ai2
(i = 1, ..., n)
(13)
This corresponds to model IV in Table I. In eq 13 f p is the pure-component reference fugacity of component i, which is correlated directly as a function of temperature in our data base. With this local model the activity coefficient is given by In yi = ailxi ai2 (14) and the parameters ail and ai2 are evaluated from two
+
612
Ind. Eng. Chem. Process Des. Dev., Vol. 23,No. 3, 1984 S C H E M A T I C O F THE P - D E L T A * A L G O R I T H M
Table 111. Average Number of Iterations to Determine the Dew Point Temperature o f a Mixture (10 Points for Each Mixture)
Pmixture
BUDET DELTA
A. six components, ideal, wide boiling; nbutane (l), isobutane (2), n-pentane (3), nhexane (4), n-heptane (5), n-hexadecane (6); pressure = 15 bar B. three components, mildly nonideal; ethanol (l), 1propanol (2), cyclohexane ( 3 ) ; pressure = 10 bar C. three components, strongly nonideal; acetone (l), water (2), acetonitrile ( 3 ) ;pressure = 1 bar
4.5
cnoccas I N T E R Y I Or LOCAL YODELS
5.0
t IUlTlALlZE DATA-BAS€
7.9
5.7
10.8
5.7
1I
*PP PROCESS RoYIYATE
PARAMETERS
YODEL
rigorous K value evaluations done on the phase envelope as described in part 1. This model is then used to solve the dew-point problem which is given as the following set of equations. n
-
I
*PROCESS
cxi- 1=0
NO
CONVERGENCE
I
DESIGN BY L I M I T I N G T H E R M O D Y N A M I C A P P R O X I M A T I O N
Figure 2. Schematic of the P-DELTA Method.
i=l
A Newton-Raphson iteration applied to this set of equations is straightforward. The following partial derivatives are required and can be determined analytically
A comparison with a conventional dew-point algorithm (BUDET) taken from the text by Prausnitz et al. (1980) is shown in Table 111. Three mixtures of increasing nonideality are compared. The conventional algorithm solves eq 11 using only temperature derivatives of the K values obtained by perturbation. While no significant difference is noticed for the most ideal case, the use of composition dependent models requires up to 50% fewer iterations for the most nonideal case studied. For these examples, the dewpoint calculations were made by first establishingthe composition range in advance. Parameters were fit once over this range; it was not necessary to update the parameters during the calculation procedure since a good accuracy was achieved over the entire range. This illustrates an advantage of composition-dependent local models.
If the composition range is not known in advance, a slightly different procedure is adopted. Coefficients ail and ai2in the activity coefficient expressions are set to zero at the start in the outside loop. This corresponds to an ideal mixture approximation. In subsequent iterations a single equilibrium-ratio evaluation permits calculating the coefficient ui2. Good convergences have been achieved in all problems with this method. Note that this procedure differs from the assumption of constant relative volatility. The performance of our method is again compared to BUDET, the algorithm given by Prausnitz et al. (1980). We have also compared our procedure with the inside-out method of Boston and Britt (1978), with the quasi-linear (QLIN) method of Hutchinson and Shewchuk (1974) and with a standard Newton-Raphson iteration, where all derivatives are computed by perturbation. Results are shown in Table IV. Tests were carried out for a number of mixture compositions. The same starting points and convergence tolerances were used for all methods. For both the P-DELTA and the inside-out methods, direct substitution was used in updating the local model parameters. Our method and the inside-out one were generally comparable. The BUDET routine required on the average twice their number of K-value evaluations; the quasilinear method fared the worst, because it requires a complete
Table IV. Average Number of Rigorous K Evaluations To Determine Bubble and Dew Point Temperature of a Mixture (10 Points for Each Mixture) bubble point BUDET
C(P= 1 bar) C(P= 10 bar)
7.8 8.4 8.0 7.6
5.9 5 .I
overall average
7.9
6.0
A B
a
insideout
mixture
Does not converge in 30 iterations.
5.8
6.8
dew point
P-DELTA BUDET 3.1 3.1 3.0
insideout
P-DELTA
4.0
9.8 13.2 24.8 15.0
10.5 8.1 13.0 8.0
4.8 8 .O 11.2 6.3
3.6
15.7
9.9
1.6
QLIN
NR (perturbation)
38.1
25.0 23.5 21.5 21.0
49.3
22.75
a
49.6
60.1
Ind. Eng. Chem. Process Des. Dev., Vol. 23, No. 3, 1984 813 S C H E M A T I C OF P-DELTA A L G O R I T H M ( I ) DEFINE OBJECTIVE FUNCTION FOR FLASH CALCULATION
~(uI,! -- = (2) GENERATE
(x,y,T,P,a)
Vu gi(IJ, Vuh,(g)
V,$(g),
FROM LOCAL MODELS
(3) INITIALIZE MODEL PARAMETERS IN Vug,(lJ), INITIALIZE
:$
AK
K = iteration counter
hK= W4) SOLVE OPTIMIZATION
Model
MINIMIZE
SUBJECT TO
Parameters
$ (2)
g,(y):o
,m)
(I=I,
h , & I ~ o (]=I,Z,
gKI
&
,P)
CONVERGENCE PRINT RESULTS
- ( 6 ) UPDATE MODEL PARAMETERS BASED UPON SOLUTION IN ( 4 ) '
AKtI
I
-
I I I
I
XA
XD
Figure 4. Outline of the P-DELTA algorithm applied to the optimization of a single-stage flash.
XF
XB
LIQUID PHASE COMPOSITION
Figure 3. Two-phase envelope for a binary system.
bubble point evaluation at every iteration. In fact, for the wide-boiling hydrocarbon mixture (150 K between bubble and dew point), the quasilinear method failed to converge in 30 iterations. Application to Vapor-Liquid Flash Calculations For the next problem we examined the optimal design of a single-stage multicomponent, equilibrium flash unit. If the single-stage flash problem is solved by an optimization method, both the usual adiabatic and isothermal flash as well as the underspecified design problem can be solved with the same program. Only a simple rearrangement of the objective function and constraints is required to treat the various types. The underspecified problem is the more interesting one. The objective is usually to find the value of the design variable (often temperature) that best meets recovery and/or purity targets on any of the components in the feed. A complete definition of the equations used for this problem is given in the Appendix. A qualitative insight into this problem can be gained by examining the phase diagram for a typical binary system shown in Figure 3. For a feed composition xF with bubble point temperature T B and dew point temperatured TD, the objective of the calculations is to locate the temperature lying along the phase envelope between TBand T D that best meets purity and/or recovery objectives on any of the components in the exit streams from the unit. For light components, high purity specifications in the exit vapor stream are favored by working at temperatures close to TB, while high recovery specifications force the temperature closer to TD.Intermediate to these two temperatures is the optimal one for the imposed design objectives. The local K-value models are used to give an accurate mathematical representation of the phase envelope in a region including the composition range spanned by X F and xD. For this problem, the range over which the local model must apply is well-known. Figure 3 shows that for a feed of known composition, it is possible to bracket the region in which the solution of the flash problem must lie. Multicomponent mixtures can be treated analogously, the two-phase region lying between the bubble and dew points of the feed. By similar reasoning, establishing the problem solution region for multiple interconnected flash units is straightforward. Hence, for this class of problem, model initialization is an easy task. The composition extremes of the feasible region are first established, and more
points may be fixed by linear interpolation at equally spaced intervals. Bubble point calculations are required at each point to provide the equilibrium data for the parameter regression. We have generally found four points to be satisfactory for all of the models. Figure 3 can again be used to illustrate this procedure; the two end points are identified as xF and xD, and two intermediate points are located along the vector connecting xp and xD. The K values determined at these compositions are used 'to estimate the local-model parameters by simple linear regression. For the type of problems considered here, parameter updating is a straightforward task. After the parameters have been initialized, the first set of inside-tier iterations can be performed, providing a first solution to the problem. In the outside tier, this solution is then used to update the model parameters in a neighborhood about this point. This is accomplished by repeating the regression, with equilibrium data at the current solution replacing the point (from the previous set) most distant from it. Once new parameters are established, the inside calculations are repeated, etc. For the following examples we have always retained the two extreme points and rejected an interior one. This tended to stabilize the change in the parameters from one iteration to the next. However, this aspect of the procedure is very flexible and other approaches are likely to work well. This aspect will be covered in more detail in a subsequent publication (part 3, Chimowitz et al., 1984). A schematic of the algorithm for solving the flash problem is shown in Figure 4. represents the objective function, girepresents the equality constraints (mass, energy balances, phase equilibrium relationships), and hi represents the inequality constraints (e.g., nonnegativity of the mole fractions). The inside loop consists of the calculations performed in step 4. In the inside loop, a nonlinear program is solved. For this calculation we have used the Han-Powell method (Han, 1975; Powell, 1978), which requires partial derivatives of the objective function and constraints. The outside loop calculations (model updating) are done in step 6; then step 4 is repeated until convergence is achieved. Convergence is assumed when two successive passes through step 4 produce results that agree to within a prespecified tolerance. The first flash calculation was carried out on a moderately nonideal ternary system consisting of ethanol (11, benzene (2), and n-hexane (3) at a pressure of 10 bar. The design objective was to locate the temperature that gave 80% recovery of the most volatile component, n-hexane, in the exit vapor stream. Table V shows the results of the calculations using two different local models to represent
614
Ind. Eng. Chem. Process Des. Dev., Vol. 23,No. 3, 1984
Table V. Comparison of P-DELTA Method with Various h a 1 Models for the Ethanol ( 1 ), Benzene (2), n-Hexane ( 3 ) System at a Pressure of 10 bar. Objective Function Is That Giving 80%Recovery of n-Hexane in Exit Vapor Stream. (Initial Conditions: xi = 0.5, yi = 0.5, i = 1, ..., 3 ; a = 0.5, T = 450. Feed Composition: z(l)= 0.2, 2(21= 0.3. z ( 3 ) = 0.51 model I NGRID V/F T,K
a
Xl XZ x3
Y1 Yz Y3
no, of calls to K-value routine no. of outside loops f(objective function)
4 0.83 424.82 0.032 0.38 0.59 0.24 0.28 0.48 37
model V
FLASH solution" iteration iteration iteration at T = 1 2 3 343.5 K
4 0.83 424.7 0.032 0.38 0.59 0.23 0.28 0.48 37
oi
Xl x3
Y1 Yz Y3
2 0.45 X 10.'
model I
XI X2 x3 Y1
Yl
Ye
no. of calls to K-value routine no. of outside loops f(objective function)
V/F
XI
2 0.22 X
T,K
Table VI. Comparison of Local Models I (Composition and Temperature Dependence) and V (Temperature Only Dependence) for Example 1 (Initial Conditions: x i = 0.5, y i = 0.5, i = 1, ..., 3, a = 0.5, T = 360. Feed Composition: z(l)= 0 . 2 , 2 ( 2 ) = 0.3, z ( 3 ) = 0.5) update type NGRID a = V/F T,K
Table VII. Iteration History of Inside Loop for P-DELTA Method with Model I1 for Highly Non-Ideal System Consisting of Acetone(l), Acetonitrile( 2), Water(3) at a Pressure of 1 Bar (Objective Function: Minimize f = (ayammne - 0.24)* t (Yaqbne- 0.65)2; Initial Conditions: x i = 0.5, y i 7 0.5, z = 1, ..., 3, a = 0.5, T = 360. Feed Composition: z ( l ) = 0.3, 2 ( 2 ) = 0.1, 2 ( 3 ) = 0.6)
single point 4 0.83 424.8 0.032 0.38 0.59 0.24 0.28 0.48 37 2 0.45 X
model V single point 4 0.83 424.6 0.032 0.38 0.58 0.23 0.28 0.48 55 5 0.24 X
lo-"
the K values. Both the temperature-only model (model V) and the composition dependent model (model I) gave accurate results with only two model updates in the outside loop. Each pass through the outside loop requires a single bubble-point calculation at the current composition conditions. It should be noted that while the model parameters were initialized over the bubble-dew point range of the feed stream, the initial conditions for the optimization algorithm (see Table V) could lie well outside of the twophase region. This indicates a subtle advantage of the quasi-Newton optimization algorithm, namely that the calculations can proceed along an infeasible path during the intermediate stages of the computations without interruption, in contrast to many other feasible-path methods. This point is further illustrated in a later example dealing with two interconnected flash units. Table VI shows the results for the same calculations with different initial conditions. Here an initial temperature of 360 K was used,which is approximately 60 K below the bubble-point temperature of the mixture. The calculations with the composition-dependent local model (model I) converge just as rapidly as in the previous w e . However, model V (temperature dependence only) required more than double the number of outside loop iterations. The strongly nonideal system of acetone (l),acetonitrile (2), and water (3) at a pressure of 1bar was used for our second example. Here the design objective was to take a feed with a concentration of 30 mol 5% acetone and produce a vapor with a 65 mol % acetone concentration. In addition, an 80% recovery of acetone in the exit vapor stream
0.3324 0.1503 0.0861 0.7636 0.6007 0.1279 0.2714 342.6
0.3319 0.1533 0.0880 0.7587 0.5954 0.1241 0.2806 343.5
0.3341 0.1521 0.0878 0.7602 0.5948 0.1244 0.2808 343.5
0.3500 0.1424 0.0849 0.7727 0.5873 0.1275 0.2852 343.5
FLASH solution is that obtained using isothermal flash routine from Prausnitz et al. (1980) at temperature given in iteration 3.
is required. This represents a true optimization problem since both design objectives cannot be met simultaneously. By suitable weighting of the objective function, the relative importance of the competing objectives can be varied. Table VI1 shows the iteration history for three successive inside loop traverses. Model 11, a composition dependent local model, was used for the mixture K values. Three outside loop iterations are required to achieve convergence. After the first iteration the results are very close to the fmal solution; the remaining iterations serve only to refiie the answer. The results predicted by the P-DELTA method are in close agreement with a rigorous solution at the final temperature (last column in Table VII). To compare how efficiently different methods utilize the thermophysical properties data base, we have solved a number of optimization problems with our algorithm and with a commonly used alternative procedure. The alternative method consists of interfacing a direct search optimization method with subroutine modules that model the process. The search method then repetitively calls the process modules changing a prescribed set of decision variables that affect the state of the process. This is repeated until a set of design variables is found that minimizes some prespecified objective function. This procedure is computationally expensive (Westerberg, 1980), but in the absence of other proven techniques for process optimization, it is still frequently used (Martin and Gaddy, 1981; Parker and Hughes, 1981). Table VI11 shows the results of this comparison for a number of optimization problems ranging from a wideW i g , six-component hydrocarbon system to the strongly nonideal ternary system of the previous example. The design problem here reduces to a single variable optimization problem with the temperature of the flash unit being the single decision variable. For this optimization the Golden Section method was used to adjust the temperature for each isothermal flash calculation. The flash algorithm used was that described in the book by Prausnitz et al. (1980). The results show that the P-DELTA algorithm required approximately 80 to 90% fewer evaluations of the rigorous equilibrium ratios when compared to the direct search procedure. The advantage becomes more favorable as the system becomes more nonideal. Finally, we considered a two-stage flash process (Figure 5 ) with recycle of the liquid stream from the second stage. The system was a wide-boiling mixture consisting of n-
Ind. Eng. Chem. Process Des. Dev., Vol. 23, No. 3, 1984 615
Table VIII. Comparison of P-DELTA and Golden Section Methods for Solving Recovery-Purity Roblems in a Multicomponent Single-Stage Flash Unit
feed, mole fraction
problem description wide-boiling hydrocarbon n-butane isobutane n-pentane n-hexane heptane decane pressure = 15 bar, 96% recovery decane in liquid exit moderately nonideal ethanol benzene n-hexane pressure = 10 bar, 80% recovery n-hexane in vapor exit highly nonideal acetone acetonitrile water pressure = 1 bar, 80% recovery acetone, 65% purity in vapor exit highly nonideal acetone acetonitrite water pressure = 1 bar, 80% recovery acetone
no. of K-value evaluations for Golden Section
P-DELTA -% Golden Section
0.15 0.15 0.15 0.15 0.15 0.25
V
64
337
19
0.2 0.3 0.5
I1
37
235
15.7
0.3 0.1 0.6
I11
34
336
10.1
0.3 0.1 0.6
I11
34
303
11.2
VAPOR PRODUCT
I
no. of local model K-value (P-DELTA evaluations only) for P-DELTA
I
7 &*T2,h
s
LIQUID LIT * PRODUCT
Figure 5. Flowsheet of a two stage flash process.
pentane (l),n-hexane (2), and n-octane (3). The feed composition to the first unit was: z(1) = 0.4,z ( 2 ) = 0.3, and z ( 3 ) = 0.3. The design objective was to produce a vapor stream from unit 2 that yielded a 60% recovery of the most volatile component (n-pentane) at a purity of at least 75% by adjusting temperatures in both stages. Pressures in both stages were taken to be 1bar. There are two degrees of freedom to the problem formulated in these terms. For the P-DELTA algorithm the inside loop consisted of the, solution of a nonlinear program with 15 equality constraints and 17 variables. In this problem the two-phase region in the two stages can be significantly different. Hence, a wide-boiling feed can give rise to a
vapor product from unit 1 (i.e., feed to unit 2) with a very narrow boiling range and vice versa. This provides a challenge for the local models which have to be capable of representing these widely varying two-phase regions for the same mixture constituents. An important aspect of this problem is the manner in which the local models are initialized and updated. Mass balance restrictions dictate that the liquid-phase compositions in the two stages must lie between the two extreme points (see Figure 3) found by a single equilibrium step up and down the phase envelope from the feed stream condition. The temperature range between these two extrema is approximately 60 K. However, local models for unit 1are initialized between points 1and 2 and for unit 2 they are initializsd between points 2 and 3. Furthermore, the local model updating procedure occurs around two distinct points in these two mutually exclusive regions, one for each flash unit. This clearly illustrates one of the essential ideas behind the local model concept, namely that the local models can be used to accurately represent discrete regions along the two-phase envelope in a manner coneistent with the process requirements. The results of solving this problem by a number of different algorithms are shown in Table IX. This problem was also solved using two direct search optimization methods, Pattern Search (Hooke an Jeeves, 1986) and Rosenbrock‘s (Rosenbrock, 19601, with the temperatures in the two stages as decision variables. With both methods care must be taken to ensure that the temperatures lie within the two-phase region of the respective feeds. However, the feed compositions are not known in advance because of the presence of the recycle stream. Hence, the recycle stream must be “torn” and converged during each step of the optimization, which is computationally expensive.
Ind. Eng. Chem. Process Des. Dev., Vol. 23, No. 3, 1984
616
3 ACETIC ACID
Table IX. Optimization of Two-Stage Flash Process. Feed n-Pentane 40 mol %, n-Hexane 30 mol %, n-Octane 30 mol %. Objective Function, 60%Recovery, 75% Purity of n-Pentane from Second Column. Pressure both Columns 1 bar pattern variable P-DELTA search Rosenbrock 0.51 342.4 0.39
a1
Tl, K
R(recycle/ mole feed)
0.57 343.34 0.55
I
---
FEED COMPOSITIONS PRODUCT COMPOSITIONS TIE-LINES
0.53 342.7 0.435 1
0.24 0.33 0.43 0.59 0.34 0.07 0.45 324.6 0.47 0.41 0.12 0.75 0.24
Xl x2
x3
Yl Yz Y3
x4
x5 x6
Y4
Y5
0.01
Y6
no. of calls to K-value routine no. of isothermal flash calculations objective function value
77
__
0.22 0.33 0.45 0.58 0.35 0.07 0.38 324.45 0.47 0.42 0.12 0.748 0.24 0.0084 3453 164
0.17 X
10-6
0.23 0.33 0.44 0.59 0.34 0.07 0.43 324.55 0.47 0.41 0.12 0.748 0.243 0.0087 4830
0
0.28 X
10-3
0 100 2 WATER
80
60
40
MOLE PERCENT
Figure 6. Optimal design of a single stage liquid-liquid flash.
23 2
0.17 X
20
I BENZENE
10-4
The results in Table IX demonstrate that the design objectives used for this particular problem can be met simultaneously (the value of the objective function at the solution is zero). In general this would not be the case, requiring a compromise between the various requirements. The superiority of the P-DELTA method in terms of utilization of the thermodynamic data base is clearly evident; rigorous K values were only required 77 times. The direct search optimization methods required several thousand rigorous evaluations of equilibrium ratios.
Application to Liquid-Liquid Flash Calculations As a final example local models were used for simple liquid-liquid equilibrium calculations. In part 1 we presented a method for developing local models for ternary liquid-liquid systems. Here the models are used to solve a simple design problem involving such systems. The problem is similar to the vapor-liquid flash problem examined earlier. Given a feed of known composition, it is desired to determine the amount of solvent required to extract a specified amount of one of the componenh from the feed. For 3 components this problem can be illustrated with a ternary diagram. Figure 6 shows the system of
benzene (l),water (2), and acetic acid (3). Water is the solvent and is used for extracting acetic acid from the benzene-acetic acid feed. For a given amount of feed and solvent the corresponding overall mixture composition lies on a straight line connecting the water vertex and feed point. The compositions of the two phases in equilibrium are given by the intersection of the tie-line through this point and the binodal curve. Increased extraction requirements tend to shift the overall composition point toward the solvent (water) vertex (Le., points A and B in Figure 6). The equations modeling this process are similar to those for the vapor-liquid case, except the liquid-liquid local models are required for representing the phase equilibrium. The local models used in the following examples are a slight modification of those presented in part 1 and are shown in Table X. Parameters for the models were determined following the procedure described in part l. For this ternary liquid-liquid system the model is quite accurate over the entire composition region, and no parameter updating was required. A number of design problems were solved using 1mole of feed at various compositions and requiring different recoveries of acetic acid. The results are shown in Table X. In all cases convergence was achieved without difficulty. The results were checked against the equilibrium distribution estimated by a flash algorithm using rigorous models; agreement was generally within 1%. For a ternary system these particular calculations could also be done graphically, although this is a tedious task. For four or more components a graphical solution becomes impossible
Table X. Results for Singlestage Liquid-Liquid Recovery Problema (Initial Conditions: x,I = 0.8, xzl = 0.1,x31= 0.1, x," = 0.1, x , =~ 0.8, x,== 0.1; I Refers to Benzene-Rich Phase; I1 t o Water-Rich Phase) feed composition specified (1 mol) recovery mol of acetic of acetic solvent benzene acid acid % required 0.67 0.67 0.5 0.5 0.4
0.33 0.33 0.5 0.5 0.6
57.6 72.3 80.0 98.0 98.0
0.41 0.74 0.85 8.12 6.56
phase I x,
0.81 0.87 0.81 0.97 0.96
x*
0.018 0.01 0.017 0.005 0.005
phase I1 x3
x,
0.017 0.12 0.168 0.022 0.032
0.014 0.007 0.014 0.001 0.0015
x2
0.67 0.25 0.67 0.94 0.917
Liquid-liquid local models: In K , = 0.2497 - 4 . 4 8 8 ~ ~- ' 3~ . 0 5 9 (benzene); ~ ~ ~ ~ ~ In K , = -0.2406 + 0 . 9 3 9 2 +~ ~0 ~. 1~ 2 9 3 (acetic ~ ~ ~ ~acid).
0 . 3 8 5 5 ~ ~(water); "~ In K , = -0.0506 In K , = In ( x t l / x i l ) .
x,
iterations
0.32 0.246 0.32 0.056 0.082
+ 6.147x,12 -
13 40 14 50 41
Ind. Eng. Chem. Process Des. Dev., Vol. 23,No. 3, 1984 817
and a numerical procedure similar to the one described above would be required.
Discussion and Conclusions In this paper we have presented an approach for using thermodynamic properties in equation-oriented solution methods for solving process design problems. The method relies on the use of local models to approximate rigorous thermodynamic properties. We have implemented these concepts in conjunction with a quasi-Newton, nonlinear optimization algorithm and solved a number of problems from an optimization perspective. The objective of using local models is basically twofold: (i) to reduce the large computational overhead associated with the alternative use of rigorous thermodynamic evaluations, and (ii) to enable the process model to be formulated in a manner which allows the use of powerful equation-oriented algorithms for its solution. The results presented here indicate that progress toward these two objectives has been made. In comparison with alternative procedures, the new method presented is both robust and reduces the use of the thermodynamic data base by at least two orders of magnitude on a number of example problems. The use of both composition dependent and temperature only local models for vapor-liquid equilibria has been examined. However, for strongly nonideal systems, use of the composition-dependent local models is favored both for the increase in accuracy and better convergence that follows from their use. In the examples given her the local models were used a t pressures up to 15 bar with reliable results. However a t higher pressures, and certainly near the gas-liquid critical region, it is likely that the local models presented here would have to be modified in some manner. Although we have not yet addressed this problem, we feel that the modification would probably be to include an additional pressure-dependent term. The local model concept has also been successfully applied to simple liquid-liquid equilibrium computations. In this case composition dependency is essential for correctly representing phase equilibria. While the method has proven advantageous in the examples presented here, the greatest potential benefits lie in its use with large-scale process simulation. Current research indicates that these benefits are significant and the results will be reported in future work. Acknowledgment The authors wish to acknowledge Dr. Mordechai Shacham for his helpful advice. We are also grateful to the Control Data Corporation which supported this work through fellowship grants and computer facilities. Formulation of Flash Optimization Problem. Here, the set of equations required to optimize the single-stage flash of a nonideal mixture is presented. The objective will be to fmd a solution that best meets imposed purity and/or recovery targets on a key component in the feed. objective function minimize
=
Wl(a!Ykey
- Rkey)’
+ W2(Ykey - Pkey)’ (A.1)
where a! is vapor/feed fraction, ykeris vapor stream purity of the key component, Rkeyis reqwed recovery of the key, P k e y is required purity of the key, and W,and W2 are weighting fadors that allow the relative importance of the factors in (A.l) to be varied.
mass balances a(yi - x i ) x i - zi = 0 phase equilibrium In yi - In Ki - In x i = 0
+
(i = 1,..., n)
(A.2)
(i = 1, ..., n)
(A.3) A linear local model for Ki of the following form is used In
(T)
= Ai,l + Ai,’xi
(i = 1, ..., n) (A.4)
where Ai,l and Ai,’ are model parameters and f/’ is the standard-state fugacity of component i, a function of temperature only. Equation A.3 then becomes In yi - In xi - Ai,l - Ai,2xi+ In P - In f:
(i = 1, ..., n) (A.5)
The constitutive equation is n
n
i=l
i=l
cxi - c y ; = 0
(A.6)
At a given pressure, there are 2n + 2 equations and 2n Furthermore, the Jacobian matrix of eq A.2-A.6 c8n be analytically developed for the optimization algorithm. There are 2n model parameters to be updated in the outside loop of the algorithm. Nomenclature Ai,l, Ai,’ = parameters in vapor-liquid local model for component i f i = fugacity of component i , bar F = feed flow rate, mol/s H = molar enthalpy, J/mol HL = molar enthalpy of liquid stream, J/mol Hv = molar enthalpy of vapor stream, J/mol HF = molar enthalpy of feed stream, J/mol Ki = equilibrium ratio component i; Ki yJxi Li,l, Li,2, Li,3 = parameters in liquid-liquid local model for component i L = liquid flow rate, mol/s P = pressure, bar P F = pressure of feed, bar Pkey = purity specification on key component Q = external heat source, J Rkey = recovery specification on key component T = temperature, K TB = bubble point temperature, K TF = temperature of feed, K TD = dew point temperature, K V = vapor flow rate, mol/s x i = mole fraction component i in liquid phase yi = mole fraction component i in vapor phase zi = mole fraction component i in feed stream
+ 1 variables.
Greek Letters = vapor to feed fraction 4 = Rachford-Rice function # = enthalpy function CP = objective function a!
Subscripts i = component i F = feed stream L = liquid stream V = vapor stream Superscripts 0 = standard state pressure, 1 bar
Literature Cited Barrett, A.; Walsh. J. J. Compuf. Chem. Eng. 1979. 3. 397. Boston, J. F.; BrM, H. I. Comput. Chem. Eng. 1978. 2, 109. Boston, J. F. In ”Computer Appiicatkms to Chemical Englmrlng”, ACS Symposium Serbs 124; Squkes, R. G.; Reklaltls. 0.V.. Ed.; American Chemical Society: Washington, DC, 1979; p 135.
Ind. Eng. Chem. Process Des. Dev. 1904, 23,618-620
610
Broyden, C. 0. Math. Comp. 1965. 19. 577. Chlmowitz.~E. H.; Anderson, T. F.; Macchletto, S.; Stutzman, L. G. Ind. Eng. Chem. Process Des. D e v . 1983, 22, 217 (Part I). Chlmowitz, E. H.; Anderson, T. F.; Macchletto, S.; Stutzman, L. F. part 3, submmed for publlcation in Ind. Eng Chem prrrcess Des. D e v . Hen, S. P. Math. Prog. 1975, 7 7 . 263. M e , R.; Jeeves, T. A. J . Assoc. Comp. Mach. 1462, 8, 212. Hutchison, H. P.; Shewchuk, C. F. Trans. Inst. Chem. Eng. 1974. 52, 325. Leesley, M. F.; Heyen, 0. Comput. Chem..€ng.1977, 7 , 109. Martin, D. L.; Gaddy, J. L. Paper presented at 74th Annual Meetlng of the AIChE, Nov 8-12, New Orleans, LA, l 9 S l . Parker, A. L.; Hughes, R. R. Comput. Chem. €ng. 1981, 5 , 123. Powell, M. J. D. Math. Rug. 1978, 74, 224. Prausnltr. J.; Anderson, T.; Greus, E.; Eckert, C.; Hsleh, R.; O'Conneii, J. "Computer Calculations for Multicomponent Vapor-LlquM and Liquld-Liq-
.
.
uld Equlllbrla". 1st ed; Prentlce-Hall: Englewood Cliffs, NJ, 1980. Rachford. H. H.. Jr.; Rice, J. D. J . Pet. Techno/. 1952, 4 . 19. Rosenbrock. H. H. Comput. J . 1960, 3 , 175. Westerberg, A. W.; Hutchison, H. P.; Motard, R. L.; Winter, P. "Process Flowsheetlng", 1st ed; Cambrldge University Press: Cambridge, England, 1979. Westerberg, A. W. In "Foundations of Computer-Alded Chemlcal Process Deslgn"; Rocaedhgs of an Intematknal Conference held at New England College, Hennlker, NH, July 8-11, I 9 8 0 Mah, R. S. H., Selder, W. D., Ed.; p 149.
Received for review August 2, 1982 Revised manuscript received M a y 16, 1983 Accepted August 22, 1983
COMMUNICATIONS Oxygen Solubility in Brlnes The soktblity of oxygen in water and brlnes, expressed as the Henry's law constant, is reported for temperatures from 0 to 300 "C. A smoothing and Interpolation formula based upon a widely used heat capacity equation fit the solubility data well over this temperature range and may be extrapolated to somewhat higher temperatures. The soiubllity of oxygen in water and brines exposed to air saturated wRh water vapor is given for temperatures from 0 "C to the normal boiling point.
As part of Federal Bureau of Mines research on corrosion in geothermal and mineral processing environments, the solubility of oxygen in sodium chloride brines and in mixed-chloride geothermal brines was determined by Cramer (1980) over the temperature range 0 to 300 "C. The smoothing and interpolation formula used to describe the effect of temperature (T) on the Henry's law constant ( k ) was an empirical five-term power series in reciprocal temperature (l/T). While this approach is useful, the smoothing and interpolation formula is more often developed from heat capacity equations of the form AC, = Act"
+ 2APOT + 6AX"P
(1)
and the application of thermodynamic principles (Gokcen, 1975; Gokcen and Chang, 1975; Hougen et al., 1959; Smith and Van Ness, 1959). This latter approach gives the smoothing and interpolation formula a thermodynamic basis. Thermal effects are more readily evaluated from an equation of this form since heat capacity data for a wide variety of materials have been expressed in this way (Hougen et al., 1959; Smith and Van Ness, 1959). Finally, analytic expressions for the partial molal enthalpy, entropy, and Gibbs energy of solution are readily obtained from eq 1 (Gokcen, 1975; Hougen et al., 1959; Smith and Van Ness, 1959). The mathematical development of the smoothing and interpolation formula from eq 1is well established (Gokcen, 1975; Hougen et al., 1959; Smith and Van Ness, 1959), and only the final equation will be given here Ink" = (Act" - A S o , ) M o o Acto In T A@"T A X O P R RT R R R (2) where Moo and Moo are constants of integration and R is the gas constant, 1.9872 cal/K g-mol. The coefficients Aa", A@", and AX" are the differences between corre-
+--
sponding terms that describe the state of the gas molecules in the vapor and the solution. AC, represents the difference between the heat capacity of the gas in solution and that in the vapor and can be merged with known thermodynamic functions for the vapor-phase gas behavior to obtain a description of the gas properties in solution. The Henry's law constant in eq 2 is ko = 4 P / y x (3) where 4 is the fugacity coefficient for the gas in the vapor phase evaluated at the temperature and total pressure of the system, P is the partial pressure of the gas in the vapor phase, x is the mole fraction of the gas dissolved in solution, and y is the activity coefficient for the dissolved gas. Since y was not known except for the case of pure water where y N 1, k" was not measured in the study but rather a "practical" or modified Henry's law constant k = yko = ~ P / x (4) which is composed of readily determined quantities. It defines the solubility of gases in nonideal, as well as ideal, solutions (Cramer, 1980). In the limit of dilute solutions, for example, oxygen dissolved in water or in very dilute salt solutions at oxygen partial pressures of 100 atm or less, y 1 and k" and k are equal. However, the salt solutions considered in this study are not dilute. The reference state for the dissolved axygen is oxygen dissolved in water at infinite dilution. The activity coefficient for the oxygen in the salt solutions is, thus, not unity and k" and k are not equal. A generalized smoothing and interpolation formula for the modified Henry's law constant, having the same form as eq 2, can be obtained by combining eq 2 and 4
This article not subject to U.S. Copyright. Published 1984 by the American Chemical Society
