New Fitting Method for Vapour-Liquid Equilibrium Data and its Application for Ethanol-Water Distillation Process Modeling

For modeling and simulation of distillation processes, the main information required is the vapour-liquid equilibrium (VLE) data. These are available as discrete values, but for modeling purpose, an analytical function would be more suitable. We developed a simple calculation protocol to obtain a single continuous function for the VLE curve on the entire concentration domain, without using thermodynamic functions. The fitting parameters of binary ethanol-water mixture VLE curve were determined. The methodology is suitable for fitting continuous functions on VLE data of any binary or multicomponent, non-ideal (even azeotropic) mixtures.


Introduction
The modeling of separation processes in chemical engineering is indispensable for the modern design of chemical equipment. The process simulation and modeling in the industry are limited many times by the lack of thermodynamic knowledge and/or data. Poor thermophysical property or equilibrium data lead to significant calculation errors and conduct to erroneous equipment design [1]. Any improvements in this direction increase the accuracy of the predictions. All design procedures of separation processes are based on mass balances that can be calculated from equilibrium data. In particular, for distillation, these are the vapour-liquid equilibrium curves. These data are available in tabular form from databases, as discrete values, or when needed experimentally determined.ation is needed. To obtain continuous curves, a fitting method is needed, but in many cases (mainly for highly non-ideal, azeotropic mixtures), which also presents a concave region, this could be a challenging task. For this purpose, equations of state [2] or Gibbs excess energy method with different activity coefficient (NRTL, Wilson, UNIQUAC) or thermodynamic models (Margules or van Laar equations) are commonly used [3].
Distillation is undoubtedly one of the most important separation unit operation in the chemical processes [4,5], and alongside with the solid-liquid extraction even in the food industry. The distillation of fermented mashes for the production of food degree or bio-ethanol and for obtaining the distilled spirits with high alcohol content [6,7] is the most important industrial distillation process from the agri-food sector. One of the reasons for choosing the water-ethanol system as subject to our study was the practical importance of the mixture, another, that this binary system is highly non-ideal. The behavior of azeotropic mixtures, with inflexion point and a concave region, makes the choice of fitting function more challenging. In this condition, inasmuch as a reliable fitting method could be developed, this will likely be applicable for further difficult systems too.
Our method is based on the concept of relative volatility (α) of the two components that connect the molar fraction of the volatile component in the vapour (y) and the liquid phase (x). In simple models, the relative volatility is considered constant, but this, even in the case of ideal mixtures is an oversimplification, giving considerable errors. In technical batch distillation calculation is usual practice to consider the relative volatility as constant, their value being the average of initial and final values or equal to the value at the average temperature [8].
A better fitting with wider applicability is obtained by using the inverse linear function (eq. 1), mainly for the temperature dependence of relative volatility of near-ideal mixtures [9]. A more general approximation for VLE is the generalized constant relative volatility model proposed by Gmehling (eq. 2) [10]. This is considered the best of the mentioned models, but with some limitations, namely, their validity is restrained only to minimum boiling azeotropic binary mixtures [11].
As the mass balance equation of batch distillation (Rayleigh integral) contains the difference of the volatile component concentration in the two phases (the driving force of the process), the simplicity of the relative volatility function increase the possibility that the Rayleigh formula (eq. 3) is analytically integrable.
Based on this information, we try to find simple functions for relative volatility with a minimal error that makes the value of Rayleigh integral (IR) more precise or even ensure their analytical integrability.

