Modeling Binary and Multicomponent Systems Containing Supercritical CO2 with Polyethylene Glycols and Compounds Relevant to the Biodiesel Production

The CPA equation of state is applied to model binary, ternary, and multicomponent mixtures that contain CO2 with polyethylene glycols or compounds relevant to biodiesel production, such as glycerol and various triglycerides. Effort has been made to evaluate the model performance on correlating both the liquid and the vapor phase compositions, which is a demanding task, revealing the model’s and parameters’ limitations, due to the rather low concentrations of heavy compounds in the vapor phase. Initially the model’s binary parameters, which in all cases were temperature independent, were estimated using experimental data for binary systems. Those parameters were used to predict the phase behavior of supercritical CO2 containing ternary and multicomponent mixtures. Since no parameter was adjusted to ternary or multicomponent systems’ data, the reported CPA results for such mixtures are considered as pure predictions. This is the final part of a series of studies [Tsivintzelis et al. Fluid Phase Equilibria 430 (2016) 75–92 and 504 (2020) 112337] that complete the parameterization of the CPA equation of state for systems relevant to the biodiesel production, which allows the application of the model to multicomponent mixtures of the relevant processes.


Introduction
The processing of vegetable seeds and oils with supercritical CO 2 has attracted significant attention mainly because the relevant extraction or purification processes require lower temperatures, while CO 2 is nontoxic, and it is easily removed from the products simply by depressurizing the system. Such pressure-induced separation diminishes the need for adequate thermal separations, as in the traditional processes with liquid and often toxic solvents. Another important advantage is that the solubility of other substances in such a supercritical solvent is strongly affected by moderate changes in temperature and pressure, rendering the solvent power of supercritical CO 2 as a tunable parameter that can be adjusted to the needs of specific processes. This allows its successful use in various demanding separations of high-value compounds from mixtures of similar fluids, such as the separation of antioxidants from natural products. Despite the disadvantage of higher capital and operational costs imposed by the use of high pressures, the aforementioned advantages are often crucial in the food, cosmetic, pharmaceutical, and nutraceutical industries [1][2][3].
However, vegetable oils are also used to produce biodiesel, whose production volume increases following the need for replacing traditional fossil fuels with fuels from renewable sources. Biodiesel is produced through the transesterification reaction of various vegetable oils with methanol (or ethanol), which results in multicomponent mixtures containing methyl (or ethyl) esters of fatty acids and glycerol as a byproduct. Recently, the prediction of their values in cases of lack of experimental data. Some multicomponent biodiesel systems were recently modeled [30].
In this study, the CPA equation of state was applied for first time to polyethylene glycols of rather low polymer molecular weight, and correlations of the pure fluid parameters are presented allowing the prediction of their values for various molecular weights. Then, the model was applied to binary systems that contain CO 2 with polyethylene glycols, glycerol, or various glycerides, and the binary parameters were obtained. Using binary parameters optimized solely to experimental data for binary mixtures, the predicting ability of the model was evaluated for ternary and multicomponent systems.

The CPA Equation of State
The cubic-plus-association (CPA) equation of state (EoS) is a combination of the SRK EoS and the association term of the SAFT type models [9][10][11]. Thus, this model is capable of describing systems with hydrogen-bonding fluids such as water, alcohols, and glycols. The EoS, in terms of pressure, is given by the following equation: where x i is the mole fraction of component i, and X Ai is the fraction of the free (not associated with other sites) sites of type A on molecule i, related to the association strength ∆ AiBj as follows: where g(ρ) = 1/(1 − 1.9n) is the radial distribution function, n = (1/4)bρ the reduced density, ρ the molar density, ε AiBj the association energy, β AiBj the association volume, and b ij = b i + b j /2. The co-volume parameter, b i , of component i is considered temperature independent. A Soave-type relationship is used for the interaction energy of the physical term, as follows: where a 0 and c 1 are pure fluid parameters, and T r is the reduced temperature. The next mixing and combining rules are used for the SRK physical term [9], whereas the one temperature-independent binary interaction parameter, k ij , was used according to the following relationships: For interactions between two self-associating fluids, i and j, the CR-1 rule is used through the following relationships: For interactions between one self-associating with one non self-associating fluid, the modified CR-1 (mCR-1) combining rule is used [31]:

