Evaluation of the Peng–Robinson and the Cubic-Plus-Association Equations of State in Modeling VLE of Carboxylic Acids with Water

The performance of the classic Peng–Robinson (PR) and the modern Cubic-Plus-Association (CPA) equations of state were evaluated in modeling isobaric and isothermal vapor–liquid equilibria (VLE) of binary mixtures of carboxylic acids (formic, acetic, propanoic or butanoic) + water. Two functionalities of the alpha term were tested in PR, the original term proposed by Soave and the Matthias–Copeman term specially developed for modeling polar compounds. Within the Soave functionality, two generalized forms of the acentric factor were studied, the original general form and the Robinson and Peng modification for values of the acentric factor larger than 0.491. In addition, the case of PR with fitted parameters from saturated properties (as commonly obtained for modern equations of state) was also evaluated. VLE calculations without the use of a binary interaction parameter are in general more accurate with the modern CPA due to the association term; however, when a binary interaction parameter is used, the performance of the PR versions studied here are on average similar to those of CPA, and in some cases even superior. The original alpha function used in the PR equation and the original generalized form of the acentric factor are the best options for modeling organic acids + water systems when the binary interaction parameter is available. Temperature-dependent binary interaction parameters are provided as a database in modeling these complex systems.


Introduction
Carboxylic acids are important commodity chemicals due to their versatile applications. They are used in polymer synthesis, as food additives, pharmaceutical intermediates, paint and coating additives, plasticizers, among many other applications [1][2][3][4]. Many of the processes involving carboxylic acids will likely have the presence of water at some stage. It is therefore of interest for engineering design to have a reliable thermodynamic model applicable to a wide range of conditions of temperature, pressure and composition, and with the lowest computational cost. Thermodynamic modeling of carboxylic acids with water is nevertheless challenging since the systems exhibit strong non-ideal behavior due the presence of self-and crossassociation as well as polar interactions [1][2][3][4][5][6].
Classic cubic Equations of State (EoS) such as Soave-Redlich-Kwong (SRK) [7] and Peng-Robinson (PR) [8] have been used for decades in the chemical industry due to their reliability and simple mathematical form, characteristics that make them easy to implement and to modify [9,10]. In addition, the method to obtain their pure component parameters is universal, requiring only inputs of the experimental critical temperature ( T c ), critical pressure ( P c ) and acentric factor ( ) [11]. If experimental values of T c , P c and are not available, these can be estimated by group-contribution methods [12]. For the case of mixtures, the binary interaction parameter can also be estimated from group-contribution methods [13,14]. The disadvantages of the cubic equations are, however, well known, especially when modeling complex mixtures [9], and although PR is in general accepted as an improvement over SRK, several modifications have been made to the equation to enhance its performance [10]. A first approach to achieve this is to modify the temperature-dependent alpha function in the attractive term to obtain better predictions of saturated pressures, since it is expected that by doing this, the properties of mixtures are consequently improved. The Mathias-Copeman (MC) alpha function [15], specifically developed for PR, is one of the most popular functions found in the literature aimed to enhance the modeling of highly polar compounds [6]. Young et al. [16] have shown in their comparison of 20 alpha functions that MC is in fact one of the best for modeling pure component properties (vapor pressures) of water and carboxylic acids. Nevertheless, there are no studies comparing the performance of the MC function with that of the original alpha function for modeling mixtures of water and carboxylic acids. It is important to point out that although a common selection, the Soave [7] and the MC alpha functions do not satisfy the consistency test for alpha functions proposed by Le Guennec [6,[17][18][19], this is a requirement to get accurate and physically meaningful results, particularly in the supercritical region.
It is well established that the applicability of cubic EoS to complex mixtures is limited since they cannot handle complex intermolecular interactions explicitly due to their empirical nature. For this reason, in more recent years, theoretically based EoS such as the Statistical Associating Fluid Theory (SAFT), developed from statistical mechanical theories, have been formulated to model interactions explicitly such as strong association (hydrogen bonding). The model has been the subject of different modifications to increase the theory and applicability giving rise to the SAFT-type family of equations [20]. An interesting development was done by Kontogeorgis et al. [21] who combined a cubic equation (SRK) with the association term of SAFT producing the Cubic-Plus-Association (CPA) EoS. Despite its name, CPA is not cubic with respect to volume.
Carboxylic acid + water mixtures using CPA have been studied before but are mainly focused on the acetic acid + water system, and for low pressures and temperatures. For instance, Kontogeorgis et al. [22] modeled the acetic acid + water system at sub-atmospheric pressure. Later, and in order to improve the representation of this system, Muro-Suné et al. [23] introduced the Huron-Vidal mixing rule coupled with an NRTL expression. The model was further studied by Breil et al. [24] and Tsivintzelis and Kontogeorgis [25], who correlated enthalpies of vaporization and compressibility factors of pure acetic acid and in mixture with water. Although an improvement in the phase equilibria was obtained with the mixing rule, it was at the expense of increasing the number of adjustable parameters. Kontogeorgis et al. [26] modeled the formic acid and the propanoic acid + water systems at sub-and atmospheric pressures requiring a large (> 2) negative binary interaction parameter to get a satisfactory correlation. Román-Ramírez et al. [27] modeled the propanoic acid + water system at temperatures up to 483 K and pressures up to 1.9 MPa with CPA and the CR1 combining rule. More recently, Ribeiro et al. [28] evaluated the performance of CPA in modeling saturated pressure, density, speed of sound, second virial coefficient, compressibility factor, enthalpy of vaporization, and isobaric heat capacity of acetic acid. The authors concluded that it is not possible to obtain a single set of pure component parameters that can describe accurately all the properties of pure acetic acid and its phase equilibria.
In CPA, the SRK model was adopted; however, a different cubic equation, such as PR, can also be used, as shown by Pfohl et al. [29] and Perakis et al. [30]. The latter authors modeled acetic acid + water mixtures concluding that the non-association contribution of the SAFT and simplified SAFT equations did not offer any advantage over the simpler cubic EoS.
On the other hand, Voutsas et al. [31] have shown that it is possible to improve the predictive capability of PR if instead of using experimental critical properties and acentric factor, these are regressed from experimental vapor pressures and saturated liquid density data to calculate the attractive and repulsive parameters. The authors showed that PR with fitted parameters provided better results compared with SAFT in modeling alkane systems.
In the present paper, a comparison of the performance of the classic PR and the modern CPA EoS is made on their ability to predict and to correlate isobaric and isothermal vapor-liquid equilibria (VLE) of binary systems of carboxylic acids (from formic acid to butanoic acid) + water over a wider range of temperature and pressure conditions than previous studies. The results with the original alpha function and the literature recommended MC function for polar compounds are compared to determine the effect of using this purpose-built function. In addition, the performance of PR with pure component parameters fitted to saturated pressures and volumes (as commonly obtained for modern equations of state) is evaluated to discern whether this approach can improve the predictive capability of PR. The CPA model as proposed by Kontogeorgis et al. [21] is employed in this work with both compounds modeled as 2B based on former studies [22,27,28].

Peng-Robinson
The PR EoS in pressure-explicit form is [8]: where parameters a and b are computed from van der Waals one-fluid mixing rules according to and the following combining rules: k ij in Eq. 4 is the binary interaction parameter introduced to correct the attractive forces. Parameters a i and b i are commonly obtained from correlations of T c , P c and of the pure compounds: The alpha function in Eq. 6 is the one proposed by Soave for the Redlich-Kwong equation [7]: where the reduced temperature ( T r ) is defined as T r = T∕T c , and m is a modified function of , originally expressed as and later modified by Robinson and Peng [32]: is applicable when ≤ 0.491 while Eq. 10 when > 0.491 . When Eq. 9 is used in this work it is referred as PR (i.e., the original PR formulation), and when Eq. 10 is used it is referred as PR-m. In the current study, PR-m is only applicable to propanoic acid and butanoic acid (Table 1).

Mathias-Copeman Alpha Function
The MC function is [15]: where MC 1 , MC 2 and MC 3 are compound specific adjustable parameters determined from VLE data. The PR equation with the MC function is referred here as PR-MC.

PR Fitted
The same Soave-type expression of that used in CPA [21] for the energy parameter is used in this work but applied to PR, this is Parameters a 0,i and c 1,i , together with b i , in Eq. 3 can be found through fitting of experimental vapor pressure ( P V ) and liquid density ( L ) data, as commonly done in regressing pure component parameters in SAFT-type models. PR with fitted parameters is referred in this work as PR-f.

CPA
CPA can be written as [33]: where the first two terms are the SRK model to account for the physical interactions, whereas the third term accounts for the association interactions. The site monomer fraction ( X A i ) is computed from  CPA propanoic acid and water parameters taken from Román-Ramírez et al. [27] Compound where association strength ( with the radial distribution function ( g): and where A i B j and A i B j are the association energy and association volume, respectively. For mixtures, these can be calculated from the CR1 combining rules as follows: The energy parameter is calculated from Eq. 12 and the combining rules for a ij and b ij are those given by Eqs. 4 and 5, respectively. Five pure component parameters are required to define a compound in CPA:

Pure Component and Binary Interaction Parameters
Critical properties and acentric factors for the compounds are given in Table 1. Pure component parameters for PR-f and CPA were obtained by fitting experimental P V and L as described elsewhere [27]. For the case of PR-MC only, P V data were used in the fitting since the model is only focused on improving this property. The temperature range considered was from the triple point up to 0.99T r , with data taken from DIPPR [34].
For the mixtures, the optimum k ij was regressed from experimental bubble point pressures ( P ) and vapor compositions ( y ) of the organic acid using the following objective function: where N is the number of data points used in the fitting and superscripts exp and calc stand for an experimental and a calculated property, respectively. Due to the limited availability of isothermal data for the butanoic acid + water system, k ij was fitted to isobaric data using the bubble temperature ( T ) instead of P in Eq. 20.

Thermodynamic Evaluation
The thermodynamic models were evaluated by comparing the average deviations between experimental and calculated values of P V and L , for the case of pure components, and of P and y for mixtures, according to Eqs. 21 and 22. Two different modes were considered for mixtures: a predictive ( k ij = 0 ) and a correlative ( k ij ≠ 0 ) mode. For the case of the butanoic acid system, the average deviations in T are also included (Eq. 21).
where stands for either P V , L , P or T.

Pure Compounds
Fitted parameters for the models are presented in Table 1 together with the estimated deviations in P V and L . The parameters for PR-MC are comparable to those reported by Young et al. [16], and the differences can be attributed to the different temperature range used in the optimization and the search algorithm. CPA parameters for the acids are also between the ranges of previously reported values by Derawi et al. [35]. However, water parameters for the same association scheme differ from previous publications in which monomer fraction data were also used in the fitting procedure; nevertheless, the deviations in P V reported in this work are lower than the values found the literature [36]. In general, for the organic acids, the smallest deviations in P V are obtained with PR-MC, whereas the largest are given by PR. Particularly, large deviations are obtained when PR-m is used for the case of propanoic and butanoic acid. PR-MC improves considerably the P V calculation compared with PR since it results in more than 90 % lower deviations. The highest improvement can be seen in particular for formic acid with a 98 % decrease. The performance of PR-MC and CPA is quite similar for this property, although for formic acid the correlations are considerably better with the former. As expected, the PR-MC model does not provide any advantage in terms of L calculation, it gave practically similar errors as those obtained with PR and PR-m. Looking at both properties, CPA gives overall the best correlation. Figure 1 shows a graphical comparison of the models for ΔP V and Δ L calculations. For the case of water, the best fitting in P V is given by PR-MC whereas the best fitting in L by PR-f. The second best fitting for both properties is obtained with CPA. It is worth nothing that the value of the repulsive parameter is practically the same between PR-f and CPA for the individual compounds. The major shift in value is observed for the a 0,i parameter when the association term is introduced, clearly showing how the attractive term is being affected by the association contribution. Experimental data from DIPPR [34] As observed for the organic acids and water, the deviations in density for PR-f are more than 90 % lower compared with PR and PR-m but these are in general higher than those of CPA. The improved density correlation with PR-f is the result of obtaining the pure component parameters from saturated properties; but as a consequence, and similarly to the SAFT-type models, the critical pressure is now overestimated [11]. There is no benefit in the L calculation using a different equation among PR, PR-m or PR-MC. An improvement in density calculation can only be achieved by applying a volume translation technique (either temperature independent or temperature dependent) as demonstrated by Jaubert et al. [37].

Predictions
The deviation in P and y ( T and y for the case of isobaric calculations for the butanoic acid system) are presented in Table 2. The best predictions are undoubtedly given by CPA with errors in pressure more than 50 % lower compared with the rest of the equations and with average deviations of 0.14 % in vapor composition. As examples of the CPA performance, Figs. 2 and 3 show a graphical comparison of CPA and that of the PR equations. Although the predictions with CPA are an improvement over the PR equations studied, the deviations can still be considered high (overall average of 50 %). In addition, like the PR equations, CPA also fails to represent the observed experimental phase behavior adequately. For instance, at 423.2 K for the propanoic acid + water system (Fig. 3), CPA wrongly predicts a heterogeneous azeotrope instead of the measured homogeneous azeotrope at the low propanoic acid concentration region. Similarly, and as shown by other studies [26], a phase diagram opposite to the experimental data is predicted with CPA for the formic acid + water at all the conditions studied, as exemplified in Fig. 4. These results show, however, that including the association term in the cubic equation does improve its predictive capabilities.
Looking only at the PR versions, the PR-m model although resulted in higher deviations in P V compared with PR, its performance in P is in this instance slightly better. Nevertheless, for practical purposes, the results can be considered similar. The use of PR-MC is in this case detrimental since larger deviations in P are found compared with the original PR, particularly for the formic acid and acetic acid systems. This is unexpected considering the PR-MC gave the lowest errors in P V . PR-f provides considerably lower deviations in P for the formic acid and acetic acid systems but for propanoic acid and butanoic acid, the deviations are the highest of all PR versions. This trend can also be observed in the isobaric calculations for butanoic acid. This shows that the CPA predictive capability can be attributed to the association term and not necessarily to the pure component parameter found through fitting of saturated properties. In general, PR-f gives the best predictions for formic acid, PR for acetic acid, and PR-m for propanoic acid and butanoic acid. The performance in predicting vapor composition is practically the same for all the PR models with average values of Δy = 20 %. Considering all the systems, the best predictive capabilities are given by PR-f whereas the poorest by PR-MC.

Correlations
The k ij for each temperature and system can be found in Table 3. The parameters were then adjusted to a linear temperature dependency with the form k ij = A ij + B ij T . Values for A ij and B ij are presented in Table 4 while the calculated deviations for the correlations using this functionality are presented in Table 2. Unlike the predictions, there is no real difference between the correlations with PR, PR-MC, PR-f and CPA since all models give rounded overall average ΔP of 7 %. Depending on the system under consideration, the correlation performance can differ with one or another model. For instance, for the case of the acetic acid system, CPA gives the best correlations, especially compared with PR, but for propanoic acid + water, CPA gives in fact the highest errors against the rest of the models. No obvious trend was found to explain this.
CPA is better at modeling the formic acid system compared with the PR equations, but its performance decreases as one moves higher in the carboxylic acid chain, in such a way that for the propanoic acid and the butanoic acid systems, CPA results in higher deviations than any of the PR models. This contradicts other studies in which a poorer CPA performance for the low chain carboxylic acids (formic and acetic acid) has been reported [52]. This can be attributed to the different values of the pure component parameters employed, due to differences in the search algorithm and temperature range covered in the fitting procedure. The CPA performance for formic acid can be appreciated in Fig. 5, showing the VLE at 13.33 kPa, where the linear temperature-dependent k ij found from isothermal data is used in the calculation.
The PR models studied show a similar trend as the results observed for CPA. With formic acid, the deviations are larger with PR followed by PR-MC and PR-f, but the order shifts when moving toward butanoic acid, for which the best performance is then given by PR, followed by PR-m, PR-MC and finally PR-f. Contrasting the findings in predictive mode, the use of Eq. 10 instead of Eq. 9 introduces slightly higher errors in the modeing (PR-m: ΔP = 10.13 % vs. PR: ΔP = 9.83 %). Recently, Pina-Martinez et al. [54] have published updated expressions of the generalized Soave alpha-function for SRK and PR. We performed calculations for the formic acid + water and propanoic acid + water systems and the results were similar to the ones obtained for PR (as defined in this work), and consequently better than the results given by PR-m (for the case of propanoic acid). As mentioned by the authors, a noticeable improvement over the original generalized function will be mainly given for molecules with larger than 0.9 [54].
A negative k ij with a relatively high magnitude (average values ranging between − 0.15 to − 0.25) was obtained with all models, indicating a large correction of the attractive forces. PR-f resulted in general in the lowest magnitudes, whereas PR and CPA gave practically the same values.
In agreement with previous works comparing cubic equations and CPA [55], it was found that the description of isobaric data is better than the isothermal case, as   Table 2 (continued)  Fig. 6 for the butanoic acid system. In this case, analogous to the isothermal calculations, PR provides better correlations than CPA; however, both models are not able to represent the experimental homogeneous azeotrope, showing a heterogeneous azeotrope instead.  All models seem to provide better predictions and correlations as the temperature increases, probably due to the weakening of the polar and association molecular interactions, with no clear advantage with one model over another at the highest temperatures. As an example, Figs. 7 and 8 show the VLE diagram for the propanoic acid system at 343.2 K and 453.2 K, respectively. At 343.2 K, only PR-f and CPA are able to correlate the homogeneous azeotrope, with PR-f giving closer correlations to the experimental data. While at 343.2 K there are clear differences in the performance of the models, at 453.2 K these are less obvious.
It is apparent that while the association term does improve the predictive capabilities of the cubic equation, it also increases the deviations for some of the systems in correlative mode. It is possible that a different association scheme, other than the 2B model used in this work, could lead to better results, if not for all, at least for some of the systems. This point is exemplified by Kontogeorgis and Folas [52] who reported a satisfactory isobaric modeling for the propanoic acid + water system with the rigorous 1A scheme for the acid and the 4C for water. Nevertheless, isothermal calculations in the present work at 303.2 K with the 1A-4C combination for propanoic acid + water resulted in higher deviations in pressure (13.46 %) than the 2B-2B case (10.75 %, Table 2). Similarly, for acetic acid + water, the 2B-2B combination resulted in better correlations (4.89 %) than the 1A-4C (5.54 %, Table 2) at 293.2 K.

Conclusions
The PR and the CPA EoS can be found nowadays in almost any process simulation software. Literature suggestions for modeling polar compounds such as carboxylic acids and water with the PR EoS include the use of a particular generalized function  of the acentric factor for values greater than 0.491, as is the case for propanoic acid and butanoic acid, and the use of the Mathias-Copeman alpha function. However, the results presented in this work show that the original generalized function of the acentric factor, in theory applicable only to values lower than 0.491, is a better option than the modified function. Alternatively, the updated expression of the   [54] can be used. It is also shown that an improvement on saturated properties, especially vapor pressure, due to the use of a particular alpha function, will not necessarily lead to a better performance on mixture calculations.
The predictive capability of CPA outperforms that of all the PR versions studied, although in most of the cases the VLE representations are only qualitative. The   [27]. Lines: EoS correlations association term in CPA does in fact improve the predictive capabilities of the physical term, but when a binary interaction parameter is used there is no real benefit in the use of the term since the performance of the cubic model is on average similar to CPA, and in some cases even superior. For engineering calculations at a low computational cost, PR with parameters obtained from correlations of critical properties and acentric factor, the original generalized function of the acentric factor and the binary interaction parameters reported in this work, is a suitable alternative to model organic acids + water systems.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.