Materials and methods
As our work is related to process modeling, it was made on a computer, using dedicated software. For data analysis and graphical presentation, we used Statistica 8.0 (StatSoft, Inc., Oklahoma, USA). The VLE data for ethanol-water binary mixtures at atmospheric pressure was taken from Banu (2006) [6], since this dataset contains the largest number of values. To find an appropriate fitting function for the relative volatility (RV), firstly, we obtained their values from the VLE data, using the defining equation (eq. 4) and making their graphical representation in function of ethanol molar concentration in the liquid phase (x). Based on the graphical representation, the possible adequate function was selected.
The second strategy was to try to choose a transformation for the denominator of the IR (the driving force), for the resulting fitted form to become analytically integrable. However, to achieve this goal is still necessary for the appropriate choice of the fitting function. For this reason, besides the relative volatility (), two new functions (and) were defined (eq. 5 and eq. 6): when for the -function a rational polynomial with unit numerator (reverse polynomial function) (eq. 7) are set, can be easily demonstrated, that the IR is integrable (eq. 8), where the xF and xW represent the initial batch composition and the waste composition, respectively.
The differences between the observed data and calculated values are expressed as relative errors (eq. 9): obs pred obs yy e y   (9) where: erelative error (deviation), yobsthe observed data, ypredpredicted value.
For the fitting procedure, we used two nonlinear estimation mode, namely the least square (with Gauss-Newton estimation) and custom loss function (minimization of the square of the difference of the observed and predicted values).

