Predicting the adsorption of n -perfluorohexane in BAM-P109 standard activated carbon by molecular simulation using SAFT- γ Mie coarse-grained force fields

Abstract

This work is framed within the Eighth Industrial Fluid Properties Simulation Challenge, with the aim of assessing the capability of molecular simulation methods and force fields to accurately predict adsorption in porous media for systems of relevant practical interest. The current challenge focuses on predicting adsorption isotherms of n-perfluorohexane in the certified reference material BAM-P109 standard activated carbon. A temperature of $T=273$  K and pressures of $p/p_0=0.1$ , 0.3, and 0.6 relative to the bulk saturation pressure p₀ (as predicted by the model) are the conditions selected in this challenge. In our methodology we use coarse-grained intermolecular models and a top-down technique where an accurate equation of state is used to link the experimental macroscopic properties of a fluid to the force-field parameters. The state-of-the-art version of the statistical associating fluid theory (SAFT) for potentials of variable range as reformulated in the Mie group contribution incarnation (SAFT-γ Mie) is employed here. The parameters of the SAFT-γ Mie force field are estimated directly from the vapour pressure and saturated liquid density data of the pure fluids using the equation of state, and further validated by molecular dynamic simulations. The coarse-grained intermolecular potential models are then used to obtain the adsorption isotherm kernels for argon, carbon dioxide, and n-perfluorohexane in graphite slit pores of various widths using Grand Canonical Monte Carlo simulations. A unique and fluid-independent pore size distribution curve with total micropore volume of 0.5802 cm³/g is proposed for the BAM-P109. The pore size distribution is obtained by applying a non-linear regression procedure over the adsorption integral equation to minimise the quadratic error between the available experimental adsorption isotherms for argon and carbon dioxide and purpose-built Grand Canonical Monte Carlo kernels. The predicted adsorption levels of n-perfluorohexane at 273 K in BAM-P109 are 72.75 ± 0.01, 73.82 ± 0.01, and 75.44 ± 0.05 cm³/g at Standard Temperature and Pressure (STP) conditions for $p/p_0=0.1$ , 0.3, and 0.6, respectively.

Introduction

The aim of the Eighth Industrial Fluid Properties Simulation Challenge (IFPSC) is to assess the capability of force fields and molecular simulation methods to accurately predict adsorption in porous media for systems of relevant practical interest. Specifically, the adsorption of n-perfluorohexane in the certified reference material BAM-P109 standard activated carbon is the purpose of the current IFPSC.

Perfluorinated molecules are known to be challenging to describe. Their structural and thermophysical properties are different to those of their analogous hydrocarbon compounds, in part due to the size, polarisation, and highly electronegative nature of the fluorinated group, which gives rise to the particularly steep nature of the repulsive interactions between perfluorinated groups (Maitland et al., 1981; Yee et al., 1992). Force fields based on a generalised Lennard-Jones (Jones, 1924a, 1924b; Lennard-Jones, 1931) or Mie potential (Grüneisen, 1912; Mie, 1903), where the repulsive and attractive exponent can be varied, allow a fine tuning of the steepness of the potential well and the range of the interaction, which have proven to be especially useful to describe the fluid phase behaviour (Potoff and Bernard-Brunel, 2009; Rubio et al., 1985) and second derivative thermodynamic properties of perfluoroalkanes (Lafitte et al., 2013).

As a further complication, the use of force fields developed for bulk fluids to describe fluid–solid interactions is not a straightforward task. In most cases empirical combining rules (e.g. the traditional Lorentz–Berthelot) are used to obtain the unlike-interaction parameters. The validity of such combining rules is however restricted to systems with similar type of molecular interactions (Desgranges and Delhommelle, 2014; Haslam et al., 2008). Although not the focus of this work, it is imperative to find robust and reliable approaches to obtain the fluid–fluid and fluid–solid interactions particularly in the case of systems with limited experimental data. It is also important to recognise that coarse-grained (CG) solid–fluid interactions will depend in a complex way on the morphology of the surface (Forte et al., 2014).