CPA Pure Fluid Parameters
For non-associating pure fluids, the CPA EoS is reduced to the SRK model and requires the knowledge of three pure fluid parameters, which, in contrast with the SRK EoS, are not estimated form the critical properties but are adjusted to pure fluid experimental vapor pressures and liquid densities. In cases of self-associating fluids, the application of the model requires the knowledge of two additional parameters, namely the association energy and the association volume, which may be simultaneously adjusted with the three physical parameters, but they can also be adopted form experimental values or ab initio calculations [31].
All the pure fluid parameters needed for calculations of this study were adopted from the literature except for the parameters of glycols, which were estimated by adjusting the model predictions to DIPPR data for vapor pressures and liquid molar volumes [32]. All parameters used in this study are presented in Table 1. For glycerol, the 3 × 2B association scheme is used to describe the self-associating behavior since glycerol contains three hydroxyl groups in each molecule [7]. In addition, here, it is worth mentioning that the parameters of triolein were predicted using the correlations presented by Tsivintzelis et al. [7] due to lack of (confirmed) experimental data.Since experimental data for vapor pressures do not exist for various polyethylene glycols, the estimation of the CPA pure fluid parameters becomes difficult. Many studies showed that the estimation of pure fluid parameters of various models, including the SAFT type ones and lattice fluid theories, by using only volumetric properties does not result in optimum parameters for phase equilibrium calculations of binary systems [33]. For this reason, various approaches to predict such parameters were developed, including group contribution methods [33][34][35]. Using the CPA model for compounds of the same family, it was shown that the co-volume parameter, b, is a linear function of molecular weight (or the van der Waals volume of the fluid), the c 1 parameter reaches a plateau for relatively high molecular weights, while the energetic parameter, a 0 , of the physical term shows a linear or a mild polynomial dependence from molecular weight (or the van der Waals volume of the fluid) [7,17,28]. However, it was shown that appropriate correlations of the three physical CPA parameters are facilitated when the association parameters are kept constant for all fluids of the same family since the two energetic parameters, i.e., the a 0 and ε parameters, are highly correlated to each other [7]. Having in mind such findings for the CPA model, correlations of the pure fluid parameters for glycols were estimated as presented in Figure 1. For a 0 , both a linear and a second-order polynomial dependence was used. Such correlations allow for the prediction of pure fluid parameters for polyethylene glycols of various molecular weights, as shown in Table 2. The critical temperature reported in Table 2 was predicted using the Joback's group contribution method for the closest ethylene glycol [36].

CO2-Polyethylene Glycols VLE
Next, the model was applied to describe the solubility of CO2 in polyethylene glycols of various molecular weights. In all cases, the a0 was adopted from the linear correlation (see Figure 1 and Table 2) since it was shown that it gives slightly more accurate results, especially for higher molecular weights.
Tsivintzelis et al. modeled the CO2-glycol vapor-liquid equilibrium (in systems containing ethylene glycol, diethylene glycol, and triethylene glycol) using two major approaches, i.e., assuming CO2 as inert compound and assuming that CO2 has one site that is able to cross-associate with glycols [15]. Using the first approach, the model satisfactorily describes the liquid phase compositions; while only using the second approach (and, consequently, using one extra binary parameter) the model satisfactorily describes the vapor phase compositions. However, since in systems with polyethylene glycols, the liquid phase is of interest; the first approach was used in this study for simplicity, and CO2 was modeled as inert compound.
The results using one temperature-independent binary interaction parameter, kij, are

CO 2 -Polyethylene Glycols VLE
Next, the model was applied to describe the solubility of CO 2 in polyethylene glycols of various molecular weights. In all cases, the a 0 was adopted from the linear correlation (see Figure 1 and Table 2) since it was shown that it gives slightly more accurate results, especially for higher molecular weights.
Tsivintzelis et al. modeled the CO 2 -glycol vapor-liquid equilibrium (in systems containing ethylene glycol, diethylene glycol, and triethylene glycol) using two major approaches, i.e., assuming CO 2 as inert compound and assuming that CO 2 has one site that is able to cross-associate with glycols [15]. Using the first approach, the model satisfactorily describes the liquid phase compositions; while only using the second approach (and, consequently, using one extra binary parameter) the model satisfactorily describes the vapor phase compositions. However, since in systems with polyethylene glycols, the liquid phase is of interest; the first approach was used in this study for simplicity, and CO 2 was modeled as inert compound.
The results using one temperature-independent binary interaction parameter, k ij , are presented in Table 3 and in Figure 2. In all cases, the CPA correlations show a satisfactory agreement with the experimental data. The values of the binary interaction parameters are low and positive for low molecular weights, while they become negative for higher molecular weights. Since increasing the molecular weight, the strong intermolecular interactions are expected to decrease, the negative values of k ij for higher molecular weights mainly reflect the inaccuracies of PEG's pure fluid parameters due to the approximate way of their estimation, i.e., prediction using the correlations of Figure 1.

