VAPOR-LIQUID EQUILIBRIUM CALCULATION FOR SIMULATION OF BIOETHANOL CONCENTRATION FROM SUGARCANE

The robustness of the simulation of bioethanol concentration from sugarcane faces two major challenges: the presence of several minor components and the nonlinear behavior of vapor-liquid equilibrium (VLE) calculations. This work assesses the effect of simplifications to overcome these difficulties. From a set of seventeen substances, methanol, n-propanol, isobutanol, 2-methyl-1-butanol and 3-methyl-1-butanol were selected through the examination of the influence of each minor component on vapor-liquid equilibrium calculations of ethanol-water-third component systems. The selection procedure was based on Txy diagrams built using the modified Raoult’s law. The influence of the ratio between the vapor phase fugacity coefficients and of the Poynting correction factor were verified. The accuracy of four correlations for vapor pressure was evaluated, and two functional-group activity coefficient models were scrutinized: the recent Functional-Segment Activity Coefficient (F-SAC) and the UNIFAC-Do model.


INTRODUCTION
In ethanol production from sugarcane, the concentration process usually comprehends two sets of distillation columns, which are responsible for most of the energy demand.Reboilers consume about 35% of the total heating energy of an autonomous first generation distillery (Dias et al., 2011) to concentrate an ethanolic mixture from around 0.5% to a requested specification of 92.6-93.8% in mass fraction of ethanol (Marquini et al., 2007;Dias et al., 2011;Furlan, et al., 2012).The feed stream is sugarcane wine, a multicomponent complex mixture having ethanol and water as major components.Tzeng et al. (2010) reported over 20 contaminants in sugarcane wine during fermentation at laboratorial scale, including alcohols, Brazilian Journal of Chemical Engineering esters and organic acids.Batista et al. (2012) reported compositions of the main process streams in industry, identifying 17 minor components, and a fusel oil stream with 38.7% of impurities in mass fraction.
The use of a large number of compounds in process simulations brings several difficulties.First, the thermodynamic models for VLE calculations will have an important increase in the number of required interaction parameters, even if group contribution approaches are used.Additionally, low concentrations make equilibrium calculations less robust.Finally, the more compounds in the simulation the larger becomes the size of the plant model, once again penalizing the robustness of the numeric methods.Batista et al. (2012) performed a simulation of a typical ethanol concentration section of an industrial plant for a sugarcane wine containing 17 minor components.Deviations from industrial data were high for minor components, and parameters of the NRTL model (Non-Random Two-Liquids from Chen et al., 1982) for several pairs of substances had to be estimated with the UNIFAC-Do model (UNIQUAC Functional-group Activity Coefficients -Dortmund), using parameters of functional groups from the Aspen Plus® simulator.
In order to simplify this problem, the sugarcane wine is often assumed to be an ethanol-water binary mixture when the overall process is simulated or optimized (Marquini et al., 2007;Dias et al., 2011;Furlan, et al.,2012).For simulation of ethanol dehydration this assumption is also common (Ravagnani et al., 2010;Figueiredo et al., 2011;Benyahia et al., 2014;Dai et al., 2014;Soares et al., 2015).
Considering the system as an ethanol-water binary mixture may not have a substantial impact for the ethanol dehydration process, since its feed stream is hydrous ethanol close to the azeotropic point with less than 0.004 mass fraction of contaminants (Batista et al., 2012).However, in the ethanol concentration step, the composition of minor components through the distillation columns may be significant.Thus, a test to select which substances are essential for a reliable simulation of the process is demanded (Brignole and Pereda, 2013).
Moreover, the vapor-liquid equilibrium in each distillation stage may be sensitive to the presence of minor components, and the successful design and operation of a distillation column for complex systems heavily relies on thermodynamic models (Poling et al., 2000).Valderrama et al. (2012) presented a review about VLE thermodynamic modeling for the alcoholic food industry.The work analyzed binary and ternary systems and concluded that better results were obtained using Raoult's modified law rather than equations of state.The NRTL model was more accurate for the systems with available experimental data for estimation of the parameters.The UNIFAC model was indicated for cases without experimental data.
Group-contributions models such as the classical UNIFAC from Fredenslund et al. (1977) or its modified form, the UNIFAC-Do from Weidlich and Gmehling (1987), are indicated for multicomponent mixtures, thus overcoming the requirement of a large set of experimental data for estimation of the model parameters.Since the number of functional groups is much smaller than the number of possible molecules, the quantity of required experimental data can be sharply decreased.Recently, Soares and Gerber (2013) created a new group-contribution model, F-SAC (Functional-Segment Activity Coefficient), combining the theoretical characteristics of the interaction energy between charge segments of COSMO-SAC (Conductor-like Screening Model -Segment Activity Coefficient) from Lin and Sandler (2002) and the empirical approach of the UNIFAC.
The present work shows a simple methodology for selection of components to be taken into account in the feed stream of the ethanol concentration process.Further, a systematic assessment is accomplished, to verify the sensitivity of the VLE calculations with respect to usual simplifying assumptions, vapor pressure correlations and thermodynamic models for activity coefficients, highlighting the performance of the F-SAC model.