In our current work, we use CG intermolecular potential models to represent the fluid molecules and a collective of atomistically detailed graphite slit pores (kernel) to represent the microporosity of the reference solid. The CG force fields employed in our molecular simulations are developed using a top-down approach (Müller and Jackson, 2014), as this provides robustness in the predictions of the key properties for adsorption such as condensed liquid densities and vapour pressures. An essential element for top-down development of force fields is an accurate equation of state (EoS) characterised by a well-defined Hamiltonian enabling one to establish a formal link between target macroscopic thermodynamic properties and the CG force field. The latest version of the statistical associating fluid theory for potentials of variable range based on the Mie interaction, SAFT-VR Mie (Lafitte et al., 2013), and its group-contribution variant SAFT-γ Mie (Papaioannou et al., 2014) has been shown to be a reliable and robust method to develop CG force fields for direct use in molecular simulation (Avendaño et al., 2011, 2013; Herdes et al., 2015; Lafitte et al., 2012; Lobanova et al., 2015; Müller and Jackson, 2014; Theodorakis et al., 2015). The same approach is also used here to develop SAFT-γ Mie CG force-field parameters that are then assessed by molecular dynamic (MD) simulations and used within Grand Canonical Monte Carlo (GCMC) simulations to predict the adsorption of the system of interest.

Modelling approaches and force fields

SAFT-VR Mie EoS

Here we use an extension of the SAFT EoS (Chapman et al., 1989, 1990) referred to as SAFT-VR Mie (Lafitte et al., 2013). The generic SAFT approach stems from the thermodynamic perturbation theory of Wertheim (1984a, 1984b, 1986a, 1986b, 1986c, 1987) to explicitly take into account the non-sphericity and directional interactions of molecules. The SAFT-VR Mie (Lafitte et al., 2013) methodology is a reformulation of the SAFT-VR EoS (Galindo et al., 1998; Gil-Villegas et al., 1997) for Mie potentials with enhanced accuracy in the description of the fluid phase behaviour (particularly in the near-critical region) and second derivative properties of fluids. Full details of the development of the theory are given in Lafitte et al. (2013).

Within the SAFT-VR Mie formalism molecules are described as chains of m tangent segments that interact through a Mie potential (Grüneisen, 1912; Mie, 1903)

$u^{\text{Mie}}(r)=C{\varepsilon }_{\text{ff}}[(\frac{{\sigma }_{\text{ff}}}{r}{)}^{{\lambda }_{\text{r},\text{ff}}}-(\frac{{\sigma }_{\text{ff}}}{r}{)}^{{\lambda }_{\text{a},\text{ff}}}]$

(1)

where r is the distance between the spherical segments, ${\sigma }_{\text{ff}}$ is the segment diameter (at which the potential is zero), ${\varepsilon }_{\text{ff}}$ is the potential depth, ${\lambda }_{\text{r},\text{ff}}$ and ${\lambda }_{\text{a},\text{ff}}$ are the repulsive and attractive exponent, respectively, and the constant

$C=\frac{{\lambda }_{\text{r},\text{ff}}}{{\lambda }_{\text{r},\text{ff}}-{\lambda }_{\text{a},\text{ff}}}(\frac{{\lambda }_{\text{r},\text{ff}}}{{\lambda }_{\text{a},\text{ff}}}{)}^{\frac{{\lambda }_{\text{a},\text{ff}}}{({\lambda }_{\text{r},\text{ff}}-{\lambda }_{\text{a},\text{ff}})}}$

(2)

ensures that the potential minimum is at $-{\varepsilon }_{\text{ff}}$ . The subindex f (or s in the case of solid particles) refers to the individual nature of the species that interact exclusively via pair-wise forces.

The equation of state for a non-associating Mie fluid is written in a closed form in terms of the Helmholtz free energy A as a sum of the following contributions

$A=A^{\text{ideal}}+A^{\text{mono}}+A^{\text{chain}}$

(3)

