Prediction of Excess Enthalpy Using Volume-Translated Peng–Robinson Equation of State

In this work volume-translated Peng–Robinson group contribution equation of state was used to calculate excess enthalpies. Four model systems were selected with the purpose to compare experimental and predicted enthalpy values at different temperatures. After the calculations were performed in Matlab software, results were verified with free software tool of Dortmunder Datenbank (DDB). In a next step, the mixing process and interaction forces were described on the basis of the sign and course of enthalpy values. The endothermic behavior of three systems could be well predicted, while for the most polar system, predictions were less precise. Furthermore, the discrepancy between experimental data from the literature and predicted values was discussed to evaluate the accuracy of the selected model. Lowest mean deviations (<75 J/mol and <15% at all temperatures) could be stated for alkane/benzene mixtures, while highest deviations could be again observed for the most polar mixture. Although the magnitude of deviations was in agreement with the literature, it could be shown that the selected temperature is of major importance for the quality of predictions. Furthermore, a review of different literature values for the n-hexane/benzene system could reveal that the reliability of experimental data has to be carefully checked.


Introduction
In order to design thermal separation processes, the knowledge of fundamental thermodynamic properties as well as equilibria conditions is necessary. Even if only binary mixtures are regarded, it is laborious to collect experimental data for every possible combination [1]. Therefore, predictive models are required that allow the description of binary and multiple component systems based on few adjustable parameters. One of the most used concepts is the so-called group contribution concept, which assumes that every mixture can be sufficiently described if the interaction between functional groups of the molecules inside the mixture are considered [2].
Originally, the group contribution concept was introduced for activity coefficient models in 1969 [3]. Shortly after the ASOG GC model, the UNIFAC model was presented and in the next years developed further by the University of Dortmund [3,4]. While in the first years a focus was set on a good description of vapor liquid equilibria, the model was continuously improved in order to describe excess enthalpies and activity coefficients with a satisfying accuracy as well [3]. However, the treatment of supercritical compounds and the prediction of other thermodynamic properties (gas solubilities, enthalpies, densities, heat capacities) is not possible with activity coefficient models [2,3,5]. For this purpose, the group contribution concept was extended to equations of state, and in 1991 the predictive Soave-Redlich-Kwong group contribution equation of state was presented [6]. Ten years later the volume-translated Peng-Robinson GC EOS was introduced by [7], mainly to address the weaknesses of PSRK EOS [5].

•
One weakness of SRK EOS is that often inaccurate results are obtained for liquid densities. Therefore, the same phenomenon can be stated for PSRK EOS, too. • Unsatisfactory results can be observed in many cases for activity coefficients at infinite dilution, heats of mixing, and for strongly asymmetric systems. This can be explained by the fact that the UNIFAC model, used in the g E mixing rule, is not able to predict those properties accurately.
In this article a focus will be set on the prediction of excess enthalpies. In general, excess properties are important tools to describe deviations from ideal behavior. In the case of liquid mixtures, they are defined as the difference between real and ideal solutions. Those differences can be directly attributed to the fact that interactions between the molecules in a mixture differ from those observed for the pure components. Excess enthalpies, also known as heats of mixing, are particularly suited to illustrate those interactions. A positive sign of excess enthalpy means that the mixing process is an endothermic event. The interactions between molecules in the mixture are therefore weaker compared with those which can be stated for pure components. In contrast, a negative value of excess enthalpy indicates an exothermic mixing process and stronger interactions. For some mixtures even both phenomena can be observed. The discriminating factor in those cases is the molar ratio of the components in the mixture [8,9].
Beside the fact that excess enthalpies give important hints for the caloric behavior of the mixing process, it is important to note that those values describe the temperature dependence of activity coefficients. Both properties are connected via the Gibbs-Helmholtz equation. Experimental data of excess enthalpies at different temperatures are important information used during the fitting procedure of parameters for group contribution models [3]. Nevertheless, investigations on how well group contribution models can predict excess enthalpies at different temperatures can rarely be found. A short overview about results reported in the literature for VTPR EOS is given in Table 1. On the basis of four selected model systems, the present paper gives a more detailed quantitative insight on how good predictions with VTPR EOS are.  [14] acetone/n-heptane small deviations at lowest temperature (283 K); bigger deviations at highest temperature (343 K)

