How to Optimize the Electrostatic Interaction between a Solid

Article pubs.acs.org/JPCC

How to Optimize the Electrostatic Interaction between a Solid Adsorbent and CO2 Edder J. García,† Javier Pérez-Pellitero,† Christian Jallut,‡ and Gerhard D. Pirngruber*,† †

IFP Energies Nouvelles, Rond Point de l’échangeur de Solaize, 69360 Solaize, France Université de Lyon, Université Lyon 1, Laboratoire d’Automatique et de Génie des Procédés, UMR 5007, CNRSESCPE, 43, Bd du 11 Novembre 1918, 69622 Villeurbanne Cedex, France

‡

S Supporting Information *

ABSTRACT: Pressure swing adsorption (PSA) is a promising method for an economically viable separation of CO2 from gas mixtures. One important step in developing efficient PSA separations is to choose an adsorbent with a good compromise between selectivity and working capacity. The problem can be seen as an optimization process where one searches the properties of the adsorbent that give the best trade-off between selectivity and working capacity. In a recent work we have developed a method to optimize the properties of nonpolar adsorbents for PSA separation of CO2. In this work we have extended the analysis to polar adsorbents. A simplified model of the electric field generated by the solid is used to gain a better understanding of the effect of the geometry and the magnitude of the field on the adsorption behavior of the system. General rules about optimal adsorbents for CO2 adsorption in porous solid can be drawn. These rules allow defining a strategy for the intelligent design of porous materials for CO2 separation applications.

1. INTRODUCTION Because of the rising concern about global warming, which is mainly caused by excessive anthropogenic CO2 emissions, an intensive research activity is going on to develop new, improved adsorbent materials for CO2 capture. The enormous progress that has been made in synthesis and design of porous solids has led to an explosion of the number of materials that might potentially be interesting for CO2 separations. Yet, the development of new materials and their testing for CO2 adsorption remains a highly empirical procedure. Although efforts toward an intelligent adsorbent design are ongoing,1−3 we are still lacking clear guidelines of what the optimal adsorbent for a given CO2 capture application should look like, for example, in terms of pore geometry and chemical composition. In recent work we have made a first step in this direction: we proposed a method to reduce complex pore systems to simplified pore geometries (spheres or cylinders) with an equivalent pore diameter that exerts the same confinement effect as the real system.4 Then analytical equations were developed that allowed calculating the adsorption isotherms of small (shapeless) molecules as a function of a few, easily available descriptors of the solid adsorbent and the adsorbate molecule (the fluid):5 the equivalent pore radius, the surface density of the solid (atoms per surface area), the volume density of the solid (atoms per unit cell volume), the depth of the solid−fluid potential well, the depth of the fluid’s intermolecular Lennard-Jones potential well, and the fluid’s excluded volume. From the adsorption isotherms we can then deduce the performance of the material in a pressure swing adsorption (PSA) process for CO2 by calculating sorbent © 2014 American Chemical Society

selection criteria, i.e., the working capacity and selectivity. The attractiveness of this approach lies in the fact that it can be easily inverted. For given operating conditions of the PSA process, we can determine the values of the solid’s descriptors that will lead to the optimal values of the sorbent selection criteria. Thus, we can draw a sketch of the ideal adsorbent material, in terms of porosity and composition. However, the method that we have recently presented only works for nonpolar or slightly polar adsorbents: the analytical model only takes into account Lennard-Jones (i.e., dispersion− repulsion) interactions between adsorbent and adsorbate, but not electrostatic forces. Hence, it is only applicable for materials where dispersion−repulsion forces are governing the adsorption process. This is a serious shortcoming since slightly polar adsorbents are not very selective for the adsorption of CO2.6−8 High CO2 selectivities are achieved by an interaction of the electrostatic field of the adsorbent, created by charged atoms of the solid, with the quadrupole moment of CO2. Therefore, the objective of the present contribution is to develop, in line with our previous work, a simplified method for calculating the electrostatic interaction between the adsorbent and CO2. The idea is to represent the charge distribution of the real solid by a simplified system that can be characterized by a few descriptors, which can then be optimized with respect to the abovementioned sorbent selection criteria. Two problems need to be solved. First, we have to find a simplified representation of the charge distribution of the solid that satisfactorily reproduces the Received: January 8, 2014 Revised: April 9, 2014 Published: April 10, 2014 9458

dx.doi.org/10.1021/jp500209v | J. Phys. Chem. C 2014, 118, 9458−9467

The Journal of Physical Chemistry C

Article

electrostatic field of the real system. Second, we have to establish a link between the electrostatic field and the adsorption isotherm, to be able to calculate sorbent selection criteria. This paper is organized as follows: We first describe the simplified representation of the real solid and compare the electrostatic fields of the real and of the simplified system. We then discuss the effect of modifying amplitude and geometrical distribution of the charges in the above-defined model systems. Finally, we link the electrostatic field of the model system to the adsorption isotherm, via a calculation of the Henry constant, which is then injected in an adsorption model derived from statistical thermodynamics (as in our previous work on nonpolar systems). In this way, we are able to calculate the working capacity of the model solids in a pressure swing adsorption process as a function of the charge distribution. We can, hence, determine the optimal charge distribution that will maximize the working capacity of the model solid. In contrast to our previous work, which only treated dispersion−repulsion forces, the present paper exclusively addresses the issue of electrostatic forces. The integration of dispersion−repulsion forces and electrostatic interactions into a unified model is beyond the scope of this contribution and will be treated in a follow-up manuscript.

