Molecular Thermodynamics of Adsorption using Discrete-Potential Systems

Molecular Thermodynamics of Adsorption using Discrete-Potential Systems — A molecular thermodynamics approach has been developed in order to describe the adsorption of fluids onto solid surfaces based on the use of discrete-potential fluid models. Using perturbation theories for fluids such as the Statistical Associating Fluid Theory (SAFT) and the Discrete Potential Theory (DPT), in combination with molecular simulation, we have formulated a two-dimensional approach to describe systems of interest for the oil industry, such as adsorption isotherms of carbon dioxide and asphaltenes. 330 Oil & Gas Science and Technology – Rev. IFP, Vol. 63 (2008), No. 3


INTRODUCTION
Over the last decade, the Statistical Associating Fluid Theory (SAFT) [1,2], based on the thermodynamic perturbation theory for fluids with highly anisotropic interactions developed in References [3][4][5][6][7][8], has been extensively used to calculate the phase behaviour of associating and nonassociating fluid systems.The SAFT approach has been applied to a wide variety of complex and industrially relevant systems, and due to its predictive accuracy and versatility, this approach is superseding well-established chemical engineering equations of state.Over the years several variations of the theory have been proposed that have improved the predicting power of the method, including a SAFT-DFT approach to study the liquid-vapour interface of associating fluids [9].See References [10,11] for reviews on the SAFT approach.
On the other hand, systems of spherical particles interacting via Discrete Potentials (DP) have been useful as primitive models of liquids and colloids.Arbitrary DP can be built via a sequence of Square-Well (SW) and Square-Shoulder (SS) potentials, and can be used to obtain by computer simulation static and dynamic properties of fluids interacting via a discretised Lennard-Jones potential, that reproduce the properties of a proper continuous LJ-potential system [12], and also the properties of DP systems under confinement [13,14].A perturbation theory for 3D DP systems has been developed [15][16][17] and applied to the study of the phase diagram of systems with liquid-liquid equilibrium [18,19].The case of two-dimensional DP systems has its specific interest.The behaviour of colloidal particles located at the air-water interface has been succesfully modelled using two-dimensional discrete-potential systems [20,21].
Charge-stabilised colloidal particles can be trapped at the air-water interface by surface tension effects, where part of the particle is submerged in the water but most of it remains exposed to the air.Stabilisation of colloidal particles both in soap-froths and clusters at interparticle distances of the order of microns cannot be interpreted solely in terms of a pair potential with short-ranged repulsion and van der Waals attraction, as in the case of molecular fluids, but also a secondary attractive well plus a long-range repulsive barrier are required [22].By using two-dimensional particles interacting via a discrete potential as shown in Figure 1, it has been possible to reproduce by Monte Carlo simulation several pattern formations observed in real colloids at the airwater interface [20].Even more, at high densities an interesting phase diagram has been reported with a liquid-solid transition involving two different types of crystal phases.Instead of a liquid-liquid transition that has been evidenced in 3D DP systems, a metastable hexactic transition appears between a liquid and the two types of crystals.This suggests that a liquid-hexactic instead of a liquid-liquid transition could be present in 2D DP systems [21].
In the case of molecules adsorbed on a substrate, it is also known from quantum mechanical calculations [23] that the substrate modifies the particle-particle potential.Instead of the LJ type potential that exists between particles from the bulk, now we have for adsorbed particles a LJ potential with an extra repulsive barrier, very similar to what happens in colloidal particles at the air-water interface.
The Statistical Associating Fluid Theory for potentials of variable range (SAFT-VR) [24,25] has been extended to treat adsorption of chain molecules, based on a quasitwo-dimensional approach to describe the adsorbed fluid [26,27].The theory developed is as general as the SAFT-VR approach, in the sense that different types of particle-particle and particle-wall potential models can be used.In this paper we consider the adsorption of non-associating as well associating chain molecules formed by spherical segments interacting via a square-well potential.These molecules adsorb onto a uniform structureless wall, and the segment-wall potential is also described by a SW interaction.Theory is then applied to model adsorption isotherms of carbon dioxide on dry activated carbon, and also to predict the behaviour of asphaltene adsorption.Although the theory is based on describing the particle-particle and particle-wall interactions as square-well potentials, the theory can be extended to general discrete potentials, as we also discuss in this paper.

