4978
Ind. Eng. Chem. Res. 2003, 42, 4978-4992
A New Correlation Method for Vapor-Liquid Equilibria and Excess Enthalpies for Nonideal Solutions Using a Genetic Algorithm. Application to Ethanol + an n-Alkane Mixtures Juan Ortega* and Fernando Espiau Laboratorio de Termodina´ mica y Fisicoquı´mica, Escuela Superior de Ingenieros Industriales, Universidad de Las Palmas de Gran Canaria, 35071 Las Palmas de Gran Canaria, Islas Canarias, Spain
This study has used a new polynomial expression with temperature-dependent coefficients as a function based on the so-called active fraction of the mixture components, z(xi), which enabled the thermodynamic vapor-liquid equilibrium quantities, the Gibbs excess energy function and corresponding activity coefficients, and the excess enthalpies all to be correlated simultaneously. Correlation is carried out by using a least-squares procedure implemented in a genetic algorithm that optimized an objective function selected for that purpose. The model is applied to a set of thermodynamic data for binary systems composed of ethanol + one of five n-alkanes (hexane to decane) collected as part of the study. Good results were obtained using the model, which yielded good estimates of the mixing properties, gE and hE, which in turn were used for comparison with the results obtained using certain more conventional models. The expression that yielded the best correlation of the excess Gibbs energy function for the binary systems considered as a function of the concentration of one of the components and temperature was (gE/RT)(T,x1) ) z1z2[(A01T + A02/T + A03) + (A21T + A22/T + A23)z12]. Introduction Experimental values for thermodynamic properties and/or correlations for those same properties hold out the possibility of a number of different uses. One very interesting use is in calculations to design separation operations at chemical plants, either directly or by means of a theoretical model. In the case of models, the concern is to formulate theories whose aim is to estimate the properties of mixtures from the properties of the pure components, and to this end, a fundamental issue is appropriate processing of experimental values using expressions that yield precise estimates of mixture properties. The first step is to prepare a simple model with a certain theoretical basis addressing the problem that is sufficiently flexible to adapt to the different situations that may arise. The following step is to process the data by inserting the model in an algorithm that will converge on an appropriate solution. Taken together, this constitutes a procedure in itself, and our purpose here is to establish such a procedure including both of the steps just mentioned above, a proposed simple theoretical model, and implementation of the model using a genetic algorithm,1-4 a type of algorithm commonly employed in engineering that has the advantage of converging on a unique solution. Genetic algorithms are not commonly used in chemical engineering to produce estimates of complex theoretical models. A literature search has disclosed two contributions putting forward quite different approaches, one by Espiau et al.5 in a regional journal and another by Rangaiah.6 The purpose was to prepare a useful method for studying systems of fluids that would provide precise estimates of values for the largest * To whom correspondence should be addressed. E-mail:
[email protected].
possible number of thermodynamic properties, in particular operations involving the phase equilibria dealt with here. The verification of the procedure proposed here was carried out on a set of actual experimental values for a series of binary systems that have been extensively studied, consisting of ethanol and one of five n-alkanes (hexane to decane). The literature has published a considerable body of values for these mixtures. Published values for the properties of interest to us here have been listed below, together with the working pressure (kPa) for the isobaric vapor-liquid equilibria (VLE) and the temperature (K) for the excess properties. (i) Ethanol + n-hexane: VLE (101.32);7,8 hE (283.15),9 (298.15),9,10 (303.15, 318.15);9,11 vE (283.15),12 (298.15),12-14 (308.15).12 (ii)Ethanol+n-heptane: VLE(24.0,53.3),15 (101.32);15-17 hE (283.15, 303.15, 318.15, 333.15, 348.15),11 (303.15, 323.15),18 (293.15),19 (298.15);20 vE (298.15),21-23 (313.15).24 (iii) Ethanol + n-octane: VLE (101.32);25,26 hE (293.15, 298.15, 303.15);27 vE (298.15).14 (iv) Ethanol + n-nonane: VLE (101.32);26 hE (293.15),29 (303.15, 318.15);11 vE (298.15).22 (v) Ethanol + n-decane: VLE (101.32);26,28 hE (293.15),29 (298.15).10,29 This data set was chosen for application of the procedure precisely because of the interest these mixtures hold out as tests in the studies on solutions. The mixtures containing ethanol have been used as tests despite the scarcity of VLE values at low ethanol concentrations in mixtures with high molecular weight hydrocarbons. VLE measurements were taken at three different pressures, and vE and hE measurements were made at two different temperatures to supplement the existing experimental values available, with a view to broadening the scope of this study as it relates to the values and processing of thermodynamic quantities.
10.1021/ie030327j CCC: $25.00 © 2003 American Chemical Society Published on Web 08/28/2003
Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4979
Thermodynamic-Mathematical Approach Under practical and theoretical points of view, the Gibbs energy function and mixing enthalpy are two of the most important excess thermodynamic properties for analyzing the behavior of fluid mixtures. VLE studies carried out at constant pressure have been accorded considerable attention in the solution thermodynamic existing databases for binary and multicomponent systems because isobaric VLE values offer very good precision and provide an excellent source for study when the excess enthalpies are also known. Other useful thermodynamic expressions include equations relating the Gibbs function and the activity coefficients. Processing of VLE values needs to take the dependence of the excess properties on temperature and pressure into account. Nevertheless, this aspect has not received all of the attention it deserves, even though bringing this consideration to bear yields more exact thermodynamic relations that appreciably reduce the amount of experimentation needed. The following are important basic relations useful in equilibrium thermodynamics.
[
]
∂(gE/RT) -hE )T RT ∂T
p,x,...
[
]
cEp ∂(hE/R) ) R ∂T
p,x,...
(2)
Substituting eq 2 in the second of the two expressions in eq 1 and integrating yield E
c h b ) aT + T 2 + T 3 + I1 R 2 3
(3)
where I1 is an integration constant. Repeating this same substitution of eq 2 in the first of the two expressions in eq 1 yields
I1 gE b c ) -a ln T - T - T 2 + + I2 RT 2 6 T
cEp ) a + bT R
cEp )a R
hE b ) aT + T 2 + I1 R 2 I1 gE b ) -a ln T - T + + I2 (5) RT 2 T
hE ) aT + I1 R I1 gE ) -a ln T + + I2 (6) RT T
(1)
For the expressions in eq 1 to be of use in the model approach proposed here, certain further considerations have to be made to make them easy to integrate. The dependence of the thermal capacity on temperature can be expressed through a linear function. So, in a binary mixture at a given concentration, may be properly well represented over a broad temperature range by the following simple equation.30
cEp ) a + bT + cT 2 R
confirmed because, of the five constants needed for gE/ RT, only four are necessary for the enthalpies hE/R and three for the excess thermal capacities cEp /R. Another negative aspect of this procedure is that other quantities are needed to be able to calculate the values of the function gE/RT from equilibrium data, for instance, the activity coefficients γi, whose calculus could be inexact despite the apparent quality of the experimental isobaric VLE values obtained. A series of additional simplifications of eqs 2-4 that result in cEp /R being a linear function of T for c ) 0 or even independent of T for b ) c ) 0 are also possible. After these two special assumptions, the previous equations respectively become
(4)
where I2 is another constant arising from the integration process. In principle, eqs 2-4 allow all of the constants to be determined when experimental values for the thermodynamic quantities in these equations are available. One route would be to use the VLE values at various temperatures to furnish gE/RT values as a function of T, giving rise to various situations in which eq 4 can be used to get the constants a, b, c, I1, and I2. Thus, eqs 2 and 3 are also defined for cEp and hE as functions of T. One shortcoming of this procedure is that such primary quantities as hE and cEp would be obtained from a firstorder differentiation and a second-order differentiation of the excess Gibbs function, respectively, losing 1 order of magnitude in each differential process. This is
For short temperature intervals, cEp /R may, in fact, be independent of T in practice, which is equivalent to taking hE/R to be linear on T, as in the second expression in eq 6. The parameters may be calculated using the set of eqs 5 and 6 from appropriate values of hE/R and gE/RT as functions of T, as already explained. The most extreme case occurs when a ) b ) c ) 0 in eq 2 or when a ) 0 in the first expression in eq 6, likewise causing the enthalpies to be temperatureindependent (a practice also used by researchers in the absence of real values to establish that dependence). In this case the second and third expressions in eq 6 are respectively reduced to
hE ) I1 R
I1 gE ) + I2 RT T
(7)
The cases considered in eqs 5 and 6 may be merged into eq 7, bearing in mind that certain quantities in those equations, such as the quantity related to ln T, may be expressed in another form, as will be discussed below. The Gibbs function is also dependent on other variables, such as the pressure and composition, xi, of the solution. In this last case, the Gibbs function gE/RT for a binary system may be well represented at a given temperature by a polynomial equation derived from a simplified form of Wohl’s equation,32 previously employed by Ortega and Toledo,31 to correlate the excess quantities and the composition. That equation is
gE/RT ) z1z2(A0 + A1z1 + A2z21 + ...)
(8)
where z2 ) 1 - z1 and z1 ) x1/(x1 + kgx2) are the socalled active fractions of each of the mixture components. The most direct way to integrate using eq 8, in view of the thermodynamic approach discussed above, is to express each coefficient Ai as a function of temperature, T, considering the generalized form in eq 4, which for this case can be written as
4980 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003
Ai ) Ai1T 2 + Ai2T + Ai3 ln T +
Ai4 + Ai5 T
(9)
and then rewrite eq 8 including the dependence of the Gibbs function on T and xi
gE
(T,x1) )
RT
m
∑ i)0
z 1 z2
(
Ai1T 2 + Ai2T + Ai3 ln T +
Ai4 T
)
+ Ai5 zi1 ) m
z1(1 - z1)
Aizi1 ∑ i)0
(10)
This general expression, in which the terms Aij are constants, provides a real approximation of the mixing Gibbs function for any range of temperatures containing the equilibrium values and excess enthalpies. The procedure may be repeated by making use of simpler expressions for the thermal capacity on temperature, thereby yielding a model with fewer parameters. The generalized form for the activity coefficients furnished by model (10) is
{
ln γj ) z1z2
[xj + (-1)j-1(1 - 2z1)] m
Aizi1 + ∑ i)0
xj (-1)
z1z2 m
j-1
xj
iAizi-1 ∑ 1 i)1
}
(11)
The first relation from eq 1 can be used to reproduce the excess enthalpies, and finding the derivative of eq 11 yields
hE
-
RT
[
( ) ( )∑ ] [∑( ) ( ) ( ) ∑ ]
) -T (1 - 2z1)
Tz1(1 - z1)
m
i)0
dA0 dT
dz1
dk
dk
dT i)0 dz1 dk
zi1 +
dk
m
Aizi1 -
dT
m
iAizi-1 (12) 1
i)0
In the model considered here, the value of k has a greater influence on curve shape than on curve size and therefore k could be considered temperature-independent. Accordingly, truncating eq 9 after the third term yields the simplified form
[( ) ( ) ( ) ]
dA0 dA1 dA2 2 hE ) z (1 z ) + z + z 1 1 1 dT dT dT 1 RT 2
(13)
Substituting the derivatives of Ai using the definition of coefficients established in eq 9:
h
E
RT
2
) z1(1 - z1)
(
2Ai1T 2 + Ai2T + Ai3 ∑ i)0
)
Ai4 T
zi1 (14)
yielding global coefficients zi1 in eq 14 identical with those obtained by direct correlation of enthalpies, hE/ RT vs z1, using an equation similar to eq 8 but adapted
to the enthalpies:
hE ) z1(1 - z1)(A10 + A11z1 + A12 z21 + ...) RT
(15)
where the terms A1i are equivalent to the coefficients zi1 in eq 14. One of the procedures used here to calculate the model parameters in eq 10 was, first, to use eq 14 to estimate the values of Aij (j ) 1-4) using the enthalpy values at different temperatures and then to calculate the values of Ai5 and k by correlating eq 10 only with the set of values for gE/RT. This two-step procedure for estimating these parameters is explained in more detail below and will be compared with the results obtained employing a least-squares procedure using a genetic algorithm to optimize the resulting objective function. Processing of Equilibrium Data A model that affords good estimates is required to achieve a good representation of thermodynamic data under conditions the same as or different from those of direct experimentation. Construction of such a model involves the following steps: (i) Proposal of a model furnishing a reasonable explanation of the phenomenon observed. (ii) Compilation of as many experimental and precise measurements of the phenomenon as possible. (iii) Application and validation of the model proposed in step i. This paper proposes use of eq 10 as a model for the Gibbs function using a least-squares procedure.33 This procedure is optimal when the errors inherent to data collection are normally distributed and independent and have the same variance. In using the model (10) to process the VLE data, it needs to be kept in mind that T and p values were taken directly, whereas xi values were obtained indirectly. These values were used to calculate γi, and then the values of the Gibbs function gE/RT. The hE values were measured directly and independently by calorimetry and were related to the gE/RT values by means of eq 1. In short, the data sets that define the isobaric equilibrium states of a binary system are, on the one hand, {(Tj, x1j, lnγ1j, ln γ2j); j ) 1, ..., n} and, on the other hand, {(x1i, (hE/RT)i); i ) 1, ..., m} where the γi values are obtained at a temperature Tj for a given concentration of component 1, x1j, while the nondimensional quantities (hE/RT)i are measured at a concentration of reference component 1, x1i, and at a temperature, T, which is generally different from that of equilibrium data. Use of a least-squares approach involves an objective function (OF) that will reflect quadratic minimization of the discrepancies between the experimental values and those estimated by the model for a given equilibrium state. The OF considered was q
OF )
[
hE
(T,x1i) ∑ i)1 RT
( )] hE
2
n
+
{[ln γ1(Tj,x1j) ∑ j)1
RT i ln γ1j] + [ln γ2(Tj,x1j) - ln γ2j]2} (16) 2
The variables of this function are the coefficients of the model for the Gibbs function, with the optimum values being the values that minimize the OF. For this case, it may not be appropriate to use conventional numerical methods for optimizing the OF
Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003 4981 Table 1. Physical Properties of the Pure Components Used in This Study compound
purity
ethanol
>99%
Tb,exptlo/K
Tb,lit.o/K
351.46
351.44a
n-hexane
>99%
341.93
n-heptane
>99%
371.18
n-octane
>99%
398.65
351.47c 341.89a 341.88c 371.57a 371.60c 398.82a,c
n-nonane
>99%
423.66
423.97a,c
n-decane
>99%
446.90
447.30a,c
a
Fexptl/(kg ‚m-3)
Flit./(kg ‚m-3)
T/K
nD,exptl
nD,lit.
785.14 767.51 654.59 636.39 679.27 661.47 698.39 682.09 713.85 698.06 726.19 711.14
784.93a
298.15 318.15 298.15 318.15 298.15 318.15 298.15 318.15 298.15 318.15 298.15 318.15
1.3595 1.3553 1.3724 1.3614 1.3853 1.3748 1.3951 1.3856 1.4030 1.3939 1.4098 1.4010
1.35941a 1.3551b 1.37226a 1.3615b 1.38511a 1.3750b 1.39505a 1.3855b 1.40311a 1.3939b 1.40967a 1.4008b
767.62b 654.84a 636.67b 679.46a 662.32b 698.62a 682.09b 713.75a 698.06b 726.35a 711.43b
Reference 34. b Reference 35. c Reference 36.
for the following reasons: (i) the OF has a large number of parameters; (ii) analytical expressions of three different functions are employed in the OF, namely, (hE/ RT)i as defined by eq 14 and ln γ1j and ln γ2j as defined by eq 11; (iii) the OF may be nonconvex and as a result may have several local optimums; and (iv) the shape of the OF may involve complex systems of nonlinear equations, making estimation of the parameters difficult. On the basis of all of these reasons, it follows that a genetic algorithm can be used to optimize the OF defined in eq 16; see the appendix. Another advantage of the procedure for estimating the coefficients of the model for the Gibbs function, eq 10, proposed here is that it can be performed in a single step instead of the two steps needed to complete the procedure described by eqs 8-14, in which four Aij coefficients (j ) 1-4) and k are estimated initially based solely on hE/RT vs x1, which are considered as fixed values in eq 10, and the values of Ai5 and a new value of k are then estimated by means of a second correlation to optimize the fit of eq 10 to the gE/RT values versus composition and temperature. In contrast, optimization of the OF using the genetic algorithm proposed here yields valid overall estimates of the three quantities mentioned in item ii above. The OF as formulated in eq 16 does not take into account the values of the nondimensional Gibbs function because that function’s values are not obtained directly but rather are calculated from the activity coefficients, and hence inclusion of these values in eq 16 would not contribute any statistical independent information. Experimental Section Components. Ethanol and n-alkanes employed in this work were of the highest commercial grade (Fluka AG). All components were degassed by ultrasounds and dried on a molecular sieve (0.3 nm, Fluka) prior to use. The component quality was tested by gas chromatography using a HP-6890 system with a flame ionization detector. The degree of purity obtained is shown in Table 1 and in all cases was consistent with the manufacturer’s specifications. For the purposes of comparison, the component quality was also checked by measuring such properties as the normal boiling point Tb,io, density F, and refractive index nD. These values also appear in Table 1 and in all cases were in good agreement with the literature values. Equipment and Procedure. The experimental VLE still, the auxiliary equipment, and the experimental procedure used were described previously by Ortega et al.37
Table 2. Coefficients Ai and k, and Standard Deviations s(YE), Obtained Using Equation 8 To Correlate the Excess Properties of Binary Systems Ethanol (1) + n-Alkane (2) YE ) 109vE/(m3‚mol-1) mixtures ethanol (1) + hexane (2) heptane (2) octane (2) nonate (2) decane (2) ethanol (1) + hexane (2) heptane (2) octane (2) nonate (2) decane (2)
kv
A1
A2
s(YE)
-3134 -4007 -4275 -4912 -4628
3234 4139 4357 5201 4414
9 7 5 6 8
5234 7272 5413 6598 8206
8 10 9 8 8
A0 T ) 298.15 K
0.733 0.892 1.031 0.985 1.282
2388 2844 3064 3185 3249
T ) 318.15 K 0.370 0.603 0.845 0.945 0.974
2535 3948 3851 4425 5194
-3664 -6538 -5320 -6790 -8875
YE ) hE/(J ‚mol-1) mixtures ethanol (1) + hexane (2) heptane (2) octane (2) nonane (2) decane (2) ethanol (1) + hexane (2) heptane (2) octane (2) nonane (2) decane (2)
kh
A1
A2
s(YE)
-4798.8 -6056.9 -5945.1 -9136.0 -8169.8
4377.3 6146.3 6211.5 7160.9 6822.8
8.0 9.9 7.0 4.1 10.1
11123.2 12533.2 14299.2 13758.0 15721.6
22.9 16.7 18.2 22.7 26.3
A0 T ) 298.15 K
0.682 0.607 0.604 0.907 0.941
3652.9 3998.5 4133.5 5652.6 5287.1
T ) 318.15 K 0.526 0.576 0.601 0.480 0.571
6185.2 6994.1 7706.9 7255.3 8231.7
-11241.2 -13205.4 -15069.5 -13499.5 -16159.2
The concentration values for the binaries ethanol (1) + n-alkane (2) at equilibrium were determined from a standard curve of density on concentration (see Table 2§ in the Supporting Information) previously fit by a second-degree equation for each system at temperatures of 298.15 and 318.15 K using samples of known composition. Densities were measured using an Anton Paar model DMA60/602 digital densimeter. The relationship F ) F(x1) was validated by verifying the quality of the results for vE vs x1; using this method, the concentrations of the mixture components at equilibrium were back-calculated from the density values of the condensed vapor and liquid phases, yielding estimates better than (0.002 ethanol mole fraction units. The excess enthalpies were measured isothermally at 298.15 and 318.15 K using a Calvet (MS80D) calorimeter calibrated electrically at regular intervals. The
4982 Ind. Eng. Chem. Res., Vol. 42, No. 20, 2003
Figure 1. Excess molar volumes (b) and correlation curves for binary mixtures C2H5OH (1) + n-CvH2v+2 (2) at (a) 298.15 K and (b) 318.15 K. Labels indicate the v values. Inset figures plot the variation of equimolar volumes with v and a comparison with literature values: (b) experimental; (4) Ruel;14 (g) Marsh and Burfitt;12 (