pore geometries (spheres and cylinders) and idealized charge distributions. This will provide a better understanding of the effect of the geometry and the magnitude of the charge on the behavior of the adsorption system. We illustrate our line of thinking by a few concrete examples: In zeolites without extraframework cations, only oxygen atoms are in contact with the adsorbate molecules. The oxygen atoms carry a negative partial charge. The silicon (or aluminum and phosphorus atoms in AlPOs) atoms are in a second layer behind the oxygen atoms. They carry a positive partial charge. One could think of representing this charge distribution by means of an inner ring of negative charges and an outer ring of positive charges, i.e., two concentric cylinders for channel-like pore systems or two concentric spheres for cage-like pore systems. To avoid attributing atomic positions, the total negative charge could be evenly distributed over the inner cylinder or sphere and the positive charge evenly distributed over the outer cylinder or sphere. However, it is easy to show using Gauss’s law that for these highly symmetric surfaces the electric field is zero in the volume embedded by the surface,12 i.e., inside the volume of the pore. The interaction potential energy of a linear quadrupolar molecule with an electric field (ξ) depends on the directional gradient of the field13,14 EΘz = 1/2Θ

2. SIMPLIFIED REPRESENTATION OF THE SOLID’S CHARGE DISTRIBUTION In force-field-based Grand Canonical Monte Carlo (GCMC) simulations of adsorption, electrostatic forces are usually dealt with in the following way: a charge is attributed to each atom of the solid. These point charges are supposed to reproduce the electrostatic field generated by the atoms in the solid. The adsorbate molecule is generally also represented by an ensemble of point charges. The overall solid−molecule electrostatic interaction can be calculated by pairwise addition of individual charge−charge Coulombic interactions. The longrange nature of this interaction makes its convergence with distance slow. Therefore, to calculate the overall electrostatic interaction of the solid, the Coulombic interaction between charges located at large distances cannot be neglected. Hence, the direct summation of all the charge−charge interactions becomes a very costly and slow method. For systems where boundary periodic conditions can be applied, the Ewald method allows overcoming this problem. The idea is to replace a single divergent summation with two convergent summations. The first summation has a convergent summation in the form of its Fourier transform, and the second has a convergent direct summation. The splitting of the summation is achieved by craftily adding and subtracting to the system a Gaussian distribution of charge. The method has been widely adopted for estimating the electrostatic contribution in atomistic simulations. The charges of the solid atoms are often attributed based on the results of DFT calculations.9−11 If the charges are correct, the above-mentioned approach leads to an accurate description of the adsorption potential. However, since each adsorbent has its proper charge distribution, it is nearly impossible to deduce general rules of the behavior of adsorption systems from this approach. To be able to describe the electrostatic forces by a few descriptors, we need a simpler representation. Instead of testing real solids, we wanted to evaluate the effect of changing the amplitude and density of charges in model systems, with simple

dξ dv

(1)

where v is the vector defining the direction of the molecules axis of symmetry and Θ is the quadrupolar moment. If the field is zero at any point in the pore volume, so is the gradient. Therefore, the interaction of CO2 with these uniformly charged surfaces is zero inside the pore. It is the deviation from a perfect spherical or cylindrical symmetry, generated by concentrating the positive charges at discrete positions that generates an electric field. This explains, by the way, why purely siliceous zeolites can be treated as nonpolar solids although they do have positively charged silicon atoms and negatively charged oxygen atoms. The regular and quite homogeneous distribution of the positive and negative charges in the solid leads to an overall electric field that is quite weak. This is not the case any more for zeolites with extraframework cations. In these solids, the positive charge is concentrated in the cations, while the negative charge is still quite evenly distributed over the surrounding framework oxygen atoms. Moreover, the cation positions in zeolite cages often correspond to highly regular geometric structures. For example, the most preferred cation sites in the supercage of zeolite A are the centers of the sixmembered ring windows delimitating the supercage. The centers of these six-ring windows form a cube. Also in the faujasite structure, the six-ring windows in the supercage are the preferred cation sites, but the geometrical arrangement is different: they form a tetrahedron. We can, therefore, use a cube and a tetrahedron as geometrical structures for placing the positive point charges (the cations). Since the negative charges are more delocalized, we distribute them homogeneously on the surface of a sphere that represents the pore surface, i.e., the oxygen atoms of the zeolite cage. Depending on the position of the cation in the zeolite cage, the diameter of the negatively charged pore surface may be equal to or larger than the diameter of the cube/tetrahedron defining the cation positions. In this way, it is possible to account for the fact that larger cations are not located in the plane of the six-ring window but shifted toward the center of the cage because of steric repulsion with the oxygen atoms. In any case, following our reasoning 9459



Article

distributions, keeping in mind that the homogeneously negative charged surface is always present.

above, the negatively charged sphere hardly contributes to the electrostatic field. The field is mainly created by the cations. This simplified representation of the charge distribution can also be applied to zeolite-like metal−organic framework (ZMOF) structures.15 It is probably less well adapted for MOFs with coordinatively unsaturated sites (CUS). Still, perfectly reproducing the electrostatic field of the real solid is not necessarily the main objective of our work. The model systems provide us with a tool to play with amplitude of the charges and their mutual distance (which is linked to the pore size) and to observe the effect on CO2 adsorption. For our case study we used two geometrical charge distributions: (1) a cubic arrangement of positive point charges in a spherical pore and (2) a hexagonal arrangement of positive point charges in a cylindrical pore. The first one is meant to simulate the charge distribution of zeolite A with 64 cations per unit cell, where all the centers of six-ring windows are occupied with a cation (see Figure 1b). The pore surface was represented

3. CALCULATION OF THE ELECTROSTATIC POTENTIAL ENERGY BETWEEN SOLID AND CO2 The carbon dioxide molecule has a permanent linear quadrupole moment. The quadrupole moment (Θ) is a second-order tensor and its components are defined by19 Θij =

∫ ρc rri jdr

(i , j = x , y , z )

(2)

where the ρc is the molecule charge density distribution and ri and rj are the components of the distance vector (r). Usually, the origin of the system is placed at the molecule center of mass, and the integral is carried out over the whole volume of the molecule. If the electric charge distribution can be represented by point charges (z), we obtain the well-known definition of the primitive quadrupole

Θij =

∑ zsrsirsj s

(3)

For linear charge distributions only diagonal elements have a nonzero value. Therefore, only one scalar can characterize the tensor. This scalar is commonly referred to as the quadrupole moment (Θ).20 It is formally defined as Θ = Θzz − Θxx

