Prediction of Adsorption and Diffusion Behaviors of CO₂ and CH₄ in All-Silica Zeolites Using Molecular Simulation

Abstract

Computational chemistry is a promising technique for the prediction of material properties. Adsorption and diffusion behaviors in zeolite micropores are important for zeolite membranes. In this study, we investigated novel non-bonding interaction parameters of all-silica zeolites for the prediction of the adsorption and diffusion behaviors by focusing on the Si atom of zeolite frameworks. Our parameters (σ = 0.421 nm, ε = 0.954 kJ mol⁻¹, and q = +1.10 e) were close to theoretically derived values, and the adsorption isotherms of CO₂ and CH₄ on several zeolites could be predicted with high accuracy. Furthermore, the parameters gave the suitable self-diffusivities of CO₂ and CH₄ within MFI-type zeolite micropores through molecular dynamics simulation. Those suggest that our derived parameters are useful for selecting zeolite structure as the membrane material.

1. Introduction

Zeolites are aluminosilicate compounds with regular-sized micropore structures, and they show unique adsorption behaviors due to their micropore structures. Therefore, they are often used as adsorbents and catalysts. It has been possible to form polycrystalline zeolite layers on porous substrates since the 1990s [1,2,3,4,5,6,7], and zeolite membranes have been attracting much attention as an energy-saving separation technology since then. Furthermore, new zeolite structures and their adsorption and diffusion behaviors have been reported. However, it has been difficult to select suitable zeolite structures and compositions for membrane separation.

There are already many studies about the prediction of the adsorption and diffusion behaviors in zeolites using grand canonical Monte Carlo (GCMC) and molecular dynamics (MC) [8,9,10,11,12,13,14]. However, valid parameters for many zeolites have not yet been reported. Talu et al. [9] assumed that the adsorbed molecule in zeolite micropores interacts with the oxygen atoms of zeolites, and the adsorption isotherms of CO₂ and CH₄ on all-silica MFI-type zeolite (silicalite-1) were simulated. Although the adsorption properties for the MFI-type zeolite could be predicted, the parameters were not valid for other types of zeolites.

The parameters for adsorbates and zeolites are necessary to predict these behaviors with high accuracy. For adsorbate molecules, Sun et al. [15,16] and Siepmann et al. [17,18,19,20] investigated the suitable parameters for Lennard-Jones 9-6 and 12-6 potentials, respectively. Their parameters can describe vapor–liquid equilibria. Rahmati et al. [12] tried to predict the adsorption isotherms on all-silica zeolites using parameters reported by Sun et al. However, the simulated amounts of adsorbed CO₂ and CH₄ on MFI and DDR-type zeolites were several times higher than those in experimental results. Although Bai et al. [11] and Vujic et al. [14] reported on the simulation parameters of all-silica zeolites for Lennard Jones 12-6 potential, their parameters could not be applied to many zeolites.

In our study, we have found the suitable parameters for many kinds of all-silica zeolites by focusing on the silicon atom of all-silica zeolites, and the simulated adsorption isotherms were compared to experimental results to check the validation of our parameters in this paper. Furthermore, we discussed the simulated self-diffusivities in MFI-type zeolite micropores.

2. Theory

The interaction between adsorbate and adsorbent atoms is described as the sum of interactions between bonded and non-bonded atoms as follows:

$$E_{\mathrm{t}}=E_{\mathrm{bond}}+E_{\mathrm{non}\text{-}\mathrm{bond}}.$$

(1)

The interaction between bonded atoms is calculated as the sum of bond stretching and angle bending as follows:

$$E_{\mathrm{bond}}=E_{\mathrm{bond}\text{-}\mathrm{stretch}}+E_{\mathrm{angle}\text{-}\mathrm{bend}},$$

(2)

$$E_{\mathrm{bond}\text{-}\mathrm{strech}}=\frac{1}{2}{k}_{\mathrm{b}}{(r-r_0)}^2,$$

(3)

$$E_{\mathrm{angle}\text{-}\mathrm{bend}}=\frac{1}{2}{k}_{\theta }{(\theta -{\theta }_0)}^2,$$

