Article pubs.acs.org/JPCC
Interplay of Montmorillonite−H2O−scCO2 System between Mechanical Behavior and Adsorption: Molecular Dynamics Weina Zhang,*,† Haixiang Hu,† Xiaochun Li,† and Zhiming Fang† †
Institute of Rock and Soil Mechanics, The Chinese Academy of Science Xiaohongshan, Wuchang, Wuhan 430071, P. R. China ABSTRACT: Montmorillonite (MMT) constitutes most of fine-grained sedimentary rock, such as shale, mudstone, and siltstone. It is also applied in underground storage of carbon dioxide, petroleum drilling engineering, and material engineering. Although adsorption and mechanical behavior in MMT−H2O−CO2 is widely studied, little is known about interplay between them. In this work, we have investigated the interplay between mechanical behavior and adsorption in MMT with different proportion of H2O and CO2 content under super critical condition by molecular dynamic (MD) method. With respect to the calculated results, adsorption of H2O and CO2 can diminish stiffness of MMT. Meanwhile, diffusion of H2O and CO2 weakens because of stronger constrains deriving from geometry structure and hydrogen bond as the system undergoing compression test. Moreover, motion of H2O is restrained by CO2 molecules as self-diffusion coefficient decrease with the increase of CO2. Hydrated MMT is more sensitive to deformation vertical to clay sheet which is confirmed by comparing changing tend of potential energy and stress. Process of adsorption can yield the stiffness of MMT structures while compression loading can slow down diffusion motions of adsobates. This impact we have revealed here is meaningful to understand mechanical properties of swelling MMT.
1. INTRODUCTION Layered clay minerals compose a large proportion of clay stone and soil and acts an important role in many fields such as underground storage of carbon dioxide,1 petroleum drilling engineering2,3 and material engineering,4,5 etc. As one of the most typical clay minerals swelling upon exposure to water,6,7 Wyoming and Texas−originated montmorillonite (MMT) samples are popular in the literature because of their abundance in nature and broad characterization. Therefore, it is meaningful to understand mechanical properties of MMT which helps us access changes caused by active chemicals, namely CO2 in underground. MMT belongs to phyllosilicate minerals containing variable amounts of water trapped in the interlayer. The interlayer is space between two tetrahedral silicate layers bonded to octahedral hydroxide layers. Theses layers are normally negative charged due to substitution of some metal ions such as Mg2+ Al3+. The whole clay phase keeps electrical neutrality by interlayer cations interlayer cations (e.g., Na+, K+, Ca2+) which are easily hydrated. The extent of clay expansion depends on the relative humidity and relate to the kind of the interlayer cations.8−10 In addition, MMT interacting with CO2 can also result in changes of basal spacing11−13 which has been verified by experimental method.11,13−17 Most of the recent molecular simulation research on the MMT−H2O−CO2 system mainly focuses on the swelling mechanism18−20 but seldom concerning about the mechanical properties. It is difficult to measure mechanical properties of MMT with wanted water or CO2 because of their nanometer scale and anisotropic behavior. However, MD can easily control proportion of water and CO2 and be able to reproduce simple experimental tests such as simple shear test, uniaxial strain test, simple tension test and so on. Recently, MD has been © XXXX American Chemical Society
Figure 1. Simulation cell of the MMT−H2O−CO2 system.
extensively applied for understanding the mechanical properties of some layered silicates including dry MMT. For example, Received: May 21, 2015 Revised: August 31, 2015
A
DOI: 10.1021/acs.jpcc.5b04873 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C
sensitive to water interlayer. Other studies of MMT−H2O, such as effect of temperature and increasing water content on elastic tensors, were also commonly investigated.21,23,25−29 To our best knowledge, however, mechanical properties of MMT with a mixture of CO2 and H2O are rarely involved. Accordingly, our work focused on the relation between varying content of CO2 and H2O and mechanical properties of MMT under supercritical condition of CO2. The paper is organized as follows: we described model construction and introduced simulation methods and details in section 2. In section 3, we organized results of uniaxial compression test. In this part, discussion of the stress−strain relationship and energy and structural parameters as a function of proportion of CO2 were also presented. In section 4, we concluded some main results based on this investigation.
Table 1. Relations between Number of CO2/H2O per Unit Cell and Density of System number of CO2 per unit cell
number of H2O per unit cell
density (g/cm3)
0 0 0.5 1 0 0.4
0 4 3.5 3 8 7.6
2.8016 2.3569 2.3868 2.4040 2.1723 2.2055
Table 2. Summary of Calculated and Experimental Structural Parameters of Models m(CO2)−n(H2O) 0−0 0−4 0.5−3.5 1−3 0−8 0.4−7.6
d001 (Å)a
a (Å)
b (Å)
V (Å3)
9.35 9.839 12.21 12.1,39 12.640 12.25 12.35 14.52 14.47
20.81
17.94
6981.98
20.76
17.97
9111.86
20.76 20.77 20.73 20.69
17.97 17.97 17.88 17.90
9142.46 9220.57 10767.4 10722.8
2. MODELS AND COMPUTATIONAL DETAILS The clay model used in the present work was atomic structure of sodium montmorillonite (MMT) proposed by Skipper et al.30 The chemical formula of the unit cell is expressed as Na0.75(Si7.75Al0.25)(Al3.50Mg0.50)O20(OH)4·nH2O with negative net charge 0.75e. The octahedral layer contributes the 67% of total charge and the tetrahedral layer contributed the remaining charge. This distribution ratio was the same as previous work.31,32 The simulated supercell consists of 4 × 2 × 2 cells as shown in Figure 1. We constructed two hydrated state MMT, namely monolayer and bilayer hydrated MMT. According to the model of Carrier et al.25 and Mazo et al.,24 there are 4 and 8 molecules per unit monolayer cell and bilayer cell. So there are 64 and 128 molecules in monolayer MMT and bilayer MMT, respectively. The clay layer is parallel to the x−y plane and perpendicular to the z-axis. The substitution sites in clay sheet are randomly chosen consistent with the Loewenstein’s rule33 to ensure that the nearest sites were not substituted at one time. Adsorbates CO2/H2O locating in interlayers are also randomly distributed according to the
a
Experimental data from refs 39 and 40, as indicated. The rest are calculated as described in the text.
Manevitch et al.21 obtained the critical stress of single lamella MMT for nonlinear compression. Teich-McGoldrick et al. detected failure mechanism of muscovite under variable pressure and temperature.22 Moreover, Zartman23 summarized failure mechanism of silicates by comparing effect of cation exchange capacity on linear elastic properties of three different minerals, namely pyrophyllite, montmorillonite, and mica. As for MMT− H2O systems, simulation of thermomechanical properties of hydrated MMT with water was presented by Mazo et al.24 Their results showed that Young’s moduli and shear moduli are
Table 3. Summary of Elastic Stiffness Tensor of MMT Containing One Layer Water Carrier et al.25
this work C11 C12 C13 C22 C23 C33 C44 C55 C66
Mazo et al.24
monolayer
bilayer
monolayer
bilayer
monolayer
bilayer
UPV42
222.35 96.97 12.35 263.18 10.34 25.62 10.47 11.59 58.66
204.73 83.91 10.92 221.52 10.38 24.95 9.76 10.54 45.48
217.0 101.1 3.9 208.7 3.9 12.2 1.7 1.9 61.5
186.8 86.5 3.4 177.2 3.6 7.8 1.6 0.3 52.5
231 105 17 229 15 80 21 22 68
189 84 6 191 5 62 2 4 56
13.8−46.1 6.9−17.8 6.1−24.3 N.A. N.A. 11.5−30.4 2.9−8.9 N.A N.A.
Figure 2. Out-of-plane Cij as a function of fraction of CO2 and H2O composition. B
DOI: 10.1021/acs.jpcc.5b04873 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C
Table 4. Self-Diffusion Coefficient of Water in MMT without CO2 (Unit in 10−10 m2/s) state
this work
expt
other simulations results
monolayer bilayer
2.0 5.46
2.043 5.1046
2.22;44 1.5, 1.845 7.37;44 1046
Figure 3. In-plane Cij as a function of fraction of CO2 and H2O composition.
Monte Carlo method. We used Metropolis MC method to determine positions and orientations of the adsorbate molecules. The MC simulation was performed for 200000 steps. The first half was used to equilibrate the system and the rest steps were used for sampled analysis. Then we chose the lowest energy frame as our initial frame for the MD calculation. In order to investigate the effect of sorption on MMT properties, each hydrated model has a fixed number of CO2 molecules per unit cell. In the monolayer MMT, we set the number of CO2 molecules equaling to 0 and 1 the same as Myshakin’s et al.18 models. Additionally, we added a case with 0.5 CO2 per unit cell in the monolayer MMT. In the bilayer hydrated MMT models, we made the number of CO2 equal to 0 and 0.4 which is close to 0.45 in the Botan’s et al.19 bilayer model. For the comparison, we constructed a dry MMT model adsorbing none CO2 or H2O. Density of the four models is listed in Table 1. The apparent density of hydrated systems increases with the number of CO2 molecules per unit cell. The clay force field (CLAYFF) developed by Cygan et al.34,35 was used to compute the energies and forces between atoms. The total energy of the system is given by ⎡⎛ ⎞12 ⎛ ⎞6 ⎤ σij σij Eij = + 4εij⎢⎢⎜⎜ ⎟⎟ − ⎜⎜ ⎟⎟ ⎥⎥ + kb(R i − R 0)2 4πε0rij r ⎝ rij ⎠ ⎦ ⎣⎝ ij ⎠ qiqj
+ kθ(θi − θ0)2
Figure 5. MSD vs time in MMT with 0 CO2−4 H2O. (1)
equilibrated state. The corresponding force field parameters are directly obtained from the literature.34,35 Relaxation of MMT with different proportion of adsorbates was performed in two stages. First, each system was relaxed using the NPT (constant pressure P, number of particles N, and temperature T) ensemble at supercritical CO2 condition (T = 304.45 K and P = 7.38 MPa). The system was equilibrated and sampled separately for 500 ps. The average lattice parameters
The first two parts describe the Coulombic and van der Waals terms. qi and qj are the charges of particles i and j separated by the distance of rij. ε0 is vacuum permittivity of Coulombic potential. εij and σij are Lennard-Jones (L-J) energetic and distance parameters. The last two parts evaluate the bonded terms. kb and qθ are the force constants of bond and angle, respectively. Ri and θi are the bond length and bond angle. R0 and θ0 are the bond length and bond angle at the
Figure 4. Uniaxial compression test. (a) Monolayer versus dry MMT. (b) Bilayer versus dry MMT. C
DOI: 10.1021/acs.jpcc.5b04873 J. Phys. Chem. C XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry C
Figure 6. MSD vs time of MMT with 0.5 CO2−3.5 H2O per unit cell. (a) H2O; (b) CO2.
were obtained from the last 500 ps and given to the final configuration. This configuration was chosen as the first state of straining the system when strain equals to 0. Then, each system was deformed at a constant straining rate. Each straining cycle was conducted in NVT ensemble (constant number of particles N, cell volume V and temperature T = 304.45 K). This NVT simulation was run for 1000 ps. Configurations in the last 500 ps were considered to achieve balance and used to calculate the mechanical properties. All configurations were sampled every 1000 steps. The time step of the two relaxation stages was set to 0.5 fs. Temperature was controlled by a Nóse−Hoover−Langevin (NHL) thermostat36 with a parameter of 0.5 ps and pressure was controlled by a Parrinello−Rahman37 barostat38 with a parameter of 1 ps. Periodic boundary conditions were used to avoid interface effects. Coulomb term was determined by Ewald summation method in an accuracy of 0.0001 kcal/mol. The cutoff distance of 12.0 Å was applied for van der Waals term. The calculation and visualization was realized by the Materials Studio Forcite and Visualizer program. Figure 1 presents one of simulation models used in this study. Variable basal space d001 includes one MMT sheet and the space between two sheets. Some of our calculated structural parameters are close to available experimental data.39,40 Other parameters are unable to compare with experimental data which is not aviliable in reference. We attempted to perform uniaxial compression test at several strain rates: 107/s and 0.2 × 108/s. Considering the accuracy and simulation time, we chose 107/s as the final strain rate. Structural data are given in Table 2. In the atomistic calculation, the internal stress tensors were obtained by the virial expression
where index i runs over all particles from 1 to N; mi, vi, and fij denote the mass, velocity, and force acting on particle i, respectively. In the present work, we used a constant strain method41 to calculate the stiffness constants.40 The constant strain method obtains the tensors via minimizing the energy of the system and deforming the cell parameters in 12 directions in order. It is possible to estimate the stiffness constants from eq 3:
Cijkl =
1 T T [(∑ mi(vv i i )) + (∑ rijfij )] V0 i = 1 i