Results and discussions
The selection of an adequate fitting functions is a somewhat challenging task, therefore no universal model for VLE of non-ideal mixtures exists. Taken into account the above considerations, we have the following results. Primarily, for fitting the VLE curve to data of ethanol-water mixtures at atmospheric pressure, we tried to find a probe function. Firstly we used the function proposed by Said and Al-Ameeri (1987) [8] for calculations in chemical engineering, especially for distillation processes (eq. 10):  As can be seen in Figure 1, the function does not fit well on the data (a large difference between the observed and predicted values), but the trend (asymptotic decrease) is correct. The predicted values are mainly below the observed data, especially at lower x values. The reason for this behavior is given by the presence of x in the numerator, which decreases the function value at lower molar concentations. To compensate this effect, we probe, functions with unitary numerators. The tested functions, the resulted fitting parameters, and the correlation coefficients (R-values) are summarized in Table 1. For all the proposed functions, the fitting is tight, the value of the correlation coefficients are very high (closed to 1). Based on the R-values one cannot be evaluating the best function, therefore we represented the relative deviation (error) between the observed and predicted values for all the three functions on the same diagram (Figure 2.). In the large interval of x, the relative deviation remains below 5% (0.05). Near the azeotropic point, ( 0.9) x  the deviation increase drastically. The lowest value of relative deviation can be observed for the α3-function, named in the following as power-law function.
The observed and the predicted value of the relative volatility by the power-law function are given in Figure 3. https://doi.org /10.37358/Rev. Chim.1949  Based on the relative deviation values (Figure 2., curve eα3), it can be concluded that our proposed power-law function is adequate as an empirical model for the relative volatility curve. Thereafter, the VLE values given by the three empirical models, namely the reverse linear [8], the Gmehling [10][11][12], and the power-law functions were compared (Figure 4). The parameters of the three empirical models are given in Table 2. and were obtained by function fitting for the reverse linear and power-law function from the literature [12]. Moreover, the Gmehling model is somehow different from the other two models, i.e., use constant relative volatility (α, RV) and give directly the VLE curve (y=f(x, α)).  Still, to compare the fitting quality of the models, we used as previously the relative deviation values represented graphically in Figure 5. The best fitting is given by the power-law model proposed by us, furthermore, this model has three fitting parameters, too. The second strategy was to choose an adequate model for the VLE curve to obtain an analytically integrable form of Rayleigh equation. As demonstrated in eq. 7., that choosing a reverse polynomial form for the ( ) / y x x  fraction (hereinafter referred to as β-function, in conformity with eq. 5.), fulfills the proposed requirement. However, to achieve a tight fitting, it is still necessary the appropriate choice of the polynomial degree. Besides the relative volatility (α), two new functions (β and γ) were defined (eq. 5 and eq. 6).
By gradual increase the polynomial degree, it was demonstrated that the best fitting is given by the 5 th degree polynomials. The further increase of the degree leads to oscillatory behavior of the fitting function, moreover, no further improvement in the relative deviation was observed (results not shown). https://doi.org /10.37358/Rev. Chim.1949 Rev. Chim., 71 (7) In the denominator of the argument of the Rayleigh integral (eq. 3), the difference of the phase composition appears, which is represented in Figure 6. It would be very expedient to be able to select the fitting function based on the visual examination of Figure 6. Despite that it can be established that the function has a probability distribution shape, we were unable to define a proper function type that fits well to these values. To simplify the problem, we introduced the intermediary functions: the fractions of β and γ (eq. 5 and 6). After that, we tried to find fitting functions for them, which can be used further to calculate the VLE curve. This procedure basically is the same that the method using α-function, but the equilibrium value of the y has more simplified form (eq. 11-13): The value of the fitting parameters for all three intermediary functions, as well as for the direct fitting on the VLE curve (ydir) are summarized in Table 3. From the equations on the three intermediary functions (eq. 11-13) could be calculate the fitted models of the VLE curves using the fitted parameters from the table.3, In addition, using eq. 14, the direct fitting function for the same curve can be obtained (Figure 8). Visually, all the predicted values fit tight with the observed VLE values taken from the literature [6].
In order to compare the fitting quality of the four empirical models, we proceed similarly as previously, i.e., we represent graphically on the same diagram ( Figure 9) the empirical model deviation from the tabulated data, and evaluate the graphic. The evaluation of the VLE empirical models shows that at a very low liquid molar fraction (x<0.01) all the models give very poor results, but from an engineering point of view, this domain does not present any real interest. Strictly, taking into account only the value of the deviation in the whole concentration domain, the ranking is as follows: the best model is the ydir function, the second is the traditional yα, the third is the yγ, and the worst result is offered by the yβ model. This is not very encouraging for using the yβ model in the batch distillation calculus, but the situation is not so bad as it seems to be since the deviation problems appear only at higher molar fractions (about x>0.85). Below, x<0.85, the relative deviation band is narrow (± 2%), suitable for precise engineering estimations. https://doi.org /10.37358/Rev. Chim.1949 Moreover, the zone above x>0.85, do not present practical importance as well as the initial concentration, xF, of the batch, very rarely exceeds the value of x= 0.2 (considering the fermentative technology). This study shows that the yα-based calculus are eligible. The direct fitting models can give excellent results, but in general, is difficult to guess directly the corresponding function form.
Finally, we compared to the best results obtained from the two strategies, namely the power-law function obtained in the first part and the results obtained with the direct fitting in the second part. The results are presented in Figure 10. Though the analysis of the figure, it can be observed that at high concentrations and after the azeotropic point the direct fitting gives a better result in comparison with the power-law function. It can be concluded, that the best model among the studied ones was the ydir based direct model, using as fitting function the 5 th degree reverse polynomial. The proposed power-law function can be used if the concentration in the liquid phase does not exceed x=0.80 (the relative deviation, in this case, is under 4%). The advantage of the power-law function consists in their simplicity, forasmuch had only 3 fitting parameters, as well as their mathematical handling could be more obvious, does not present major difficulties.
The difference of the molar fractions   yx  in the two phases may be expressed in terms of a power series (eq. 15), named the Ocon active fraction (ZT) [13][14][15].
Beyond mole related quantities, the Ocon active fraction contains another constant RV, which should be determined by fitting [14], or from thermodynamic considerations [15]. We propose a similar, although a slightly different, more direct method, which makes possible a relatively straightforward way for the fitting procedure.
We also tried to achieve the integrability of the Rayleigh integral, that criteria is not fulfilled by the power series method proposed. However, the advantage of the above method is their applicability in the fitting of other thermodynamic functions and properties [13]. https://doi.org/10.37358/RC.20.7.8219

Conclusions
The modeling of innovative distillation methods (cyclic distillation) with general nonlinear VLE behavior requires more accurate non-iterative design algorithms [16], and the presented fast, and accurate method could be a part of the modeling arsenal. Using our method, even the most recent programming algorithm based on the classical McCabe-Thiele graphical method [17] could be improved and simplified. This would be possible as the above-mentioned algorithm for VLE estimation uses continuous piecewise linear approximation, with inherent inaccuracy, when the number or position of the data points are chosen inadequately. Even using professional modeling software packages, like Aspen Hysys and Pro/II great difficulty could meet to obtain a convergent solution when the relative volatility presents large variations and the VLE data are given in discrete form [18]. Due to very tight fitting of the proposed functions on the observed data, with minor deviation, the new VLE empirical models are suitable for modeling and simulation purposes, resolving the major source of a calculation error in aforementioned modeling cases.