CO2-Glycerol Binary System
In contrast to polyethylene glycols, glycerol is a low-molecular-weight compound. It has three hydroxyls in each molecule, and consequently, it presents considerable hydrogen bonding. Furthermore, Tsivintzelis et al. showed that systems containing CO2 and alcohols are better described by the CPA model if strong specific interactions between the positively charged carbon atom of CO2 and the oxygen of the hydroxyl group are taken into consideration, as indicated by various ab initio studies, assuming one site that is able Figure 2. Solubility of CO 2 in polyethylene glycols. Experimental data (points [38]) and CPA calculations (lines) using the k ij values of Table 3. Table 3. Binary interaction parameters (k ij ) and deviations from experimental data for CO 2 -polyethylene glycols (the average molecular weight is shown in parenthesis).

System
Temp.

CO 2 -Glycerol Binary System
In contrast to polyethylene glycols, glycerol is a low-molecular-weight compound. It has three hydroxyls in each molecule, and consequently, it presents considerable hydrogen bonding. Furthermore, Tsivintzelis et al. showed that systems containing CO 2 and alcohols are better described by the CPA model if strong specific interactions between the positively charged carbon atom of CO 2 and the oxygen of the hydroxyl group are taken into consideration, as indicated by various ab initio studies, assuming one site that is able to cross-associate with alcohols on every CO 2 molecule [29,31]. In addition, only using the same approach for CO 2 -low-molecular glycol systems, the very small concentration of glycols in the vapor phase can be adequately described [15].
In this work, the CO 2 -glycerol VLE was modeled using two approaches, i.e., assuming CO 2 as inert compound and, secondly, assuming one positive site on every CO 2 molecule that is only able to cross-associate with glycerol's negative sites (the 3 × 2B association scheme was used for glycerol). Thus, only the binary interaction parameter (k ij ) was adjusted in the first case, while both the k ij and the cross-association volume (β cross ) were optimized in the second one (using the mCR-1 rule). The obtained binary parameters and the deviations of model calculations from experimental data are presented in Table 4. Some representative results are presented in Figure 3, while more results are shown in Figures S1-S6 of the supplementary information file. Similar conclusions to CO 2 -glycol mixtures [15] were obtained; i.e., using a single binary interaction parameter, the model satisfactorily describes the liquid phase, but only considering CO 2 as a cross-associating compound and, consequently, using two adjustable binary parameters, the model satisfactorily describes the vapor phase compositions. Thus, such a mixture that presents various self-and cross-hydrogen-bonding interactions is better described if, in the CPA model, all possible association interactions are explicitly accounted for. Here, it is worth mentioning that, as shown by the CPA correlations, the glycerol mole fraction in the vapor phase presents a minimum when plotted against pressure ( Figure 3). This resembles the minimum that is observed for the CO 2 -water system and reflects the existence of vapor-liquid-liquid equilibrium in relatively low temperatures. Similarly to the minimum of CO 2 -water system, such behavior of the CO 2 -glycerol mixture, predicted by the model, may be of interest in various separation processes.
mixtures [15] were obtained; i.e., using a single binary interaction parameter, the model satisfactorily describes the liquid phase, but only considering CO2 as a cross-associating compound and, consequently, using two adjustable binary parameters, the model satisfactorily describes the vapor phase compositions. Thus, such a mixture that presents various self-and cross-hydrogen-bonding interactions is better described if, in the CPA model, all possible association interactions are explicitly accounted for.  [39][40][41][42] and CPA correlations using the interaction parameters of Table 4.  [39][40][41][42] and CPA correlations using the interaction parameters of Table 4.