where $A^{\text{ideal}}$ corresponds to the free energy of an ideal gas, $A^{\text{mono}}$ is the change in free energy due to the segment–segment repulsion and dispersion interactions, and $A^{\text{chain}}$ is the contribution due to the possibility of bonding between the segments to form a chain of m segments (Lafitte et al., 2013). As associating molecules are not considered in this work, the relevant term is omitted here. For a pure fluid, the parameters of the SAFT-VR Mie EoS correspond to those that define the potential, i.e. m, ${\sigma }_{\text{ff}},{\varepsilon }_{\text{ff}},{\lambda }_{\text{r},\text{ff}}$ , and ${\lambda }_{\text{a},\text{ff}}$ . Hence, by estimating the EoS parameters from macroscopic experimental data, one is simultaneously obtaining the corresponding parameters of the underlying intermolecular potential. This procedure guarantees the accuracy of the potential in representing the macroscopic properties with a reliable effective (average) molecular model of the fluid and represents a unique bridge between the microscopic and macroscopic scales.

Molecular models

The SAFT-VR Mie (Lafitte et al., 2013) (and SAFT-γ Mie (Papaioannou et al., 2014)) EoS is based on the assumption that the basic building blocks of a molecule are spherical segments. This is particularly useful to describe molecules in a CG manner: we dispense with the atomistic detail and represent molecules in terms of a tangent chain of beads. In our current study, the carbon dioxide molecule is modelled as a tangent dimer, n-perfluorohexane as a rigid trimer, and argon as a sphere. Note that potentials developed with this methodology are referred to as SAFT-γ Mie force fields in recognition of the more generic group-contribution nature of the approach (Müller and Jackson, 2014).

There is a corresponding-states relationship between the exponents in a Mie potential in such a way that an infinite number of pairs ( ${\lambda }_{\text{r}}$ and ${\lambda }_{\text{a}}$ ) provide essentially the same description of the fluid phase properties (Ramrattan et al., 2015). For simplicity we fix the value of the attractive exponent to the London value of ${\lambda }_{\text{a},\text{ff}}=6$ throughout this work.

The estimation of the three remaining parameters ( ${\sigma }_{\text{ff}},{\varepsilon }_{\text{ff}}$ , and ${\lambda }_{\text{r},\text{ff}}$ ) of equation (1) for the fluids of interest is undertaken through a comparison with experimental data. In particular, based on the experimental information available for fluids adsorbed in BAM-P109, we focused on argon, carbon dioxide, and perfluorohexane. SAFT-VR Mie parameters for argon are already available (Dufal, 2013; Dufal et al., 2014). Argon is represented as a single Mie sphere model (i.e. m = 1) and the corresponding parameters are obtained using experimental data for the vapour pressures and saturated liquid densities in a range of temperatures up to 95% the critical point. Previously published models for carbon dioxide (Dufal, 2013; Dufal et al., 2014; Lafitte et al., 2013) and n-perfluorohexane (Lafitte et al., 2013) are not used in our current work because non-integer SAFT chain models cause some complications (e.g. the specification of the bond length) when used directly in molecular simulation. Instead we develop a two-segment m = 2 model for carbon dioxide and three-segment m = 3 model for n-perfluorohexane using vapour pressure and saturated liquid density data in a range of temperatures of $\pm 20$  K in the vicinity of the temperature of current interest ( $T=273$  K). The resulting SAFT-γ Mie intermolecular model parameters obtained are collected in Table 1 together with the deviations of the SAFT-VR Mie description from the available experimental data (Dunlap et al., 1958; Khairulin et al., 2007; Lemmon et al., 2014). The parameters in Table 1 are obtained using the EoS with a standard parameter estimation procedure (see Lafitte et al. (2013) for details). The comparison of the description of the vapour–liquid coexistence envelop and vapour pressure curve with the corresponding experimental data can be seen in Figures 1 to 3.

Figure 1.

Coexistence of vapour–liquid densities, vapour pressures, and surface tensions for argon. The dashed red curve corresponds to the SAFT-VR Mie EoS calculations, the continuous blue curve corresponds to correlated experimental data (Lemmon et al., 2014), and the orange circles correspond to molecular dynamics simulations, using the same parameters as the theory (color online).

Figure 2.

Coexistence of vapour–liquid densities, vapour pressures, and surface tensions for carbon dioxide. Symbols as in Figure 1.

Figure 3.

Coexistence of vapour–liquid densities, vapour pressures, and surface tensions for n-perfluorohexane. Symbols as in Figure 1.

Table 1.