PREDICTING ADSORPTION ISOTHERMS WITH A TWO-DIMENSIONAL SAFT-VR APPROACH
The model that we have studied is a single-component fluid in the presence of a uniform wall.The fluid consists of N spherical particles with a diameter σ.We represent by u pw (z; w , λ w ) the interaction exterted by the wall on a particle, where z is the perpendicular distance of the particles from the wall, w is the energy depth of the potential, and λ w σ defines the range of the potential.We can describe this sytem as composed of two subsystems: a) the fluid whose particles are near to the wall, i.e., when z ≤ λσ, and that from now on we will describe as the adsorbed fluid, and b) the fluid whose particles are far from the wall, i.e., when z > λσ, that we call the bulk fluid.The length scale that characterises the adsorbed fluid is given by λ w σ.This approach is formally valid if we do not take into account the interface between the adsorbed and bulk fluids.
Due to the presence of the wall, the pair interaction between particles is different for the adsorbed and bulk phases [23].We will denote by u pp (r; , λ) and u ads pp (r; ads , λ ads ) the pair potential for particles in the bulk and adsorbed phases, respectively, where and λ (or ads and λ ads ) are parameters that describe the energy depth and the range of the potential for the bulk and adsorbed particles, respectively.
Introducing the density of particles of the fluid, ρ, as a function of z, we can describe the amount of adsorbed particles using the coverage Γ, defined as where ρ b is the density of bulk particles, i.e., ρ(z → ∞).
Once we have defined the length scale of the adsorbed fluid as λ w σ, then we can rewrite Equation 1as The integral of ρ(z) in the right-hand term of Equation 2 can be identified as a surface or adsorbed density, Equation 3 can be identified with a two-dimensional density, in such a way that the properties of the adsorbed fluid are identified with a 2D fluid.Since in thermodynamic equilibrium, the chemical potentials of the adsorbed (μ ads ) and bulk (μ b ) phases are the same, for a given temperature T and bulk pressure P, then ρ ads can be obtained by solving the previous equation.The formal identification of the adsorbed fluid as a 2D system can be described as a decoupling of the pair of coordinates x, y from the coordinate z for each adsorbed particle.In a similar way to the definition of the adsorbed density, we can define a 2D potential φ(x, y; 2D , λ 2D ) from the actual pair potential u ads pp (r; ads , λ ads ) by integrating over the Z-coordinate, φ(x, y; 2D , λ 2D ) = dzu ads pp (r; ads , λ ads ) According to the definitions introduced previously, the partition function of the adsorbed fluid is given as where V ads is the volume containing the adsorbed fluid, Λ is the de Broglie thermal wavelength Λ = h/(2πmkT ) 1/2 in terms of the Planck's constant h and the mass m of the particles, and Q ads is the configurational partition function given by In a previous communication Z ads was written incorrectly (Equation 15 in [27]), and Equation 6 is the right expression.By introducing the adsorption area S , we have that V ads = S λ w σ, and Q ads can be rewritten as where Q 1D and Q 2D are the one-and two-dimensional configurational partition functions, respectively, given by and In Equation 9we have applied the Mean Value Theorem with z * as the specific value of the distance to the wall implied by this theorem.
The Helmoltz free energy of the adsorbed fluid is given by In this equation, A 2D is the Helmholtz free energy of a twodimensional fluid interacting via the potential φ(x, y), and that can be modelled by perturbation theory, where A HD is the excess Helmholtz free energy for a harddisk fluid, β = 1/kT , and a 2D 1 and a 2D 2 are the first two terms of the perturbation expansion in 2D.It is also relevant to stress that the expression given for Q 1D is exact, and this information appears in Equation 11 as an external field that essentially can be used to define the energy parameter of the wall, w .For the simple case where the wall-particle interaction is given by a square-well interaction of range λ w σ and energy depth w , then we have that u pw (z * ) = − w .
For the bulk fluid, we consider an analogous perturbation expression, where A HS is the excess Helmholtz free energy for a hardsphere fluid.Obtaining the chemical potentials μ ads and μ b from Equations 11 and 13, we can rewrite Equation 4 in the following way, where μ 3D and μ 2D are the chemical potentials for 3D and 2D fluids, respectively, and μ w is the contribution to the chemical potential due to the wall,