(4)

where k_b and k_θ are the force constants for bond stretching and angle bending, respectively. The interaction between the non-bonded atom pair is calculated as the sum of van der Waals and coulomb interactions as follows:

$$E_{\mathrm{non}\text{-}\mathrm{bond}}=E_{\mathrm{vdW}}+E_{\mathrm{coulomb}},$$

(5)

$$E_{\mathrm{vdW}}=4{\epsilon }_{ij}\left[{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^{12}-{\left(\frac{{\sigma }_{ij}}{r_{ij}}\right)}^6\right],$$

(6)

$$E_{\mathrm{coulomb}}=\frac{1}{4\pi {\epsilon }_0}\cdot \frac{q_i{q}_j}{r_{ij}},$$

(7)

where the depth of interaction ε_ij and zero-interaction distance σ_ij for the pair of different atoms are calculated as follows:

$${\epsilon }_{ij}=\sqrt{{\epsilon }_i{\epsilon }_j},$$

(8)

$${\sigma }_{ij}=\frac{1}{2}({\sigma }_i+{\sigma }_j).$$

(9)

When two atoms are in the same structure and separated by three covalent bonds (so-called 1–4 interaction), the interaction is treated by non-bonded interaction with a scaling factor of 0.5. Non-bonded interaction is ignored for the directly bonded atoms (1–2 interaction) and two atoms separated by two bonds (1–3 interaction), since they are included in the bond-stretching and angle-bending interactions.

3. Materials and Methods

3.1. Synthesis of All-Silica Zeolites

3.1.1. BEA-Type Zeolite

All-silica BEA-type zeolite was synthesized via a hydrothermal process [13]. Tetraethyl orthosilicate (Tokyo Chemicals Industry, Tokyo, Japan) was added to tetraethylammonium hydroxide solution (TEAOH, 40 wt.%, Aldrich, St Louis, MO, USA), and the mixture was stirred overnight at 313 K for hydrolysis. After the by-produced ethanol was vaporized at 373 K, water and hydrogen fluoride solution (46%, Wako Pure Chemicals Industry, Tokyo, Japan) was added. The mixture was mixed using a spatula to obtain a homogeneous white paste. The molar composition of the paste was 1 SiO₂: 0.5 TEAOH: 0.5 HF: 7 H₂O. The paste was transferred into an autoclave, and a hydrothermal reaction was carried out at 423 K for 14 days under a rotation of 60 rpm. After the reaction, solids were recovered by filtration, washed with distilled water, dried at 383 K overnight, and calcined in air at 773 K for 10 h to remove the structure-directing agent. Thus, the all-silica BEA-type zeolite particles were obtained.

3.1.2. CDO-Type Zeolite

CDO-type zeolite particles were prepared via the structure conversion of layered silicate PLS-1 [21]. Potassium hydroxide, tetramethylammonium hydroxide solution (TMAOH, Aldrich), and 1,4-dioxane (Dx, Wako Pure Chemical industry) were dissolved into distilled water. Fumed silica particles (Cabot Corp., Cab-o-sil M5) were added to the solution, and the gel was transferred into the autoclave after mixing at room temperature for 6 h. The molar composition of the gel was 1 SiO₂: 0.007K₂O: 0.22 TMAOH: 3.4 Dx: 16 H₂O. A hydrothermal reaction was carried out at 423 K for 14 days. Solids were recovered by centrifugal separation, washed with distilled water, and dried in air at 383 K overnight to obtain layered silicate PLS-1 particles. Then, the PLS-1 particles were ion-exchanged in a 0.1 mol/L-HCl solution at 333 K for 3 h. The particles were recovered by filtration, washed with distilled water, and dried in air at 383 K overnight. Finally, the PLS-1 particles were calcined in air at 873 K for 2 h to convert the structure from PLS-1 to CDO-type zeolite.

3.1.3. CHA-Type Zeolite

All-silica CHA-type zeolite was synthesized hydrothermally [22]. Tetraethyl orthosilicate and N,N,N-trimethyl-1-ammonium hydroxide solution (TMAAOH, Sachem, Osaka, Japan) were mixed and stirred for 2 h at 373 K for hydrolysis. A hydrofluoric acid solution (Wako Pure Chemicals Industry) was added to the solution, after the solution was concentrated via evaporation. The molar composition of the mixture was 1SiO₂: 0.5TMAAOH: 0.5HF: 3H₂O. The mixture was transferred to the autoclave, and a hydrothermal reaction was carried out at 423 K for 24 h. Solids were recovered by filtration and washed with water until the pH of the washing water became neutral. After drying at 383 K overnight, the particles were calcined in air at 823 K for 10 h to burn out the structure-directing agent.

3.1.4. DDR-Type Zeolite

All-silica DDR-type zeolite was synthesized via a hydrothermal reaction [23]. Adamantamine (ADA, Tokyo Chemicals Industry) was dissolved into ethylenediamine (EDA, Tokyo Chemicals Industry), and distilled water was added to the solution. The mixture was stirred for 1 h at 373 K after stirring at room temperature for 1 h, and then, the solution was cooled in an ice bath. Tetramethyl orthosilicate (Tokyo Chemicals Industry) was added to the solution, and the mixture was stirred at 373 K overnight to obtain a clear solution. The molar composition of the solution was 1SiO₂: 4.1EDA: 0.47ADA: 112H₂O. The solution was added to the autoclave, and a hydrothermal process was carried out at 423 K for 25 days. Solids were recovered by filtration, washed with water, and dried in air at 383 K overnight. Finally, the particles were calcined in air at 973 K for 6 h to remove the structure-directing agent.

3.1.5. MFI-Type Zeolite

All-silica MFI-type zeolite was synthesized via a hydrothermal process [24]. Triethylamine (TEA, Tokyo Chemicals Industry) and tetraethyl orthosilicate (TEOS, Tokyo Chemicals Industry) were added to distilled water, and the mixture was stirred for 3 h at room temperature for hydrolysis. Tetrapropylammonium bromide (TPABr, Tokyo Chemicals Industry) was dissolved into the solution, and the solution was stirred for 4 h at 363 K. The molar composition of the solution was 1SiO₂: 2TEA: 1.6TPABr: 200H₂O. Then, the solution was poured into an autoclave, and a hydrothermal reaction was carried out at 403 K for 24 h. After cooling to room temperature, solids were recovered by filtration and washed with distilled water. The washing and filtration were repeated until the pH of the washing water became neutral, and the particles were dried in air at 383 K overnight. Finally, the particles were calcined in air at 823 K for 10 h to remove a structure-directing agent, in order to obtain all-silica MFI-type zeolite particles.

3.1.6. STT-Type Zeolite

All-silica STT-type zeolite was synthesized via the same procedure as for the CHA-type zeolite [22]. The composition of a synthesis gel was 1SiO₂: 0.5TMAAOH: 0.5HF: 10H₂O. A hydrothermal reaction was performed at 423 K for 11 days. After solids were filtered, washed with distilled water, and dried in air at 383 K overnight, the particles were calcined in air at 823 K for 10 h for the removal of the structure-directing agent.

3.2. Characterization

The crystal structures of zeolites were identified via X-ray diffraction (Smart Lab., Rigaku, Tokyo, Japan), and the morphology was observed using a scanning electron microscope (TM-1000, Hitachi High Technology, Tokyo, Japan).

3.3. Adsorption Experiment

The amount of adsorbed CO₂ and CH₄ on zeolites was determined using an adsorption unit equipped with a constant volume cell (Belsorp Max, Japan BEL, Osaka, Japan) at 298 K and 343 K. Approximately 1 g of zeolite particles was added to the cell and evacuated at 673 K overnight to remove adsorbed water. The pretreated particles were used for the adsorption experiment without exposure to atmosphere. The pretreatment was carried out before each measurement. Either CO₂ or CH₄ were added to the cell, and the pressure was decreased because of the adsorption of these gases on the zeolite. When the pressure reached a constant level, the amount of adsorption was determined by the reduction of pressure. The difference in the amount of adsorption for each measurement was less than 1%.

3.4. Simulation

The adsorption isotherms of CO₂ and CH₄ on the all-silica BEA, CDO, CHA, DDR, MFI, and STT-type zeolites were simulated via a grand canonical Monte Carlo (GCMC) technique using a software program (Biovia, Materials Studio 2017 Sorption). For GCMC simulation, fugacity was given to the canonical ensemble, and the number and location of molecules with the lowest potential energy were calculated probabilistically. The cutoff distance of the van der Waals interaction was 1.25 nm, and the Ewald summation method was used for the integration of the coulomb interaction. The total number of Monte Carlo cycles was 10⁶, and the average of the final 10⁵ steps was used as the simulation result. The fugacity was assumed to be equal to the pressure in this study, since the difference between fugacity and pressure was less than 5% below 1 MPa.

Table 1. Crystal structures of all-silica zeolites.

Zeolite	Space Group	a (nm)	b (nm)	c (nm)	α (°)	β (°)	γ (°)	Composition
BEA	P4122 (No.91)	1.266	1.266	2.641	90.00	90.00	90.00	Si64O128
CDO	Pnma (No.62)	1.836	1.378	7.367	90.00	90.00	90.00	Si36O72
CHA	R3m (No.166)	0.942	0.942	0.942	94.07	94.07	94.07	Si12O24
DDR	R3m (No.166)	1.386	1.386	4.089	90.00	90.00	120.00	Si120O240
MFI	Pnma (No.62)	2.002	1.990	1.338	90.00	90.00	90.00	Si96O192
STT	P21/N (No.14)	1.309	2.167	1.373	90.00	102.58	90.00	Si64O128

shows the crystal structure data of all-silica zeolites used in our calculation. The crystal structures were taken from the IZA zeolite database [21]. The bond interactions in the zeolite framework were ignored in this study.

Figure 1 shows the atomistic models of CO₂, CH₄ and MFI-type zeolite. The model of the CO₂ molecule reported by Harris et al. [25] was used in this study. This model can describe the gas–liquid coexistence curve including the critical point region. The carbon atom was connected to two oxygen atoms by chemical bonds which are 0.1149 nm long, and bond stretching was ignored (k_b = 0). The original angle of O=C=O was 180°, and the force constant k_θ = 1236 kJ mol⁻¹ rad⁻². For CH₄, the model reported by Siepman et al. [17] was used. The carbon atom was connected to four hydrogen atoms with the bond length of 0.11 nm, and each H-C-H angle was 109.5°. Although the bond stretching and angle bending are ignored in this model (k_b = k_θ = 0), the gas–liquid coexistence curve can be expressed.

Table 2. Non-bonding interaction parameters of CH₄ and CO₂.

Molecule	Element	σ (nm)	ε/k (K)	q (e)	Ref.
CH4	C	0.3730	148.0	0	[17]
	H	---	---	0
CO2	C	0.2757	28.1	0.6512	[25]
	O	0.3033	80.5	−0.3256

lists the non-bonding interaction parameters for CO₂ and CH₄.

The self-diffusivities of CH₄ and CO₂ in the MFI-type zeolite channel were also simulated via a molecular dynamic technique (Biovia, Materials Studio 2017 Forcite Plus). Several molecules of CH₄ and CO₂ were adsorbed by GCMC, and the molecular dynamic simulation was carried with the time step of 1 fs. The total simulation time was 1 ns, and the mean square displacements of every 1 ps were plotted against the simulation time. The self-diffusivity was calculated by the Einstein equation as the slope of the mean square displacement with respect to time.

4. Results and Discussion

4.1. Characterization

Figure 2 shows the SEM images of the all-silica zeolite particles synthesized in this study. The shape of the BEA-type zeolite was a truncated bipyramid with the size of 2–10 μm, and the CHA-type zeolite was cubic crystals with the size of ca. 2 μm. The CDO and STT-type zeolite were plate crystallites, and the sizes were circa 2 × 2 × 0.2 μm³ (CDO) and 5 × 10 × 0.1 μm³ (STT). Figure 3 shows the XRD patterns of the zeolite particles synthesized in this study. The XRD patterns were identical to those of corresponding zeolites, and no peaks of impurities were observed. The BEA-type zeolite has two kinds of polymorphs, A and B, and the containing ratio is estimated by the peak patterns at 2θ = 6−8° [26]. The BEA-type zeolite particles synthesized in this study were estimated to be A/B = 1/1. These results suggest that the zeolite crystallites with no impurities can be synthesized.

4.2. Adsorption Isotherm

Figure 4 and Figure 5 compare the simulated adsorption isotherms of CO₂ and CH₄ on the MFI and DDR-type zeolites to experimental results, respectively. Sun et al. [27] and Abdul-Rehman et al. [28] reported the adsorption isotherms of CO₂ and CH₄ on MFI-type zeolites at high pressures, and those for the DDR-type zeolite were reported by Himeno et al. [29]. Our experimental results agreed well with these reported data. The parameters were then optimized, so that the simulation results matched the experimental data. As a result, the parameters listed in

Table 3. Non-bonding interaction parameters of all-silica zeolites.

Elements	σ (nm)	ε (kJ mol−1)	q (e)
Si	0.421	0.954	1.10
O	---	---	−0.55

were obtained. The simulation results were more consistent with the experimental data by giving the van der Waals parameters to only the Si atom of zeolites, as compared to cases where the parameters were given to both O and Si atoms. The distance between the Si atoms at zero interaction potential σ was the same as the van der Waals diameter of Si (0.42 nm [30]), and the potential depth of ε was close to the interaction potential calculated by the Kirkwood–Muller equation (1.0 kJ mol⁻¹) [31]. Furthermore, assuming that the bond length of Si-O is 0.15 nm, the partial charge of the Si atom could be calculated to be 1.02 e from the dipole moment of the Si-O bond (1.23 × 10⁻²⁹ C m [32]). These suggest that our parameters are theoretically reasonable.

Figure 6 compares the simulated adsorption isotherms of CH₄ and CO₂ on BEA, CDO, CHA, and STT-type zeolites to the experimental data. The simulated values were almost the same as the experimental results, and the differences were below 0.09 mol kg⁻¹ for CH₄, and 0.30 mol kg⁻¹ for CO₂. The same amounts of adsorbed CO₂ and CH₄ at different temperatures gave the same heats of adsorption between experiments and simulation. As described below, the isotherms on the BEA-type zeolite had the largest differences in previous studies. The heat of adsorption on the BEA-type zeolite was 18.0 kJ mol⁻¹ for CO₂ and 15.3 kJ mol⁻¹ for CH₄ experimentally, and the simulated values were 18.5 and 14.6 kJ mol⁻¹ for CO₂ and CH4, respectively. Pham et al. [11], Bai et al. [13], and Vujic et al. [14] reported the non-bonding parameters for zeolites, as shown in

Table 4. Non-bonding interaction parameters of all-silica zeolites, averages of differences between the amounts of adsorbed obtained by the experiment and simulation, and their standard deviations.

Element	Parameter	Bai et al. [11]	Pham et al. [13]	Vujic et al. [14]	This Work
Si	σ (nm)	0.231	---	0.297	0.421
	ε (kJ mol−1)	0.185	---	0.266	0.954
	q (e)	1.5	2.0	1.413	1.10
O	σ (nm)	0.3304	0.2806	0.3011	---
	ε (kJ mol−1)	0.442	0.744	0.432	---
	q (e)	−0.75	−1.0	−0.7065	−0.55
Average of difference (mol kg−1)		0.09	0.13	0.09	0.06
Standard deviation (mol kg−1)		0.09	0.18	0.10	0.06

, and they also simulated the adsorption isotherms on several kinds of zeolites. The adsorption isotherms of CO₂ and CH₄ on the BEA, CDO, CHA, DDR, MFI, and STT-type zeolites were also simulated using their parameters in this study. Their parameters could thus predict the adsorption isotherms on some zeolites. However, the accuracies for the adsorption on the BEA and CDO-type zeolites were significantly low. Figure 7 compares the simulated adsorption amount using their parameters for the experimental data, and Table 4 summarizes the average difference between the experimental and simulated values and its standard deviation. The average of difference and standard derivation were calculated as follows:

$$\mathrm{Average}\mathrm{of}\mathrm{difference}(\%)=\frac{100}{N}{\displaystyle \sum _{}^N\frac{|q_{\exp }-q_{\mathrm{sim}}|}{q_{\exp }}},$$

(10)

$$\mathrm{Standard}\mathrm{derivation}(\%)=\frac{100}{N}{\displaystyle \sum _{}^N|\frac{|q_{\exp }-q_{\mathrm{sim}}|}{q_{\exp }}-\frac{1}{N}{\displaystyle \sum _{}^N\frac{|q_{\exp }-q_{\mathrm{sim}}|}{q_{\exp }}}|},$$

(11)

where q_exp and q_sim are the amount of adsorbed obtained by the experiment and simulation, respectively. The average of the difference and standard deviation of our parameters were 0.06 mol kg⁻¹ and 0.06 mol kg⁻¹, respectively. These were lower than those for the other parameters. Therefore, we conclude that our parameters can predict the amounts of adsorbed CO₂ and CH₄ on all-silica zeolites with high accuracy.

4.3. Self-Diffusivity Adsorption Isotherm

Next, the self-diffusivities were simulated via molecular dynamics to confirm the effectiveness of our parameters. The MFI-type zeolite was the selected zeolite, since many experimental data of self-diffusivities are reported. Figure 8 shows the influences of the amount of adsorbed CO₂ and CH₄ on the self-diffusivities of CO₂ and CH₄ in the MFI-type zeolite at 298 K. When three CH₄ molecules were adsorbed to the unit cell (0.52 mol kg⁻¹), the diffusivity was estimated to be 1.9 × 10⁻⁸ m² s⁻¹. The diffusivity decreased with an increase in the number of adsorbed molecules, and it was 3.8 × 10⁻⁹ m² s⁻¹ per twenty-one CH₄ molecules adsorbed (3.64 mol kg⁻¹). The self-diffusivity of CO₂ showed the same trend. Caro et al. [33] and Karger et al. [34] used the PFR-NMR method to measure the self-diffusivities of CH₄ and CO₂ in the all-silica, MFI-type zeolite. Their data are also plotted in Figure 8. Our simulation results were in good agreement with their experimental results.

Finally, the effect of temperatures on the self-diffusivities of CH₄ and CO₂ in MFI-type zeolite was simulated. Figure 9 shows the temperature dependencies of the self-diffusivities of CO₂ and CH₄ in the MFI-type zeolite. For the comparison with the PFG-NMR experimental data of Caro et al. [33] and Karger et al. [34], nineteen CH₄ or fifteen CO₂ molecules were adsorbed, and the molecular dynamic simulations were carried out. These number of molecules correspond to the amount of adsorbed 3.29 or 2.60 mol kg⁻¹, respectively. Millot et al. [35] and Talu et al. [36] determined the diffusivities from the permeability of CH₄ through a polycrystalline membrane and single crystal membrane, respectively. The self-diffusivities of CH₄ simulated using our parameters were similar to the experimental results, as shown in Figure 9a. Further, Vujic et al. [14] simulated the CO₂ diffusivity in the MFI-type zeolite by their non-binding parameters, as listed in Table 4. However, their simulation results are 2–3 times higher than the experimental results of Karger [34], as shown in Figure 9b. By contrast, the CO₂ diffusivities using our parameters were close to the experimental results. This means that our parameters can be used for the prediction of the self-diffusivities in MFI-type zeolite with a higher accuracy than previously reported parameters.

5. Conclusions

In this study, all-silica BEA, CDO, CHA, DDR, MFI, and STT-type zeolite particles were synthesized, and the amount of adsorbed CO₂ and CH₄ was determined at 298 K and 343 K. When the non-bonding parameters (σ_Si = 0.421 nm, ε_Si = 0.954 kJ mol⁻¹, and q_Si = 1.10) were used, the simulated adsorption isotherms showed better agreement with the experimental results than those calculated using previously reported parameters. Those parameters also provided self-diffusivities close to experimental data, compared to the other parameters. These results suggest that our derived parameters are useful for selecting zeolite structure as the membrane material.