(4)

Considering that the carbon atom is placed at the origin and the oxygen atoms are placed along the z-axis and using eqs 3 and 4 the quadrupole moment of CO2 can be calculated by Θ = −2zOrC2 = O

(5)

where z0 is the partial charge on the oxygen atoms and r is the distance between the oxygen and carbon atoms. The experimental CO2 quadrupole moment is −14.31 ± 0.74 × 10−40 Cm.2,21,22 This value is larger than those of other simple molecules such as N2, O2, H2, COS, and ethane.13,21,22 3.1. Evaluation of the Electric Field Gradient. In the literature many approaches have been used to study electrostatic interactions between quadrupolar molecules and adsorbent materials.13,14,23−27 As mentioned before, the most common method consists of summing up the Coulombic potential between all possible pairs of point charges between adsorbate and adsorbent, as is usually done in the framework of atomistic Monte Carlo simulations.23−27 Alternatively, one can calculate the electric field gradient (EFG) generated by the solid and directly estimate its interaction with the adsorbate quadrupole moment.13,14 We initially followed the latter approach and attempted to calculate the EFG−CO2 interaction by evaluating eq 1. The differential term in this equation is the rate of change of the electric field along the molecule’s axis of symmetry. The rate of change of the electric field is usually known as the electric field gradient. Therefore, we need to compute the EFG in the direction of the molecule’s axis of symmetry. By going through some tensor algebra (see Supporting Information), one can show that the electrostatic energy is proportional to the norm of the directional derivative of the electric field

Figure 1. (a) Charges forming a cube, (b) unit cell of LTA with 64 cations per unit cell, (c) charges forming a hexagon, and (d) pore of CPO-27-Ni. Atoms: Ni, blue; O, red; C, gray; H, white, Na, purple.

by a sphere running through the corner of the cube that defines the cation positions. The negative charge was distributed over the surface of the sphere by placing small point charges with a mesh density of 0.2 point charges per Å2. The magnitude of the negative point charges was appropriately chosen to compensate the total positive charge, thus the overall system is neutral. The second system considered is a one-dimensional cylindrical pore forming a hexagonal lattice. The cations were placed to form a hexagon. The radial distance between the center and the vertices of the hexagon is the same as the radius of the cylindrical negative charged surface. This system is similar to honeycomb MOFs with coordinatively unsaturated sites, e.g., CPO-27 (also known as MOF-74), STA-12, and STA-16.16−18 To study a system similar to those structures the distance between successive rings of positive charges was fixed at 2.9 Å, which is approximately the distance between rings of CUS in CPO-27. For the sake of simplicity hereafter we will refer to these systems as the cubic and the hexagonal charge

EΘz = 1/2Θ || Dvξ ||

(6)

The electrostatic potential is a function of the position of the quadrupole and the orientation of its axis of symmetry. A 9460



Article

The CO2−solid electrostatic interaction depends on the position in the pore, r = (x,y,z), and the orientation of CO2 molecules defined by the polar and azimuth angles II′ and ω′. The average potential energy at each position, E̅ (x,y,z)elec, can be found by averaging the solid−CO2 potential over all the orientations with a Boltzmann factor according to eq 11.

rigorous treatment of this interaction should consider all possible directions of the molecule axis. We made the simplification that if the electric field is large enough the molecules will stay aligned with the direction that produces the strongest interaction. For a weak electric field or at high temperatures, molecules can freely change their directions, hence this approximation is no longer applicable, but these are not the situations that we are most interested in. The perfect alignment of the CO2 axis with the EFG corresponds to the maximal electrostatic energy that can be reached at each point. It can be calculated without knowing the direction of the EFG because by definition the norm of a matrix is the maximum value that the norm of a vector obtained as the product of a matrix and a vector can reach28

EΘz,max ≈ 1/2Θ ||∇ξ ||

2π

E ̅ (r)elec =

|| v̂ ||s = 1

π

The Cartesian coordinates defining the position of CO2 molecules can be converted into spherical coordinates (r, Π, ω). Then, for a spherical pore, the average potential energy at a given radial distance (r) is given by

(7)

E (r )elec,sph = 2π

π

∫0 ∫0 Eelec ̅ (r , Π, ω)exp( − Eelec ̅ (r , Π, ω)/kT )sin ωdωdΠ

(8)

2π

π

∫0 ∫0 exp(−Eelec ̅ (r , Π, ω)/kT )sin ωdωdΠ

where the index s defines the type of matrix norm. The norm that satisfies eq 8 is the 2-norm (s = 2) also known as the spectral norm or the induced norm.29 The 2-norm can be defined as the largest singular value of the matrix. To obtain this value, the characteristic polynomial of the matrix must be solved. A 3 × 3-matrix leads to a third degree polynomial. Hence, the roots can be found by the Cardon’s method or any other iterative approach. The main advantage of this method is that the EFG converges faster than the electric potential (Coulomb’s law) calculation. Also, it gives directly the orientation of the CO2 molecule that maximizes the electrostatic interaction. However, very complex equations are obtained even for simple charge distributions. Moreover, numerical methods are involved in the calculation of the EFG norm. Since a purely analytic solution of the mathematical problem is not accessible, the calculation of the EFG presents little advantage compared to a direct numerical Coulombic calculation. This method is detailed in the next section. 3.2. Direct Coulombic Calculations. The solid−CO2 electrostatic interactions can be evaluated by Coulomb’s law through the Ewald summation technique. The quadrupole moment of CO2 is depicted by three point charges, according to the model proposed by Harris and Yung.30 The solid is represented by a set of point charges, applying the simplifications discussed in Section 2. The electrostatic potential generated by the solid at any point r is, thus, calculated from z 1 V (r) = ∑ i 4πε0 i |r − ri| (9)

(12)

By using cylindrical coordinates (r, ω, z), the cylindrical pore can be treated as E (r )elec,cyl 2π

=

l

∫0 ∫0 Eelec ̅ (r , ω , z)exp( − Eelec ̅ (r , ω , z)/kT )dz dω l

2π