Adsorption of Associating Fluids
The theory presented in the last section has been extended to the general case of adsorption of homonuclear and heteronuclear associating chain molecules, within the Statistical Associating Fluid Theory for potentials of variable range (SAFT-VR) framework [24,25].In the simplest case of homonuclear molecules, we consider that the fluid is formed by N molecules, and each molecule has m spherical segments of diameter σ.The particle-particle interaction is given now by two different contributions, an isotropic square-well interaction and a set of anisotropic site-site interactions that are also described by square-well potentials.We assume that the particle-wall interaction is described by a square-well potential.The thermodynamic equilibrium between the bulk and adsorbed fluid is obtained from the equality of the chemical potentials, as given by Equation 14; however, the chemical potentials are derived from the SAFT-VR expressions for the Helmholtz free energies for 3D and 2D chain-molecule fluids.The Helmholtz free energy for chain molecules in the SAFT-VR approach is written as separate terms, which account for the ideal, monomer, chain and association contributions to the free energy.The expressions for the 3D system can be found in Reference [24].For the 2D case we have that The expressions for the non-associating systems have been reported recently elsewhere [27] and we summarise them in Appendix I, including the association free energy term.

Extension to DP Systems
The 2D approach presented in this section can be easily extended to the case when particles interact via discrete potentials, using the DPT approach presented in References [15,16].Here we describe the main details related to this extension.We first consider the case of a system of N disk particles of diameter σ contained in an area S .Particles interact via an arbitrary discrete potential given by Discrete pair interaction used to simulate the behaviour of colloidal particles at the air/water interface.This model contains a secondary minimum and a secondary maximum.
where u HD is the hard-disks repulsive contribution, q is the number of discrete steps, and the potential φ i is defined by where i could be positive or negative, and λ 0 = 1.The set of steps φ i defines the shape of a general potential u d , and by a proper selection of repulsive and attractive steps it is possible to obtain approximated versions of continuous model potentials, such as the Lennard-Jones interaction, or more complex potentials in the case of colloidal particles, i.e., the interaction for colloidal particles at the air-water interface depicted in Figure 1.In the case of adsorption phenomena, we know that the particle-particle interaction is modified by the presence of the substrate and that the LJ potential is modified near the wall [23].The actual pair potential has a repulsive barrier, similar to what happens in colloidal particles.This effect can be modelled by a pair potential for the adsorbed particles similar to the potential given in Figure 1.
The Helmholtz free energy is obtained from a hightemperature expansion to second order, similarly to the 3D DP theory [15,16]: where A HD is the free energy for the hard-disk reference fluid, and a S 1 and a S 2 are the first-and second-order perturbation terms for a square-well (ε i < 0) or square-shoulder (ε i > 0) fluid, respectively.Then, the thermodynamic properties of a general discrete-potential system are given by the properties of elementary systems interacting via the SW and SS potentials.
Since the properties of the SS fluid can be given in terms of the SW fluid, by taking care of the sign of ε i , i.e., and then the 2D DPT approach, as in the case of the 3D case, depends only on the properties of a 2D SW fluid, and the 2D equation of state presented in Appendix I can be used for this purpose.Equation 20 is an exact relation, whereas Equation 21 is valid within the local compressibility approximation [28].
The extension of the theory to chain molecules is straightforward, since the SAFT approach requires knowing the Helmholtz free energy of the segments that form the chain molecules, and the contact value of the corresponding radial distribution function for segments.This last quantity can be obtained from perturbation theory using Equation 26 and a self-consistent determination of the pressure, as explained in Reference [24].The final expression for g 1 is where λ + i and λ − i indicate the right and left limits, respectively.