CO 2 -Glycerides
The phase behavior of five CO 2 -triglyceride mixtures containing tricaprylin, tributyrin, trilaurin, trimyristin, and triolein was investigated. In all cases, both CO 2 and triglycerides were modeled as inert compounds, and only one temperature-independent binary parameter, k ij , was used. The obtained binary parameter and the deviations of model calculations from experimental data are presented in Table 5. Overall, it was observed that the modeling of such mixtures present significant difficulty, mainly due to the very low concentration of the heavy compound (triglyceride) in the vapor phase. Forcing the model to satisfactorily describe the vapor phase increases the deviations for the liquid phase compositions. Consequently, in all cases the model correlations for the CO 2 mole fraction of the liquid phase range between 5-11%. In more detail, considering the CO 2 -tricaprylin system, representative calculations are presented in Figure 4 and in Figures S7 and S8 of the supplementary information file. The model represents very satisfactorily the very low vapor phase content of tricaprylin but underestimates to a low extent the CO 2 mole fraction of the liquid phase.   [43] and CPA correlations using the binary interaction parameter of Table 5.
Results for the CO2-tributyrin are shown in Figure 5 and Figures S9-S11 of the supplementary information file. For this system, the binary interaction parameter was adjusted solely to liquid phase data since data for the vapor phase are not available. The model shows satisfactory agreement with the experimental data. Figure 5. CO2-tributyrin VLE at 298 K. Experimental data [44,45] and CPA correlations using the binary interaction parameter of Table 5.
Results for the CO2-trilaurin system are illustrated in Figure 6 and in Figures S12-S14 of the supporting information file. For this system, only one set of data for liquid phase  [43] and CPA correlations using the binary interaction parameter of Table 5.
Results for the CO 2 -tributyrin are shown in Figure 5 and Figures S9-S11 of the supplementary information file. For this system, the binary interaction parameter was adjusted solely to liquid phase data since data for the vapor phase are not available. The model shows satisfactory agreement with the experimental data.
Results for the CO 2 -trilaurin system are illustrated in Figure 6 and in Figures S12-S14 of the supporting information file. For this system, only one set of data for liquid phase compositions is available (see Figure S14 of the supporting information file), and consequently, the binary parameter was regressed using mainly vapor phase compositions, which, however, are available for high pressures. It is shown that the model satisfactorily describes the solubility of the glyceride in the vapor or the supercritical fluid phase up to approximately 200 bar, while higher deviations are observed for higher pressures.  [43] and CPA correlations using the bi-nary interaction parameter of Table 5.
Results for the CO2-tributyrin are shown in Figure 5 and Figures S9-S11 of the supplementary information file. For this system, the binary interaction parameter was adjusted solely to liquid phase data since data for the vapor phase are not available. The model shows satisfactory agreement with the experimental data.  [44,45] and CPA correlations using the binary interaction parameter of Table 5.
Results for the CO2-trilaurin system are illustrated in Figure 6 and in Figures S12-S14 of the supporting information file. For this system, only one set of data for liquid phase compositions is available (see Figure S14 of the supporting information file), and consequently, the binary parameter was regressed using mainly vapor phase compositions, which, however, are available for high pressures. It is shown that the model satisfactorily describes the solubility of the glyceride in the vapor or the supercritical fluid phase up to approximately 200 bar, while higher deviations are observed for higher pressures.  [44,45] and CPA correlations using the binary interaction parameter of Table 5.  [46,48,49] and CPA correlations using the binary interaction parameter of Table 5.
Representative calculations for the CO2-trimyristin system are presented in Figure 7 and in Figures S15 and S16 of the supporting information file. No liquid phase compositions are available for this system, and consequently, the binary interaction parameter was adjusted only to vapor phase compositions. Similarly, for the CO2-trilaurin system, the model satisfactorily describes the very small solubility of trymyristin in the fluid phase up to approximately 200 bar, while higher deviations are observed at higher pressures.  [46,48,49] and CPA correlations using the binary interaction parameter of Table 5.
Representative calculations for the CO 2 -trimyristin system are presented in Figure 7 and in Figures S15 and S16 of the supporting information file. No liquid phase compositions are available for this system, and consequently, the binary interaction parameter was adjusted only to vapor phase compositions. Similarly, for the CO 2 -trilaurin system, the model satisfactorily describes the very small solubility of trymyristin in the fluid phase up to approximately 200 bar, while higher deviations are observed at higher pressures.
Representative calculations for the CO2-trimyristin system are presented in Figure 7 and in Figures S15 and S16 of the supporting information file. No liquid phase compositions are available for this system, and consequently, the binary interaction parameter was adjusted only to vapor phase compositions. Similarly, for the CO2-trilaurin system, the model satisfactorily describes the very small solubility of trymyristin in the fluid phase up to approximately 200 bar, while higher deviations are observed at higher pressures.  [46,48] and CPA correlations using the binary interaction parameter of Table 5.
Results for the CO2-triolein mixture are presented in Figure 8 and in Figures S17-S20 of the supplementary information file. As presented in these Figures, the experimental data are scattered, while for the glyceride content of the vapor phase, data from different sources may differ up to two orders of magnitude. As presented in Table 5, the model correlations present an average deviation form experimental glyceride mole fraction of the vapor phase equal to 139%, but to a large extent, such value is attributed to the  [46,48] and CPA correlations using the binary interaction parameter of Table 5.
Results for the CO 2 -triolein mixture are presented in Figure 8 and in Figures S17-S20 of the supplementary information file. As presented in these Figures, the experimental data are scattered, while for the glyceride content of the vapor phase, data from different sources may differ up to two orders of magnitude. As presented in Table 5, the model correlations present an average deviation form experimental glyceride mole fraction of the vapor phase equal to 139%, but to a large extent, such value is attributed to the scattering of the experimental data, and Figure 8 and Figures S17-S20 present the capabilities of the model in a more realistic way. Furthermore, the reader should have in mind that, as mentioned above, the triolein pure fluid parameters were predicted using appropriate correlations in absence of (confirmed) experimental data. scattering of the experimental data, and Figures 8 and S17-S20 present the capabilities of the model in a more realistic way. Furthermore, the reader should have in mind that, as mentioned above, the triolein pure fluid parameters were predicted using appropriate correlations in absence of (confirmed) experimental data.  [51,53] and CPA correlations using the binary interaction parameter of Table 5.

