conditions on steady-state activity, the present calculations demonstrate that physical factors may sometimes be the cause instead. Acknowledgment
The machine computations were performed a t the Massachusetts Institute of Technology Computation Center. The National Science Foundation provided financial support to George W. Roberts in the form of a fellowship during this study. Nomenclature
D
=
E
= = = =
k k’
KA, KB Ki K
=
=
effective diffusivity of reactant in catalyst, sq. cm./sec. defined by Equation 2 reaction rat’: constant, g. moles/(cc.) (sec.) reaction rate constant, g. moles/(cc.) (sec.)(atm.)2 adsorption coefficients for the reactants, (atm.)-l adsorption c,oefficient for any product, (atm.)-l an over-all adsorption coefficient, (atm.) K = [ K A -. DA Z (Kiui/Di)]/w i
L P
R TA
= = = =
thickness of catalyst mass, cm. partial pressure of reactant or product, atm. gas law constant, (atm.)(cc.)/(g. moles) (OK.) reaction rate of component A, g. moles/(cc.) (sec.)
t T
=
X
=
time, sec.
= temperature,
OK.
distance in from exposed catalyst face, cm.
GREEKLETTERS 11 = effectiveness factor = void fraction in catalyst e = stoichiometric coefficient, negative for a reactant V @L = dimensionless modulus, = (L2/D~Ca,,)(observed reaction rate/gross catalyst volume) = dimensionless modulus, = L(k‘RT/ 1 K I DA)”~ +.w w = 1 FKi[pi,a ( p ~ , s ~ i D ~ / D i ) l
+
+
I
SUBSCRIPTS A, B = reactants z = any product species = initial conditions inside catalyst pellet 0 = conditions at outside (exposed) catalyst surface S literature Cited
(1) Hartman, J. S., S. M. thesis, Department of Chemical Engi-
neerinc. Massachusetts Institute of Technolow. 1965. G. W., Satterfield, C. N., IND.Kid. CHEM.FUNDAMENTALS 4, 288 (1965). ( 3 ) Zbid., 5 , 317 (1966). RECEIVED for review February 18, 1966 ACCEPTED August 18; 1966 ( 2 ) Rob&,
TREATMIENT OF ISOBARIC BINARY VAPOR=1,l Q U ID EQ U I LI BR I U M DATA LUH C. TAO Department of Chemical Engineering, University of Nebraska, Lincoln, N e b ,
A new approach i s proposed to treat isobaric, binary vapor-liquid equilibrium data so that correlations will be consistent thermodynamically and will also include other physical constraints. NDUSTRIAL
distillation always operates at a constant pressure.
I Therefore, isobaric vapor-liquid equilibrium data are important for design and operation of distillation equipment. The simplest system is a binary, the treatment of which is discussed here. The most common treatment of isothermal vapor-liquid equilibrium data is to check thermodynamic consistency and to find a numerical correlation, both from log(yIly2) as a function of composition (8). This treatment is successful because the term involving volume of mixing in the isothermal Gibbs-Duhem equation ( 5 ) is usually negligible. When this same method is used to treat isobaric data, the thermodynamic consistency test and correlation may be complicated by the existence of the term involving heat of mixing in the GibbsDuhem equation. Only when this term is negligible will the method for isothermal data succeed. I n view of the scarcity of heat of mixing data, isobaric data have not been satisfactorily treated. C hao (3) proposed the only available correlation based on a modified Redlich-Kister equation, and its failure to observe some thermodynamic constraints was observed by Naphtali (7). A desirable correlation regardless of isothermal or isobaric data should satisfy thermodynamic consistency criteria, which include Q and log(yI/yz) as corresponding integrals and
differentials with respect to x1 (70). The usual treatment implies only the over-all integration criterion. For isobaric data, the heat of mixing term in d Q / d x l makes the Van Laar, Margules, and Redlich-Kister equations impossible to fit the log(yl/yz) profile since these equations are restricted to a form where the over-all integral must be zero. However, if a satisfactory treatment can be found, theoretically sound correlations will provide constants for extrapolation of homolog components and multicomponent systems, and heat of mixing of solutions may be estimated for process design and operation. Theory and Method
The dimensionless excess free energy of a binary solution is defined in Equation 1 (70) and its derivative in Equation 2 is the isobaric Gibbs-Duhem equation ( 5 ) .
Q
=
GE/2.303 R T =
dQ Y1 - = log dx 1 7 2
XI
+ x z log y2
(1)
WM dT 2.303 R T 2 dxl
(2)
log yi
- Z where Z
=
T h e Q function has three constraints: unimodality (having a single hump), Q ( x l = 0) = 0, and Q(xl = 1.0) = 0. T o satisfy the last two properties, any Q function containing a VOL. 6
NO. 1
FEBRUARY 1967
83
factor of x l x 2 is suitable. However, the analytical form of the Q function to observe unimodality is not obvious. From extensive work on vapor-liquid equilibrium data, the most common types of suitable functions appear to be the Van Laar, Margules, and Redlich-Kister equations ( 4 ) , in which the Margules equation corresponds to a two-constant Redlich-Kister equation. These equations are summarized in Table I. Since Equation 1 does not contain the heat of mixing term, apparently, a natural start for treatment of data would be to find a suitable Q function observing all three properties mentioned. This implies computation of the most probable values of constants of equations in Table I by least-squares methods or some suitable optimum search methods. Once a Q function is determined, its corresponding differential derivative can be found from Table I. Therefore, Z(x1) can be computed from Equation 2 by using calculated dQ/dxi from Table I and experimental values of log(yi/yz). T h e heat of mixing is usually a unimodal function with constraints of HM(xl = 0) = 0 and HM(xl = 1) = 0 exactly like those of the Q function. Therefore, for any binary without an azeotrope, the Z function probably can be fitted by similar functions in Table I. For an azeotropic binary, dT/dxl has different signs for solutions on two sides of the azeotrope. Therefore, the Z function would be unimodal only within 0 _< X I _< xla or x, _< X I _< 1. Also, as (dT/dxi),=O, the Z function must satisfy Z(x1,) = 0. These additional constraints for the Z function can be conveniently met by replacing ~ 1 x 2 by x1x2(x1 ~ 1 ~ )Three . probable forms of 2 functions are also tabulated in Table I. If the calculated Z values have random signs along x i , they signify the propagated experimental errors of data and further correlation of Z values is meaningless. Only when Z values are consistent does correlation of possible Z functions become significant and consequently HM may be estimated according to the definition of Z. Once a Q function and a Z function are computed, analytical functions of activity coefficients can be obtained by solving Equations 1 and 2 for log 71 and log 7 2 . This gives Equations 3 and 4. An examination of Equations 3 and 4 shows that they follow all necessary thermodynamic constraints represented by Equations 5 and 2 and other constraints of Q and Z functions.
-
Qc
-
QO
=
~~log(~1/~z)dxi
++
-
- XZ)~]
Equation QV. A = 0.245341 B = 0.307307 S.S. = 0.000398 Equation QM. A = 0.267030 C = 0.008543 S.S. = 0.000397 Equation QR. B = 0.261891 C = 0.008318 D = 0.029894 S.S. = 0.000363 Among the 40 data points for each equation, only two points have positive 2 values and the rest are all negative. If these Z values are fitted into functions in Table I, the coefficients of Z functions with the least sum of squares in each of the Q functions are :
QV:
ZV
U = -0.732218 V S.S. = 0.032126
QM:
ZM
U
= -0.182377 V = S.S. = 0.053681
QR:
ZV
U
= -0.291474 V = -0.117592 S.S. = 0.033353
ZR: z
04
l&EC FUNDAMENTALS
=
XlXZ(X1
- xio) term.
-0.086001 0.154898
log yi = 0.023169 xz2/(O.245341 xi
(5)
log
72
+ 0.307307 0.062971 ~ 1 ~ ~ ~ / ( 0 . 7 3 2xi2 + 1 8 0.086001 = 0.018497 xiz/(O.245341 + 0.307307 + 0.062971 ~ 1 ~ ~ ~ / ( 0 . 7 3 2 2 + 1 80.086001 XZ)'
XI
~
2
)
XI
Commonly Used Functions Derivatives of Q Function dQ/dxi = AB(Bxz2 - AxlZ)/(Axi B~2)2 dQ/'dxi = B(xz - xi) C(6xixz - 1) dQ/dxi = B(xz - xi) C(6xix2 - 1) D ( x i 2-Functions
++
- Xl,)[U
~ 2 )
~
XZ)
The other two sets give relatively close values of 71 and 72 a t infinite dilution. All three sets fit data points well, except a few points in the terminal regions.
+ +
ZV: z = uVxlxz(x1 - X l r r ) / ( U X 1 Vxz) ZM: Z = XIXZ(X~ - X I , ) [ U V(XI - X Z ) ]
For nonazeotropic systems, delete the ( x i
=
(4)
method or suitable optimum-seeking methods.
+
The vapor-liquid equilibrium data of the n-pentane-benzene system at 760 mm. of I l g (6) are used to illustrate the proposed treatment. Coefficients in Q and Z functions of the RedlichKister type were computed by the conventional least-squares method and of the Van Laar type by the method of Smith and Tao (9). The latter method for nonlinear equations would locate only an approximate solution in the second cycle but is less time-consuming than many available optimumseeking methods that locate the exact solution. The Q functions (Table I) computed have the following numerical values and corresponding sum of squares (S.S.)
Therefore, the three sets of activity coefficient functions are explicitly obtained. The simplified set, QV-ZV, for activity coefficients is shown below:
1. Compute Q functions (Table I) by the least-squares
Table 1.
Example
(3)
T h e proposed procedure for treating data thus involves the following steps :
Q Functions Q V : Q = A B x ~ x z / ( A x ~Bxz) Q M : Q = x i ~ z [ B C(xi xz)l QR: Q = X I X Z [ B C ( X I - XZ) 4-D ( x 1
2. Calculate dQ/dxl (Table I) by using Q functions developed by step l . 3. Calculate Z (Equation 2) by using dQ/dxl in step 2 and experimental activity coefficients, 4. Calculate the Zfunction (Table I) by methods mentioned in step 1. 5. Substitute the computed coefficients into Equations 3 and 4 to obtain analytical functions of activity coefficients. 6. Compute the estimated heat of mixing from the definition of Z (Equation 2).
+ V(x1 - + W(x1 - x2)7 x2)
+
+
- ~ 2 ) ( 8xixz - 1)
0
I
.04
zI -104
data equption QV 7
0
.2
0 Figure 1.
-.4
Q vs.
0
.4
,
XI
x,
,6
I
.8
1.0
for n-pentane (1 )-benzene (2)
data 0 equation QV- -equation QV-ZV .
,
.
.2 Figure 2.
.
*4
0 0.
.
XI
.6
.8
1.0
Log yI/y2 VS. x1
I equation QV--equation QV-ZV-
0 Figure 3.
.2
.4
.6
-8
1.0
Liquid phase activity coefficient
The first set, QV-217,is used to present the computation results graphically. Figures 1 and 2 show conformity toward thermodynamic consistency constraints. Figure 3 indicates that activity coefficient equations satisfactorily represent data points except a t ends where experimental errors are usually rather large (70). Finally, Figure 4 plots estimated heats of mixing a t various bubble temperatures of the solution. They are about 2.5 times larger than the heat of mixing of benzeneheptane solutions measured experimentally at 20’ C. ( 2 ) . Discussion
This proposed approach differs from the conventional treatment of vapor-liquid equilibrium data in starting a
Figure 4. perature
Estimated heat of mixing at bubble tem-
correlation of Q rather than of log(yl/y*). This makes it possible to compute Z for further correlation. Z was correlated without a temperature profile, which was used only in estimating H‘. Since the curve fitting of Q and Z is made sequentially, the best representation cannot be determined from the sum of squares of Q functions alone. If the fitting of the log(yl/yz) profile is the criterion, QV-ZV appears to be the best of the three in the example. However, if fitting y or log y is the criterion, further computation is necessary to determine the best equation. Besides the three types tabulated, it is possible to include others, such as a four-constant Redlich-Kister equation, However, the basic approach is the same. For isothermal vapor-liquid equilibrium data, Z in EquaRT)(dP/dxl). If this term is not tion 2 equals -(V‘/2.303 negligible and the binary consists of components with a very large difference in vapor pressure, this proposed method definitely will apply also. Experimental determination of heat of mixing at the bubble point is difficult, since the volatile component in the pure state will be a vapor, not a liquid, and a small relative error in the latent heat of condensation will cause a very large error in H’. Probably, extrapolation of low temperature calorimetric data at constant liquid composition will be desirable. T h e heat of mixing obtained from vapor-liquid equilibrium data is of interest in that experimental work on phase equilibria is less tedious and less costly than precision-type calorimetric work. T o ensure the validity of the computed heat of mixing, it is desirable to correct liquid phase activity coefficients for nonideal gas behavior of vapors (7). This was illustrated by Tao (70)and would make the estimated heat of mixing smaller than that of Figure 4. This shows the importance of vapor phase imperfection in the estimation of heat of mixing. Since any correction of liquid phase activity coefficients requires a set of equations of state a priori, no attempt was made at such corrections in the illustration given here. Based on the general concept of statistics, propagated errors from unavoidable experimental errors will be less if more experimental runs are made to give more data points in a computation. An increased number of data points definitely enhances the reliability of this method and the example was chosen for the large number of data points available for this binary. Finally, the present approach with suitable equations should be useful for a consistency test in general, as this approach would cover both integral and differential tests for local and over-all composition ranges. The method is better than the usual graphic approach, in that graphic differentiation inVOL. 6
NO. 1
FEBRUARY 1967
85
volves a very large operation error. For common isothermal data, consistency means that the same Q function fits both Q and the log(yI/y*) profile and for isobaric data consistency means an additional agreement between estimated and experimental H M values. Of course, the algebraic functions used imply less flexibility than the graphic means to compensate for the cited advantages.
Z y
= =
a designated term in Gibbs-Duhem equation liquid phase activity coefficient
SUBSCRIPTS 1 , 2 = component identity = azeotrope c, b = arbitrary liquid composition a
literature Cited Ac knowledgment
Sidney G. Martin and Adrian C. Snyder assisted in the computer work. Nomenclature
GE = excess free energy HJf = heat of mixing of liquid solution
P
total vapor pressure of a liquid solution dimensionless excess free energy gas constant = absolute temperature V+f= volume of mixing of liquid solution = liquid solution composition, mole fraction x
Q R T
= = =
(1) Black, C., Ind. Eng. Chem. 50, 391 (1958). (2) Brown, C. P., Mathieson, A. R., Thynne, J. C. J., J. Chem. SOC.1959, p. 4141.
(3) Chao, K. C., Ind. Eng. Chem. 51, 93 (1959). (4) Hougen, 0. A., Watson, K. M., Ragatz, R. A., “Chemical Process Principles,” Part 11, 2nd ed., pp. 904-7, Wiley, New York, 1947. (5) Ibl, N. V., Dodge, B. F., Chem. Eng. Sci. 2,120 (1953). (6) Myer, H. S., Znd. Eng. Chem. 47, 2215 (1955). ( 7 ) Naphtali, L. M., Ibid., 51, 1318 (1959). (8) Redlich, O., Kister, A. T., Zbid., 40, 345 (1948). ( 9 ) Smith, L. C., Tao, L. C.,Chem. Eng. 70,193 (1963). (10) Tao, L. C., Znd. Eng. Chem. 56, 36 (1964). RECEIVED for review June 13, 1966 ACCEPTED September 16, 1966
DEGREES OF FREEDOM FOR UNSTEADYSTATE DISTILLATION PROCESSES G . M I C H A E L HOWARD
University of Connecticut, Storrs, Conn. The usual method of determining the number of degrees of freedom for a distillation column is shown to b e incorrect when applied to unsteady-state processes. The error is in neglecting the various holdups in the system, which provide extra degrees of freedom. The additional degrees of freedom exactly equal the number of holdups which may b e specified. The design of columns for steady-state operation is not affected by these additional degrees of freedom, but they should b e considered when developing mathematical models or analyzing possible control systems in which the number of independent controllers is equal to the degrees of freedom.
THE problem
of deciding how many variables must be specified to ensure a unique solution of steady-state distillation tower design problems is a familiar one. T h e basic concept is that of the phase rule. I n a system with a large number of physical variables and a set of relationships among them, how many of the variables may be freely chosen? T h e rest will be uniquely determined by the various relationships. If too few variables are specified, the problem will have an infinite number of solutions. If too many variables are specified, the problem will have no solution unless the extra variables happen to fit the exact solution. The problem has been treated by several authors (7-6) with Smith’s summary of Kwauk’s (3) work the most lucid explanation of the general theory. Most recently, Murrill (4)applied the concepts to the selection of instrumentation schemes for a distillation tower. The author recently attempted to use this approach in determining the conditions which could be specified in solving the equations describing the unsteady-state behavior of distillation columns and found that the previously published tables of degrees of freedom do not apply. T h e holdup in the systems has not been considered. In making steady-state calculations, the amount of holdup is immaterial; but in 86
I&EC FUNDAMENTALS
transient operation, the response time is a direct function of holdup. Thus, some provision must be made for including holdup in the analysis. Determination of Degrees of Freedom with Holdup
Expressions for the degrees of freedom of the various parts of a distillation column may be derived by noting that the holdup is analogous to a stream. I t has the same variables, except that a quantity is involved instead of a flow rate. T h e development here follows the notation and procedure of Smith. A single phase stream has C 2 variables which completely specify the stream. These are C - 1 concentration variables, a rate (or quantity) and two other intensive variables (temperature and pressure). There are no restrictions among these variables other than the condition that all values be positive. An ideal plate is shown in Figure 1. I t involves five material streams (vapor holdup is neglected) and one heat stream to give a total number of variables
+
N , = 5(C
+ 2) + 1 = 5 C + 11
The relationships among the variables depend on the way in which the plate is defined, but the total number of relation-