The SAFT-γ Mie CG force fields for the fluids of interest (the attractive exponent is fixed at ${\lambda }_{\text{a},\text{ff}}=6$ in all cases) and percentage absolute average deviations (% AAD) of the theoretical description from the experimental vapour pressures and saturated liquid densities in the range of temperatures at which the parameter estimations were done.

	Model parameters						% AAD
	Ref.	m	${\sigma }_{\text{ff}}$ /Å	( ${\varepsilon }_{\text{ff}}/k_{\text{B}}$ )/K	${\lambda }_{\text{r},\text{ff}}$	${\lambda }_{\text{a},\text{ff}}$	p	${\rho }_{\text{L}}$	Exp. Ref.
Ar	Dufal (2013), Dufal et al. (2014)	1.0	3.4038	117.84	12.085	6.0	0.2	0.7	Lemmon et al. (2014)
CO2	This work	2.0	2.8514	190.14	13.77	6.0	0.07	0.4	Lemmon et al. (2014)
C6F14	This work	3.0	4.4337	298.30	20.106	6.0	0.06	0.3	Dunlap et al. (1958), Khairulin et al. (2007), Lemmon et al. (2014)

$k_{\text{B}}$ refers to Boltzmann’s constant.

Molecular simulation of fluid phase equilibria

Direct MD simulations of the vapour–liquid equilibria (VLE) of carbon dioxide, argon, and n-perfluorohexane are carried out using the GROMACS open source suite (Van Der Spoel et al., 2005). In the direct technique, a film of liquid surrounded by vapour (with two stable vapour–liquid interfaces) is simulated to obtain the VLE and interfacial properties of the system.

The simulations are performed in the canonical (NVT) ensemble, where the total density and temperature are kept constant. The parallelepiped cells have aspect ratios of $\mathrm{Lz}/\mathrm{Lx}=6$ and $\mathrm{Lx}=\mathrm{Ly}=4$ , 5, and 6 nm to accommodate N = 1740 carbon dioxide, 4394 argon, and 1152 n-perfluorohexane molecules, respectively.

The MD simulations are thermostated every 2 ps using a Nose–Hoover algorithm (Frenkel and Smit, 2002; Hoover, 1985; Nosé, 1984), and all non-bonded interactions are truncated at $r_{\text{c}}=6{\sigma }_{\text{ff}}$ . The systems are simulated with a time step of 0.01 ps for at least 20 ns. All reported properties are obtained as appropriate averages, taken over the final half of the configurations explored. The simulated vapour–liquid coexistence densities, vapour pressures, and surface tensions are also included in Figure 1 to 3 along with the experimental and SAFT-VR Mie description for completeness. The agreement between the theory, simulation, and experiments is found to be excellent. The vapour pressure of n-perfluorohexane at 273 K predicted by the SAFT-γ Mie CG force field is $p^0=0.008664$  MPa, which is in excellent agreement with the experimental value of 0.008548 MPa (Dunlap et al., 1958; Khairulin et al., 2007; Lemmon et al., 2014).

Molecular simulation of the adsorbed phase

The adsorption isotherms are determined using GCMC simulations in the $\mu \mathrm{VT}$ ensemble. A slit-like pore geometry is used to model a solid with the structure of graphite. Two solid blocks (each comprising six carbon layers) are considered. Boundary conditions are imposed in the Cartesian directions parallel to the surfaces of the solid. The pair potentials are truncated at a distance proportional to the fluid segment diameter such that $r_{\text{c}}=6{\sigma }_{\text{ff}}$ , at which the value of the potential is assumed to be negligible. No long-range corrections are employed. The solid is described based on the hexagonal honeycomb lattice typical of graphite where the space between adjacent carbons is 1.42 Å, with the parameters commonly employed for atomistically detailed carbon surfaces (Steele, 1973): ${\sigma }_{\text{ss}}=3.40$  Å, ${\varepsilon }_{\text{ss}}/k_{\text{B}}=28.0\text{K},{\lambda }_{\text{r},\text{ss}}=12.0,{\lambda }_{\text{a},\text{ss}}=6.0$ .

The unlike-interaction parameters between segments of different type (in this case fluid and solid segments) are obtained using the Lorentz–Berthelot combining rules, i.e.