Ternary and Multicomponent Systems
Next, the model was applied to predict the phase behavior of ternary and multicomponent mixtures and the results were compared to experimental data from literature [55][56][57][58][59][60]. In all cases, no parameter was regressed using ternary or multicomponent systems data, but all binary parameters were adopted from the respective sub-binary mixtures. Consequently, the calculations that are presented in this section are considered as pure predictions. All the used binary parameters are presented in Table 6.  [51,53] and CPA correlations using the binary interaction parameter of Table 5.

Ternary and Multicomponent Systems
Next, the model was applied to predict the phase behavior of ternary and multicomponent mixtures and the results were compared to experimental data from literature [55][56][57][58][59][60]. In all cases, no parameter was regressed using ternary or multicomponent systems data, but all binary parameters were adopted from the respective sub-binary mixtures. Consequently, the calculations that are presented in this section are considered as pure predictions. All the used binary parameters are presented in Table 6.  [7] a For combining rule 1 (CR1) and modified combining rule 1 (mCR-1), see ref. [31]. b Parameters adopted from methanol-methyl laurate in absence of parameters for the particular system. Initially, the model was applied to describe the vapor-liquid equilibrium of mixtures containing CO 2 with glycerol and methanol or ethanol. In all cases, one positive association site was assumed on every CO 2 molecule that is only able to cross-associate with the negative sites of glycerol, methanol, or ethanol, and consequently, the modified CR-1 (mCR-1) combining rule was used for estimating the cross-association parameters [31]. Further, as shown in Table 6, the glycerol-methanol (or ethanol) cross-association parameters were calculated using the CR-1 combining rule [31].
The obtained results are shown in Figures 9-11. Considering the CO 2 -glycerolmethanol system, the model presents an average absolute deviation (in bubble point pressure) from the experimental data shown in Figure 9 [55] equal to 8% and from experimental data shown in Figure 11a [59] equal to 4%. As shown in Figure 9, deviations increase as the mixtures become richer in CO 2 . (mCR-1) combining rule was used for estimating the cross-association parameters [31]. Further, as shown in Table 6, the glycerol-methanol (or ethanol) cross-association parameters were calculated using the CR-1 combining rule [31].
The obtained results are shown in Figures 9-11. Considering the CO2-glycerol-methanol system, the model presents an average absolute deviation (in bubble point pressure) from the experimental data shown in Figure 9 [55] equal to 8% and from experimental data shown in Figure 11a [59] equal to 4%. As shown in Figure 9, deviations increase as the mixtures become richer in CO2.   . Phase boundaries for CO2-glycerol-methanol ternary system. Experimental data (points, [55]) and CPA predictions using the binary interaction parameters of     Figure 11. CO2 mole fractions in liquid and vapor phases of (a) CO2-glycerol-methanol and (b) CO2glycerol-ethanol ternary systems containing alcohol to glycerol ratio equal to 113. Experimental data [59,60] and CPA predictions using the binary interaction parameters of Table 6.
On the other hand, considering the CO2-glycerol-ethanol system shown in Figure 10, the model describes in a more accurate way mixtures that are rich in CO2. In this case, the average absolute deviation (in bubble point pressure) from the experimental data shown in Figure 10 [56] and Figure 11b [60] is equal to 14% and 2%, respectively.
Next, CPA was applied to predict the phase behavior of the ternary CO2-methanollauric acid mixture. Such mixture is reactive since it contains both an alcohol and an organic acid, which renders the accurate measurement and modeling of the phase behavior a demanding task. Similarly to the previous systems, CO2 was modeled assuming one positive association site that can only cross-associate with methanol. The obtained results are presented in Figure 12. The model shows an average absolute deviation (in bubble point pressure) from the experimental data shown in Figure 12 [57] equal to 9%. Ferreira et al., who obtained the experimental data, used the Peng-Robinson EoS and in total six binary parameters (two binary parameters per sub-binary system), while two of them (those referring to the methanol-lauric acid mixture) were adjusted to the ternary system's data [57]. According to their results, the experimental data are satisfactorily correlated. Figure 11. CO 2 mole fractions in liquid and vapor phases of (a) CO 2 -glycerol-methanol and (b) CO 2glycerol-ethanol ternary systems containing alcohol to glycerol ratio equal to 113. Experimental data [59,60] and CPA predictions using the binary interaction parameters of Table 6.
On the other hand, considering the CO 2 -glycerol-ethanol system shown in Figure 10, the model describes in a more accurate way mixtures that are rich in CO 2 . In this case, the average absolute deviation (in bubble point pressure) from the experimental data shown in Figure 10 [56] and Figure 11b [60] is equal to 14% and 2%, respectively.
Next, CPA was applied to predict the phase behavior of the ternary CO 2 -methanol-lauric acid mixture. Such mixture is reactive since it contains both an alcohol and an organic acid, which renders the accurate measurement and modeling of the phase behavior a demanding task. Similarly to the previous systems, CO 2 was modeled assuming one positive association site that can only cross-associate with methanol. The obtained results are presented in Figure 12. The model shows an average absolute deviation (in bubble point pressure) from the experimental data shown in Figure 12 [57] equal to 9%. Ferreira et al., who obtained the experimental data, used the Peng-Robinson EoS and in total six binary parameters (two binary parameters per sub-binary system), while two of them (those referring to the methanol-lauric acid mixture) were adjusted to the ternary system's data [57]. According to their results, the experimental data are satisfactorily correlated. ganic acid, which renders the accurate measurement and modeling of the phase behavior a demanding task. Similarly to the previous systems, CO2 was modeled assuming one positive association site that can only cross-associate with methanol. The obtained results are presented in Figure 12. The model shows an average absolute deviation (in bubble point pressure) from the experimental data shown in Figure 12 [57] equal to 9%. Ferreira et al., who obtained the experimental data, used the Peng-Robinson EoS and in total six binary parameters (two binary parameters per sub-binary system), while two of them (those referring to the methanol-lauric acid mixture) were adjusted to the ternary system's data [57]. According to their results, the experimental data are satisfactorily correlated. Figure 12. Phase boundaries for CO2-methanol-lauric acid ternary system. Experimental data (points, [57]) and CPA predictions using the binary interaction parameters of Table 6   Finally, the CPA EoS was applied to predict the phase behavior of a CO 2 mixture with a real biodiesel sample ( Figure 13). The experimental data from Araujo et al. were used, and biodiesel was considered as a multicomponent mixture of C16:0, C18:0, C18:1, C18:2, and C18:3 methyl esters, according to the analysis mentioned in that study [58]. The CO 2 -methyl ester binary interaction parameters were adopted from Tsivintzelis et al. [7], while, since such compounds have similar molecular structure, the corresponding esterester binary interaction parameter was set equal to zero. The overall average absolute deviations of model predictions from experimental data in CO 2 mole fraction is 8% (7% for 323 K, 8% for 333 K, and 10% for 343 K). It is observed that the deviation increases with temperature; however, the performance of the model is considered rather satisfactory if we take into account that calculations are pure predictions and also the complexity of such multicomponent mixture. Finally, the CPA EoS was applied to predict the phase behavior of a CO2 mixture with a real biodiesel sample ( Figure 13). The experimental data from Araujo et al. were used, and biodiesel was considered as a multicomponent mixture of C16:0, C18:0, C18:1, C18:2, and C18:3 methyl esters, according to the analysis mentioned in that study [58]. The CO2-methyl ester binary interaction parameters were adopted from Tsivintzelis et al. [7], while, since such compounds have similar molecular structure, the corresponding esterester binary interaction parameter was set equal to zero. The overall average absolute deviations of model predictions from experimental data in CO2 mole fraction is 8% (7% for 323 K, 8% for 333 K, and 10% for 343 K). It is observed that the deviation increases with temperature; however, the performance of the model is considered rather satisfactory if we take into account that calculations are pure predictions and also the complexity of such multicomponent mixture. Figure 13. VLE of CO2-biodiesel system. Experimental data (points, [58]) and CPA predictions (lines) using the binary interaction parameters of Table 6.

