J. Phys. Chem. C 2008, 112, 11295–11300
11295
Computer Simulation for Adsorption of CO2, N2 and Flue Gas in a Mimetic MCM-41 Shengchi Zhuo,†,‡ Yongmin Huang,† Jun Hu,† Honglai Liu,*,† Ying Hu,† and Jianwen Jiang*,‡ Laboratory for AdVanced Material and Department of Chemistry, East China UniVersity of Science and Technology, Shanghai 200237, China, Department of Chemical and Biomolecular Engineering, National UniVersity of Singapore, Singapore 117576 ReceiVed: April 20, 2008; ReVised Manuscript ReceiVed: May 17, 2008
Adsorption of pure and mixed CO2 and N2 is simulated in a mimetic MCM-41. The full-atom MCM-41 model is constructed by caving cylindrical pores from an amorphous silica matrix and energetically optimized. Dreiding force field is used for the dispersive interaction with the atomic charges estimated from the densityfunctional theory calculations. The optimized MCM-41 maintains a hexagonal array of the mesoscopic pores as evidenced by the three characteristic peaks in the XRD pattern. The pore surface of MCM-41 is corrugated and coated with hydroxyls and defects. The pore size exhibits a Gaussian distribution with an average radius of 14.38 Å close to the experimental value. Simulated adsorption isotherms and isosteric heats of CO2 match well with the experimental data. CO2 adsorbs preferentially at the active sites near the pore surface, while N2 tends to adsorb homogeneously on the pore surface. In CO2 and N2 mixture as a flue gas, CO2 is more adsorbed than N2. The selectivity of CO2 over N2 drops rapidly with increasing temperature and depends weakly on pressure. At temperatures higher than 400 K, the selectivity approaches a constant and pressure has no discernible effect. 1. Introduction Greenhouse gas CO2 has increased by 30% in the atmosphere since the late 19th century. This has directly or indirectly caused a number of severe environmental impacts. The rapidly increasing concentration of CO2 is largely due to the vast amount of flue gases emitted by industrial and utility power generation plants using carbon-based fossil fuels. Flue gas is generally at or slightly above atmospheric pressure and typically consists of N2 and CO2. Sequestration of CO2 has been one of the most pressing issues in environmental protection. Prior to the sequestration, however, CO2 must be captured and separated from the flue gas. Extensive studies have been carried out in order to develop cost-effective techniques to capture and separate CO2, such as solvent absorption, cryogenic distillation, adsorption and bioconversion. Among these, adsorption in porous materials is regarded as cost-effective and technically feasible. Toward this end, it is crucial to select a proper sorbent. Within a variety of commonly used sorbents, mesoporous MCM-41 is potentially useful. MCM-41 was first synthesized by Beck and Kresge in 1992.1,2 Composed of silica, MCM-41 contains a hexagonal array of nearly uniform mesoscopic pores. The pore diameter can be controlled to vary from 15 Å to greater than 100 Å. As a result, MCM-41 is able to accommodate relatively large molecules, which cannot intercalate microporous molecular sieves. MCM-41 also possesses a large specific area and a pore volume of about 700 m2/g and 0.7 cm3/g, respectively. In addition, the pores of MCM-41 are geometrically simple and do not interconnect with one other. Along with experiments reported for the use of MCM-41 in catalysis, gas adsorption and separation, there have been a * Corresponding authors. E-mail:
[email protected] (H.L.); chejj@ nus.edu.sg (J.J.). † Laboratory for Advanced Material and Department of Chemistry, East China University of Science and Technology. ‡ Department of Chemical and Biomolecular Engineering, National University of Singapore.
number of theoretical and simulation studies. However, a reliable model of MCM-41 has to be built first. In most cases, MCM41 is simplified as a smooth and homogeneous cylindrical pore and a one-dimensional potential is then used to describe the interaction between the skeleton and fluid molecules. For instance, Ravikovitch et al.3 studied the adsorption of N2 in MCM-41 by a nonlocal density functional theory with the TFM model,4 which is a widely used one-dimensional potential model. Maddox and Gubbins5 also used one-dimensional model in Monte Carlo simulation to study the adsorption of Ar and N2 in MCM-41 and buckytube. One drawback in the onedimensional potential model is that the solid-fluid potential and the fluid density vary only in the direction normal to the pore wall and there is a strong layering in the adsorbed phase.6 As an improvement, Maddox and Gubbins7 developed a twodimensional potential model, which is radial and axial dependent. In their simulation work, the skeleton of MCM-41 were composed of randomly distributed O atoms and divided into eight sections. Each section exerted individual interaction with fluid molecules. Koh et al.8 built a more realistic full-atom model of MCM-41. They inserted randomly O and Si atoms into a box based on two criteria: (1) the position was not located in the free volume of the pore (2) the position was not within 3 Å from any other atoms. He and Seaton9 proposed three different full-atom models of MCM-41. Model 1 contained three concentric cylinders of O atoms with a layer of Si atoms between each layer of O atoms. This model was homogeneous at the atomic scale. Models 2 and 3 were generated by removing the O and Si atoms from R-quartz and amorphous silica matrixes, respectively. Benoit et al.10 also reported two models of MCM41 with different surface roughness. In the work, we construct a highly mimetic model of MCM41 by molecular-mechanical optimization. Adsorption of pure and mixed CO2 and N2 in the constructed MCM-41 model is subsequently investigated using computer simulation. In section 2, the mimetic model of MCM-41 is described, followed by
10.1021/jp803428n CCC: $40.75 2008 American Chemical Society Published on Web 07/04/2008
11296 J. Phys. Chem. C, Vol. 112, No. 30, 2008
Zhuo et al. TABLE 1: Lennard-Jones Potential Parameters and Atomic Charges molecule
site
σ(Å)
ε/kB(K)
q (e)
MCM-41
Si O H C O N
3.804 3.033 2.846 3.473 3.033 3.263
155.858 48.115 0.0503 47.813 48.115 38.914
0.1222 -0.06157 0.0318 0.604 -0.302 0
CO2 N2
Figure 1. (a) Amorphous silica matrix; (b) MCM-41 model before optimization; (c) MCM-41 model after optimization. Key: Si, yellow; O, red; H, white.
Figure 2. Clusters in the DFT calculations for atomic charges. Key: Si, yellow; O, red; H, white.
simulation methodology for adsorption. In section 3, we characterize the MCM-41 model and compare the simulation results with experimental results in terms of adsorption isotherm and isosteric heat. Concluding remarks are summarized in section 4. 2. Model and Methodology The mimetic MCM-41 was constructed by carving out cylindrical pores in an amorphous silica matrix as shown in Figure 1, parts a and b. These pores have radius of about 15.0 Å and are hexagonally arranged. During the construction, unsaturated Si and O atoms were generated. To maintain the original hybridization, O and H atoms were added separately to saturate the cleaved bonds. As a result, hydroxyls were introduced to the pore surface in a random manner. However, Si atoms linked with three hydroxyls were removed. The MCM41 model generated in this way possesses 7-8 hydroxyls per square nanometer.10,11 Before the constructed MCM-41 model was subject to molecular-mechanism optimization, the atomic charges were calculated using the density-functional theory (DFT) as implemented in Dmol3 package.12 The charge distribution of framework directly affects the interaction between sorbent and sorbate, therefore, plays an important role in determining adsorption properties, especially for polar sorbates. Figure 2 shows three different clusters used in the DFT to calculate atomic charges. Model a contains a number of tetrahedral SiO4 and hydroxyl groups, which can form hydrogen bonds. Thus, the effect of
hydrogen bond on atomic charge is examined. Model b bears a six-member ring and an eight-member ring, which allows to identify the influence of ring structure. Model c is an atomic cluster cut from the amorphous silica matrix. Each cluster was geometrically optimized with Perdew-Wang functional13 and DNP (double numerical plus polarization) basis set. In addition, the semicore pseudopotential was used. This allows only the softer valence electron wave functions to be explicitly treated, which is usually the portion that controls the chemistry and can significantly reduce computational cost. We would point out that no unique method is available currently to determine the atomic charges. The most commonly used methods include the electrostatic potential (ESP) fitting,14–17 Mulliken popular analysis18 and Hirshfeld popular analysis.19 Mulliken analysis largely depends on the basis set used and usually overestimates the charges, while Hirshfeld method usually underestimates the charge. In ESP method, the electrostatic potentials are fitted at grids located with equal density on different layers around the molecule. In this work, the atomic partial charges were obtained by ESP method. The charges calculated from the three different clusters were found to be approximately the same, revealing that the atomic charge is insensitive to the chosen cluster. The averaged charges of the three clusters are listed in Table 1. The dispersive interactions in MCM-41 were modeled by Dreiding force field.20 Dreiding force field includes a harmonic valence term, a cosine-Fourier expansion torsion term, and an umbrella functional form according to the Wilson out-of-plane definition. The van der Waals interactions were described by the Lennard-Jones 12-6 potential. Hydrogen bonding was represented explicitly by the Lennard-Jones 12-10 potential. Electrostatic interactions were described by the atomic monopoles and a screened (distancedependent) Coulombic term. Table 1 gives the Lennard-Jones potential parameters for the atoms in MCM-41 model. The MCM-41 model was optimized by molecular-mechanism method using Forcite module in Materials Studio.12 The simulation box was 64.2 × 149.8 × 42.8 Å and contained eight cylindrical pores. Periodic boundary conditions were exerted in all three dimensions. To achieve rapid convergence, three algorithms, namely, the steepest descent, adjusted basis set Newton-Raphson and quasi-Newtown methods were subsequently used. Both the van der Waals and hydrogen bond interactions were evaluated by the atom-based method with cutoff distance of 15.5 Å and 4.5 Å, respectively. Coulombic interactions were calculated by the Ewald sum with an accuracy of 10-4 kcal/mol. After the optimization, the MCM-41 model was characterized by the X-ray diffraction pattern. Adsorbates CO2 and N2 were mimicked by three- and twosite molecules, respectively. The dispersive interactions were obtained using the Dreiding force field with the parameters in Table 1. The lengths of C-O bond and N-N bond were 1.152 Å and 1.102 Å, respectively. CO2 molecule was optimized by Forcite package and the atomic charges were calculated by the Dmol3 module. The atomic charges of C and O atom were
Adsorption of CO2, N2 and Flue Gas
J. Phys. Chem. C, Vol. 112, No. 30, 2008 11297
Figure 3. Characterization of MCM-41 model: (a) XRD pattern of the optimized MCM-41 model. The inset is the snapshot. Key: Si, yellow; O, red; H, white. (b) Distribution of pore radius in the MCM-41 model. The curve is a Gaussian distribution fit to the data.
Figure 4. (a) Adsorption isotherms and (b) isosteric heats of CO2 in MCM-41 at 298 K from simulation and experiment.26 M41C14, M41C16, and M41C22 are three experimental samples of MCM-41 with radii of 13.50, 15.05, and 19.45 Å, respectively.
+0.604e and -0.302e, respectively, close to the charges in the elementary physical model for CO2.21 Grand canonical Monte Carlo simulations were performed to simulate gas adsorption in the MCM-41. Dreiding force field was used to describe the interactions of sorbate-sorbate and sorbate-sorbent. Electrostatic interactions were summed by the Ewald sum with an accuracy of 10-4 kcal/mol, while the van der Waals interactions were evaluated by atom-based method with cutoff length of 15.5 Å. Four types of trial moves were conducted randomly in the simulation, including deletion and insertion, rotation with the maximum amplitude of 5 degree and translation with the maximum step of 1 Å. The numbers of trial moves in a typical simulation was 1 × 107, but additional trial moves were used at high coverage. The first 5 × 106 steps were used to achieve the equilibration, while the later 5 × 106 production steps were used to obtain the ensemble averages. 3. Results and Discussion 3.1. Characterization of the MCM-41 Model. The optimized MCM-41 is shown in Figure 1c. The lattice parameters of the MCM-41 model before and after optimization are listed in Table 2. After optimization, the structure shrinks by 10% in each direction and slightly switches to be nonorthorhombic. Nevertheless, the cylindrical pore shape and the hexagonal array are fairly well maintained. Figure 3a shows the XRD pattern of the optimized MCM-41 model. There are three characteristic peaks at 2θ ) 3.04°, 5.28° and 6.08°, respectively. In
TABLE 2: Lattice Parameters of the MCM-41 Model before and after Optimization parameter
a/Å
b/Å
c/Å
R/deg
β/deg
γ/deg
before after
64.18 58.01
149.76 133.88
42.78 39.17
90 90.37
90 89.21
90 90.76
particularly, the peak at 2θ ) 3.04° is significantly larger. These indicate the diffractions on the (1 0 0), (1 1 0) and (2 0 0) surfaces of a hexagonal array,2,3,22 as schematically illustrated in the inset of Figure 3a. The pore radius ra was estimated by measuring the distance between the pore center and the H atoms on the pore surface.10 Figure 3b shows the pore radius in the optimized MCM-41 model has a Gaussian distribution. The average pore radius is 14.38 Å, slightly smaller than the initial radius before optimization (15.0 Å). This pore size in our MCM41 model is very close to 13.48 Å in an experimental sample.9 The smoothness of a pore surface at the atomic level plays a major role in adsorption. Here, a coefficient λ defined as23
λ)
Sc 2nπralz
(1)
was used to depict the smoothness of the pore surface in MCM41. In eq 1, Sc is the Connolly surface24,25 estimated with a probe of diameter 3.3 Å, which is the kinetic diameter of CO2. The number of pore in the model is n and lz is length of the pore along axis. The value of λ depicts the degree of smoothness. A
11298 J. Phys. Chem. C, Vol. 112, No. 30, 2008
Figure 5. Adsorption isotherms of CO2 in the MCM-41 at 303 K from simulation and experiment.9 The BET specific surface area and BJH average pore radius of the experimental MCM-41 sample are 1013.7 m2/g and 13.48 Å, respectively. The Connolly surface area and average pore radius of our MCM-41 model are 1047.7 m2/g and 14.38 Å, respectively.
unit value indicates a perfect smooth cylindrical pore at the atomic level, while a value greater or less than 1 indicates a degree of corrugation on the surface. In our mimetic MCM-41 model, Sc and λ was found to be 1047.7 m2/g and 1.22, which implies the MCM-41 model is not smooth. The estimated surface area is in good agreement with the experimentally measured value of 1013.78 m2/g.9 3.2. Adsorption of Pure Gas. Figure 4 shows the simulated and experimental adsorption isotherms and isosteric heats of CO2 in MCM-41 at low pressures. Experimental data were determined in three samples M41C14, M41C16, and M41C22 with different pore radii of 13.50, 15.05, and 19.45 Å, respectively.26 The three experimental samples give nearly identical loading. Fairly good agreement is observed between simulation and experiment, though the simulated loading is a bit higher than the experimental data at relatively high pressures. The experimental data of isosteric heat are rather scattered and the simulated results match well with the data measured in sample M41C16. The good agreement between simulated and experimental adsorption isotherms and heats validates the rationality of the mimetic MCM-41 model. Note that the isosteric heat is usually used to ascertain the adsorption mechanism as the isosteric heat is quite sensitive to the change of adsorption energy. Within the low pressure range in Figure 4, the isosteric heat decreases with increasing loading. We can expect that it will increase with further increasing pressure due to the corporative attraction between sorbate molecules. However, it will finally drop approaching saturation because the space available for adsorption becomes less. Figure 5 shows the adsorption isotherms of CO2 in MCM41 at high pressures from simulation and experiment. The specific surface area and BJH average pore radius of the experimental MCM-41 sample were measured to be 1013.7 m2/g and 13.48 Å,9 in good accordance with 1047.7 m2/g and 14.38 Å of the mimetic MCM-41 model. Although the simulated results are slightly greater than the experimental data, the agreement between the two is satisfactory. At low pressures (20 and 100 kPa) CO2 molecules are adsorbed at the active sites on/near the MCM-41 surface, such as hydroxyls and defects on the pore surface shown in Figure 6. With increasing pressure at 1500 kPa, the loading increases; CO2 molecules start to occupy the region next to the pore surface
Zhuo et al.
Figure 6. Snapshots of CO2 adsorbed in MCM-41 at 298 K and different pressures (20, 100 and 1500 kPa). Key: Si, yellow; O, red; H, white; C, black.
Figure 7. Density distributions of (a) CO2 and (b) N2 at 400 kPa and 298 K. The pore skeleton is represented by ball and stick. Key: Si, yellow; O, red; H, white.
and annular layers are gradually formed. At still higher pressures (not shown), the center region is expected to be occupied, as observed for gas adsorption in cylindrical carbon nanotubes.27,28 Figure 7 shows the density distributions of CO2 and N2 at 400 kPa. We can see that CO2 and N2 exhibit different adsorption regions in the MCM-41 pores. The density of CO2 at the active sites near the pore surface is greater than at other regions. However, N2 tends to adsorb homogeneously on the pore surface. This is attributed to the favorable interaction between CO2 and sorbent. CO2 is a three-site molecule with a relatively large quadrapole moment. The pore surface of MCM-41 is covered with polar hydroxyls, which can interact more strongly with CO2 than N2. 3.3. Adsorption of Flue Gas. Flue gas is mimicked here by a mixture consisting of 15% CO2 and 85% N2 in mole fraction. The adsorption isotherms and isosteric heats of pure and mixed CO2 and N2 are shown in Figure 8. Apparently, CO2 has a higher loading and isosteric heat than N2. The loadings of CO2 and N2 in the flue gas are lower than in pure gas, particularly for N2 at high pressures. Similar trend also appears in the isosteric heat, which is more distinct for N2. This indicates that CO2 is more preferentially adsorbed than N2 in MCM-41, as observed above. Consequently, CO2 adsorption is not significantly influenced by the presence of N2. The favorite sites on the pore surface are covered mostly by CO2 rather than N2, which leads to the
Adsorption of CO2, N2 and Flue Gas
J. Phys. Chem. C, Vol. 112, No. 30, 2008 11299
Figure 8. (a) Adsorption isotherms and (b) isosteric heats for the adsorption of pure and mixed CO2 and N2 at 298 K from simulation. The lines are drawn to guide the eye.
Figure 9. Selectivity of CO2 over N2 in MCM-41 from simulation as a function of (a) total pressure and (b) temperature. The lines are drawn to guide the eye.
decrease in loading and isosteric heat for N2 in the flue gas upon comparison with pure N2. Separation of CO2 and N2 in the flue gas is described by selectivity S as defined by:
S)
xCO2/xN2 yCO2/yN2
(2)
where xCO2 and xN2 are the mole fractions of CO2 and N2 in adsorbed phase, while yCO2 and yN2 are the mole fractions of CO2 and N2 in the bulk phase. Figure 9a shows the dependence of S on the total pressure of the flue gas at three different temperatures 298, 373 and 573 K. With increasing pressure, the selectivity drops at a low temperature 298 K. This is because the loadings of CO2 and N2 rise as pressure increases. At low loadings, CO2 is more preferentially adsorbed than N2, which contributes to a higher selectivity. The active sites run out rapidly with increasing loading and further adsorbed sorbate molecules have to locate at less favored sites. At these sites, CO2 is not so competitive to N2. As a result, the selectivity of CO2 over N2 declines with the increase of pressure and gradually reaches a constant value of about 7. However, at a higher temperature 373 or 573 K, the selectivity is nearly independent of pressure. The reason is that the active sites cannot strongly capture gas molecules with a larger kinetic energy at a higher temperature. As a consequence, the distinction between CO2 and N2 diminishes and temperature becomes the dominant factor for selectivity. As evidenced in Figure 9b, the selectivity decreases rapidly with increasing temperature. In general, the
selectivity at a given temperature is a weak function of total pressure. At temperatures higher than 400 K, the selectivity tends to approach a plateau and total pressure has no discernible influence on the selectivity. These reveal that temperature plays a more important role than pressure in determining the separation factor. 4. Conclusions We have developed a mimetic model of MCM-41 by carving the cylindrical pores in an amorphous silica matrix and subsequent geometric optimization. The atomic partial charges for the model are not sensitive to the cluster chosen in the DFT calculations. The XRD pattern presents three characteristic peaks at the 2θ ) 3.04°, 5.28° and 6.08°, corresponding to the diffractions on the (1 0 0), (1 1 0) and (2 0 0) surfaces of a hexagonal array. The pore surface is not smooth and coated with hydroxyl groups and defects as active sites for adsorption. The pore radius has a Gaussian-like distributed and an average value of 14.38 Å, close to experimental measurement. Good agreement is found between simulated and experimental isotherms and heats of CO2 adsorption. Due to the stronger Coulombic interaction, CO2 is more preferentially adsorbed on the active sites in the mimetic MCM-41, while N2 tends to adsorb homogeneously on the pore surface. The loading and isosteric heat of CO2 are almost the same in pure gas and mixture, however, they are smaller for N2 in mixture than in pure N2. The selectivity of CO2 over N2 drops rapidly with increasing temperature and reaches a constant at high temper-
11300 J. Phys. Chem. C, Vol. 112, No. 30, 2008 atures. Temperature is a dominant factor to govern the selectivity, though pressure also influences weakly at low temperatures. Acknowledgment. This work was supported by the National Natural Science Foundation of China (Grants 20776045 and 20736002), Program for Changjiang Scholars and Innovative Research Team in University of China (Grant IRT0721), the 111 Project (Grant B08021) of China, and the National University of Singapore (Grants R-279-000-243-123 and R-279000-198-112/133). References and Notes (1) Kresge, C. T.; Leonowicz, M. E.; Roth, W. J.; Vartuli, J. C.; Beck, J. S. Nature 1992, 359, 710. (2) Beck, J. S.; Vartuli, J. C.; Roth, W. J. J. Am. Chem. Soc. 1992, 114, 10834. ´ .; Neimark, A. V.; Schu¨th, (3) Ravikovitch, P. I.; Domhnaill, S. C.O F.; Unger, K. K. Langmuir 1995, 11, 4765. (4) Tjatjopoulos, G. J.; Feke, D. L.; Mann, J. A., Jr J. Phys. Chem. 1988, 92, 4006. (5) Maddox, M. W.; Gubbins, K. E. Int. J. Thermophys. 1994, 15, 1115. (6) Ravikovith, P. I.; Neimark, A. V. Langmuir 2006, 22, 11171. (7) Maddox, M. W.; Olivier, J. P.; Gubbins, K. E. Langmuir 1997, 13, 1737.
Zhuo et al. (8) Koh, C. A.; Montanari, T.; Nooney, R. I.; Tahir, S. F.; Westacott, R. E. Langmuir 1999, 15, 6043. (9) He, Y. F; Seaton, N. A. Langmuir 2003, 19, 10132. (10) Benoit, C.; Hung, F. R.; Roland, J. M.; Pellenq; Siperstein, F. R.; Gubbins, K. E. Langmuir 2006, 22, 194. (11) Jeanette, M. S.; Enrique, I. Chem. Eng. Sci. 2001, 56, 4205. (12) Materials Studio; Accelrys Inc.: San Diego, CA. (13) Perdew, J. P.; Wang, Y. Phys. ReV. B 1992, 45, 13244. (14) Singh, C. U.; Kollman, P. A. J. Comput. Chem. 1984, 5, 129. (15) Bakalarski, G.; Grochowski, P.; Kwiatkowski, J. S.; Lesyng, B.; Leszczynski, J. Chem. Phys. 1996, 204, 301. (16) Bayly, C. I.; Cieplak, P.; Cornell, W. D.; Kollman, P. A. J. Phys. Chem. 1993, 97, 10269. (17) Merz, K. M., Jr J. Comput. Chem. 1992, 13, 749. (18) Mulliken, R. S. J. Chem. Phys. 1955, 23, 1833. (19) Hirshfeld, F. L. Theor. Chim. Acta B 1977, 44, 129. (20) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. J. Phys. Chem. 1990, 94, 8897. (21) Harris, J. G.; Yung, K. H. J. Phys. Chem. 1995, 99, 12021. (22) Shinae, J.; Sang, H. J.; Ryong, R.; et al. J. Am. Chem. Soc. 2000, 122, 10712. (23) Vyas, S.; Dickinson, D. S.; Armstron-Poston, E. Mol. Simul. 2006, 32, 135. (24) Connolly, M. L. Science 1983, 221, 709. (25) Connolly, M. L. J. Appl. Crystallogr. 1983, 16, 548. (26) He, Y. F; Seaton, N. A. Langmuir 2006, 22, 1150. (27) Jiang, J. W.; Sandler, S. I. Phys. ReV. B 2003, 68, 245412. (28) Jiang, J. W.; Sandler, S. I. Langmuir 2004, 20, 10910.
JP803428N