SAFT-VR-2D
The prediction for the a 1 and a 2 terms has been compared with Monte Carlo values derived according to the procedure described in Reference [29], applied to a 2D SW fluid.Very good agreement has been obtained for SW systems with different ranges.In Figure 2, we present the case of the fluctuation term a 2 when λ = 1.5.Notice that the local compressibility approximation gives a very good prediction of a 2 , compared with MC values.The agreement between theoretical and simulated values is higher in this case than for the corresponding prediction in 3D systems, suggesting that the LCA is a very good approximation in two-dimensional SW systems.As an overall test of the 2D perturbation theory, we present in Figure 3 the theoretical predictions for the internal energy as a function of the density, for four different values of the SW range, and compared with NVT-MC data.A similar agreement between theory and computer simulation values is observed for the pressure.
In a previous article [27], the adsorption theory presented in the previous section was tested against computer The second-order perturbation term a 2 as a function of the density ρ * for a two-dimensional fluid composed by hard disks of diameter σ with an attractive pair potential described by a square well with an energy depth and a range λ = 1.5.
The solid line corresponds to the theory, using the local compressibility approximation, and circles correspond to Monte Carlo (MC) simulation data.The reduced density is defined as ρ * = Nσ 2 /S , where N is the number of disks contained within an area S .The MC values were obtained using 108 particles.simulation using the Gibbs Ensemble Monte Carlo method for inhomogeneous fluids [30].Excellent agreement was obtained between theory and simulation data.It is important to stress here that the simulation results were obtained without using the quasi-two-dimensional hypothesis, so the agreement observed is a direct validation of the quality of Adsorption isotherms for a monomeric square-well (SW) fluid confined between planar walls, for a temperature T * = 1.5.The wall-particle interaction is a square-well potential, with a depth energy w and range λ w .Solid circles are 2D packing fractions obtained from the integration of 3D density profiles simulated with the Gibbs Ensemble Monte Carlo technique for inhomogenous fluids (GEMC-IF) [30].The particle-particle SW range is λ = 1.5, whereas the particle-wall SW range is λ w = 0.2453.Different isotherms are reported, according to the value of * = w / .Lines correspond to the theoretical predictions obtained with the SAFT-VR approach modelling of the bulk (3D) and adsorbed (2D) fluids, for non-associating (solid lines) and associating particles (dotted lines).The associating fluid was modelled using two SW sites with energy depth assoc / = 5 and a bonding volume K = 0.05.
the 2D approach used in the theory.The effect of the association in the adsorption isotherms is presented in Figure 4 for a system of monomeric SW particles with two association sites, and for three different values of the particlewall interaction ( * = w / = 2, 4, 8).The results are compared with theoretical and simulation predictions for the non-associating SW system, for the same range, λ = 1.5.The differences introduced by the association depend on the value of the particle-wall energy.For high values of * the association increases the amount of adsorbed fluid with respect to the non-associating case.However, by reducing the value of * , a competing effect arises between wall and association interactions.In the case of * = 4, adsorption is enhanced by the association at low and moderate values of the bulk packing fraction η, but this tendency changes at high values of η, where adsorption decreases with respect to the non-associating case.At lower values of the wallparticle energy, i.e. * = 2, the reduction of the adsorbed amount is even greater, with the appearance of negative adsorption (Γ < 0).