Overview of Multicomponent Systems Modeled with the CPA EoS
This work completes an extensive work on modeling with the CPA EoS the phase behavior of mixtures relevant to the biodiesel process, which may contain methanol, ethanol, glycerol, fatty acids, esters of fatty acids, glycerides, and CO2. In this series of articles Figure 13. VLE of CO 2 -biodiesel system. Experimental data (points, [58]) and CPA predictions (lines) using the binary interaction parameters of Table 6.

Overview of Multicomponent Systems Modeled with the CPA EoS
This work completes an extensive work on modeling with the CPA EoS the phase behavior of mixtures relevant to the biodiesel process, which may contain methanol, ethanol, glycerol, fatty acids, esters of fatty acids, glycerides, and CO 2 . In this series of articles (this work and [7,8,30]), pure fluid parameters for such compounds are reported, along with correlations against the molecular weight or the van der Waals volume, which allow the prediction in cases of lack of experimental data. Moreover, the binary parameters for all corresponding binary systems are provided, and in many cases, correlations of such parameters with the van der Waals volume are reported. The capabilities and the limitations of the model are revealed through the prediction of the phase behavior of ternary and multicomponent systems without adjusting any parameter to such experimental data, but only using the sub-binary mixture parameters that were estimated solely from data for binary mixtures. Thus, no "fine-tuning" of those interaction parameters was performed by considering ternary or multicomponent systems data. In this sense, all calculations for ternary or multicomponent mixtures are considered as pure predictions.
In this work and in all our related previous works [7,8,30], we modeled alcohols using the 2B association scheme, water with the 4C association scheme, and glycerol with the 3 × 2B association scheme. In all cases, cross-association interactions among those hydrogen-bonding fluids were taken into consideration. Furthermore, esters of fatty acids, glycerides, and CO 2 were modeled assuming that their molecules have sites that are only able to cross-associate with other hydrogen-bonding fluids such as alcohols, water, and glycerol. Using one or two binary parameters per binary system, the CPA EoS was able to satisfactorily correlate the phase behavior of binary mixtures. The same approach to account for cross-association interactions was used in ternary or multicomponent mixtures, and it was shown that the model presents a rather satisfactory prediction ability. Numerous ternary or multicomponent mixtures were investigated. Table 7 provides an overview of multicomponent mixtures considered in this work and our previous studies [7,8,30,61] but also in other literature studies with CPA calculations [18][19][20][21][22][23][24][25].