∫0 ∫0 exp(−Eelec ̅ (r , ω , z)/kT )dz dω (13)

The length of the cylinder (l) was chosen to be equal to the lattice parameter in the direction of the cylinder axis of symmetry. For cylindrical and spherical pores the average electrostatic energy in the pore (⟨Eelec⟩) is obtained by integrating over the radial position ⟨Eelec,sph⟩ =

∫0

Racc

Eelec,sph(r )exp( −Eelec,sph(r )/kT )r 2dr

∫0

Racc

exp( −Eelec,sph(r )/kT )r 2dr (14)

⟨Eelec,cyl⟩ =

∫0

Racc

Eelec,cyl(r )exp( −Eelec,cyl(r )/kT )r dr

∫0

Racc

exp( −Eelec,cyl(r )/kT )r dr (15)

All these average potential energies are temperature dependent. The average electrostatic energy ⟨Eelec⟩ can be calculated by integrating the electrostatic potential up to the accessible pore radius (Racc). The accessible pore radius is obtained by subtracting the solid−fluid collision distance σsf from the pore radius R

where ri and zi are, respectively, the position and the charge value of the ith particle in the solid and ε0 is the permittivity in vacuum. The sum is carried out over all the point charges in the solid. The electrostatic energy between CO2 at a given position in the solid is then calculated by summing up the interaction between the point charges (zC, zO) placed at the atomic positions of both the carbon (r) and the oxygen (r + d and r − d) atoms of the CO2 molecule and those representing the solid. This calculation is done according to Eelec(r) = zC·V (r) + zO·V (r + d) + zO·V (r − d)

2π

∫0 ∫0 exp(−Eelec(x , y , z , Π′, ω′)/kT )dΠ′dω′ (11)

where we have used the fact that max ||∇ξ·v̂ ||s = ||∇ξ ||s

π

∫0 ∫0 Eelec(Π′, ω′)exp(−Eelec(Π′, ω′)/kT )dΠ′dω′

R acc = R − σsf

(16)

In our model solids, R corresponds to the radius of the sphere or cylinder carrying the negative charges. σsf was set to 3.36 Å, which is valid for the interaction of CO2 with zeolites.5 To estimate the effectiveness of the electrostatic interactions in the cavity, let us define the following function ψ by

(10) 9461


The Journal of Physical Chemistry C ψ=

∫0

Racc

∫0

(exp( −Eelec,sph(r )/kT ) − 1)r 2dr Racc

exp( −Eelec,sph(r )/kT )r 2dr

Article

The fact that the simplified charge distribution reproduces quite accurately the electrostatic field of the real solid allows us to explore the effects of both pore size and charge on the adsorption of CO2 by means of the simplified picture. This is addressed in the next section. 4.2. CO2−Solid Electrostatic Interactions for Cubic and Hexagonal Charge Distributions. 4.2.1. Cubic Charge Distribution. We first discuss the shape of the electrostatic field in the model solid that represents zeolite Na64LTA. Figure 3(a) and (b) shows, respectively, the potential energy profile over a plane parallel to the faces intersecting at the center of the cubic lattice and at one of the lattice faces. The electrostatic potential at the center of the lattice is zero. It is also close to zero in the middle of the edges. In these highly symmetric zones of the unit cell (inversion points), the electric field gradient is zero. Thus, the electrostatic interaction of CO2 with the solid vanishes. On the other hand, the CO2−solid electrostatic interaction is maximized along the diagonal between the center of the unit cell and the cationic sites. This direction has the strongest EFG. To increase its interaction with the EFG, CO2 molecules must be aligned in the direction of the diagonal of the cube. Figure 4 shows the electrostatic interaction of CO2 as a function of the pore radius and the charge of the cations. Concerning the latter parameter |⟨Eelec⟩| steadily increases as the charge is increased. When looking to the effect of the pore radius, a region of maximal interaction at a pore radius of 6.8 Å is identified. This can be explained as follows: for small cavities the high symmetry of the charge distribution causes the CO2− solid interaction to vanish. As the pore radius is increased, the region of the pore where the EFG is strong is increased. However, for large pores, a region where CO2 electrostatic interactions are ineffective is added. Therefore, the average electrostatic interaction decreases. The trend is confirmed by Figure 5, which shows the evolution of the ψ function, i.e., the fraction of the pore volume where the electrostatic potential is nonzero (see eq 17). 4.2.2. Hexagonal Charge Distribution. For comparison, we have studied the effect of the pore radius and the magnitude of the charge on |⟨Eelec⟩| for a hexagonal charge distribution. The cation−cation distance along the axis of the cylinder was kept constant in this exercise. For reasons of symmetry, the electric field components cancel out along the central axis of the pore and in the middle of the hexagonal edges. Figure 6 shows the pore radius dependence of ⟨Eelec⟩ in the hexagonal system. We observe a maximum in the electrostatic potential at R = 4 Å. Below this pore radius, the central symmetry (CO2 is obliged to be close to the center of the pore where the EFG is zero) as well as the hard-sphere repulsion of the pore wall reduce the interaction potential. The decrease of the electrostatic potential between 4 and 8 Å is also attributed to symmetry effects (see Figure 7, which shows the fraction of the pore where the potential is nonzero). At large pore radii, the hexagonal system converges toward a constant ⟨Eelec⟩. Up to 15 Å we do not observe the drop of the electrostatic potential at high pore radii that was observed in the cubic system. The reason for this difference arises from the conception of both models. While in the case of the cubic distribution the distance between electrostatic charges is homogeneously varied in the three dimensions of the space, in the hexagonal case we choose to keep the cation−cation distance constant along the cylindrical symmetry axis. There-

(sphere) (17)

ψ=

∫0

Racc

∫0

(exp( −Eelec,cyl(r )/kT ) − 1)r dr Racc

(cylinder)

exp( −Eelec,cyl(r )/kT )r dr (18)