Adsorption Isotherms for Carbon Dioxide
In Reference [27], we presented the application of the 2D SAFT-VR theory described in the previous section to Adsorption isotherms for carbon dioxide adsorbed onto dry activated carbon at temperature T = 318.2K. Gibbs and absolute adsorptions are reported scaled with their respective values at pressure P = 0.40 MPa.Theoretical predictions are compared with experimental data.Dotted lines correspond to the predictions using the SAFT-VR approach, and solid lines correspond to the square-well quadrupolar fluid of Reference [35].Molecular parameters for the 3D SAFT-VR prediction are taken from Reference [32].Circles and triangles are the experimental values for Gibbs and absolute adsorptions, respectively, taken from Reference [31].Molecular parameters are reported in Table 1.
modelling adsorption isotherms of nitrogen and methane adsorbed onto dry activated carbon, and we found excellent agreement for pressures up to 13 MPa.In Figure 5, we present the predictions for carbon dioxide adsorbed onto dry activated carbon.Gibbs and absolute adsorption isotherms were obtained according to Equations 2 and 3, respectively, and compared with experimental data reported in Reference [31].Bulk and adsorbed fluids were modelled using the same parameters reported in previous studies for carbon dioxide, considered as a chain molecule with two segments, within the SAFT-VR approach [32].In Table 1, the molecular parameters used are reported.Carbon dioxide is represented by a non-spherical molecule (m = 2) in both the bulk and adsorbed phases.The diameter of the spherical segments, σ, is also the same for bulk and adsorbed phases, whereas the parameters of the SW attractive potentials are different in order to account for the influence of the substrate on the pair potential.Following Sinanoglu and Pitzer's work on the adsorption of a LJ fluid [23], the energy well depth of the particle-particle potential in the adsorbed phase was obtained from its bulk phase value using a reduction factor of 20%, i.e. ads = 0.8 .The range λ ads is determined from the ratio between the critical temperatures of the bulk (T bulk c ) and a monolayer adsorbed phase (T ads c ), R c = T ads c /T bulk c , as we explain now.This ratio is a function of the SW parameters of the bulk and adsorbed phases, R c = R c (m, σ, , λ, ads , λ ads ).Since all the parameters except λ ads are known, this last quantity can be tuned in order to reproduce the experimental value of R c .For noble gases and methane adsorbed onto graphite, it is known that 0.39 ≤ R c ≤ 0.43 [26].According to the experimental determination of the phase diagram of a carbon dioxide monolayer adsorbed onto graphite, T ads c = 127.5K [33].Since T bulk c = 304.136K [34], then R c = 0.42.In Table 1, we also report the parameter values for the particle-wall potential.We can interpret λ w σ as the mean vibration amplitude of the adsorbed particles [26], that can be modelled by a harmonic vibration.For the case of Kripton adsorbed onto graphite, this criteria gives λ w > 0.1305.On the other hand, from geometrical considerations of the adsorption of spherical particles, it is found that the upper limit that can be used to describe monolayer adsorption is λ w = 0.8165, that is the value that was considered in this work.Finally, the energy parameter w was fitted in order to reproduce the adsorption isotherm data.
We can see that the prediction obtained in both cases is very good for the absolute adsorption isotherm up to 13 MPa, whereas in the case of the Gibbs adsorption there is a strong deviation above 6MPa, very different to the behaviour of methane and nitrogen [27].These differences can be in part traced back to the strong quadrupolar moment of carbon dioxide.To confirm this hypothesis, we also present in Figure 5 the prediction obtained when we describe carbon dioxide as a monomeric quadrupolar SW fluid in the bulk, using the equation of state developed in [35].We can see that the inclusion of the quadrupolar moment at the level of the bulk fluid improves the theoretical prediction.A complete modelling of adsorption of quadrupolar fluids requires the development of a 2D quadrupolar equation of state.