${\sigma }_{\text{fs}}=\frac{{\sigma }_{\text{ff}}+{\sigma }_{\text{ss}}}{2}$

(4)

and

${\varepsilon }_{\text{fs}}=\sqrt{{\varepsilon }_{\text{ff}}{\varepsilon }_{\text{ss}}}$

(5)

while the following relationship is used for the unlike attractive and repulsive exponents to ensure that the geometric mean is retained for the unlike van der Waals attractive interaction (Lafitte et al., 2013)

${\lambda }_{\text{k},\text{fs}}-3=\sqrt{({\lambda }_{\text{k},\text{ff}}-3)({\lambda }_{\text{k},\text{ss}}-3)},\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}}\text{k}=\text{a},\text{r}$

(6)

and no attempt was made to improve the description at low pressures (Henry’s law region).

The GCMC simulations are carried out using a standard procedure starting with an empty pore which is filled until equilibrium is attained (Allen and Tildesley, 1989; Frenkel and Smit, 2002; Nicholson and Parsonage, 1982). For each Monte Carlo cycle, both the displacement/reorientation of a randomly chosen fluid molecule and a random creation/destruction of a fluid molecule are attempted. The systems are equilibrated at each pressure level and then averages are determined from a further two million cycles. Block averages are also obtained every 10,000 cycles and are used to compute standard deviations to calculate uncertainties.

In the GCMC simulations the temperature and the activity (Müller, 2010) (which is directly related to the chemical potential) are specified. Simulations are also performed for bulk fluids that would hypothetically be in equilibrium with the adsorbed fluid (i.e. these are performed at the same temperature and activity) so that the density of the bulk fluid can be calculated as a simple average. This density is then used to determine the pressure at a given temperature using the SAFT-VR Mie EoS. Since the EoS reproduces faithfully the volumetric properties of the fluids (see Figures 1 to 3), this procedure is in essence equivalent to generating these pressures from pure component simulations, albeit more efficiently.

Development of a kernel of adsorption isotherms

In order to characterise the experimental data provided in the challenge documentation, we built a GCMC kernel of isotherms which are used to represent the known experimental adsorption isotherms of argon and carbon dioxide in BAM-P109. This procedure guarantees that we obtain a pore size distribution (PSD) compatible with the molecular modelling strategy avoiding any uncertainties in the underlying theoretical treatment as would be the case if the PSD provided in the challenge were used directly.

The pore widths that constitute our kernel are chosen to be consistent (or commensurate) with the expected PSD of BAM-P109 as suggested in the challenge documentation. We employ 11 values: H = 0.320, 0.488, 0.656, 0.824, 0.992, 1.160, 1.328, 1.496, 1.664, 1.832, and 2.0 nm. The distances between the solid surfaces, measured from the centre of the carbon atoms, correspond to $H+{\sigma }_{\text{ss}}$ . For each one of the pore widths and fluids, six simulations are performed to span the range of pressures of interest for the adsorption of argon and carbon dioxide.

Unique PSD curve and estimation of the pore volume for BAM-P109

An important step towards the prediction of the adsorption isotherm is the determination of a reliable PSD of the BAM-P109 activated carbon. The porosity analysis, provided as part of the challenge documentation, comprises three differential pore volume distributions, two from argon (using quenched solid density functional theory, QSDFT, and non-local density functional theory, NLDFT, at 87 K) and another from carbon dioxide (using NLDFT at 273 K), which we reproduce in Figure 4. Several problems are evident from an inspection of Figure 4. One of the most obvious is the non-uniqueness of the PSD, as the analysis for each fluid (argon and carbon dioxide) and each temperature differs considerably. Another evident problem is the suggestion of unrealistically small pore sizes by means of the more refined theory (QSDFT).

Figure 4.

Differential pore volume distribution of BAM-P109. The blue curve (triangles) corresponds to the quenched solid density functional theory analysis performed with argon at T = 87 K, the green curve (squares) to the non-local density functional theory (NLDFT) analysis performed with argon at T = 87 K, the red curve (diamonds) to the NLDFT analysis performed with carbon dioxide at T = 273 K. Information taken from the benchmark data of the Challenge. The black curve (filled circles) is the PSD determined for use in our current work (color online).