The function ψ provides the fraction of the pore volume where CO2 molecules are affected by the electrostatic potential. This is the section of the pore where Eelec ≫ kT. For cavities where molecules have an important electrostatic interaction in the whole pore volume ψ is equal to 1. Conversely, for cavities where the electrostatic potential is significantly lower than kT, ψ is zero. Further technical details of the calculation method (Ewald summation, etc.) are provided in the Supporting Information.

4. RESULTS 4.1. Validation of the Simplified Model. To validate the methodology used to depict the electrostatic field generated by the solid, the electrostatic potential calculated using a full atomistic model and the simplified model were compared. We compared the electrostatic field of an LTA zeolite with 64 Na+ atoms per unit cell with our model system, i.e., a spherical cavity with point charges placed at the corner of a cube and negative charges distributed over the surface of the sphere (see Figure 1). The diagonal distance between the center of the cube and the corner was fixed to 6.21 Å. The radius of the negatively charged sphere had the same value. This distance matches the pore radius of the super cage in LTA. The charges used in the full atomistic representation of the Na64LTA zeolite were determined from DFT quantum mechanics calculations. Details about these calculations can be found in the Supporting Information. Figure 2 shows that there is a good agreement

Figure 2. Comparison of the electrostatic potential between CO2 in a cubic charge distribution (cationic charge equal to +0.9e) and a full atomistic representation of the Na64LTA (charges determined by DFT calculations).

between both representations when the charge of the cations used in the simplified representation is equal to +0.9e. This value is close to that used for Na atoms in the full atomistic representation (+1.00e). The simplified representation of the solid can reproduce the shape and the depth of the potential well. 9462



Article

Figure 3. Average electrostatic interaction potential of a model solid with a cubic charge distribution: (a) potential profile at the face [0,0,1] passing through the middle of the lattice and (b) potential profile in one of the faces [0,0,1] passing in the edge of the lattice.

Figure 7. ψ for a hexagonal lattice as a function of the pore radius. Figure 4. Average electrostatic potential for a cubic lattice as a function of the pore radius.

faster values close to one when the charge is increased and the pore radius is reduced than in the cubic pore. In the case of the hexagonal distribution, the electrostatic interaction is more efficiently distributed in the pore volume than in the case of the cubic configuration. 4.3. Calculation of the Working Capacity of the Model Solids. The ultimate goal of our work is to determine which model charge distribution would lead to an optimal performance of the solid in a PSA process. Since PSA processes cycle between adsorption and desorption steps, the best adsorbent is not necessarily the material with the highest heat of adsorption (i.e., the highest average electrostatic potential). We have to find the right balance between a fairly strong adsorption and an easy regeneration in the desorption step. This trade-off can be quantified by calculating the so-called working capacity of the solid. The working capacity is defined as the difference between the adsorbed amount under adsorption conditions (pressure, composition, temperature) and under desorption conditions (pressure, composition, temperature). Selectivity for adsorption of CO2 vs the other components in the gas mixture is another important criterion, but it will not be treated in the present contribution since we only focus on the adsorption of CO2. The strategy for establishing a link between the model charge distribution and the working capacity is the same as used in our previous work on dispersion−repulsion forces: integration of the electrostatic potential over the pore volume yields the Henry constant of adsorption.

Figure 5. ψ for a cubic lattice as a function of the pore radius.

Figure 6. Average electrostatic potential for a hexagonal lattice as a function of the pore radius.

K sph =

fore, the overall charge density decreases more gradually when increasing pore radius. Figure 7 reports the evolution of the function ψ for the hexagonal charge distribution. It can be noticed that ψ reaches

Kcyl = 9463

3 · 3 kT ·R acc 2 · 2 kT ·R acc

∫0

∫0

R acc

R acc

exp[−Eelec,sph(r )/kT ]r 2dr exp[−Eelec,cyl(r )/kT ]r dr

(19)

(20)



Article

where K is the Henry constant normalized by the accessible pore volume. This Henry constant is then injected into the Ruthven Statistical Model (RSM),31 which allows calculating the fractional loading of the pore volume as a function of pressure, i.e., the isotherm. The form of the RSM equation, as used in the present work, is given in eq 21 θ=

1 smax s s s ⎡ ⎤ s 2εff s smax (Kp) (bsmax ) ⎢ Kp(bsmax ) + ∑s = 2 (s − 1) ! 1 − smax exp smaxkT ⎥ ⎢ ⎥ s s 2 s εff ⎢ 1 + Kp(bs ) + ∑smax (Kp) (bsmax) 1 − s ⎥ exp max s=2 s! smax smax kT ⎥ ⎣⎢ ⎦ (21)

(

(

) ( ) ) ( )

Figure 9. Henry constant calculated for the model solids, with cubic or hexagonal charge distributions, as a function of the pore radius. The dashed line represents the optimal value of the Henry constant for a PSA process operating between 5 and 1 bar CO2.

where θ is the fractional loading of the pore volume, p the pressure, b the excluded volume of CO2, smax the maximum number of CO2 molecules that fits into the pore volume, and εff the potential-well depth of the intermolecular CO2−CO2 interaction. The choice of the parameters b and εff for CO2 was discussed in our previous work.5 From the isotherm we can calculate the working capacity for an isothermal PSA process, as a function of the partial pressure of CO2 in the adsorption and desorption step. In particular, a fractional working capacity is calculated from the RSM model according to Δθ = θads(K , pads,CO2 ) − θdes(K , pdes,CO2 )

the case of the hexagonal system, for z = 1, the same qualitative trend is observed, but the maximal Henry constant is reached at a lower pore radius, and the maximum is more pronounced. Below a pore radius of 7 Å, the Henry constant in the hexagonal system is very high, but these high values of K are not optimal for the adsorption process. For both model solids, the optimal values of K are approached when the pore radius is in the range of 7−9 Å, if the cation charge is z = 1. In the cubic system, lowering the cation charge moves us away from the target K value. Lowering the cation charge to 0.5 already reduces the electrostatic field too much and leads to suboptimal Henry constants, irrespective of the pore radius of the charge symmetry. Figure 10 translates the result of Figure 9 into a fractional working capacity. For the two systems with a cation charge of 1,

