Activity coefficient models for accurate prediction of adsorption azeotropes

In this study seven adsorption azeotropes involving binary systems and zeolite-based adsorbents were systematically investigated. Pure component isotherms and mixed-gas adsorption data were taken from published literature except for the benzene–propene system on silicalite, which is newly presented in this work using molecular simulations. Experimental adsorbed phase composition and total amount adsorbed of the azeotropic systems were compared with the predictions of several models including: the ideal adsorbed solution theory (IAST), the heterogeneous ideal adsorbed solution theory (HIAST) and the real adsorbed solution theory (RAST) coupled with the 1-parameter Margules (1-Margules) and the van Laar equations. In the latter two models an additional loading parameter was incorporated in the expression of the excess Gibbs energy to account for the reduced grand potential dependency of the activity coefficients in the adsorbed phase. It was found that the HIAST and RAST–1-Margules models were able to predict the azeotropic behaviour of some systems with good accuracy. However, only the RAST–van Laar model consistently showed an average relative deviation below 3% compared to experimental data for both the adsorbed phase composition and the total amount adsorbed across the systems. This modified van Laar equation is therefore preferable in those engineering applications when the location of adsorption azeotropes is required with great accuracy and when there is lack of detailed characterization of the adsorbent that is needed to carry out molecular simulations.