METHODOLOGY
Vapor-liquid equilibrium calculation Valderrama et al. (2012) used the gamma-phi approach (Poling et al., 2000) for their equilibrium calculations to study alcoholic distillations of musts made from fermented grapes. (1) The exponential term is the Poynting correction factor, here calculated using the DIPPR equation (Design Institute for Physical Properties) for molar volume.The fugacity coefficients ratio was estimated here using the Soave-Redlich-Kwong equation of state, with its original mixing rule (Soave, 1972).However, Vapor-liquid equilibrium calculation for simulation of bioethanol concentration from sugarcane for the initial stage of selection of minor components, the term i ℑ was considered equal to one and equa- tions 1 and 2 reduce to the modified Raoult's law In the present work, four models for vapor pressure were evaluated: the Antoine equation (Antoine, 1888, cited by Thomson, 1946), the extended Antoine equation (from the Design Institute for Physical Properties -DIPPR), the Wagner equation (Wagner, 1973) and the correlation of Chemical Engineering and Materials Research Information Center (CHERIC).

Functional-Segment Activity Coefficient (F-SAC) model
The activity coefficient γ i for the molecule i in solution may be computed as the sum of a combinatorial or entropic term plus a residual or enthalpic term: The F-SAC model employs the combinatorial term of the UNIFAC-Do model (Weidlich and Gmehling, 1987;and Gmehling et al., 1993).The residual contribution comes from the difference between the free energy for charge restoration of a solute molecule i in solution and for charge restoration in a pure liquid i (Gerber and Soares, 2010): where the Gibbs free energy to restore the charge of the solute molecule i in solution is the summation of (5) The activity coefficient of each charge segment σ m is calculated iteratively: where the energetic misfit constant αʹ is 8544.6 kcal.Å 4 /mol and the standard contact radius is assumed as 1.07 Å.The term E HB accounts for hydrogen bond effects and can be estimated from experimental data.
The σ-profile for each molecule is a summation of the property of the functional groups and these from several charge segments as follows: , then the σ-profile for each functional group is determined by three parameters Q k + (absolute area with positive charge), Q k − (absolute area with negative charge) and σ k + (charge density of positive segment).

Simulation and assessment
Equilibrium calculations were performed through the algorithm introduced by Luyben ( 2007) with temperature updating by the Newton-Raphson method.The program was written in C++ language and run in Codeblocks, a cross-platform IDE.The parameters of the Antoine equation and CHERIC correlation were provided by Poling et al. (2000) and by the CHERIC website, F-SAC parameters by Soares and Gerber (2013) and Soares et al. (2013) and the other parameters were taken from the APV82 PURE28 databank of Aspen Plus ® .The sources of experimental data were reported throughout the work.
The comparative analysis was based on the variation between the calculated and experimental values of the property ξ: Absolute deviation and relative absolute deviations were also used:

Minor compounds selection
The minor compounds selection procedure started by defining an initial set of substances.The set was chosen based on the data obtained by Batista et al. (2012), since they represent the actual composition of an industrial process.The highest concentration of each compound, among all the streams of the distillation process, was considered (Table 1).
Next, Txy diagrams were calculated by Raoult's modified law for each ethanol-water-third component (fixed at its highest mole fraction value indicated in Table1) system, and compared with the calculated curve and experimental data of the ethanol-water binary system.The simulations and the experimental data are at 1.013 bar.
The third components that displayed a significant detachment between the binary and ternary systems were selected.As shown in Figure 1, n-propanol, isobutanol, 3-methyl-1-butanol and 2-methyl-1-butanol complied with this condition.Methanol does not follow the condition, but it was also chosen because it is more volatile than ethanol and it is the main contaminant in the hydrous ethanol stream.Although Hayden O'Connell's model for organic acids and Henry's law for uncondensed gases would be more appropriate to the vapor-liquid equilibrium calculations for acetic and propionic acids and gas carbonic systems, these substances did not have a noticeable influence.

Assesment of VLE calculation
The pressure in the real process is not atmospheric.Through the columns, the pressure drops from bottom to top.According to Batista et al. (2012), the pressure ranges from 0.9 to 1.6 bar.Temperature varies from 350 to 390 K.These ranges were used to analyze each term of equations 1 and 2, to validate simplifications of the modified Raoult's law, and to compare models.

Vapor pressure
Table 2 presents vapor pressures calculated by the extended Antoine equation for the minimum and maximum temperature values.The Table also shows the average values of absolute deviations between the experimental and calculated vapor pressures for the four evaluated models.Average values near to zero must be considered with caution, because they may indicate that the used parameters were estimated from the same data set.If that is the case, an unwanted bias will be present, because the literature does not always clarify which experimental data were used in its estimation.Nevertheless, the extended Antoine equation and Wagner equation were more accurate in spite of significant deviations for the less volatile substances, a trend that was noticed for all equations.

Poynting correction factor
Figure 2 shows a similar behavior of the Poynting correction factor for both minimum and maximum pressures and for all substances.The highest deviations from unity were -0.0007 for methanol at 390 K and 0.9 bar, and around 0.006 for 2-methyl-1-butanol and 3-methyl-1-butanol at 350 K at 1.6 bar.In the distillation column, pressure and temperature increase from bottom to top.That means a better representation of the concentration process by the left side of

Ratio between fugacity coefficients
The fugacity coefficients of the saturated pure substances were calculated using the Soave-Redlich-Kwong cubic equation of state in the considered range of temperature.They deviated from unity when the temperature increased (Figure 3).The farthest values were for methanol and ethanol components, but they are more concentrated in superior stages of the columns, where temperatures are lower.In Figure 4, three substances (water and two whose fugacity coefficient of the saturated pure substance are farther from one) were chosen to illustrate the behavior of the ratio between the fugacity coefficients of the pure substance and of the component in the mixture (first term of Equation 2).This ratio ranges between 0.965 at 373 K for water in the water-methanol system (Figure 4.c) and its highest value is equal to 1.01 at 352 K for ethanol in the methanol-ethanol system (Figure 4.b).
Figure 4.a and the left sides of Figure 4.b and 4.c (with lower temperatures) represent the ethanol concentration process because methanol has a similar profile as ethanol within the columns.Therefore, the ratio between the fugacity coefficients within the real process will be closer to unity than the value shown in Figure 4.