Asphaltenes
Asphaltenes are a complex solubility class of compounds whose tendency to precipitate from petroleum causes serious operational problems in all facets of oil production.Asphaltenes are stabilised in crude oils by natural resins.The action of these surfactant-like agents on the asphaltene aggregation and precipitation processes is thought to be of paramount importance.Due to their molecular constitution, asphaltenes and resins have a mutual intrinsic effect on the stability of molecular self-assembling, either in the form of asphaltene-resin micellar association or asphalteneasphaltene association, that promotes precipitation.Due to the complexity of the molecular structure of asphaltenes and resins, to generate molecular-based equations of state for these substances is not an easy task.Usually, a whole series of aspects must be taken into account in order to satisfactorily obtain an equation of state for its use in industrial applications, and a compromise has to be made among all of them: the primitive model developed, the accuracy of the equation of state obtained, and, finally, the capability of this equation of state to predict the right phase diagrams.
A SAFT primitive model for asphaltene precipitation was developed in References [36,37], that describes the thermodynamic properties of an asphaltene-resin system considered as a solution, within the McMillan-Mayer approach, that is, introducing the Potential of the Mean Force (PMF) as the effective potential to take into account in the theoretical description.The aggregation behaviour of asphaltenes in asphaltene-resin mixtures within different host media has been studied on the basis of this PMF approach, using Molecular Mechanics and Molecular Dynamics [38], and according to the obtained computer simulation results, an important feature observed in the PMF of aggregated systems is the presence of repulsive barriers.
In the same spirit of the approach followed in References [36,37] to modelling asphaltenes, the SAFT-VR methodology has been used to describe the precipitation of asphaltenes in several Mexican crude oils [39].A very good prediction of the precipitation envelopes has been obtained, and we are currently working on the extension of this theory to predict the behaviour of asphaltene-resin mixtures in porous media.A first approximation followed in this study has been the modelling of adsorption of asphaltenes onto flat walls using the SAFT-VR-2D theory discussed in previous sections.Preliminary results using this approach are presented in Figures 6 and 7, where we are only modelling the effect on the asphaltene molecules and we are not considering the resins.Adsorption properties were obtained using the molecular parameters reported in Reference [39] for the bulk phase, whereas the parameters for the adsorbed phase were obtained following the same approach used for other substances [26,27] and in the previous section.Due to the lack of enough experimental data to model asphaltene adsorption, such as R c , several approximations were used, based on the results obtained for other systems and in order to have a qualitative prediction of the phenomena.Summarising the results that are known for the systems that have been studied under this theoretical approach (i.e., Ar, Kr, Xe, CH 4 , N 2 and CO 2 ) we have that 0.39 ≤ R c ≤ 0.43, 0.1305 ≤ λ w ≤ 0.8165 and 7.8 ≤ w / ≤ 10.We then selected for the asphaltene system R c = 0.40, λ w = 0.2453 and w / = 8.0, which are basically the same mean-field values used for adsorption of noble gases onto graphite [26] (see Tab. 2).
In Figure 6, the adsorption isotherm for a temperature T = 298 K is presented in terms of the 2D asphaltene Prediction for the adsorption isotherm of asphaltenes adsorbed onto a model planar wall at temperature T = 298 K. Molecular parameters were taken from the 3D SAFT-VR model for asphaltenes of Reference [39], see Table 2.The wall-particle energy was taken as * = 8.0, and γ and η correspond to the 2D and 3D packing fractions, respectively.Prediction for the fraction of non-associated molecules χ as a function of the bulk packing fraction η, corresponding to the adsorption isotherm of asphaltenes shown in Figure 6.Solid and dashed lines correspond to bulk and adsorbed fluid, respectively.packing fraction as a function of the corresponding 3D value.There are two features that are very interesting in this figure.First, the amount of adsorbed asphaltene is very high even at low bulk densities, suggesting that the associating behaviour enhances adsorption.However, the second point to stress from Figure 6 is that adsorption decreases as the bulk density increases until a relative minimum is reached before increasing again.Since we have already seen that the shape of adsorption isotherms is qualitatively modified as a consequence of a competing effect between the particle-wall and the site-site interactions (Fig. 4), then we can understand that the shape of the adsorption isotherm of Figure 6 can be modified as a consequence of the interplay between these interactions.Therefore, we can conclude that this very simple model predicts that under confinement we could expect that the asphaltene adsorption would tend to decrease at moderate densities.In Figure 7, we report the fraction of non-associated asphaltenes in the bulk and adsorbed onto the wall, χ.Whereas χ decreases as η increases for bulk asphaltenes, following the same pattern that would be expected for any associating fluid under the same conditions, in the case of the adsorbed phase a different behaviour arises, since χ has a relative maximum at low densities.This result suggests that association in the adsorbed phase occurs in four stages.In the first one, adsorption occurs without enhancing association, in such a way that dχ dη > 0. By increasing the bulk density, a second stage appears where χ gets its maximum value and after this, a bulk density region exists where dχ dη < 0. In the fourth and last stage, the adsorbed fluid tends to reach a saturated value where dχ dη = 0. Notice in Figure 7 that the whole process occurs for low values of χ, i.e., the adsorbed molecules are mainly associated, and then precipitation would be enhanced by the wall, with respect to the bulk.
Since the right development of a theoretical equation of state for very complex substances such as asphaltenes requires being able to test the assumptions made as well as the reliability of the primitive model, accomplishing this goal requires the use of computer simulation.We performed Monte Carlo simulations in order to test the approach followed in Reference [39].In Figures 8 and 9, we present the predictions for the internal energy, U/N, and heath capacity at constant volume, C v /Nk, respectively, obtained with the SAFT-VR approach description of asphaltene precipitation presented in Reference [39].The results are compared with Monte Carlo simulation data obtained for the same model.Since the 2D and 3D SAFT-VR EOS are based on secondorder perturbation theories, in both figures we include the effects that would arise if we consider higher-order perturbation terms (SAFT-VR + HOT), using the SW EOS of Reference [40].As can be observed in both figures, SAFT-VR + HOT gives a better representation of the properties, although in the case of the heath capacity the deviation is higher.
As pointed out in Reference [41], when the SAFT-VR molecular parameters for an associating fluid, obtained by fitting to experimental data, are used to study the same molecular model by computer simulation, then multiple bonding effects are observed in the simulated configurations The internal energy U * as a function of the density ρ * for an associating square-well fluid according to the primitive model of asphaltenes presented in Reference [39].The particles have two associating sites.The reduced internal energy is defined as U * = U/ , where is the energy depth of the SW fluid, and the reduced density is given by ρ * = Nσ 3 /V, where N is the number of particles contained in a volume V and σ is the diameter of the particles.The solid line corresponds to the prediction given by the SAFT-VR approach, according to the molecular model of Reference [39], and the dotted line is the SAFT-VR approach with higher-order perturbation terms (SAFT-VR + HOT), according to the expressions given in Reference [40].
Circles correspond to Monte Carlo simulation data using 108 particles in the NVT ensemble.The heath capacity at constant volume C v * as a function of the density ρ * for an associating square-well fluid according to the primitive model of asphaltenes presented in Reference [39].The particles have two associating sites.The reduced heath capacity is defined as C v * = C v /Nk, and the reduced density is given by ρ * = Nσ 3 /V, where N is the number of particles contained in a volume V and σ is the diameter of the particles.The solid line corresponds to the prediction given by the SAFT-VR approach, according to the molecular model of Reference [39], and the dotted line is the SAFT-VR approach with higher-order perturbation terms (SAFT-VR + HOT), according to the expressions given in Reference [40].Circles correspond to Monte Carlo simulation data using 108 particles in the NVT ensemble.that are not considered in the theory.We have observed the same effect in the simulation of the primitive model of asphaltenes that we are considering.The computer simulation data used in Figures 6 to 9 were obtained using the same molecular parameters reported in Reference [39].Multiple bonding arises in the simulated system as a consequence of the size of the sites determined when the theoretical parameters are fitted to experimental data, as shown in Figure 10.The same problem arises when we study confined asphaltenes and we have to bear in mind this effect when theoretical values are compared against simulation data, as in Figures 8 and 9.
Standard techniques to rotate molecules used for simple non-associating molecular fluids must be used with caution when applied to systems with associating sites, since associated clusters could significantly reduce the efficiency of Monte Carlo movements.In Appendix II, we report a rotation procedure that has been useful in our studies of bulk associating systems as well as under confinement, based on a two-vector representation of the rotation matrix introduced in Reference [42].We have adapted this two-vector scheme to be used in molecular simulations of associating systems, and we have found that this algorithm enables us to control the rotation parameters in a more robust way in comparison with other algorithms (like quaternions), particularly for non-linear molecules and without reducing the speed of the simulation.