Introduction
Adsorption-based processes have been recently expanding their presence in many industrial applications, due to their high selectivity, high throughput and high energy efficiency, coupled with a low maintenance and ease of operation. Particularly, they play a fundamental role in catalytic reaction processes and in several separations and purifications of gases and liquids [1][2][3]. However, the complexities of solid surfaces and our inability to characterize exactly their interactions with adsorbed molecules limits our understanding of the adsorption processes. Unlike vapour-liquid equilibrium, gas adsorption is characterized by an additional degree of freedom because of the presence of the solid adsorbent, which imposes a potential field upon the adsorbed phase [4]. Besides, it has been extensively reported that experimental and simulated data for mixed-gas adsorption on energetically heterogeneous surfaces such as activated carbons and zeolites can exhibit negative deviations from ideality [5][6][7]. These deviations become sometimes large enough to form an adsorption azeotrope, which determines phenomena of composition crossovers and selectivity reversals. Several works in the literature have shown experimental data of adsorption azeotropes for mixtures of polar and non-polar components adsorbed into heterogeneous adsorbents including 13X, H-modernite, ZSM-5 and silicalite [6,[8][9][10][11]. Among the causes of non-ideal and azeotropic behaviour energetic heterogeneity of adsorption sites, steric exclusion effects and interactions between molecules in the adsorbed phase have been proposed [6,12]. Whatever the cause of non-idealities, engineers require accurate thermodynamic models to correlate and predict multicomponent adsorption equilibria in support of rigorous process simulation tools. This is particularly relevant in the design of adsorbers when it is required to locate the adsorption azeotropes with great accuracy and when there is lack of detailed characterization of the adsorbent that is needed to carry out molecular simulations. By and large there are two different cohorts of thermodynamic models used to describe multicomponent adsorption equilibria, broadly classified as extended models and solution theories. Given that the conventional extended Langmuir model is unable to predict non-idealities, a variety of multisite extended Langmuir models have been proposed to include adsorbent heterogeneity and overcome this issue. They include: the heterogeneous extended Langmuir model, the multiregion extended Langmuir model, the extended dualsite Langmuir model, the extended triple-site Langmuir model, and the multiregion multisite extended Langmuir model [13][14][15][16][17][18]. All the previous models have the clear advantage of using only pure component adsorption parameters and have shown some degrees of improvement in the prediction of non-ideal mixtures. However, the results were generally poor, with azeotropic behaviours predicted only qualitatively, with respect to both the adsorbed phase composition and the total amount adsorbed. Among the purely predictive solution theories, the ideal adsorbed solution theory (IAST) suffers from the intrinsic inability to predict non-idealities [19] while the heterogeneous ideal adsorbed solution theory (HIAST) and the non-ideal adsorbed solution theory (NIAST) often improved the predictions of nonideal and azeotropic mixtures to a limited extent [20][21][22].
A substantial improvement in the prediction of non-idealities came with the vacancy solution theory (VST), which included the vacancies as a component [23], and the real adsorbed solution theory (RAST), which introduced activity coefficients in the adsorbed phase, thus mimicking the γ-φ approach used in VLE. Unfortunately, early and recent attempts rather blindly applied the RAST coupling it with conventional VLE activity coefficient models (Wilson, UNI-QUAC, NRTL) but neglected the variation of the excess Gibbs energy with the reduced grand potential [24,25]. Accounting for this dependency is fundamental because experimental data are usually obtained at constant pressure, while the reduced grand potential varies with composition. For instance, Talu and Zwiebel [11] have shown that neglecting this contribution introduces errors larger than 30% in the prediction of the total adsorbed amount of adsorption azeotropes. Subsequent works codified the model requirements for thermodynamic consistency and reduction to an ideal adsorbed solution at the limit of zero loading [4,7]. In particular, Valenzuela and Myers [26] introduced for the first time the one-parameter Margules equation with an exponential dependency upon the reduced grand potential, which agrees with experiments and molecular simulations from zero loading up to saturation. As the regular-solution theory failed to agree with experiments [5], Siperstein and Myers [6] proposed the ABC equation in which the quadratic excess Gibbs free energy had a linear dependence with temperature. The authors claimed that the three parameters accounting for excess enthalpy, excess entropy and excess reciprocal loading were necessary to describe non-idealities of mixtures as a function of temperature, composition and loading. After regressing the parameters based on the experimental binary system data, the average relative deviation in the calculated composition was 5% for the six investigated systems, with larger errors for the highly non-ideal and azeotropic mixtures.
It is clear from the previous studies that none of the thermodynamic models proposed so far was able to predict the azeotropic behaviour of adsorption mixtures with great accuracy. The results have also been presented in a fragmentary way as the authors often assessed their models with respect to a single specific variable, such as the adsorbed phase composition, the total amount adsorbed or the total pressure. This work is aimed to systematically investigate seven adsorption azeotropes involving binary systems and zeolitebased adsorbents in order to predict both the adsorbed phase composition and the total amount adsorbed with an accuracy below 3%, which would be required for the correct design of separation units [27]. Table 1 lists the investigated adsorption azeotropes along with their operating conditions and relevant industrial applications. Pure component isotherms and mixed-gas adsorption data were taken from published literature except for the benzene-propene system on silicalite, which is newly presented in this work using molecular simulations. The results of binary systems were compared and fully validated with the data reported by Ban et al. [8]. It should be noted that molecular simulations might be able to predict in-silico many more adsorption azeotropes than the one presented in this work. However, this study would limit the analysis to experimental azeotropes and the simulated azeotrope benzene-propene on silicalite, given its industrial relevance [28]. Details and results of the molecular simulations carried out to generate pure component isotherms and mixed-gas adsorption data for the system benzene-propene on silicalite can be found in the Appendix.
Aside the known IAST, HIAST and RAST-1-Margules models, a RAST model coupled with the van Laar equation has been proposed to tackle the stringent accuracy target in the prediction of adsorption azeotropes. The resulting threeparameter expression for the excess Gibbs energy is asymmetrical in composition and has an exponential dependency upon the reduced grand potential, thus being thermodynamically consistent [4,7]. Another contribution of the current investigation is to discuss the behaviour of adsorption azeotropes in relation to the infinite dilution activity coefficients and their visualization on the reduced grand potentialcomposition diagram at constant temperature and pressure, therefore confirming and extending the conditions for azeotrope existence defined by Siperstein [12].

Thermodynamic models
As previously mentioned, in this study experimental adsorbed phase composition and total amount adsorbed of the azeotropic systems were compared with the predictions of several models including the IAST, the HIAST and the RAST coupled with the 1-Margules and the van Laar equations. In all the models the pure component isotherms were described by the dual-site Langmuir model, for which a heterogeneous adsorbent consists of two different sites of different energies and the overall isotherm of a gas is obtained by adding the contributions of each site [29]. The mathematical expression is the following: Adsorption isotherm parameters were fitted using Origin software [30] with values summarized in Table 2. For 13X and H-modernite model parameters were fitted to experimental isotherm data while for NaX they were fitted to model isotherm data obtained from a modified virial equation [6], given the complexity to extract experimental data from a log-log plot. For silicalite model parameters
In the IAST the gas phase is ideal and in equilibrium with the adsorbed phase, which is assumed an ideal solution following the Raoult's law: where P i 0 is the surface pressure of component i at the temperature and reduced grand potential of the mixture. The IAST is based on the iso-reduced-grand-potential condition stating that each component in the adsorbed phase has the same reduced grand potential at equilibrium. This last condition is expressed as follows: Finally, the total amount adsorbed is calculated by: where q i 0 is the amount adsorbed of component i calculated at P i 0 . Considering the dual-site Langmuir model the reduced grand potential is analytically obtained solving the integral in Eq. (3) as: The large success of the IAST relies on its simple thermodynamic framework and its ability to predict multicomponent adsorption equilibria from pure component adsorption isotherms.

Heterogeneous ideal adsorbed solution theory
In order to predict the behaviour of non-ideal adsorption mixtures, Valenzuela et al. [22] proposed a heterogeneous ideal adsorbed solution theory (HIAST) by incorporating the effect of energetic heterogeneity into the IAST. The resulting theory behaves ideally on a particular adsorption site, but energetic heterogeneity causes a segregation in the composition of the adsorbed phase. The authors used the Langmuir model and the discrete binomial energy distribution with a perfect positive correlation of sites. As a result, the HIAST always improved the predictions, but they were still not in quantitative agreement with the experiments. However, it was also reported that there are cases in which the predictability of the HIAST can be inferior to that of the IAST if the energy distribution parameters are not properly chosen [39].
To avoid the uncertainties in the energy distributions, in this study an approach similar to that proposed by Myers [20] was used to define the HIAST on the basis of the dualsite Langmuir model. Since the two energetic adsorption sites are independent of one another, the IAST is applied to each single site, where the equality of the reduced grand potential for each pure component requires that: As the saturation capacities are different for the components of the binary systems, both perfect positive and perfect negative correlations were tested on each of the two Langmuir sites, and the best correlation was chosen by comparing the predictions with experimental mixed-gas adsorption data [18]. Similarly to the IAST, the HIAST has the great advantage of requiring only pure component adsorption isotherms to predict multicomponent adsorption equilibria. However, to cope with the poor model predictions, the same Valenzuela et al. [22] proposed to incorporate into the HIAST a binary adsorption parameter whose value needs to be determined from experimental data.

Real adsorbed solution theory
The real adsorbed solution theory (RAST) is derived from the IAST by introducing the activity coefficients in the adsorbed phase with the purpose of predicting the behaviour of non-ideal adsorption mixtures. As a consequence, the equation of equilibrium for multicomponent adsorption becomes [19]: The activity coefficients are related to the excess Gibbs energy of the adsorbed solution as follows: The total amount adsorbed is obtained by the ideal contribution plus the excess contribution: where the excess contribution is defined as: Therefore, given the behaviour of pure adsorbates, multicomponent mixtures can be fully characterized knowing an expression for the excess Gibbs energy as function of temperature, composition and reduced grand potential.

1-parameter Margules equation
The simplest thermodynamically consistent model proposed to describe the excess Gibbs energy is the 1-parameter or 2-suffix Margules equation with an exponential dependency upon the reduced grand potential [5,7]. The resulting twoparameter expression of the excess Gibbs energy provides a complete description of isothermal, mixed-gas adsorption equilibrium over the entire range of surface coverage. Even the three-parameter ABC equation introduced by Siperstein and Myers [6] reduces to the 1-parameter Margules equation in case of isothermal conditions. Moreover, this equation has been successfully applied to the common tangent plane approach to validate the presence of adsorption azeotropes [40].
The mathematical expression of the excess Gibbs energy for a binary system is the following: where A and B are two equation parameters. The excess contribution to the total amount adsorbed is obtained by inserting Eq. (16) into Eq. (15): The activity coefficients of the components are given by substituting Eq. (16) into Eq. (12):

van Laar equation
In order to refine the model predictions of adsorption azeotropes, the RAST can be coupled with a modified threeparameter van Laar equation for which the excess Gibbs energy has an asymmetrical composition dependency and the same exponential dependency upon the reduced grand potential: Similar to the 1-Margules equation, the excess contribution to the total amount adsorbed and the activity coefficients can be evaluated as follows: It should be noted that the van Laar equation reduces to the 1-Margules equation in the special case of A 12 = A 21 .
In the previous 1-Margules and van Laar equations the activity coefficients satisfy the requirements of being unity at the limit of composition approaching unity (x i → 1) or loading approaching zero (ψ → 0). At high loading (ψ → ∞) the activity coefficients approach a constant value corresponding to saturation. From the thermodynamics of liquid mixtures, the parameters A, A 12 and A 21 are empirical constants with units of energy and could be assumed, at constant composition, temperature independent (athermal solutions) or proportional to the reciprocal of absolute temperature (regular solutions) [41].

Estimation of activity coefficient model parameters
In this study the activity coefficient model parameters, namely A and B for the 1-Margules equation and A 12 , A 21 and B for the van Laar equation, were estimated based on experimental data of binary systems at fixed temperature and pressure ( Table 1). The parameter estimations were carried out using gPROMS software [42] by minimizing the sum of relative errors in the calculated adsorbed phase composition and total amount adsorbed, resulting in the following objective function: A constant variance of 10 -5 was set for x 1 and q TOT whereas absolute and relative tolerances were both fixed to 10 -6 . Values of these parameters derived from experimental data for the adsorption azeotropes are reported in Table 3 along with their standard deviations.
For the system iC 4 H 10 -C 2 H 4 -13X the parameter B was considered temperature independent, being the loading parameter of excess Gibbs energy. With reference to the system CO 2 -C 3 H 8 -NaX, Siperstein and Myers [6] obtained comparable values for the regressed parameters A and B of the 1-Margules equation by minimizing the error in the total pressure and selectivity. From Table 3 it should be noted that parameters A, A 12 and A 21 showed all negative values, thus confirming the negative deviations from Raoult's law, while parameter B ranged 0.05-0.93 across all the systems. Although all the model parameters exhibited low standard deviations in the interval 10 -7 -10 -4 , it is not possible to specify a unique set of parameters as there will always be an uncertainty, which ultimately depends on the accuracy of the experimental data. In addition, the parameters correlate in the neighbourhood of the optimal point. Figure 1 shows the confidence ellipses of the van Laar parameters A 12 and A 21 for the system C 6 H 6 -C 3 H 6 on silicalite. These are obtained from the eigenvalues and eigenvectors of the covariance matrix of the parameters. In particular, the direction vectors of the ellipse axes are the eigenvectors, whereas the length of the vectors corresponds to the eigenvalues [43]. The elliptical form of the deviation plot indicates that the two parameters are partially coupled. For instance, the smallest ellipse identifies a region in the A 12 -A 21 plane wherein any point gives a set of parameters that reproduces at least 90% of the experimental data within the set variances. It is expected that the ellipses become smaller as the quality and/or quantity of experimental data rises and as the suitability of the model for the excess Gibbs energy increases [41].
To validate the effectiveness of the parameter estimations, Fig. 2 shows a comparison of the experimental and calculated adsorbed phase composition and total amount adsorbed using the van Laar equation for the system CO 2 (1)-C 3 H 8 (2) on NaX. As the points lie very close to the 45° reference line, the model is able to predict the multicomponent adsorption equilibria with great accuracy.

Model predictions of adsorption azeotropes
After estimating the activity coefficient model parameters, the RAST coupled with the 1-Margules and van Laar equations as well as the IAST and the HIAST were assessed on the prediction of adsorption equilibria of the seven azeotropic binary systems listed in Table 1. Model predictions were compared with experimental adsorbed phase composition and total amount adsorbed by means of average relative deviations (ARD), defined as: where ARDs are expressed in percentage (%) and n is the number of experimental data points. Table 4 reports the ARDs for the adsorbed phase composition and total amount adsorbed between experimental data and predictions from the four models assessed in this study.
Overall, the thermodynamic models can be ranked in order of increasing accuracy as IAST < HIAST < RAST-1-Margules < RAST-van Laar, but there are some exceptions in which the HIAST was worse than the IAST or better than the RAST in the predictions of the adsorbed phase composition, the total amount adsorbed or both. However, with the exclusion of x 1 for the C 3 H 8 (1)-H 2 S(2)-H-modernite system, only the RAST-van Laar model consistently showed an average relative deviation below 3% compared to experimental data for both the adsorbed phase composition and the total amount adsorbed across the systems. Figures 3-6 show isothermal, isobaric x-y diagrams and q TOT -y diagrams including experimental data and different model predictions for selected adsorption azeotropes.
In the system iC 4 H 10 (1)-C 2 H 4 (2) on 13X the pure component isotherms cross each other at low surface coverage leading to an azeotropic behaviour as seen in the binary system experimental data [10]. In Fig. 3 all models apart from the IAST captured the existence of the azeotrope at 298 K, but the selectivity reversal was predicted at different compositions, with the prediction of the RAST-van Laar model being the closest to the experimental data. Since the results of the HIAST and RAST-1-Margules models were not satisfactory, the asymmetrical three-parameter van Laar equation is required to accurately predict both the adsorbed phase composition and the total amount adsorbed (Fig. 3).
A similar trend was observed for the system C 2 H 4 (1)-CO 2 (2) on 13X, as shown in Fig. 4. Also in this case the azeotrope was accurately located only by the RAST-van Laar model while HIAST and RAST-1-Margules models reported a large deviation in the adsorbed phase composition in the range of 15-22%. It should be noted that the deviation in the total amount adsorbed predicted by the HIAST was even larger than the IAST (Fig. 4b). This discrepancy, also found in the system iC 4 H 10 (1)-C 2 H 4 (2) on 13X at 323 K (Table 4), highlights the fact that energetic heterogeneity is insufficient to describe these adsorption azeotropes. In complex heterogeneous adsorbents such as 13X zeolite, in fact, the effect of partial steric exclusion of larger molecules from pores accessible to smaller molecules can become prominent, thus increasing the negative deviations from Raoult's law and leading to adsorption azeotropes [6,22].
The system C 3 H 8 (1)-CO 2 (2) on H-modernite analysed in Fig. 5 was extensively studied due to its highly non-ideal behaviour [11]. Model predictions provided by the HIAST greatly improved those of the IAST with ARDs in the magnitude of 5-8%. The discrepancies can be explained by Table 4 Average relative deviations for the adsorbed phase composition and total amount adsorbed between binary system data and predictions from different models Exp. experimental data, Mol. sim. molecular simulation data obtained in this study a Outlier data point with y 1 = 0.1168 was excluded from the calculations b Outlier data points with y 1 < 0.05 were excluded from the calculations specific electrostatic interactions as well as steric exclusion effects [13,22]. Once again, the lack of detailed understanding of adsorption azeotropic behaviour can be compensated with activity coefficient models able to reduce the average relative deviations down to 1.5-5.8% and 1.3-1.6% using the 1-Margules and the van Laar equations, respectively. In reference to the same H-modernite adsorbent, H 2 S is adsorbed much more strongly than C 3 H 8 due to its permanent dipole that generates strong electric field-dipole interactions. This leads to large ARDs in the adsorbed phase composition for all models, including the RAST-van Laar model (Table 4).
For the system CO 2 (1)-C 3 H 8 (2) on NaX reported by Siperstein and Myers [6], it is worth noticing that it is the least non-ideal system investigated in this study. In fact, the IAST exhibited ARDs of around 8%, which were reduced to less than 3% using the HIAST and the RAST models, as reported in Table 4. In particular, the parameter estimations of the 1-Margules and van Laar equations (Table 3) produced the same value for B while the values of A 12 and A 21 were very close to the value of A. Therefore, the two RAST models showed virtually identical ARDs in both the adsorbed phase composition and the total amount adsorbed.  The last azeotropic system investigated in this study was C 6 H 6 (1)-C 3 H 6 (2) on silicalite with results shown in Fig. 6. Details of the molecular simulations carried out to generate pure component isotherms and mixed-gas adsorption data are reported in the Appendix. In this case, the model predictions of the HIAST resulted even marginally superior to those of the RAST models, with ARDs below 3%. This occurred because silicalite adsorbent is characterized by three adsorbing sites: straight channels, sinusoidal channels and their intersections. However, for this system the use of the dual-site Langmuir model is justified as the straight and sinusoidal channels have equal adsorption strengths, so that molecules located at these two channel interior positions compete all together with the molecules located at the channel intersections [44]. A graphical representation of the pore space of silicalite is reported in the Appendix. The benzene has a strong affinity with the intersections, given the higher available space compared to the channels. The strong affinity of benzene with the intersections leads to its favourable adsorption over propene up to y 1 ≈ 0.5. However, at higher compositions, the fugacity of benzene at 373 K and 100 kPa is not  Fig. 6 (a) x-y diagram and (b) q TOT -y diagram showing molecular simulation data and different model predictions for the system C 6 H 6 (1)-C 3 H 6 (2) on silicalite at 373 K and 100.0 kPa high enough to force the benzene molecules in the channels, which are free to be filled by propene, thus causing the reversal of selectivity for y 1 > 0.5. It can be concluded that energetic heterogeneity is the main responsible for the formation of this adsorption azeotrope.
Eventually, Table 5 reports the average relative deviations between experimental data and model predictions from different works in the literature in terms of both adsorbed phase composition and total amount adsorbed. Compared to the extended models and the adsorbed solution theories, the RAST model coupled with the modified van Laar equation proposed in this study showed by far the lowest deviations across all the azeotropic systems. However, it should be remembered that the UNIQUACbased spreading-pressure-dependent (SPD) model utilizing the binary data to evaluate the necessary parameters was also able to predict adsorption azeotropes on H-modernite with great accuracy [11]. Up to now, only adsorbed solution theory models utilizing experimental data to evaluate the necessary parameters and accounting for the reduced grand potential dependency of the activity coefficients have been successful to accurately predict binary adsorption azeotropes.

Implications of azeotropic behaviour
Deviations of binary adsorbed mixtures from ideality can be visualized by plotting activity coefficients and excess functions in the adsorbed phase composition domain. Isothermal, isobaric activity coefficients and excess Gibbs energy are shown in Fig. 7 using the 1-parameter Margules and van Laar equations for the system C 2 H 4 (1)-CO 2 (2) on 13X. As expected, the excess Gibbs energy is largely negative and the system exhibits negative deviations from Raoult's law, leading to azeotropic behaviour. The activity coefficients are less than unity in the entire composition domain, with the activity coefficient of C 2 H 4 at infinite dilution (γ 1 ∞ ) being equal to 0.12 using the more accurate van Laar equation (Fig. 7a). Given the composition dependency at constant reduced grand potential, the excess Gibbs energy is symmetrical with the 1-Margules equation and asymmetrical with the van Laar equation. However, under the isobaric conditions imposed on Fig. 7b, both curves show their minima displaced towards CO 2 , the component with the higher loading.
With the aim of determining the conditions where azeotropic behaviour is observed, Siperstein [12] showed that the existence of an adsorption azeotrope can be defined from relations between the infinite dilution activity coefficients and the surface pressure of the pure components at the same temperature and reduced grand potential of the mixture. In particular, these relationships can be obtained from the slope of the bubble point pressure curve while it approaches the pure component. However, the azeotropic behaviour can also be determined under the more common experimental conditions of constant temperature and pressure. As the ψ vs P function has a trend similar to the T boil vs P function typical of VLE, negative deviations from Raoult's law imply a maximum-reduced-grand-potential azeotrope at constant temperature and pressure. Figure 8 shows the ψ-x-y diagram using the IAST and RAST-van Laar models for the system C 2 H 4 (1)-CO 2 (2) on 13X.
The van Laar equation correctly predicts the existence of the adsorption azeotrope at the mole fraction of C 2 H 4 of around 0.1 and the maximum reduced grand potential of around 13.4 mol kg -1 . Since CO 2 (2) is the strongly adsorbed component on 13X adsorbent, it follows that P 2 0 < P 1 0 . Focusing on the slope of the "equivalent" bubble point curve as it approaches pure CO 2 , a negative azeotrope is observed when: In addition, according to the findings reported by Siperstein [12], the condition for observing a negative azeotrope at a given temperature and reduced grand potential is given by: Therefore, the system analysed in Fig. 8 confirms the condition expressed by Eq. (28) as γ 1 ∞ = 0.12 < P 2 0 / P 1 0 = 0.41. This condition is also confirmed by all the other adsorption azeotropes investigated in this study. Fig. 7 (a) Activity coefficients and (b) excess Gibbs energy using the 1-parameter Margules and van Laar equations for the system C 2 H 4 (1)-CO 2 (2) on 13X at 298 K and 137.8 kPa Fig. 8 ψ-x-y diagram using the IAST and RAST -van Laar models for the system C 2 H 4 (1)-CO 2 (2) on 13X at 298 K and 137.8 kPa

Conclusions
The prediction of adsorption equilibria of multicomponent non-ideal mixtures is a challenging task for the design of gas separation adsorbers, especially in case of azeotropic behaviour where phenomena of composition crossovers and selectivity reversals occur. In this work, seven adsorption azeotropes involving binary systems and zeolite-based adsorbents were systematically investigated with the aim of predicting both adsorbed phase compositions and total amounts adsorbed with great accuracy. Pure component isotherms and mixed-gas adsorption data were taken from published literature except for the benzene-propene system on silicalite, which was newly presented in this work using molecular simulations. Several models were tested including the ideal adsorbed solution theory (IAST), the heterogeneous ideal adsorbed solution theory (HIAST) and the real adsorbed solution theory (RAST) coupled with the 1-parameter Margules and the van Laar equations. The activity coefficient models were thermodynamically consistent as an exponential dependency upon the reduced grand potential was incorporated in the expression of the excess Gibbs energy of the adsorbed solution.
After carrying out appropriate parameter estimations on experimental binary system data, it was found that the RAST-van Laar model consistently showed an average relative deviation below 3% compared to the experimental adsorbed phase composition and total amount adsorbed across the systems. Up to now, the modified van Laar equation of this study is the most accurate thermodynamic framework to comprehensively predict binary adsorption azeotropes among the extended models and the adsorbed solution theories proposed in the literature. Therefore, it can serve to correlate and predict non-ideal multicomponent adsorption equilibria in support of rigorous process simulation tools. This is particularly relevant in the design of adsorbers when the location of adsorption azeotropes is required with great accuracy and when there is lack of detailed characterization of the adsorbent that is needed to carry out molecular simulations.
The conditions for observing azeotropic behaviour proposed by Siperstein [12] and based on relations between the infinite dilution activity coefficients and the surface pressure of the pure components were confirmed in this study. Moreover, under the conditions of constant temperature and pressure it was shown that the existence of a negative adsorption azeotrope can be visualized by a maximum in the reduced grand potential function along the adsorbed phase and gas phase composition domains.

Appendix
In this work, Grand Canonical Monte Carlo (GCMC) simulations were carried out to investigate the adsorption of pure benzene and propene on MFI-type silicalite at 373 K, and their binary system which exhibits an azeotrope at 373 K and 100.0 kPa, as shown by Ban et al. [8]. The ortho-MFI configuration was used to model the rigid solid framework, whose atom coordinates have been reported by Van Koningsveld et al. [45]. The MFI structure consists of straight and sinusoidal channels which cross at intersections. The schematic representation of the ortho-MFI framework is shown in Fig. 9. The pore network representation was derived using PoreBlazer [46,47].
The rigid framework was modelled using chargeless silicon atoms and charged oxygen atoms. The interaction between adsorbates and MFI framework was concentrated on the oxygen atoms of the MFI, leaving the silicon atoms as inert [48,49]. The selected force field used a united atom approach as that presented by Ban et al. [8]. The Fig. 9 Highlight of the accessible pore network of the ortho-MFI: in yellow the accessible pore network; in black the solid framework; in dashed red the sinusoidal channel; in solid light blue line the straight channel adsorbate-zeolite atoms interactions were modelled with Lennard-Jones (L-J) potentials where the cut-off length of the intermolecular interactions was set to 12 Å. Jorgensen mixing rules were used for non-identical united atoms interactions. The Coulomb interactions between the charged sites of the benzene and the oxygen atoms of the MFI were modelled by means of Ewald summation. The benzene was modelled with the nine-site model described by Siepmann and co-workers [50,51], which considers six CH chargeless sites of the benzene ring and three extra charged sites to correctly represent the quadrupole moment generated by the π bonding system. Lennard-Jones parameters for guest-guest and guest-host interactions along with details of the MFI framework are given in Table 6.
RASPA software was used to run the molecular simulations [52]. The number of cycles was changed according to each simulation point for both pure component isotherms at 373 K and binary system at 373 K and 100.0 kPa to ensure that the simulations reached a reliable statistical average. This was monitored by assessing the number of molecules present in each simulation cycle. The highest loading of 9.8 molecules per unit cell was achieved for propene adsorption at 373 K and 250.0 kPa. The inaccessible pockets of the MFI framework were blocked with inert rigid spheres to avoid unnecessary attempts of insertion and improve the convergence performance. The results of the molecular simulations for pure component isotherms and mixed-gas adsorption data are reported in Table 7. It should be highlighted that the adsorption of the C 6 H 6 (1)-C 3 H 6 (2) system at 373 K and 100.0 kPa was fully validated with the results presented by Ban et al. [8] on the same silicalite.  Table 7 Pure component isotherms and mixed-gas adsorption data for the system C 6 H 6 (1)-C 3 H 6 (2) on silicalite at 373 K. Binary system is at 100.0 kPa C 6 H 6 C 3 H 6 Binary system P (kPa) q (mol kg -1 ) P (kPa) q (mol kg -1 ) y 1 (-) x 1 (-) q TOT (mol kg -1 ) 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/.