Ind. Eng. Chem. Res. 2007, 46, 5437-5445
5437
Modeling of Simultaneous Chemical and Phase Equilibria in Esterification of Acetic Acid with Ethanol in High-Pressure Carbon Dioxide Andrzej Wyczesany* Institute of Organic Chemistry and Technology, UniVersity of Technology, 31-155 Krako´ w, Poland
This paper presents modeling of simultaneous chemical and phase equilibria (CPE) in esterification of acetic acid with ethanol in carbon dioxide at high pressure. Nine models based on the cubic equations of state (EOSs) have been applied. Four of them use classical mixing rules with or without the binary interaction coefficients kij, one model treats this coefficient as a function of temperature and composition, and the last four methods apply the mixing rules combining EOS and excess Gibbs free energy models. Two procedures had to be applied to obtain the agreement between the calculated and experimental results. First, the thermodynamic data referring to the acetic acid Gibbs free energies of formation were corrected to describe precisely the literature experimental chemical equilibrium constants of the considered esterification reaction in the gas phase. Second, since the literature VLE data of the system CO2-ethyl acetate turned out to be inaccurate, the binary interaction coefficients of this mixture were fitted from the ternary VLE data of the system CO2-ethanol-ethyl acetate. 1. Introduction The esterification of acetic acid with ethanol in CO2 at high pressure was studied experimentally by Blanchard and Brennecke.1 They stated that equilibrium conversion to ethyl acetate and water at 60 °C and 58.6 bar reached 72% and was 9% greater than the conversion obtained in the liquid phase without CO2. In many reactions, conversions are limited by chemical equilibrium. In some of these reactions, application of highpressure CO2 as a solvent can increase the conversion. The experimental investigations of such processes should be preceded by suitable calculations that would indicate the best ranges of the reaction parameters. However, calculations will be helpful only in the cases for which the reliable thermodynamic models are available. The present state of knowledge enables us to model precisely the chemical equilibrium in the gas phase only. The equilibrium composition calculated for an ideal model usually does not deviate too much from the experimental values even at elevated pressures and for strongly nonideal systems. These deviations can be described quite well by the equation of state (EOS) methods. The accuracy of the chemical equilibria predictions in the liquid phase is considerably lower. This question is also much less examined. The problem of chemical equilibrium calculations in diisopropyl ether synthesis from water and propene in the liquid phase at elevated pressure was investigated by Wyczesany.2 Comparison of the calculated results with the experimental data revealed the significant influence of the model applied. For the investigated four-component mixture, concentrations calculated by some methods differed by more than 10 mol % from the measured values. The prediction of chemical and phase equilibria (CPE) is still a less examined problem. In such cases, the model used for calculations affects the results much more than in the case of independent vapor-liquid equilibria (VLE) or chemical equilibria in one phase. This problem illustrated by example of esterification processes at atmospheric pressure was studied by * To whom correspondence should be addressed. E-mail: awyczes@ chemia.pk.edu.pl. Tel.: 48-12-6282764. Fax: 48-12-6282037.
Wyczesany.3 As it turned out, the equilibrium compositions calculated with the use of methods usually considered as reliable for VLE differed significantly. The ASOG method (neglecting association of acetic acid in the gas phase) predicted equilibrium in the vapor phase only, while the original UNIFAC method predicted that the equilibrium product formed the liquid phase below the bubble point. The considered esterification process performed in carbon dioxide at high pressure was not investigated in the respect of modeling. This reaction was used only as an example for testing of the computational algorithm presented by Burgos-Solo´rzano et al.4 The numerical values of some parameters kij in the three used models were chosen arbitrarily to show that the program can calculate even the liquid-liquid-vapor (LLV) + chemical equilibria problems. The authors stated that none of the three models could describe the experiment presented by Blanchard and Brennecke.1 The aim of this paper is to find a proper model for this process. Nine models based on the cubic EOS have been chosen. Four of them use the classical mixing rule with kij ) 0.0 or kij fitted from the experimental VLE data. In one model, kij ) f(T, xi). The last four methods use the mixing rule from combining EOSs and excess free energy models. (1) Soave-Redlich-Kwong EOS with the coefficient kij equal to 0.05 (denoted below as SRK) (2) SRK EOS with kij fitted from the experimental VLE data (SRK-KIJ) (3) Peng-Robinson EOS with kij equal to 0.06 (PR) (4) PR EOS with kij fitted from the experimental VLE data (PR-KIJ) (5) Peng-Robinson-Stryjek-Vera EOS with kij being a function of temperature as well as composition (PRSV)7-8 (6) Predictive SRK EOS9-10 (PSRK) (7) MHV2 model11 (MHV2) (8) Method based on the PRSV EOS with the modified Wong-Sandler mixing rules12 (WS-M) (9) Model based on the SRK EOS with the mixing rules worked out by Twu and Coon13 (TC) All the above methods have been incorporated into the VCS algorithm14 calculating CPE. The accuracy of the individual
10.1021/ie070260q CCC: $37.00 © 2007 American Chemical Society Published on Web 07/04/2007
5438
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007
approaches has been evaluated by comparison of the calculated results with the experimental data.1 2. Description of the Algorithm and the Methods Applied The VCS algorithm is a useful tool for calculation of the equilibrium compositions of multicomponent systems. Its comprehensive description and the program code have been presented by Smith and Missen.14 The original algorithm VCS refers to systems consisting of many pure condensed phases and at most two multicomponent ideal phases. In this paper, a modification allowing the treatment of the multicomponent phases as nonideal ones is introduced. The original algorithm calculates the Gibbs free energy change ∆G(m) of each indej pendent reaction j in each iteration (m) due to eq 1. N
∆G(m) j )
νij µ(m) ∑ i i)1
j ) 1, 2, ..., RN
(2)
where T is the temperature, P is the total pressure, R is the gas constant, µ°i is the standard chemical potential of species i referred to the ideal gas state, ni is the mole numbers of species i, nT is the total number of moles including inerts, and φVi is the fugacity coefficient of species i in the gas-phase calculated from the EOS. For the liquid phase, the equation expressing the chemical potential depends on the approach. In the activity coefficients models, the equation has the following form:
µi ) µ°i + RT ln
psi
+ RT ln
φsi
+
(4)
where φLi is the fugacity coefficient of liquid component i calculated from the EOS. The models SRK, SRK-KIJ, PR, and PR-KIJ use classical mixing rules:
amix )
bixi
∑i ∑j xixj(aiaj)0.5(1 - kij)
(5) (6)
The method PRSV uses eqs 7 and 8:
amix )
∑i ∑j xixj(aiaj)0.5(1 - xikij - xjkji)
(7)
(8)
0.42478R2Tc2 [f(Tr)]2 Pc
(9)
f(Tr) ) 1 + (0.48 + 1.574ω - 0.176ω2)(1 - Tr0.5)
(10)
f(Tr) ) 1 + C1(1 - Tr0.5) + C2(1 - Tr0.5)2 + C3(1 - Tr0.5)3 for Tr < 1 (11) f(Tr) ) 1 + C1(1 - Tr0.5) for Tr > 1
(12)
where constants C1-C3 are the parameters of pure components. The mixing rules are expressed by eqs 5 and 13 and the excess molar Gibbs free energy Gex γ is computed from the activity coefficients γi calculated from the UNIFAC method.15 amix ) bmix
(
Gex γ
-0.64663
+
ai
RT
∑x b + -0.64663∑x ln i
i
bi
i
i
)
bmix
i
(13)
The method MHV2 uses the SRK EOS with eqs 11 and 12. The mixing rules are defined by eqs 5 and 14, and the molar excess Gibbs free energy Gex γ is calculated with the modified version of the UNIFAC method.16
g1
(
amix
-
bmixRT
ai
∑i xib RT i
) [( ) ( ) ] amix
+ g2
2
-
bmixRT
Gex γ
psi )
µi ) µ°i + RT ln P + RT ln ni - RT ln nT + RT ln φLi
∑i
a)
+ RT ln ni RT ln nT + RT ln γi (3)
VLi (P
where psi is the saturated vapor pressure of pure component i, φsi is the fugacity coefficient of saturated vapor of pure component i, VLi is the liquid molar volume of pure component i, and γi is the activity coefficient of species i. In the case of the EOS approach, the chemical potentials of the liquid components are defined as:
bmix )
k(2) k(2) ij ji , kji ) k(1) + ji T T
For every binary mixture, the following constants should be (2) (1) (2) estimated: k(1) ij , kij , kji , and kji . The model PSRK uses the SRK EOS, but f(Tr) in eq 9 is calculated from eqs 11 and/or 12 not from the original eq 10.
(1)
The modification is based on introduction of the expressions describing the chemical potential µi of real components in an explicit form to eq 1. For the gas phase, the chemical potential is defined as follows:
µi ) µ°i + RT ln P + RT ln ni - RT ln nT + RT ln φVi
kij ) k(1) ij +
RT
ai
∑i xi b RT
2
)
i
+
∑i xi ln
bmix (14)
bi
The recommended values of constants g1 and g2 are -0.478 and -0.0047, respectively. The model WS-M uses the PRSV EOS and the mixing rules expressed by eqs 5 and 15. According to eq 16, the excess molar Helmholtz energy Aex,∞ is equal to the excess γ 12 The last value is com. molar Gibbs free energy Gex γ puted from the activity coefficients γi calculated with the UNIQUAC method having two adjustable parameters Aij and Aji. amix ) bmix
(
Aex,∞ γ
-0.62323
+
ai
RT
∑x b + -0.62323∑x ln i
i
i
i
i
)
bmix bi
(15)
Gex(T, x, P ) low) ) Aex(T, x, P ) low) ) Aex(T, x, P ) ∞) (16) The version of model TC applied in this work uses the SRK EOS, and the parameter a of the individual components is calculated from eq 17.
a)
0.42478R2Tc2 N(M-1) Tr exp(L(1 - TNM r )) Pc
(17)
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007 5439
mixing rules with one or two adjustable parameters. It is applicable to the systems containing only nonpolar compounds or for the mixtures with light gases. 3. Thermodynamic Data Referring to the Chemical Equilibrium
Figure 1. Plot of log(Ka) versus 1000/T for values obtained from eq 24 and those calculated from various thermodynamic data.
where L, M, and N are the parameters of pure components. The methods PSRK and MHV2 use an ideal solution as a reference fluid to develop the EOS mixing rules. The model TC uses the van der Waals fluid as a reference state where avdW and bvdW have the following definitions:
avdW )
bvdW )
∑i ∑j xixj(aiaj)0.5(1 - kij)
( )
∑j ∑i xixj
b i + bj 2
(18)
(1 - lij)
avdW bvdW RT bmix ) avdW AE 1RTbvdW RT ln(2) amix ) bmix
(
avdW AE bvdW ln(2)
(19)
)
)
(20)
(21)
N
AE RT
N
)
xj ∑ j)1
∑xkτkjGkj k)1 N
(22)
∑xkGkj
k)1
where N is the number of components and the remaining parameters have the following definitions:
τij )
log(Ka) )
Aij , Gij ) exp(-Rijτij), Aij * Aji, Rij ) Rji T
(23)
The mixing rules worked out by Twu and Coon are more flexible than those in the other models used in this paper. When the system contains polar substances, five adjustable parameters can be applied. Three of them describe the excess molar Helmholtz energy and the remaining two represent the van der Waals interaction coefficients (kij and lij). However, when AE is set to zero, eqs 20 and 21 are reduced to the van der Waals
649 + 0.042 T
(24)
The numerical values of eq 24 were calculated with the assumption that the chemical equilibrium constant of acetic acid dimerization has the following definition:
log(Kdim) )
where kij and lij are the binary interaction coefficients. The mixing rules for the parameters a and b are defined by eqs 20 and 21 and the excess molar Helmholtz energy AE is calculated from the NRTL equation.
(
In the case of CPE, the calculated composition depends not only on the model coefficients fitted from the experimental data and describing the phase equilibrium but also on standard chemical potentials µ°i responsible for the chemical equilibrium. These potentials are identified with the standard Gibbs free energies of formation. Russel and Kabel17 performed a critical analysis of the published experimental data referring to the chemical equilibrium constant for the esterification of acetic acid with ethanol in the gas phase. As a result of the analysis, the temperature function of the gas-phase equilibrium constant Ka in form 24 was found:
3000 - 10.149 T
(25)
The numerical values of Ka can also be calculated from thermodynamic tables if they contain the data needed for computation of the standard Gibbs free energy of reaction at T (∆G°r,T).
Ka ) exp
(
)
-∆G°r,T RT
(26)
Figure 1 shows the lines representing temperature functions of Ka calculated from the following: eq 24, the thermodynamic tables of Stull et al.18 and Reid et al.,19 and the ChemCAD databank.20 It can be seen that values of Ka obtained from the thermodynamic tables or the ChemCAD databank differ significantly from the values calculated from eq 24 determined on the basis of experiments. These discrepancies are caused by the inaccurate table data of acetic acid. The temperature function of standard chemical potential can be precisely described by eq 27. Table 1 presents the coefficients A1-A5 of this equation determined from tables18 for the needed components. In the case of acetic acid, two sets of these coefficients are presented. The first set was obtained from the original data, and the second one was fitted to the values calculated from eq 24. The fitting was performed in the following way: - At first, the values of Ka were calculated for several T values from eq 24. - Next these values were put into eq 26 to obtain several µ°HAc values for the considered temperatures (µ°i of the remaining components were calculated from the coefficients presented in Table 1). - In the last step, the “fitted” coefficients A1-A5 were computed by the least-squares method.
µ°i ) A1 + A2T + A3T2 + A4T3 + A5T ln(T)
(27)
Blanchard and Brennecke1 stated that equilibrium conversion of the considered esterification process in the liquid
5440
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007
phase without CO2 reached 63% at 60 °C. Table 2 presents results of equilibrium compositions calculations for both sets of coefficients A1-A5 referring to acetic acid. The activity coefficients were calculated by the modified UNIFAC method,16 and coefficients φsi in eq 3 were calculated by Hayden and O’Connel’s version21 of the virial equation. The corrected values of acetic acid chemical potentials allow the prediction of the equilibrium conversion with much better precision. Kang et al.22 published the data referring to CPE in the esterification of acetic acid with ethanol. The data contain the following: mole fractions of the individual components in liquid and vapor phases and temperature and pressure for each experimental point. The individual phases fractions and the initial feed compositions are not given. These data have also been used for verification of the corrected potentials of acetic acid. Transformation of eq 3 allows the definition of the chemical equilibrium constant in the liquid phase in two manners:
Ktab ) exp
[∑ ( -νi
i
µ°i
RT
+ ln
Kexp )
psi
+ ln
φsi
+
)]
VLi (P - psi ) RT
∏i (xiγi)ν ) ∏i (ai)ν i
i
(28) (29)
The constant Ktab is calculated basing on the individual component table data and the coefficients φsi independent of the composition. The constant Kexp defined by eq 29 represents the value obtained from the experimental activities calculated on the base of the experimental mole fractions and the activity coefficients computed by a chosen method. If the experimental results and the table data were perfect, the quotient Kexp/Ktab should be equal to 1. To check the table data accuracy, such quotients were calculated for all the experimental points presented in ref 22. Next, these values were used to calculate the arithmetic mean due to eq 30. Four methods were chosen to compute the activity coefficients: the ASOG method,23 the Wilson equation with the coefficients taken from ref 24, and two modified versions of the UNIFAC method (the UNIF-87 from ref 16 and the UNIF-93 from ref 25). Table 3 presents the AVQK values calculated for every nonideality description method and for two ways of acetic acid chemical potential calculation (from the original tables, AVQK-ORG, and from the corrected values due to eq 24, AVQK-NEW). For the original acetic acid data, the ratio Kexp/Ktab often exceeds 2 independently of the nonideality description method. On the other hand, the AVQK-NEW values are close to 1 and this fact confirms that the acetic acid chemical potential data should be corrected. The UNIF-87 method gives the AVQK-NEW value the closest to 1. Table 3 presents also the values ∆y and ∆T (defined by eqs 31) representing mean deviations in phase equilibrium prediction. These values prove that the UNIF-87 method describes the mixture nonideality with the best precision.
AVQK )
∆y )
1
1
NP
∑
Kexp
NP i)1 Ktab
(30)
NP
∑ |yi,eksp - yi,obl|100,
NP i)1
∆T ) 1
NP
∑ |Ti,eksp - Ti,obl|
NP i)1
(31)
∆P )
1
Pi,eksp - Pi,obl
∑| NP
NP i)1
|
100
Pi,eksp
(32)
The results from Tables 2 and 3 clearly prove that the corrected values of acetic acid chemical potentials accurately describe various experimental data referring to chemical equilibrium. Therefore, these values were used in the next parts of this work. 4. Estimation of Parameters for the Individual Models The models SRK, PR, and PSRK are entirely predictive. One parameter kij should be fitted for the SRK-KIJ and PR-KIJ methods, and two UNIQUAC parameters (Aij and Aji), for the method WS-M. The model PRSV needs at most four parameters (2) (1) (2) (k(1) ij , kij , kji , and kji ), while the method TC can use five adjustable coefficients (Aij, Aji, Rij, kij, and lij). The model MHV2 using the UNIF-87 method needs the group interaction parameters for the main groups COOH-CO2. The estimation of the parameters was carried out by minimization of the objective function 33 using the procedure BSOLVE.26 Depending on the requirements, the weight factor WT enables better pressure or composition fitting. The critical parameters and the acentric factors used in the individual models were taken from ref 19.
OF ) WT
∑i
[
] [ ] [ ]
Pexp - Pobl i i Piexp
2
+
∑i
exp obl y1,i - y1,i
2
+
yexp 1,i
∑i
exp obl y2,i - y2,i
yexp 2,i
2
(33)
The system considered consists of 5 components forming 10 binary mixtures. The following data were used for estimation of the model parameters for the individual two-component systems: EtOH-HAc,27,28 EtOH-EtAc,29 EtOH-H2O,30,31 EtOH-CO2,32,33 HAc-EtAc,28,34 HAc-H2O,28,35 HAc-CO2,36,37 EtAc-H2O,38-40 EtAc-CO2,41-43 and H2O-CO2.37 The highpressure VLE data were used for the binary mixtures containing CO2 and the system EtOH-H2O30 whereas for the remaining mixtures the atmospheric pressures data were applied. The binary interaction parameters estimation for the mixture EtAc-CO2 caused many troubles. The parameter kij values of this system computed from the data,41-43 with the models used differed significantly. For the model SRK-KIJ, they were equal -0.0271, -0.2050, and -0.5020, respectively. The same refers to the models PR-KIJ, PRSV, WS-M, and TC. The interaction coefficients estimated in this way used for calculations of VLE of the ternary mixture CO2-EtOH-EtAc44 gave very inaccurate prediction for all the considered models. This prediction with the methods PSRK and MHV2 was also very poor. Apparently, the authors of these methods have estimated the interaction parameters of the main groups COOC-CO2 from the data of refs 41 or/and 42 which are not reliable. Because of these difficulties, the binary interaction parameters of the mixture EtAc-CO2 were estimated from the data of the ternary mixture CO2-EtOH-EtAc for all the considered models. In this estimation, the coefficients of the remaining two binary mixtures (CO2-EtOH and EtOH-EtAc) were treated as the known valuessevaluated previously from the suitable binary data. For the two methods (PSRK and MHV2) using two versions of the UNIFAC model, the ternary data mentioned above were used
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007 5441
Figure 2. Mean deviations ∆y and ∆P predicted by the methods PSRK and MHV2 and obtained during the parameter estimations for the models SRK-KIJ, PR-KIJ, PRSV, WS-M, and TC. Table 1. Coefficients A1-A5 of Eq 27a
a
compound
A1
A2 × 102
A3 × 106
A4 × 1010
A5 × 102
EtOH HAc (values obtained from the original data) HAc (values fitted to eq 24) EtAc H2O
-216.8542 -420.4794
-23.9665 -18.7803
-40.5774 -46.3813
31.2574 63.8043
7.27096 6.10801
-425.6300 -416.5507 -238.3579
-11.3186 -31.1652 -4.20032
-48.6888 -82.8030 -6.46282
88.9941 105.6811 14.5700
5.20598 11.13450 1.34245
The term µ°i is in kilojoules per mole, and T is in kelvin.
Table 2. Calculated Equilibrium Compositions for Both Sets of Coefficients A1-A5 in the Esterification of Acetic Acid with Ethanol for Initial Feed Composition of 1 mol EtOH and 1 mol HAc, T ) 333 K, and P ) 1 bara A1-A5 of HAc obtained from the original data18 xi γi EtOH HAc EtAc H2O conversion (%) a
0.2219 0.2219 0.2781 0.2781
1.1320 0.8410 1.6714 2.2381 55.62
A1-A5 of HAc fitted to eq 24 xi γi 0,1758 1,1746 0,1758 0,8188 0,3242 1,6735 0,3242 2,2201 64.83
The mixture after the reaction is the liquid phase only.
for estimation of the new values of the interaction parameters of COOC-CO2 groups. Table 4 presents the accuracy of VLE prediction for the system CO2-EtOH-EtAc with the parameters estimated from the various data using the PRSV model. The predictions were very inaccurate for the parameters estimated from the binary data of the mixture EtAc-CO2. The mean deviations ∆P obtained during the estimation of the EtAc-CO2 binary interaction coefficients from the ternary data are still significant but considerably lower than in the case of VLE prediction for the parameters estimated from the binary data. The last four columns of Table 4 present the mean deviations ∆y and ∆P obtained with the PSRK model for the ternary mixture CO2EtOH-EtAc using the original values of COOC-CO2 group interaction parameters and the new values estimated from these ternary data. The deviations are much lower for the new values. The calculated group interaction coefficients for the main groups COOC-CO2 have the following numerical values: (336.479, -44.0561) and (2621.66, -84.9289) for the models PSRK and MHV2, respectively. The last model did not contain these coefficients for the main groups COOH-CO2. They were
estimated using data from refs 36 and 37, and their numerical values are as follows: (196.918, 154.993). Tables 5 and 6 present the binary interaction parameters for the following models: SRK-KIJ, PR-KIJ, PRSV, WS-M, and TC. Since the parameters of binary mixture EtAc-CO2 were evaluated from the ternary data, the temperature function of the coefficient kij (eq 8) was not considered. Therefore, the (2) parameters k(2) ij and kji for this mixture were equal to zero. For the model TC in the case of binary mixtures with carbon dioxide, the fitted parameters represent the van der Waals interaction coefficients only. Figure 2 presents the mean deviations ∆y and ∆P obtained during estimations of the binary interaction parameters for the models SRK-KIJ, PR-KIJ, PRSV, WS-M, and TC as well as these deviations illustrating the accuracy of the VLE prediction for the methods PSRK and MHV2. The analysis of Figure 2 proves that the model PRSV with the parameter kij being the function of temperature and composition is definitively the best correlation method for the system containing EtOH, HAc, EtAc, H2O, and CO2. The model WS-M takes second place, and the model TC is the third one in the ranking. Both models are inaccurate in the case of the H2OCO2 mixture. For the majority of the binary mixtures, the predictive methods PSRK and MHV2 using the UNIFAC model are less accurate than the TC model. It should be emphasized however that both predictive methods work quite well in the case of the systems H2O-CO2 and EtAc-H2O. The last mixture is difficult for modeling since its components do not mix in the entire range of concentrations. The models SRK-KIJ and PRKIJ are the least accurate. In the case of the mixtures EtAcH2O and H2O-CO2, both models are very imprecise. All the investigated models, even those with four or five adjusted parameters, have a problem with very accurate cor-
5442
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007
Table 3. Comparison of the Values AVQK-ORG, AVQK-NEW, ∆y, and ∆T for the Considered Methods Describing the Mixture Nonideality AVQK-ORG AVQK-NEW ∆y (mol %) EtOH HAc EtAc H2O
ASOG
Wilson
UNIF-87
UNIF-93
2.478 1.108
2.600 1.165
2.065 0.925
1.936 0.867
∆T (K)
4.90 1.38 3.27 1.20
∆y (mol %)
∆T (K)
0.73 1.02 2.10 1.51
3.41
∆y (mol %)
∆T (K)
0.32 0.41 1.41 1.74
1.42
∆y (mol %)
∆T (K)
0.58 1.04 2.13 1.92
0.51
1.63
Table 4. Accuracy of VLE Prediction of the Models PRSV and PSRK for the System CO2-EtOH-EtAc PRSV 1a ∆y CO2 EtOH EtAc
0.97 0.55 0.53
PSRK
2a ∆P
∆y
20.98
1.17 0.65 0.64
3a ∆P
∆y
34.18
0.60 0.13 0.63
4a ∆P
∆y
41.50
0.88 0.20 0.87
5b ∆P
∆y
10.98
0.98 0.55 0.52
6c ∆P
∆y
∆P
23.27
0.42 0.19 0.29
13.72
a The parameters were estimated from the binary data41-43 and the ternary system,44 respectively. b Results obtained for the original COOC-CO group 2 interaction coefficients. c Results obtained for the new COOC-CO2 group interaction coefficients estimated from the ternary data.
Table 5. Parameters of the Models SRK-KIJ, PR-KIJ, and PRSV Fitted for All the Considered Binary Mixtures SRK-KIJ
EtOH-HAc EtOH-EtAc EtOH-H2O EtOH-CO2 HAc-EtAc HAc-H2O HAc-CO2 EtAc-H2O EtAc-CO2 H2O-CO2
PR-KIJ
PRSV
kij
kij
k(1) ij
-0.03052 0.03284 -0.09430 0.09327 -0.02918 -0.1372 0.03654 -0.2844 0.1740 -0.1093
-0.03387 0.02550 -0.09980 0.09490 -0.02910 -0.1543 0.04291 -0.2812 0.1632 -0.08958
-0.03220 0.6410 -0.00318 0.3527 -0.02621 0.06132 0.2101 1.009 0.6491 0.04616
k(1) ji
k(2) ij
k(2) ji
-0.08657 -0.9147 -0.05528 -0.005543 -0.01637 -0.07621 -0.08087 -0.6665 0.1161 6.679
-0.2835 -217.8 -21.96 -80.34 -4.358 -70.91 -46.80 -359.3 0 -40.69
16.63 335.7 -20.12 31.92 1.034 -38.89 31.57 127.8 0 -2228.0
Table 6. Parameters of the Models WS-M and TC Fitted for All the Considered Binary Mixtures WS-M EtOH-HAc EtOH-EtAc EtOH-H2O EtOH-CO2 HAc-EtAc HAc-H2O HAc-CO2 EtAc-H2O EtAc-CO2 H2O-CO2
TC
Aij
Aji
Aij
Aji
aij
kij
lij
-91.781 -53.070 -39.255 169.72 -236.53 730.12 436.34 747.36 2480.2 963.81
-8.1818 250.40 286.28 56.736 422.66 -370.01 -243.36 -40.423 -130.79 -147.73
25.99 211.9 -586.1 0 -260.8 803.4 0 871.0 0 0
-391.8 365.4 333.1 0 1628.0 -62.45 0 204.7 0 0
0.6003 0.7353 0.7777 0 0.6195 0.8000 0 0.3086 0 0
0.04012 -0.02237 0.02386 0.08738 0.01984 -0.1694 -0.007281 -0.2249 0.1450 -0.05143
0 0 0 -0.008887 0 0 -0.03285 0.06735 0 0.02471
relation of the data for two mixtures: EtOH-HAc and HAcEtAc. These two systems were also not very precisely described by the methods based on the activity coefficients γi.3 These facts indicate that the VLE data accuracy of the considered two mixtures is not very high. 5. Results of CPE Calculations and their Discussion The esterification of acetic acid with ethanol in CO2 was studied experimentally1 at 333 K and 58.6 bar. The initial feed consisted of equimolar mixtures of reactants and products (with the ratio nreactants:nproducts ) 3:1) and carbon dioxide (with an initial mole fraction equal to 0.3). The equilibrium product contained two phasessgaseous and liquidstheir concentrations are given in Table 7. The authors of ref 1 stated that equilibrium conversion of this esterification process can be shifted from 63% without carbon dioxide to 72% with CO2. The fractions of both phases in equilibrium were not given.1 To obtain these data,
they were calculated on the basis of the carbon dioxide concentration in the initial feed and the final equilibrium in both phases. Assuming that the mole fraction of inert CO2 (I) in the initial feed is equal to 0.3 and the amounts of the remaining components result from the proportions given above, the initial elemental balance can be described as follows: C ) 1.4, H ) 3.5, O ) 1.05, and I ) 0.3. The final balance for the experimental equilibrium product presented in Table 7 is described as follows: C ) 1.3634, H ) 3.4834, O ) 1.0217, and I ) 0.3. So, the amounts of elements before and after reaction are not exactly the same (some amounts of C, H, and O disappeared), but the balance accuracy is fairly good. Table 7 contains the experimental data and results of CPE calculations for all the considered models. The influence of the nonideality modeling method on the calculated equilibrium composition is significant. The simplest models SRK and PR fail completely. They predict that the system in equilibrium is
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007 5443 Table 7. Comparison of the Calculated Equilibrium Results (mole percent) with the Experimental Data exp data EtOH Hac EtAc H2O CO2 fraction, %
3.00 4.00 5.00 2.00 86.00 3.45
EtOH Hac EtAc H2O CO2 fraction, % conv, %
12.00 8.00 25.00 27.00 28.00 96.55 72.00
SRK
10.43 10.43 24.57 24.57 30.00 100.0 70.20
SRK-KIJ
PR
PR-KIJ
PRSV
PSRK
MHV2
WS-M
TC
0.22 0.04 1.03 0.48 98.23 7.74
Gas Phase 0.23 0.04 1.10 0.50 98.13 7.66
0.30 0.08 2.03 0.66 96.93 4.38
0.26 0.05 1.38 0.58 97.73 6.03
0.28 0.07 1.75 0.46 97.44 15.80
0.28 0.05 1.78 0.55 97.34 6.19
0.16 0.06 1.21 0.43 98.14 8.80
5.97 5.99 31.86 31.90 24.28 92.26 84.21
Liquid Phase 10.39 6.22 10.39 6.24 24.61 31.57 24.61 31.62 30.00 24.35 100.0 92.34 70.30 83.53
10.57 10.58 25.93 25.99 26.93 95.62 71.08
8.93 8.94 28.21 28.26 25.66 93.97 75.98
10.06 10.11 31.12 31.36 17.35 84.20 75.66
10.39 10.44 26.78 26.86 25.56 93.81 72.10
10.38 10.38 27.87 27.95 23.42 91.20 72.92
a liquid phase only. If one parameter, kij, is fitted from the experimental VLE data, these models are able to predict the existence of a two-phase system in equilibrium. However, both models SRK-KIJ and PR-KIJ predict too high conversion (8384%) and the liquid-phase compositions differ considerably from the experimental ones. The PRSV model gives the best description, predicting as well the liquid-phase composition as the fraction of both phases very similar to the experimental data. The calculated equilibrium conversion is equal to 71.08% and differs by less than 1% from the measured value. Similarly as in the previous paragraph, the WS-M method takes the second place in the ranking. This model calculates the conversion nearly equal to the experimental one but predicts the fraction of both phases and the liquid-phase composition a bit less accurately than the PRSV method. The model TC predicts accurately the equilibrium conversion but calculates a too low concentration of CO2 in the liquid phase. The models PSRK and MHV2 give also fairly good results even though both methods predict a bit too high equilibrium conversion. Model MHV2 calculates also the greatest amount of the gas phase and in consequence a too low concentration of CO2 in the liquid phase. Taking into account the accuracy of the VLE correlation, the models considered have been classified from the most to the least precise in the previous paragraph. This order is almost identical with regard to the accuracy of CPE prediction. All the investigated models predict very similar compositions of the gas phase with a CO2 concentration of 97-98%. From among the remaining components, ethyl acetate is present in the biggest amount. In the experimental data, the concentration of carbon dioxide is lower (86%) but the ester takes the second place with respect to the amount. It should be noticed that the gas fraction is small. In such a case, the inevitable disagreements in the elemental balance between the experimental data and the calculated results have much greater influence on the numerical values of the gas-phase concentrations than in the case when fractions of both phases would be comparable. The CPE calculations were also performed using the interaction parameters of the models estimated from the VLE data of the binary mixture EtAc-CO2.41-43 In all the cases, independently of the model applied, the calculated equilibrium product was the liquid phase which did not reach the bubble point. This confirms the earlier conclusion that the binary VLE data of the mixture EtAc-CO2 are unreliable. 6. Conclusions Although the reacting mixture in the esterification of acetic acid with ethanol in CO2 at high pressure is a system difficult to model, it can be described with different accuracy by models
based on the cubic equations of state. The accurate prediction of chemical equilibrium is also essential for this modeling. Reference 17 presents the temperature dependency of the chemical equilibrium constant of the considered esterification in the gas phase (eq 24) being the result of the critical review of the experimental data. Figure 1 has shown however that the line obtained from eq 24 differs significantly from the lines calculated on the basis of various thermodynamic tables. Since it was assumed that acetic acid standard chemical potentials were responsible for those differences, they were corrected to satisfy eq 24 together with the standard chemical potentials of the remaining components. The usefulness of the modified temperature dependency of acetic acid standard chemical potential was tested on the basis of other experimental data referring to the esterification of acetic acid with ethanol in the liquid phase (Table 2) and separately in a two-phase system (Table 3). Both tables clearly indicated that the modified dependency gave the results much closer to the experimental values than the temperature function obtained from the original tables. The binary interaction coefficients of the models considered were estimated using two-component experimental VLE data. The only exception was the mixture EtAc-CO2, since its experimental data seem to be unreliable. The parameters estimated from these binary data gave very inaccurate VLE prediction of the ternary mixture CO2-EtOH-EtAc independently of the model applied. These parameters predicted also the one-phase system in the CPE calculations of the esterification of acetic acid with ethanol in CO2 for all the models. Therefore, the binary parameters of the system EtAc-CO2 were estimated from the data of the ternary mixture CO2-EtOH-EtAc. For the methods PSRK and MHV2, the ternary data were used to calculate new numerical values of the COOC-CO2 group interaction parameters of the suitable versions of the UNIFAC method. Analysis of the mean deviations ∆y and ∆P obtained during the individual model parameters estimations proved the greatest accuracy of the PRSV model. The methods WS-M and TC were less precise than the PRSV, but they were more accurate than the models using the group contribution methods (PSRK and MHV2). The models SRK-KIJ and PR-KIJ were the least accurate. The same hierarchy of accuracy was found when the results of CPE calculations were compared with the experimental data.1 Only the models SRK and PR failed completely predicting the one-phase system in the equilibrium state. The most precise model PRSV predicts the equilibrium conversion and fractions of the individual phases with an accuracy of 1%, whereas the liquid-phase concentration was predicted with an accuracy of
5444
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007
2-3%. All the models (with the exception of SRK and PR) calculated very similar compositions of the gas phase that were different, however, from the experimental data. The reasons for this disagreement are as follows. The binary VLE data referred to high pressure for the mixtures with CO2 and the system EtOH-H2O only. The experimental data for the remaining twocomponent mixtures are available only for the atmospheric pressures. The data for the systems EtOH-HAc and HAcEtAc seem to be not very accurate since the mean deviations ∆y and ∆P obtained during the correlation were relatively high for all the models considered. These two systems were also not very precisely described by the methods based on the activity coefficients γi.3 The experimental amounts of the individual elements in the feed and product1 are similar but not exactly the same. Since the amount of gas phase is small, the disagreements in the elemental balance can have quite significant influence on the differences between the calculated and measured values of the individual component concentrations. Nomenclature A1-A5 ) constants in eq 27 Aij, Aji ) parameters of NRTL or UNIQUAC equations (K) AE ) molar excess Helmholtz free energy (J/mol) AVQK ) factor defined by eq 30 a ) attraction parameter in the equation of state (Pa m6/mol2) b ) covolume (m3/mol) C1-C3 ) constants in eqs 11 and 12 Gxγ ) molar excess Gibbs free energy computed from activity coefficients γi (J/mol) ∆G(m) j ) Gibbs free energy change of reaction j in iteration (m) (J/mol) g1 and g2 ) constants of the MHV2 model Ka ) chemical equilibrium constant in the gas phase Kdim ) chemical equilibrium constant of acetic acid dimerization in the gas phase (mmHg-1) kij ) binary interaction parameter (2) (1) (2) k(1) ij , kij , kji , and kji ) parameters to be estimated in model PRSV L ) constant in eq 17 lij ) binary interaction parameter M ) constant in eq 17 N ) number of species or constant in eq 17 NP ) number of experimental points ni ) number of moles of species i nT ) total number of moles in a phase, including inert species P ) total pressure (Pa or bar) psi ) saturated vapor pressure of pure component i (Pa or bar) ∆P ) absolute mean deviation between experimental and calculated equilibrium pressure, relative to the experimental pressure (%) R ) gas constant (J/(mol K)) RN ) number of independent reactions T ) temperature (K) ∆T ) absolute mean deviation between experimental and calculated equilibrium temperature (K) VLi ) molar volume of pure liquid species i (m3/mol) WT ) weight factor x, y ) mole fractions in the liquid and the vapor phase, respectively ∆y ) absolute mean deviation between experimental and calculated equilibrium composition in the vapor phase (mol %)
Greek Letters Rij ) parameter of NRTL equation γi ) activity coefficient of species i µi ) chemical potential of species i (J/mol) µ°i ) standard chemical potential of species i (referred to the ideal gas state) (J/mol) νij ) stoichiometric coefficient of species i in reaction j φVi , φLi ) fugacity coefficient of species i in the vapor and the liquid phase, respectively φsi ) fugacity coefficient of saturated vapor of pure component i (calculated for T and psi ) ω ) acentric factor Subscripts c ) critical property cal ) calculated exp ) experimental mix ) mixture r ) reduced property tab ) calculated from table values vdW ) van der Waals fluid AbbreViations EtOH ) ethanol HAc ) acetic acid EtAc ) ethyl acetate Literature Cited (1) Blanchard, L. A.; Brennecke, J. F. Esterification of acetic acid with ethanol in carbon dioxide. Green Chem. 2001, 3, 17. (2) Wyczesany, A. Application of equation of state methods for chemical equilibria calculation in strongly non-ideal liquid phase. Inz˘ ynieria Chemiczna i Procesowa 2003, 24, 411 (in polish). (3) Wyczesany, A. Nonstoichiometric Algorithm of Calculation of Simultaneous Chemical and Phase Equilibria. 1. Influence of the Method Modeling Nonideality of Systems on the Calculated Equilibrium Composition at Low Pressure. Ind. Eng. Chem. Res. 1993, 32, 3072. (4) Burgos-Solo´rzano, G. I.; Brennecke, J. F.; Stadtherr, M. A. Validated computing approach for high-pressure chemical and multiphase equilibrium. Fluid Phase Equilib. 2004, 219, 245. (5) Soave, G. Equilibrium constants from a modified Redlich-Kwong equation of state. Chem. Eng. Sci. 1972, 27, 1197. (6) Peng, D.; Robinson, D. B. A New Two-Constant Equation of State. Ind. Eng. Chem. Fundam. 1976, 15, 59. (7) Stryjek, R.; Vera, J. H. PRSV: An Improved Peng-Robinson Equation of State for Pure Compounds and Mixtures. Can. J. Chem. Eng. 1986, 64, 323. (8) Stryjek, R.; Vera, J. H. PRSV: An Improved Peng-Robinson Equation of State with New Mixing Rules for Strongly Nonideal Mixtures. Can. J. Chem. Eng. 1986, 64, 334. (9) Holderbaum, T.; Gmehling, J. PSRK: A Group Contribution Equation of State Based on UNIFAC. Fluid Phase Equilib. 1991, 70, 251. (10) Fischer, K.; Gmehling, J. Further development, status and results of the PSRK method for the prediction of vapor-liquid equilibria and gas solubilities. Fluid Phase Equilib. 1996, 121, 185. (11) Dahl, S.; Michelsen, M. L. High-Pressure Vapor-Liquid Equilibrium with a UNIFAC-Based Equation of State. AIChE J. 1990, 36, 1829. (12) Orbey, H.; Sabdler, S. I. A comparison of Huron-Vidal type mixing rules of mixtures of compounds with large size differences, and a new mixing rule. Fluid Phase Equilib. 1997, 132, 1. (13) Twu, C. H.; Coon, J. E. CEOS/AE Mixing Rules Constrained by vdW Mixing Rule and Second Virial Coefficient. AIChE J. 1996, 42, 3212. (14) Smith, W. R.; Missen, R. W. Chemical Reaction Equilibrium Analysis: Theory and Algorithms; Wiley: New York, 1982. (15) Hansen, H. K.; Rasmussen, P.; Fredenslund, Aa.; Schiller, M.; Gmehling, J. Vapor-Liquid Equilibria by UNIFAC Group Contribution. 5. Revision and Extension. Ind. Eng. Chem. Res. 1991, 30, 2352. (16) Larsen, B. L.; Rasmussen, P.; Fredenslund, Aa. A modified UNIFAC group-contribution model for prediction of phase equilibria and heats of mixing. Ind. Eng. Chem. Res. 1987, 26, 2274.
Ind. Eng. Chem. Res., Vol. 46, No. 16, 2007 5445 (17) Russel, W. H.; Kabel, R. L. Thermodynamic Equilibrium in the Vapor Phase Esterification of Acetic Acid with Ethanol. AIChE J. 1968, 14, 606. (18) Stull, R. D.; Prophet, E. W.; Sinke, G. C. The Chemical Thermodynamics of Organic Compounds; John Wiley & Sons, Inc.: New York, 1969. (19) Reid, R. C.; Prausnitz, J. M.; Poling, B. E. The Properties of Gases and Liquids, 4th ed.; McGraw-Hill: New York, 1987. (20) ChemCAD Software, version 5.6; Chemstations Inc.: Houston, 2006. (21) Hayden, J. G.; O’Connel, J. P. A Generalized Method for Predicting Second Virial Coefficients. Ind. Eng. Chem., Process Des. DeV. 1975, 14, 209. (22) Yun, W. K.; Youn, Y. L.; Won, K. L. Vapor-liquid equilibria with chemical reaction equilibrium. Systems containing acetic acid, ethyl alcohol, water and ethyl acetate. J. Chem. Eng. Jpn. 1992, 25, 649. (23) Kojima, K.; Tochigi, K. Prediction of Vapor - liquid equilibria by ASOG method; Kodansa: Tokyo, 1979. (24) Barbosa, D.; Doherty, M. F. The influence of equilibrium chemical reactions on vapor-liquid phase diagrams. Chem. Eng. Sci. 1988, 43, 529. (25) Gmehling, J.; Li, J.; Schiller, M. A Modified UNIFAC Model. 2. Present Parameter Matrix and Results for Different Thermodynamic Properties. Ind. Eng. Chem. Res. 1993, 32, 178. (26) Kuester, I. L., Mize, J. H. Optimization technique with FORTRAN; McGraw-Hill, INC: New York, 1973. (27) Rius, A.; Otero, J. L.; Macarron, A. Equilibres liquide-vapeur de me´langes binaires donnant une re´action chimique: syste`mes me´thanol-acide ace´tique; e´thanol-acide ace´tique; n-propanol-acide ace´tique; n-butanol-acide ace´tique. Chem. Eng. Sci. 1959, 10, 105. (28) Miyamoto, S.; Nakamura, S.; Iwai, Y.; Arai, Y. Measurement of Isothermal Vapor-Liquid Equilibria for Binary and Ternary Systems Containing Monocarboxylic Acid. J. Chem. Eng. Data 2001, 46, 1225. (29) Furnas, C. C.; Leighton, W. B. Vapor-liquid equilibrium data. Ethyl alcohol-ethyl acetate and acetic acid-ethyl acetate systems. Ind. Eng. Chem. 1937, 29, 709. (30) Barr-David, F.; Dodge, B. F. Vapor-Liquid Equilibrium at High Pressures. The Systems Ethanol-Water and 2-Propanol-Water. J. Chem. Eng. Data 1959, 4, 107. (31) Ma¸czyn´ski, A.; Ma¸czyn´ska, Z. Verified Vapor-liquid equilibrium data. Binary one-liquid systems of water and organic compounds; PWN: Warszawa, 1981; Vol. 5, p 57, data 40. (32) Jennings, D. W.; Leem, R. J; Teja, A. S. Vapor-Liquid Equilibria in the Carbon Dioxide + Ethanol and Carbon Dioxide + l-Butanol Systems. J. Chem. Eng. Data 1991, 36, 303. (33) Joung, S. N.; Yoo, C. W.; Shin, H. Y.; Kim, S. Y.; Yoo, K. P.;
Lee, C. S.; Huh, W. S. Measurements and correlation of high-pressure VLE of binary CO2-alcohol systems (methanol, ethanol, 2-metoxyethanol and 2-ethoxyethanol). Fluid Phase Equilib. 2001, 185, 219. (34) Macedo, E. A.; Rasmussen, P. Vapor-Liquid Equilibrium for the Binary Systems Ethyl Acetate-Acetic Acid and Ethyl Propionate-Propionic Acid. J. Chem. Eng. Data 1982, 27, 463. (35) Ito, T.; Yoshida, F. Vapor-Liquid Equilibria of Water-Lower Fatty Acid Systems: Water-Formic Acid, Water-Acetic Acid and WaterPropionic Acid. J. Chem. Eng. Data 1963, 8, 315. (36) Jonasson, A.; Persson, O.; Rasmussen, P.; Soave, G. S. Vaporliquid equilibria of systems containing acetic acid and gaseous components. Measurements and calculations by a cubic equation of state. Fluid Phase Equilib. 1998, 152, 67. (37) Bamberger, A.; Sieder, G.; Maurer, G. High pressure (vapor + liquid) equilibrium in binary mixtures of (carbon dioxide + water or acetic acid) at temperatures from 313 to 353 K. J. Supercrit. Fluids 2000, 17, 97. (38) Kato, M.; Konishi, H.; Hirata, M. New apparatus for isobaric dew and bubble point method. Methanol-water, ethyl acetate-ethanol, water1-butanol, and ethyl acetate-water systems. J. Chem. Eng. Data 1970, 15, 435. (39) Metrel, I. Vapor-liquid equilibrium. IL. Phase equilibria in the ternary system ethyl acetate-ethanol-water. Collect. Czech. Chem. Commun. 1972, 37, 366. (40) Ma¸czyn´ski, A.; Skrzecz, A. Verified Vapor-liquid equilibrium data. Binary two-liquid systems of water and organic compounds; PWN: Warszawa, 1983; Vol. 7, p 53, data 7. (41) Chrisochoou, A.; Schaber, K.; Bolz, U. Phase equilibria for enzymecatalyzed reactions in supercritical carbon dioxide. Fluid Phase Equilib. 1995, 108, 1. (42) Wagner, Z.; Pavlicˇek, J. Vapour-liquid equilibrium in the carbon dioxide-ethyl acetate system at high pressure. Fluid Phase Equilib. 1994, 97, 119. (43) Vazquez da Silva, M.; Barbosa, D.; Ferreira, P. O.; Mendonc¸ a, J. High pressure phase equilibrium data for the systems carbon dioxide/ehtyl acetate and carbon dioxide/isoamyl acetate at 295.2, 303.2 and 313.2 K. Fluid Phase Equilib. 2000, 175, 19. (44) Cheng, C.-H.; Chen, Y.-P. Vapor-liquid equilibria for the ternary system of carbon dioxide + ethanol + ethyl acetate at elevated pressures. Fluid Phase Equilib. 2006, 242, 169.
ReceiVed for reView February 19, 2007 ReVised manuscript receiVed May 9, 2007 Accepted May 30, 2007 IE070260Q