CONCLUSIONS
In this study, we have presented an extension of the SAFT-VR and DPT approaches to model thermodynamic properties of confined fluids.Using a quasi-2D approach, we have shown that the method can describe accurately computer simulation results as well as properties of real substances, such as adsorption isotherms of carbon dioxide.We have also presented the extension of the SAFT-VR modelling of asphaltenes when adsorption phenomena occurs.The theory relies on the determination of eight parameters used to describe the shape of the particles and the three different types of involved potentials: the molecular shape parameter, i.e., m, five parameters related to particle-particle interactions in the bulk and adsorbed phases, i.e., σ, , λ, ads and λ ads , and another two parameters that describe the SW particle-wall interactions, i.e., λ w and w .Assuming that m and σ are the same for the bulk and adsorbed particles, the first four parameters can be determined from bulk phase properties, whereas ads and λ w can be predicted based on theoretical descriptions of the effects of the substrate on the particle-particle interaction in the adsorbed phase, and the size of an adsorbed monolayer, respectively.With all this information, the particle-wall range potential, λ w , is determined from the critical temperatures of the adsorbed and bulk phases, R c .The last parameter to be determined, w , is obtained by fitting directly to an adsorption isotherm of the system under study.Although the theory is limited by neglecting the adsorbed-bulk interface, our results indicate that the SAFT-VR approach can give useful insights into the behaviour of associating systems under confinement.
In the last equation, g 2D is the contact value of the 2D monomer-monomer radial distribution function, K AB is the bonding volume determined by sites A and B, and f AB is the Mayer function of the A-B site-site bonding interaction, that is given by a SW potential with range λ s and energy depth The corresponding expressions for the 2D version of the SAFT-VR approach follow the same structure as the 3D case.
The HD Helmholtz free energy is obtained from the Henderson equation [43], and the corresponding contact value of the radial distribution function, g HD , is given by where γ is the two-dimensional packing fraction.The firstorder perturbation term of the Helmholtz free energy is given by a where, as in the 3D case, the radial distribution function at contact, g HD , is evaluated with an effective packing fraction, γ e f f , given by γ where and The second-order term is given by the local compressibility approximation (LCA) result [28] Finally, g 2D 1 (σ) is given by the self-consistent expression