Calculation of Excess Enthalpy Using VTPR EOS
Before the results will be presented, the equations that were used to calculate excess enthalpies are given within this section. The formula of VTPR EOS first mentioned by [7] can be seen in Equation (1). It is a so-called cubic equation of state, which connects the state variables temperature (T), pressure (P), and (molar) volume (v). The symbol R stands for the ideal gas constant.
The modifications in comparison with PSRK EOS can be seen if the three coefficients a, b, and c are regarded.

•
For the calculation of parameter a, which accounts for attractive forces, a modified α function (so-called TWU-α function) is used (pure component i; see Equations (2) and (3). A proper description of the sub-and supercritical range is possible with the help of adjustable parameters L, M, N, and the usage of critical temperature (T crit ) and pressure (P crit ). Compared with Mathias-Copeman-α function, used in PSRK EOS, results are more reliable at high reduced temperatures [2,3]. The good applicability is also shown in recent publications; the authors of [15] could successfully model catalytic hydrogenation processes by modification of original Peng-Robinson EOS with TWU-α function. For binary mixtures a new g E mixing rule leads to improved prediction accuracies (Formula (4)) [1,5,10]. The parameter z stands for the mole ratio of component i.
• In order to get a better description for asymmetric systems parameter b, which accounts for repulsive forces, is not calculated by a linear mixing rule (binary mixture). Instead, a quadratic mixing rule as shown in Equation (5) • With the introduction of parameter c a volume translation is implemented. Together with the usage of PR instead of SRK EOS, this leads to improved results for liquid densities. The volume translation parameter c for pure component i can be calculated on the basis of the difference between experimental density v exp and the density calculated with PR EOS (v PR ) at the reference temperature T re f as shown in Equation (7) (temperature-independent volume translation) [1,7]. For all the compounds mentioned in this article, c values published by [1] were used. If binary mixtures are regarded, a simple linear mixing rule can be used (Equation (8)) [10].
For the application of the g E mixing rule (parameter a) the Gibbs energy (g E ) has to be derived. The residual part of UNIFAC model was proposed as suitable GC model [10]. Formulas (9)- (11) show the calculation based on the residual part of the activity coefficient of component i (γ res,i ) and group activity coefficients Γ for each individual subgroup k in the mixture (Γ k ) or for the pure component i (Γ [4,5]. The number of the subgroup k in the compound i is considered with the parameter v To achieve better results over a large temperature and pressure range similar to the modified UNIFAC model (Dortmund), three energy parameters a m,n , b m,n , and c m,n are adjusted to account for temperature-dependent interactions (often described as difference U m,n − U n,m ) between different functional groups m and n (Formula (12)) [1,4]. Parameter Ψ is finally obtained.
Furthermore, the calculation of surface area fraction θ and group mole fraction X is given in Equations (13) and (14). The relative van der Waals group surface areas Q for the individual subgroups were either adopted from the PSRK model or adjusted for the VTPR model [5,16].
Excess enthalpies (h E ) can then be calculated with equations of state using the so-called residual enthalpies (h − h id ) of mixture and pure components (Formula (15)) [5,17].
The residual enthalpy itself can be derived on the basis of mixing parameters for pure components (a, b, c) and the mixture as shown in Equation (16) [5]. In this formula z describes the compressibility factor z = P·v R·T .
Furthermore, the temperature derivative of parameter a must be calculated. This was done in agreement with the descriptions of [5,18]. First of all, on the basis of the g E mixing rule Equation (17) can be derived [5].
Parameter a i describes the parameter a for pure component i. Therefore, the derivative can be expressed as shown in Formulas (18) and (19). The derivative of the g E part can expressed with the help of Formula (20) [5].
Finally, the derivative of the activity coefficient can be calculated in agreement with the descriptions of [18] (Formulas (21)-(25)). Two constants C and D were introduced in order to simplify the program code in Matlab.

Results
To obtain the results for excess enthalpy, a program code was written in Matlab analysis software, which is built of four important parts: 1. The main m-file that executes the program contains the basic information needed for the calculations (name of compounds, temperature, pressure) as well as experimental data of the system. The first mentioned information is passed to the excess enthalpy function. 2. Calculations of excess enthalpy are performed in a function stored as separate file.
It interacts not only with the main m-file but also with other functions described in the next point. 3. In order to improve the clarity of the Matlab program, additional functions were written to read necessary data out of an Excel file, to calculate mixing parameters for pure components and the mixture, to carry out calculations of the g E part, and to differentiate the g E part with respect to temperature.
4. The required pure component properties and the parameter for VTPR EOS were saved in an Excel file, which allows to use an extension of the program for other components in future. The data were retrieved from [1,12,13,16,19].
An illustration of the deviations between predicted results and experimental data is given in Figure 1. First of all, it is remarkable that three of four mixtures showed positive values of excess enthalpy, which means respective mixing processes are endothermic [8]. In all cases this behavior was predicted by VTPR EOS, too. Only for the mixture acetone/water heat was released during the mixing process (exothermic behavior) if the molar ratio of acetone was low (<0.3-0.5) [8]. Again, this course of enthalpy values was predicted by VTPR EOS. However, the molar ratio at which the excess enthalpy changed the sign was shifted.
In addition, it can be seen that the system n-hexadecane/benzene showed the biggest values of excess enthalpy, which indicates that only weak interaction forces between the molecules can be found in the mixture (compared with pure components). Furthermore, the asymmetric curvature of the system 1-hexene/n-BuOH indicates that interaction forces between 1-hexene molecules clearly differ from those of n-BuOH molecules [8]. The opposite can be stated for the system n-hexane/benzene which gives an almost symmetric curvature course.  [21,[25][26][27][28][29] Finally, from this figure it is visible that the accuracy of predicted results depends on the temperature, which is regarded as follows: • For the model system n-hexadecane/benzene ( Figure 1a) predictions with VTPR EOS are nearly independent of temperature, although experimental data clearly differ at the selected temperatures. At low molar ratios of benzene a good fit can be stated when experimental data at 323.15 K are regarded, while at higher molar ratios the fit is better for a temperature of 298.15 K. • For the model system 1-hexene/n-BuOH (Figure 1b) predictions are better for a temperature of 313.15 K. Although this time the predicted results show a temperature dependence at 298.15 K, a large gap between predicted and experimental data has to be stated. • For the model system acetone/water (Figure 1c) deviations between experimental data and predicted values are rather big at all temperatures. Again, a temperature dependence can be stated for experimental as well as predicted values. • For the model system n-hexane/benzene (Figure 1d) predicted results are nearly independent of temperature. Predictions are best for a temperature of 293.15 K. However, it is visible that experimental data itself show big fluctuations, especially at 298.15 K. This circumstance will be further discussed in the next section.

Discussion
Before entering a more profound discussion it seems worth to mention that the results obtained with Matlab software were compared with the predictions of DDB software tool. The values for absolute deviations of the selected model systems at different temperatures are given in Figure 2. Arithmetic mean average deviations (AMADs) were calculated as shown in Equation (26) using the values calculated with Matlab (h pred ), those given by DDB software tool (h exp ), and the total number of datapoints N. From Figure 2 the conclusion can be drawn that in all cases deviations were negligibly small (<0.05 J/mol in all cases), and therefore the Matlab code was accurate enough to return valid results.
(a) (arithmetic) mean absolute deviation in J/mol (b) maximum absolute deviation in J/mol In a next step relative and absolute errors were calculated on the basis of experimental data and predicted values. Mean absolute values were calculated as shown in Formula (26), and arithmetic mean relative deviations (AMRD) were computed according to Equation (27) using results from Matlab (h pred ), experimental data (h exp ), and the total number of datapoints N.
In the last section it was already mentioned that the system n-hexane/benzene showed the biggest fluctuations in experimental data. Therefore, this system will be discussed at first. In Figure 3 experimental data from different sources measured at 25°C are shown. On the basis of those data the following conclusions can be drawn: • Highest values for heats of mixing are measured by the newest publications [28,29]. The authors used the same method (Calvet microcalorimeter) for their measurements and state an inaccuracy of excess enthalpy measurements below 1%. • Slightly lower values for excess enthalpy are measured by Jones et al. [27]. A Brass calorimeter was used. Errors of this apparature were estimated to be below 2% for aqueous/alcohol solutions [30]. Jones et al. [27] compared their results graphically with those obtained by Prigogine et al. [31] and Schnaible et al. [32]. Good agreement between the three publications could be stated at molar ratios of benzene above 0.5, while at lower molar ratios results of Schnaible et al. [32] and Progogine et al. [31] were higher. Finally, excess enthalpies measured by Romani et al. [33] were in the same range as the results of Jones et al. [27], following descriptions of Mato et al. [29]. • Significantly lower values were measured by Ridgway et al. [26] and Baluja Santos [34]. The authors of the first mentioned publication measured their data in a Dewar vessel and state an accuracy around 1% [26]. According to the graphical presentation of Mato et al. [29], Andrade et al. [35] also obtained maximum values for excess enthalpy, which are comparable to those of Ridgway et al. [26].
The predictions of VTPR EOS show low deviations to most of the experimental data. This can be quantitatively confirmed with the help of mean absolute and relative errors, which are given in the following Figure 4.  Except for the data of Baluja Santos [34], absolute mean deviations were below 55 J/mol and relative mean errors were below 10%. Together with the other literature data mentioned in the bullet points above, the results of Baluja Santos [34] seem to be questionable. However, for a final evaluation, a reliability check would be necessary. A proper description how this can be done mathematically is seldom described in the literature. Schmid [1] briefly stated that a polynomial function might be a useful tool, but no details are given. Another possibility might be an evaluation of the measurement aperture and its error sources, which lies beyond the scope of this article. Nevertheless, the big differences illustrate the importance of verifying that reliable experimental data are chosen in order to obtain accurate model predictions.
At this point the absolute and relative deviations of the three other systems should be introduced. Those are shown, together with deviations of n-hexane/benzene system at other temperatures, in the two bar charts in Figure 5. From these charts the following observations can be made: • Again, it can be seen that the accuracy of predictions was temperature-dependent. Especially for the system 1-hexene/n-BuOH (absolute error 2.5 times higher at 298.15 K compared with the value at 313.15 K) the difference was pronounced for the temperatures which were investigated. However, it has to be noted that the deviations at 313.15 K were low (8.2% or 38 J/mol). • For the system n-hexane/benzene the same temperature dependence can be stated, although differences were low at all temperatures shown in Figure 5 (<10% or <25 J/mol). For measurements at 298.15 K literature data were already discussed with the help of Figure 4. If the dataset of Baluja Santos [34] is not regarded, relative as well as absolute mean errors are close to those obtained at 323.15 K or slightly higher. • The other alkane/benzene system (n-hexadecane/benzene) shows slightly higher error values (12-20% or 52-73 J/mol). For both systems it is striking that predictions of VTPR EOS were nearly temperature-independent (Figure 1a,d). An explanation can be given when the main group interaction parameters are regarded (see Table 2; main group 1-CH2 and 3-ACH are relevant). For those mixtures no energy parameter c was determined, which causes the temperature dependence of the derivative of g E to be missing (see also derivation of temperature derivation described in Section 2 and especially Formulas (24) and (25)). • Predictions exhibit the greatest deviations for the most polar system acetone/water, which is in agreement with the conclusion drawn based on Figure 1. The high relative and absolute errors (>50% and >200 J/mol) confirm that care must be taken if such a system is regarded. However, it can be said that with VTPR EOS still a tremendous improvement could be achieved compared to previous models. To underline this statement, the results for excess enthalpy of acetone/water system are shown for PSRK EOS in Figure 6. The above-mentioned error values raise the question whether the observed deviations are in agreement with the literature. To answer this question some sources should be named which deal with accuracy of h E predictions with GC models: • Schmid [1] gives relative deviations for different group combinations using VTPR, PSRK, and the modified UNIFAC (Dortmund) model. For the combination H 2 O/ CH 2 CO a value of ∼20% is given for the modified UNIFAC (Dortmund) and a value of ∼25% is given for VTPR EOS. However, no information on how these deviations were calculated or what data were used is described. • Hayashi et al. [36] evaluated predictions of excess enthalpy data using the ASOG GC model. For 51% of 871 binary data sets the maximum absolute deviation is below 100 J/mol, but also for 21.8% a maximum deviation >300 J/mol was stated.
For the predictions presented in this article, maximum absolute deviations were for 5 of 11 datasets (system at one temperature taken as dataset) below 100 J/mol, for 3 datasets deviations were between 100 and 150 J/mol, and 3 datasets showed even higher deviations. • Gmehling et al. [37] gives a mean deviation of 30% or 103.8 J/mol from investigations of 6000 datasets using the modified UNIFAC model. Several authors describe that accuracy of predictions with VTPR EOS is comparable with the modified UNIFAC (Dortmund) model [1][2][3]12]. In the present work only mean deviations of acetone/water system were higher than 30%.
• Vigh et al. [38] gives deviations for excess enthalpy data predicted with the modified UNIFAC models (Lungby and Dortmund). For the system n-hexane/benzene a value 7.6% is given for a temperature range between 25 and 45 • C, which is close to the predictions presented in this article if the results of Baluja Santos [34] are neglected.
For the system acetone/water a relative deviation of 30.6% (25-35 • C) was observed, but maximum experimental excess enthalpy was taken as reference value for all datapoints used for the calculation of relative errors. If this modification is applied to the datasets that were presented beforehand, a good agreement with the literature can be confirmed. • Chen et al. [39] investigated deviations for the modified UNIFAC (Dortmund) model and PSRK model. For alkane/benzene systems a mean deviation of 66 J/mol or 6.6% is given for the modified UNIFAC (Dortmund) model, while for PSRK EOS clearly higher errors are mentioned (454 J/mol or 43%). If an accuracy close to modified UNIFAC (Dortmund) model is expected, a good agreement can be stated with regard to the results presented in this article (results of Baluja Santos [34] neglected).
Overall, it seems reasonable to conclude that the results are in an expected range. Finally, it has to be stressed that investigations with additional binary systems will give a clearer picture on how well experimental data can be predicted by VTPR EOS.

Conclusions
Excess enthalpy values of four selected model systems could be predicted using the volume-translated Peng-Robinson equation of state. The necessary calculations were performed with Matlab analysis software.
In a first step the results allowed to draw conclusions about the kind and strength of interaction forces in the mixture compared with pure components. In addition, the caloric behavior of the mixing process could be assessed. In three of four cases an endothermic mixing process, which corresponds with weaker interaction forces in the mixture (compared with pure components), could be stated. Only for the acetone/water system was an exothermic behavior (and therefore stronger interaction forces in the mixture) found for low molar ratios of acetone.
From the deviations between experimental and predicted data of the selected model systems, it could be shown that predictions of strongly polar systems is challenging: mean deviations for the system acetone/water were above 50% or 200 J/mol at all temperatures. In contrast, deviations were lowest for alkane/benzene systems reaching a minimum mean value of 1.6% or 10 J/mol for the system n-hexane/benzene at 20 • C. Furthermore, especially for the systems n-hexane/benzene and 1-hexene/n-BuOH, the accuracy of predictions was strongly dependent on the temperature that was regarded. For the system 1-hexene/n-BuOH deviations at 25 • C were 2.5 times (absolute values) or almost 4 times higher (relative values) than at 40 • C.
Finally, a literature review of excess enthalpy measurements of n-hexane/benzene system at 25 • C could underline the importance of carefully checking experimental data: while eight sources showed maximum values between 800 and 950 J/mol, one source presented maximum values of 650 J/mol.

Software
Matlab analysis software (version 2020b, Update 1, 64-bit) was used for all calculations. To verify results Dortmund Data Bank Explorer Version 2015 (2.6.1.16) was used. The software can be downloaded from the following website: http://www.ddbst.com/ free-software.html#New (accessed on 8 November 2020).
Evaluation of data was carried out in Microsoft Excel ® (part of Microsoft 365 MsO package; 64-bit). constant used for calculation of temperature derivative of parameter a,g molar Gibbs energy, J · mol −1 h molar enthalpy, J · mol −1 L parameter TWU-α-function, -M parameter TWU-α-function, -N parameter TWU-α-function, -P pressure, bar Q group surface area UNIFAC model, -R universal gas constant, J · mol −1 · K −1 T absolute temperature, K U m,n /U n,m interaction parameter UNIFAC model, J · mol −1 v molar volume, dm 3 · mol −1 v i k number main group k in component i, -X group mole fraction UNIFAC model,z compressibility factor,z i , z j molar ratio of component i/j, -