Given an experimental adsorption isotherm $\Gamma (p)$ , related to the PSD (denoted by f(H) as a probability density function (Meyer and Klobes, 1999)) and the fluid adsorption isotherm kernels $\rho (H,p)$ by the adsorption integral equation

$\Gamma (p)={\int }_{H_{\text{min}}}^{H_{\text{max}}}f(H)\rho (H,p)\text{dH}$

(7)

a consistent procedure to obtain the PSD is to invert equation (7) for each fluid, which corresponds to a Fredholm integral equation of the first kind. This can be achieved by using an appropriate regularisation procedure.

We aim to determine an averaged unique single-fluid independent PSD f(H) by minimising the quadratic error between the experimental adsorption isotherms $\Gamma (p)$ for argon and carbon dioxide with their respective GCMC adsorption isotherms kernels $\rho (H,p)$ .

Results and discussion

In Figures 5 to 7 we show the results obtained by simultaneously fitting the experimental data to the kernels generated for each fluid: the averaged single fluid independent PSD f(H) (with the ancillary pore volume variable, which we estimated to be a value of 0.5802 cm³/g at $p/p_0=$ 0.6), the adsorption isotherms of argon at 87 K and carbon dioxide at 273 K. The resulting unique PSD is characterised by a smoother profile and a mean-pore size which appears to be shifted towards slightly larger pores, when compared with those provided in the challenge documentation (see Figure 4). This is consistent, however, with the trends observed by other researchers using similar GCMC kernels (Herdes et al., 2005; Miyahara et al., 2014) and provides confidence in our findings. With the average PSD and the construction of an appropriate GCMC kernel $\rho (H,p)$ for n-perfluorohexane, the evaluation of the total adsorption $\Gamma (p)$ from equation (7) is an straightforward task. Our predictions of this key target of the challenge are presented in Table 2.

Figure 5.

Calculated probability density function for the pore size distribution of BAM-P109 (color online).

Figure 6.

Result of the fitting of the experimental data (continuous curve) and GCMC kernel (circles) for argon at T = 87 K.

Figure 7.

Result of the fitting of the experimental data (continuous curve) and GCMC kernel (circles) for carbon dioxide at T = 273 K.

Table 2.

Predicted adsorption of n-perfluorohexane Γ(p) in BAM-P109 at 273 K and standard deviation, SD, for various pressures $p/p_0$ relative to the bulk saturated pressure of the fluid model (p₀ = 8.664 kPa for C₇F₁₆).

$p/p_0$	Γ(p)/(STP cm3g−1)	±SD/(cm3g−1)
0.1	72.75	0.01
0.3	73.82	0.01
0.6	75.44	0.05

Conclusions

We have demonstrated a synergetic use of molecular simulation methodologies, accurate force fields, and experimental data for the prediction of adsorption isotherms of challenging fluid–solid systems. It is expected that the determination of a unique PSD based on simultaneous consideration of the isotherms of two distinct fluids will provide a sound qualitative prediction of the unknown adsorption of a fluid, as the main goal specified in the challenge. Had there been data for the adsorption isotherm of n-perfluorohexane, these could have been considered within our methodology to provide further confidence in the calculated values. The interaction between carbon and perfluorinated moieties is non-trivial and this is most likely to be the largest source of uncertainty in the predicted results.

After the disclosure of the Benchmark results (Ross et al., 2016), we were pleased to discover that our predictions had the correct behaviour. Our prediction and the experimental data (made available at a special session sponsored by CoMSEF at the AIChE annual meeting) are compared in Figure 8. As expected, we over-predicted the adsorption (by approximately 20%) as our model surface was a perfect graphite slit pore. In hindsight, the consideration of a structure with a tuneable degree of disorder such as random carbon platelets could (Kumar et al., 2011, 2012) have improved the accuracy of the predictions as was demonstrated by the entry of Sarkisov (2016) to this challenge.

Figure 8.

Comparison of the total adsorption predicted (diamonds) for n-perfluorohexane at T = 273 K and the experimental data (circles) released after the disclosure of the benchmark results. The curves are just guides to the eye.