(22)

Figure 8 shows the fractional working capacity as a function of the Henry constant for an adsorption pressure of 5 bar CO2

Figure 8. Fractional working capacity of CO2 as a function of the Henry constant. pads = 5 or 10 bar and pdes = 1 bar.

and a desorption pressure of 1 bar CO2. For these operating conditions, the optimal Henry constant maximizing the fractional working capacity is 6 × 103 molecules·bar−1·Å−3. When repeating the calculation at a higher adsorption pressure of 10 bar (Figure 8) we surprisingly observe that the optimal value of K is almost not affected (it slightly decreases). Note, however, that the optimal K may be quite different when looking at PSA applications operating in much higher or lower pressure ranges as is the case of VSA (vacuum swing adsorption). Having defined the target value of the Henry constant, we calculated the Henry constants of our model solids as a function of the pore radius (Figure 9). The evolution of the Henry constant as a function of the pore radius follows the evolution of the average electrostatic potential (Figure 4 and Figure 6). For the cubic charge distribution, the Henry constant first increases with increasing pore radius and then decreases. In

Figure 10. Fractional working capacity calculated for the model solids, with cubic or hexagonal charge distribution, as a function of the pore radius, calculated for a PSA process operating between 5 and 1 bar CO2.

the fractional working capacity remains high in a quite large range of pore radii, i.e., from 6 to 12 Å for the cubic system and from 7 to 15 Å for the hexagonal system. At high pore radii the working capacity in the cubic system drops because the electric field does not cover the whole pore volume any more (Figure 5). This drop is less pronounced in the hexagonal system because it maintains a high charge density along the axis of the cylinder. To convert the fractional working capacity into a capacity per volume of adsorbent, θ must be multiplied with the porosity of the adsorbent φ (pore volume divided by the volume of the 9464



Article

design of adsorbents: make large pore solids containing highly charged cations. When the cation charge is less than one the electrostatic potential rapidly drops below the optimal values (for the PSA applications in a partial pressure range of 10−1 bar). In MOFs, the bond between framework cations and the ligand is not fully ionic. Therefore, the effective cation charge is quite weak.34,35 We can deduce that MOFs, which only contain framework cations, will probably not be able to generate an electrostatic field that is sufficiently strong to optimize the interaction with CO2. The situation may be different for ZMOFs, which contain extraframework cations, as in zeolites. Note also that MOFs may interact with CO2 by other mechanisms than LennardJones forces and electrostatic interactions.36,37

adsorbent). Although the exact relation between porosity and pore radius depends on the pore geometry and the wall thickness of the pore, we can safely state that φ increases with increasing pore radius. That means that within the range of optimal pore radii in Figure 10 the higher values are more favorable to optimize the working capacity on a per volume base.

5. DISCUSSION 5.1. Symmetry Effects in the Cubic and Hexagonal Charge Distributions. In an adsorbent that interacts with the adsorbate only via dispersion−repulsion forces, the optimal pore radius is defined by the trade-off between repulsion and attraction. Low pore radii lead to low Henry constants because repulsion dominates. The Henry constant is maximal, when the pore radius is close to the radius of the adsorbate and decreases beyond that value. For electrostatic interactions, the evolution of the Henry constant with the pore radius is more complex and depends on the symmetry of the charge distribution. In cubic arrangements, the high symmetry of the system leads to a mutual cancellation of the electrostatic fields of each individual cation. This effect is very marked at low pore radii and then decreases, when the fields of different cations do not interfere with each other anymore. As a result, the Henry constant increases. Further augmentation of the pore radius moves the charges too much away from each other, and their electrostatic field cannot cover the whole pore volume any more. Hence, the Henry constant drops. The above-described general trend of an initial rise and subsequent fall of the electrostatic interaction potential are perturbed by a shallow minimum of the magnitude of ⟨Eelec⟩ at 8 Å (Figure 4). This phenomenon is probably linked to additional inversion symmetry points, which are located at 1/4 of the cube’s diagonal axis. For a hexagonal charge distribution, i.e., a cylindrical pore, with constant distance between the charges along the cylinder axis, the evolution of the Henry constant as a function of the pore radius is different. The hexagonal system is less centrosymmetric than the cubic one. As a result, the cancellation of the electrostatic field components in the center of the pore is less pronounced. Low pore radii (down to 4 Å), therefore, lead to higher Henry constants than in the cubic system. On the contrary, the minimum of the average electrostatic interaction potential at R = 8 Å (Figure 6) is more pronounced. As in the cubic system, the minimum is ascribed to the presence of additional inversion points, which are located at 1/4 of the hexagon’s diagonal. Above R = 8 Å, the Henry constant increases slightly again, in spite of the dilatation of the charges in the radial direction. The drop of the Henry constant at high values of R that occurs in the cubic system is not observed because the cation−cation distance along the cylinder axis remains constant. 5.2. Strategy for Adsorbent Development. The porosity of an adsorbent increases with increasing pore size. Potentially, such large pore adsorbents may have a very high volumetric adsorption capacity. That explains the vivid interest in ultraporous MOF material with very large pore sizes.32,33 However, if CO2 interacts only via van der Waals forces with the adsorbent surface, the large pore volume is not efficiently exploited because the Lennard-Jones potential does not extend to the center of the pore. Electrostatic interactions are more far reaching. They are able to cover a large fraction of the pore volume up to high pore radii, provided that the cation charge is high (Figure 5 and Figure 7). This defines a strategy for the