Activity coefficient
In the last two sections, the two terms of Equation 2were separately analyzed.The general result of the multiplication of them ( i ℑ ) is a deviation in the second decimal place (above unity for less volatile components and below unity for more volatile components).The highest deviations occur when the compound is more dilute, in other words, when the activity coefficient of the component reaches its maximum.In consequence, at dilute concentrations, the term i ℑ has a much lower weight in equation ( 1) than the activity coefficient, and it may be considered equal to one.Thus, one follows in the modified Raoult's law.Nine binary systems were selected (Table 3) to test this statement.Calculations of the bubble pressure were carried by the modified Raoult's law using two models of activity coefficient based on group contribution methods.The results were compared to experimental data with variable temperature and pressure, whose thermodynamic consistency was reported by the reference.Table 3 also presents deviations obtained by other authors using other two models for activity coefficient (UNIFAC -UNIQUAC Functional-group Activity Coefficients; and NRTL -Non Random Two Liquids).
These results were similar to those obtained by Valderrama et al. (2012) for ethanol-water-third component ternary systems using three models based on local composition.Higher deviations were observed for more complex molecules.
Although the studied systems of other references might be different from the ones in this work, according to Table 3, the F-SAC model revealed to be promising, with a better representation of polar substance interactions.
The limitation of group contribution methods is the inability of discerning isomers such as 2-methyl-1--butanol and 3-methyl-1-butanol.On the other hand, bad fits were attained for simple systems like ethanol--methanol calculated by UNIFAC-Do and water-isobutanol by both models.
Figure 5 illustrates the calculations for the binary system of greatest interest (ethanol-water) at two temperatures.The two models estimate well the equilibrium.
Figure 6 explores the sources of the deviations.Deviations in the pressure estimation show to be more influenced by temperature while deviations in the mole fraction of vapor phase are more influenced by the mole fraction of the liquid phase.This behavior indicates a deficiency or lack of fit of the models at high temperatures and high dilutions.

CONCLUSION
In this paper, five substances (methanol, n-propanol, isobutanol, 2-methyl-1-butanol and 3-methyl-1--butanol) were selected as most important to be considered in VLE calculations, from a set of seventeen minor components present in the ethanol from the sugarcane concentration process.The used criterion was the influence of each component in the vapor-liquid equilibrium of ethanol-water-third compound system, when compared to the ethanol-water binary system.
For the elected substances, assuming typical industrial process conditions of pressure and temperature, two out of four vapor pressure correlations were indicated to be applicable based on the average absolute deviation (extended Antoine and Wagner).The Poynting correction factor and the ratio between the fugacity coefficients in the gamma-phi approach could be assumed to be equal to one.Furthermore, two models based on the group-contribution concept for activity calculation were investigated.The results reveal a slightly superior performance of the F-SAC model over

Figure 1 :
Figure 1: T-x-y diagrams for ethanol-water-selected minor component at 1.013 bar: experimental data of ethanol-water binary system at 1.013 bar from Kurihara et al. (1993) (○), calculated curves without (-) and with presence of minor component (-).

Figure 2 .
Figure2shows a similar behavior of the Poynting correction factor for both minimum and maximum pressures and for all substances.The highest deviations from unity were -0.0007 for methanol at 390 K and 0.9 bar, and around 0.006 for 2-methyl-1-butanol and 3-methyl-1-butanol at 350 K at 1.6 bar.In the distillation column, pressure and temperature increase from bottom to top.That means a better representation of the concentration process by the left side of Figure 2.a and the right side of Figure 2.b, in other words, where Poynting correction factors are close to one.

Figure 2 :
Figure 2: Poynting correction factor with molar volume by DIPPR's equation and vapor pressure by the extended Antoine equation for the considered temperature range and minimum (a) and maximum (b) pressures.

Figure 3 :
Figure 3: Fugacity coefficients of the saturated pure substance calculated using the Soave-Redlich-Kwong equation of state with vapor pressures calculated by the extended Antoine equation.

Figure 4 :
Figure 4: Fugacity coefficient of saturated pure substance ϕ i sat , fugacity coefficient of the component in the mixture in the vapor phase V i φ ˆ and the ratio

Figure 5 :
Figure 5: Vapor-liquid equilibrium of the isothermal ethanol-water binary system at two temperatures: experimental data from Gmehling et al. (1982) (○) and calculated curves by UNIFAC-Do model (-) and F-SAC model (---) with vapor pressure calculated by extended Antoine equation.

Figure 6 :
Figure 6: Deviations in: (a) pressure (bar) as a function of temperature (K), and (b) mole fraction of ethanol in the vapor phase as a function of the mole fraction of ethanol in the liquid phase for the ethanol-water system*.

Table 1 :
Highest concentrations (in mass and mole fractions) of each minor component and their stream, from experimental data obtained by Batista et al. (2012).

Table 2 :
Vapor pressures (bar) calculated by the extended Antoine equation and averages of absolute deviations (bar) between experimental data and calculated values for the four considered correlations.