Conclusions
In this study, the CPA pure fluid parameters for polyethylene glycols were obtained by estimating the trends of such parameters with the molecular weight. Using one temperatureindependent binary interaction parameter, the model accurately correlates the solubility of CO 2 in liquid polyethylene glycols of various molecular weights. Subsequently, the model was applied to describe the phase behavior of the CO 2glycerol system, and it was shown that accounting for cross-association interactions and using two temperature-independent binary parameters, the model accurately describes both phases in equilibrium and especially the vapor phase, which contains a very low amount of the heavy compound, rendering the accurate modeling of such systems difficult. In addition, CPA correlations revealed that the glycerol content of the vapor phase shows a minimum when plotted against pressure, which may be of interest in separation processes, similarly to the CO 2 -water system. Then, the model was applied to correlate the phase behavior of CO 2 -triglycerides in a wide temperature range in cases for which experimental data are available. Using one temperature-independent parameter, the CPA EoS is able to describe rather satisfactorily the phase behavior, including the very low vapor phase content of the heavy compound. In order to obtain reasonable compositions for the vapor phase, relatively increased but still satisfactory deviations for the liquid phase were obtained. Subsequently, the model was applied to predict the dew points of three ternary systems. Using no optimized parameters to such data, the model presents average absolute deviations in dew point pressures ranging from 6-14%. Finally, the CPA EoS rather satisfactorily predicted the VLE of a multicomponent mixture containing CO 2 and a real biodiesel sample, showing an average absolute deviation on CO 2 mole fractions equal to 8.2%.
This work completes a series of studies in modeling the phase behavior of mixtures relevant to the production of biodiesel [7,8,30]. In such studies, all pure fluid and binary parameters of most important systems are presented, and the predicting ability of the model is evaluated. It is shown that the CPA model presents a satisfactory prediction ability for such systems including mixtures of real biodiesel samples.