6. SUMMARY AND CONCLUSIONS The selectivity of adsorbent for adsorption of CO2 often relies on electrostatic interactions between the electrostatic field of the solid and the quadrupole moment of CO2. The electrostatic field of the solid is generated by the charge distribution of the solid’s atoms. The ambition of this work was to represent the charge distributions of real solids by simplified model systems. For these model systems, we then analyzed qualitatively and quantitatively what happens when the pore radius and the cation charge changes. Our approach for defining model charge distributions was the following: according the Gauss law, the electric field of perfectly symmetric and homogeneous charge density distributions is zero. It is the deviation from symmetry and the concentration of homogeneous charge density into localized point charges that produces a strong electric field. We, therefore, used the heuristic that the charge distribution of a zeolite with extraframework cations can be described by a set of positive point charges, corresponding to the cation positions and a negative charge that are evenly distributed over the pore surface. They correspond to the lattice oxygen atoms. The field is mainly generated by the cations. The electrostatic field generated by this extremely simplified model is close to that of the real zeolite. The analysis of the model charge distribution with cubic and hexagonal arrangement of positive point charges shows that symmetry plays a big role. At low pore radii, the high symmetry of the system leads to a mutual cancellation of the fields of the individual cations at low pore radii, i.e., when the individual fields interpenetrate. Increasing the size of the systems reduces the interpenetration and, therefore, increases the overall electrostatic potential up to a maximum value. Above this value, the cations are so far from each other that their fields do not interpenetrate but also cannot cover the whole pore volume any more. The different symmetries of the hexagonal and cubic system lead to differences in the evolution of the electric field as a function of pore radius, which proves that the geometrical arrangement of the cation positions is an extremely important parameter (as important as the overall charge density). For the two aforementioned model charge distributions, we have calculated the pore radii that will maximize the working capacity of the adsorbent in an adsorption process working between 5 and 10 bar adsorption pressure (partial pressure of CO2) and 1 bar desorption pressure. The highest working capacities are obtained when the pore radius is in the range of 6−12 Å for the cubic system and 7−15 Å for the hexagonal system. The optimal range extends to higher radii than for solids that only interact with CO2 via van der Waals 9465



Article

(13) Buckingham, A. D.; Disch, R. L.; Dunmur, D. A. Quadrupole Moments of some Simple Molecules. J. Am. Chem. Soc. 1968, 90, 3104−3107. (14) Drain, L. E. Permanent Electric Quadrupole Moments of Molecules and Heats of Adsorption. Trans. Faraday Soc. 1953, 49, 650−654. (15) Liu, Y. L.; Kravtsov, V. C.; Larsen, R.; Eddaoudi, M. Molecular Building Blocks Approach to the Assembly of Zeolite-Like MetalOrganic Frameworks (ZMOFs) with Extra-Large Cavities. Chem. Commun. 2006, 1488−1490. (16) Dietzel, P. D. C.; Morita, Y.; Blom, R.; Fjelag, H. An In Situ High-Temperature Single-Crystal Investigation of a Dehydrated Metal-Organic Framework Compound and Field-Induced Magnetization of One-Dimensional Metal-Oxygen Chains. Angew. Chem. 2005, 117, 6512−6516. (17) Wharmby, M. T.; Mowat, J. P. S.; Thompson, S. P.; Wright, P. A. Extending the Pore Size of Crystalline Metal Phosphonates toward the Mesoporous Regime by Isoreticular Synthesis. J. Am. Chem. Soc. 2011, 133, 1266−1269. (18) Rosi, N. L.; Kim, J.; Eddaoudi, M.; Chen, B. L.; O’Keeffe, M.; Yaghi, O. M. Rod Packings and Metal-Organic Frameworks Constructed from Rod-Shaped Secondary Building Units. J. Am. Chem. Soc. 2005, 127, 1504−1518. (19) Glaser, R.; Wu, Z.; Lewis, M. A Higher Level ab initio QuantumMechanical Study of the Quadrupole Moment Tensor Components of Carbon Dioxide. J. Mol. Struct. 2000, 556, 131−141. (20) Greenhow, C.; Smith, W. Molecular Quadrupole Moments of N2 and O2. J. Chem. Phys. 1951, 19, 1298−1300. (21) Chetty, N.; Couling, V. W. Measurement of the Electric Quadrupole Moments of CO2 and OCS. Mol. Phys. 2011, 109, 655− 666. (22) Graham, C. Measurement of the electric quadrupole moments of CO2, CO, N2, Cl2 and BF3. Mol. Phys. 1998, 93, 49−56. (23) Yang, Q.; Zhong, C.; Chen, J. F. Computational Study of CO2 Storage in Metal-Organic Frameworks. J. Phys. Chem. C 2008, 112, 1562−1569. (24) Yang, Q.; Zhong, C. Electrostatic-Field-Induced Enhancement of Gas Mixture Separation in Metal-Organic Frameworks: A Computational Study. ChemPhysChem 2006, 7, 1417−1421. (25) Wang, S.; Yang, Q.; Zhong, C. Molecular Simulation Study of Separation of CO2 from Alkanes Using Metal-Organic Frameworks. J. Phys. Chem. B 2006, 110, 26507−26507. (26) Pirngruber, G. D.; Raybaud, P.; Belmabkhout, Y.; Cejka, J.; Zukal, A. The Role of the Extra-Framework Cations in the Adsorption of CO2 on Faujasite Y. Phys. Chem. Chem. Phys. 2010, 12, 13534− 13546. (27) Maurin, G.; Belmabkhout, Y.; Pirngruber, G.; Gaberova, L.; Llewellyn, P. CO2 Adsorption in LiY and NaY at High Temperature: Molecular Simulations Compared to Experiments. Adsorption 2007, 13, 453−460. (28) Meyer, C. D. Matrix Analysis and Applied Linear Algebra; Society for Industrial and Applied Mathematics: Philadelphia, 2000. (29) Basar, Y.; Weichert, D. Mathematical fundamentals. In Nonlinear Continuum Mechanics of Solids: Fundamental Mathematical and Physical Concepts; Springer: Berlin, 2000. (30) Harris, J. G.; Yung, K. H. Carbon Dioxides Liquid-Vapor Coexistence Curve and Critical Properties As Predicted by A Simple Molecular-Model. J. Phys. Chem. 1995, 99, 12021−12024. (31) Ruthven, D. M.; Loughlin, K. F. Sorption of Light Paraffins in Type-A Zeolites. Analysis and Interpretation of Equilibrium Isotherms. J. Chem. Soc., Faraday Trans. 1 1972, 68, 696−708. (32) Furukawa, H.; Ko, N.; Go, Y. B.; Aratani, N.; Choi, S. B.; Choi, E.; Yazaydin, A. O.; Snurr, R. Q.; O’Keeffe, M.; Kim, J.; Yaghi, O. M. Ultrahigh Porosity in Metal-Organic Frameworks. Science 2010, 329, 424−428. (33) Farha, O. K.; Eryazici, I.; Jeong, N. C.; Hauser, B. G.; Wilmer, C. E.; Sarjeant, A. A.; Snurr, R. Q.; Nguyen, S. T.; Yazaydin, A. O.; Hupp, J. T. Metal-Organic Framework Materials with Ultrahigh Surface