APPENDIX II
In this section, we report an efficient rotation algorithm to rotate systems with sites, based on the procedure described in Reference [42].If we denote by e 0 the unit vector along a molecular axis of a particle, representing a linear or nonlinear molecule, then a new orientation e can be generated through a vector v obtained by a random orientation within the surface of a unit sphere [44] e = e 0 + λv e (39) where λ is a parameter that controls the size of the rotation and e is the norm of vector e.If we denote by R the rotation matrix that transforms e 0 into e according to the relation e = Re 0 (40) then a basic question is to determine R, if we know the initial and final vectors, e 0 and e.This problem has been solved previously [42], where a new parametrisation of matrix R was given in terms of these two vectors and an angle ζ.
Basically, matrix R is written as the product of another two matrices R(e 0 , e, ζ) = S(e, ζ)T(e 0 , e) where T gives a rotation about the axis orthogonal to both vectors, and S gives a rotation around e through an angle ζ.We now describe these rotations.Given e 0 and e, we generate their dot and cross-products: N = e 0 × e = n sin φ (43) where φ is the angle that rotates e 0 into e about the axis n, which is a unit vector orthogonal to e 0 and e.Through the use of quaternion algebra, we can simplify the description of T in terms of n and φ.The quaternion associated with T is q = q 0 + q (44) where Using this quaternion, Equations 42 and 43, and standard trigonometric identities, then the diagonal matrix elements of T are given by: T 33 = q 2 0 − q 2 1 − q 2 2 + q 3 For a linear molecule, it is not necessary to apply matrix S.However, for a non-linear molecule, once we have obtained e, it is necessary to apply T and S to have the general rotation matrix that corresponds to the rotation of the molecule once we have followed the random process selection of the new orientation e.To apply S numerically, we have to restrict the value of angle ζ to the range 0 ≤ ζ ≤ π.In this way, the values of the trigonometric functions are Using the procedure outlined in this appendix, we have performed computer simulation of different models of linear and non-linear molecules.Specifically for the case with primitive models that includes site-site interactions, we have concluded that this rotation algorithm allows one to control the efficiency of the MC moves in a better way than by applying standard procedures, such as quaternions.

Figure 3
Figure 3 The internal energy as a function of the packing fraction for a two-dimensional fluid composed by hard disks of diameter σ with an attractive pair potential described by a square well (SW) with an energy depth and a range λ = 1.5.The reduced internal energy is defined by U * = U/ .Results are presented for several values of the range, λ = 1.3 (open circles), 1.4 (solid line), 1.5 (dots) and 1.6 (triangles).Solid circles represent Monte Carlo simulation data for 108 particles in the NVT ensemble.

Figure 10 Snapshot
Figure 10Snapshot of a Monte Carlo configuration for a simulated state for the primitive model of asphaltene presented in Reference[39], for a reduced density ρ * = 0.1 and a reduced temperature T * = 1.5.Notice the presence of multiple-bonding configurations.

TABLE 1
Values of the molecular parameters used to describe carbon dioxide adsorption onto dry activated carbon