Supplementary Materials:
The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/molecules27185785/s1, Figure S1: CO 2 -glycerol VLE. Glycerol mole fraction in the vapor phase. Experimental data and CPA correlations using the interaction parameters of Table 4 of the main article; Figure S2: CO 2 -glycerol VLE. CO 2 mole fraction in the liquid phase. Experimental data [41,42] and CPA correlations using the interaction parameters of Table 4 of the main article; Figure S3: CO 2 -glycerol VLE. CO 2 mole fraction in the liquid phase. Experimental data [39,41,42] and CPA correlations using the interaction parameters of Table 4 of the main article; Figure S4: CO 2 -glycerol VLE. CO 2 mole fraction in the liquid phase. Experimental data [42] and CPA correlations using the interaction parameters of Table 4 of the main article; Figure S5: CO 2 -glycerol VLE. CO 2 mole fraction in the liquid phase. Experimental data [42] and CPA correlations using the interaction parameters of Table 4 of the main article; Figure S6: CO 2 -glycerol VLE. CO 2 mole fraction in the liquid phase. Experimental data [39] and CPA correlations using the interaction parameters of Table 4 of the main article; Figure S7: CO 2 -tricaprilyn VLE at 313 K. Experimental data [43,44] and CPA correlations using the binary parameter of Table 5 of the main article; Figure S8: CO 2 -tricaprilyn VLE at 393 K. Experimental data [43] and CPA correlations using the binary parameter of Table 5 of the main article; Figure S9: CO 2 -tributyrin VLE at 288 K. Experimental data [44,45] and CPA correlations using the binary parameter of Table 5 of the main article; Figure S10: CO 2 -tributyrin VLE at 303 K. Experimental data [44,45] and CPA correlations using the binary parameter of Table 5 of the main article; Figure S11: CO 2 -tributyrin VLE at 338 K. data [44,45] and CPA correlations using the binary parameter of Table 5 of the main article; Figure S12: CO 2 -trilaurin VLE and LLE at 308 K. Experimental data [48,49] and CPA correlations using the binary parameter of Table 5 of the main article; Figure S13: CO 2 -trilaurin VLE at 328 K. Experimental data [48] and CPA correlations using the binary parameter of Table 5 of the main article; Figure S14: CO 2 -trilaurin VLE at 353 K. Experimental data [47] and CPA correlations using the binary parameter of Table 5 of the main article; Figure S15: CO 2 -trimyristin VLE at 308 K. Experimental data [48] and CPA correlations using the binary parameter of Table 5 of the main article; Figure S16: CO 2 -trimyristin VLE at 328 K. Experimental data [48] and CPA correlations using the binary parameter of Table 5 of the main article; Figure S17: CO 2 -triolein VLE at 313 K (vapor phase compositions). Experimental data [53,54] and CPA correlations using the binary parameter of Table 5 of the main article; Figure S18: CO 2triolein VLE at 323 K (vapor phase compositions). Experimental data [50,53,54] and CPA correlations using the binary parameter of Table 5 of the main article; Figure S19: CO 2 -triolein VLE at 333 K. Experimental data [51][52][53] and CPA correlations using the binary parameter of Table 5 of the main article; Figure S20. CO 2 -triolein VLE at 343 K. Experimental data [52,53] and CPA correlations using the binary parameter of Table 5