interactions. However, high working capacities are only achieved when the cation charge is close to 1 (or higher in a cubic system with large pores). For lower charges the adsorption of CO2 is not sufficiently strong (for the pressure range considered here). This allows defining a strategy for the design of solids for CO2 capture applications: we have to aim for large pore solids containing highly charged cations.

■

ASSOCIATED CONTENT

S Supporting Information *

Details about the calculation of potential energy of CO2 via the electric field gradient or via the Coulombic method; details about the DFT calculations; and some additional figures of the potential energy of CO2 as a function of the cation charge in a cubic and hexagonal charge distribution. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Author

*Fax: (+33) 437 702 066. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■ ■

ABBREVIATIONS PSA, pressure swing adsorption; EFG, electric field gradient; CUS, coordinatively unsaturated site. REFERENCES

(1) Martin, R. L.; Willems, T. F.; Lin, L. C.; Kim, J.; Swisher, J. A.; Smit, B.; Haranczyk, M. Similarity-Driven Discovery of Zeolite Materials for Adsorption-Based Separations. ChemPhysChem 2012, 13, 3595−3597. (2) Lin, L. C.; Berger, A. H.; Martin, R. L.; Kim, J.; Swisher, J. A.; Jariwala, K.; Rycroft, C. H.; Bhown, A. S.; Deem, M.; Haranczyk, M.; Smit, B. In Silico Screening of Carbon-Capture materials. Nat. Mater. 2012, 11, 633−641. (3) Wu, D.; Yang, Q. Y.; Zhong, C. L.; Liu, D. H.; Huang, H. L.; Zhang, W. J.; Maurin, G. Revealing the Structure-Property Relationships of Metal-Organic Frameworks for CO2 Capture from Flue Gas. Langmuir 2012, 28, 12094−12099. (4) Garcia, E. J.; Perez-Pellitero, J.; Jallut, C.; Pirngruber, G. D. Quantification of the Confinement Effect in Microporous Materials. Phys. Chem. Chem. Phys. 2013, 15, 5648−5657. (5) Garcia, E. J.; Perez-Pellitero, J.; Jallut, C.; Pirngruber, G. D. Modeling Adsorption Properties on the Basis of Microscopic, Molecular, and Structural Descriptors for Nonpolar Adsorbents. Langmuir 2013, 29, 9398−9409. (6) Nicholson, D.; Gubbins, K. E. Separation of Carbon DioxideMethane Mixtures by Adsorption: Effects of Geometry and Energetics on Selectivity. J. Chem. Phys. 1996, 104, 8126−8134. (7) Ducrot-Boisgointier, C.; Parmentier, J.; Faour, A.; Patarin, J.; Pirngruber, G. D. FAU-Type Zeolite Nanocasted Carbon Replicas for CO2 Adsorption and Hydrogen Purification. Energy Fuels 2010, 24, 3595−3606. (8) Peng, X.; Cao, D.; Zhao, J. Grand Canonical Monte Carlo Simulation of Methane-Carbon Dioxide Mixtures on Ordered Mesoporous Carbon CMK-1. Sep. Purif. Technol. 2009, 68, 50−60. (9) Manz, T. A.; Sholl, D. S. Chemically Meaningful Atomic Charges That Reproduce the Electrostatic Potential in Periodic and Nonperiodic Materials. J. Chem. Theory Comput. 2010, 6, 2455−2468. (10) Campana, C.; Mussard, B.; Woo, T. K. J. Chem. Theory Comput. 2009, 5, 2866. (11) Hirshfeld, F. L. Bonded-Atom Fragments for Describing Molecular Charge densities. Theor. Chem. Acc. 1977, 44, 129−138. (12) Zahn, M. Electromagnetic Field Theory: A Problem Solving Approach, 2nd ed.; Krieger Pub Co: Malabar, FL, 2003. 9466



Article

Areas: Is the Sky the Limit? J. Am. Chem. Soc. 2012, 134, 15016− 15021. (34) Babarao, R.; Jiang, J. W.; Sandler, S. I. Molecular Simulations for Adsorptive Separation of CO2/CH4Mixture in Metal-Exposed, Catenated, and Charged Metal-Organic Frameworks. Langmuir 2009, 25, 5239−5247. (35) Zheng, C.; Liu, D.; Yang, Q.; Zhong, C.; Mi, J. Computational Study on the Influences of Framework Charges on CO2 Uptake in Metal-Organic Frameworks. Ind. Eng. Chem. Res. 2009, 48, 10479− 10484. (36) Garcia, E. J.; Mowat, J. P. S.; Wright, P. A.; Pérez-Pellitero, J.; Jallut, C.; Pirngruber, G. D. Role of Structure and Chemistry in Controlling Separations of CO2/CH4 and CO2/CH4/CO Mixtures over Honeycomb MOFs with Coordinatively Unsaturated Metal Sites. J. Phys. Chem. C 2012, 116, 26636−26648. (37) Wu, H.; Simmons, J. M.; Srinivas, G.; Zhou, W.; Yildirim, T. Adsorption Sites and Binding Nature of CO2 in Prototypical MetalOrganic Frameworks: A Combined Neutron Diffraction and FirstPrinciples Study. J. Phys. Chem. Lett. 2010, 1, 1946−1951.

9467


How to Optimize the Electrostatic Interaction between a Solid

Recommend Documents