Molecular Dynamics Simulation of Nanoscale Elastic Properties of Hydrated Na-, Cs-, and Ca-Montmorillonite

Abstract

The knowledge of nanoscale mechanical properties of montmorillonite (MMT) with various compensation cations upon hydration is essential for many environmental engineering-related applications. This paper uses a Molecular Dynamics (MD) method to simulate nanoscale elastic properties of hydrated Na-, Cs-, and Ca-MMT with unconstrained system atoms. The variation of basal spacing of MMT shows step characteristics in the initial crystalline swelling stage followed by an approximately linear change in the subsequent osmotic swelling stage as the increasing of interlayer water content. The water content of MMT in the thermodynamic stable-state conditions during hydration is determined by comparing the immersion energy and hydration energy. Under this stable hydration state, the nanoscale elastic properties are further simulated by the constant strain method. Since the non-bonding strength between MMT lamellae is much lower than the boning strength within the mineral structure, the in-plane and out-of-plane strength of MMT has strong anisotropy. Simulated results including the stiffness tensor and linear elastic constants based on the assumption of orthotropic symmetry are all in good agreement with results from the literature. Furthermore, the out-of-plane stiffness tensor components of C₃₃, C₄₄, and C₅₅ all fluctuate with the increase of interlayer water content, which is related to the formation of interlayer H-bonds and atom-free volume ratio. The in-plane stiffness tensor components C₁₁, C₂₂, and C₁₂ decrease nonlinearly with the increase of water content, and these components are mainly controlled by the bonding strength of mineral atoms and the geometry of the hydrated MMT system. Young’s modulus in all three directions exhibits a nonlinear decrease with increasing water content.

1. Introduction

Montmorillonite (MMT) consists of well-defined layers separated by interlayer spaces and has great water absorption potential [1,2,3,4]. Thus, it is often used as clay barriers or adsorption material in environmental geotechnical engineering, for example, the backfill in the nuclear waste disposal and the cushion in the landfill [5,6,7,8,9,10]. It is very challenging to directly measure the fundamental mechanical properties of MMT [11,12]. Atomic force acoustic microscopy is used to determine the elastic properties of clay mineral aggregates by measuring adhesion forces [13,14,15]. Nanoindentation and ultrasonic pulse velocity technology (UPV) are also used to measure the elastic stiffness constants of clays and shale samples [16,17,18]. However, these techniques are powerless to measure the anisotropic mechanical properties of MMT, due to the complications in sample preparation and even the testing process.

In contrast to difficulty in experimental methods, molecular dynamics (MD) simulation is a supplemental way to model and understand the properties of hydrated clay minerals on an atomic scale [1,19,20,21]. With the development of polymer–clay nanocomposites in recent years, several studies on the nanomechanical properties of MMT minerals or hydrated MMT systems have emerged. The authors in [18] reported for the first time the stress–strain response of rock-forming minerals by MD simulation. The authors in [22] performed MD simulations to study the flexibility and the mechanical behavior of a single clay layer in a completely exfoliated state and found that the critical stress of a clay layer bending fracture is 0.8 MPa in both in-plane directions. In [23], the authors presented results for the elastic properties of a single lamella of MMT by MD simulation and further pointed out that the elastic constants characteristic of a thin plate is closely related to the thickness of the nanoplate. Based on the analogous minerals and modulus-density relations, ref. [24] deduced a convergence of opinion in the range 178–265 GPa for the elastic moduli of smectite clay platelets. The authors in [25] systematically simulated the structure and the thermomechanical properties of an isolated clay nanoplate and MMT crystals intercalated by water or polyethylene oxide. In [26], the authors clarified the linear elastic properties including tensile moduli, shear moduli, and potential failure mechanisms as a function of cation density and stress for the minerals pyrophyllite, MMT, and mica. In [27], the authors examined the mechanism of bending, the stored energy, and the failure of several clay minerals. They revealed that molecular contributions to the bending energy include bond stretching and bending of bond angles in the mineral as well as rearrangements of alkali ions on the surface of the layers. The authors in [28] further simulated the elastic and structural properties of muscovite as a function of temperature, pressure, and strain. The results demonstrate that the elastic properties of muscovite depend on both temperature and pressure. MD with the CLAYFF force field was used in [29] to simulate isothermal isobaric water adsorption of interlayer MMT, and nanoscale elastic properties of the clay-interlayer water system are calculated from the potential energy of the model system. Similar simulations have also been carried out on MMT by implementing the elastic bath method [30]. The mechanical properties of Na-MMT were investigated in [31] from the view of water content by MD, and the hysteresis phenomena of elastic constant during the swelling and shrinking was found. In addition, researchers found that the mechanical behavior of MMT crystal exhibits a clear dissymmetry between compression and tension and an important dependency on mean stress [32].

The above simulation studies are generally based on the strain method, that is, the stress is calculated by changing the size of the simulation box. The stress method and the large-scale fluctuation method are also usually adopted. For example, using steered molecular dynamics (SMD) technology, the mechanical response of the interlayer of hydrated MMT was evaluated by [33]. Using grid-computing technology, ref. [34] simulated MMT clay large-scale systems containing up to approximately ten million atoms; large-scale systems exhibit emergent behavior with increasing size, and the material mechanical properties were calculated based on thermal bending fluctuations. From the surveyed status of nanoscale mechanical properties of clay minerals, it can be found that the stress–strain response of minerals can be revealed by different methods (constant strain, SMD, or the large-scale fluctuation method) based on different force fields and different platforms. However, the moisture content under most simulated hydration conditions is usually arbitrarily selected and generally restricts some movement directions of clay plates or even the interlayer water molecules which affect the mechanical response characteristics of the hydrated mineral system.

The goal of this study is to simulate the nanoscale mechanical properties of hydrated MMT minerals with three cations (Na, Ca, and Cs) at variable water contents. In the simulation, all degrees of freedom of atoms are released, and the water content in different hydration states is determined according to stable-state thermodynamic conditions. To the best of our knowledge, no such treatments have been reported to study the mechanical properties of the MMT mineral up to now. The quantitative values of the basal spacing for different compensation cationic MMTs are firstly studied. The relative stabilities of different states are determined by comparing the immersion energy and hydration energy between MMT systems, and, then, the corresponding moisture content of thermodynamics stable states is obtained. Finally, the nanoscale elastic properties of MMT are further simulated by the constant strain method under the stable state water content.

2. Model and Methods

2.1. Molecular Models

Wyoming-type MMT, which has been widely simulated, was selected as the model clay in this study, as shown in Figure 1. The initial unit cell parameters are a = 5.23 Å and b = 9.06 Å, and the c-axis increases with increasing water content, the unit cell angles are α = γ = 90°, β = 99°. The x and y axes were chosen parallel to the clay layer, while the z-axis was orthogonal to the clay layer. The unit cell of MMT was replicated (4a × 2b × 1c) in the x, y, and z directions so that the supercell contained 8-unit cells with approximate dimensions in the XY plane of 21 Å × 18 Å and 9.6 Å in the z-axis. One of every 8 octahedral Al³⁺ atoms was substituted with an Mg²⁺ atom and one of every 4 tetrahedral Si⁴⁺ atoms was substituted with an Al³⁺ atom. The isomorphic substitutions were distributed randomly obeying Lowenstein’s rule. This extent of substitution led to a layer charge of −0.75 e per unit cell and the unbalanced charge will be compensated by Na⁺, Ca²⁺, and Cs⁺ which are commonly adopted in engineering applications [12,35]. In addition, the interlayer space was randomly filled with a given number n of water molecules, the number n ranged from 0 to 160, and a total of 23 clay–water-ion systems will be evaluated. The structural formula of the simulated layer of Na- and Cs-MMT is M_0.75{Si_7.75Al_0.25}[Al_3.5Mg_0.5]O₂₀(OH)₄·nH₂O (M = Na or Cs), whereas for the Ca-MMT, it can be expressed as Ca_0.375{Si_7.75Al_0.25}[Al_3.5Mg_0.5]O₂₀(OH)₄·nH₂O.

2.2. Simulation Protocol

MD simulations were conducted in Forcite module in the Material studio 8.0 simulation package and the CLAYFF force field, developed by [36]. The latter is a flexible model which describes hydrated mineral systems through nonbonded electrostatic and Lennard–Jones terms. It can accurately reproduce structural and spectroscopic properties of clays, as well as dynamical and energetic properties of clay interlayer and aqueous interfaces. The interlayer water molecules were described with the flexible version of the simple point charge (SPC) potential. The potential energy of the system was evaluated with an 8.5 Å cutoff, a 0.5 Å spline, and buffer width for short-range van der Waals interaction. The Ewald summation for the electrostatic interaction was calculated with a precision of 10⁻³ kcal/mol. Periodic boundary conditions were imposed on three dimensions to avoid interface effects.

The initial conformation of hydrated MMT needs optimization of the geometry first by searching for the minimum potential energy of the system. A smart algorithm was chosen for the optimization with a maximum number of 5000 cycles, and the convergence thresholds for the specified maximum energy and displacement changes were set as 1.0 × 10⁻⁴ kcal/mol and 5.0 × 10⁻⁵ Å, respectively. Following this, the final configuration from optimization was used as the initial configuration for the equilibration stage simulation in the NPT (isothermal–isobaric) ensemble at P = 1 atm and T = 298 K. Each system was equilibrated for 100 ps (105 steps) with a time step of 1 fs. The temperature and the pressure were controlled by a Nosé–Hoover–Langevin thermostat of 1 ps and a Berendsen barostat of 0.1 ps, respectively. Finally, following the thermal equilibration configurations, a microcanonical NVE ensemble was then performed for a 1 ns production simulation, and the time step was still set as 1 fs (NVE: Each system in the ensemble has the same energy, and the number of particles and volume of each system are usually the same). The data from all 1 ns trajectories of the NVE simulation runs were recorded every 500 fs and were used to analyze the microstructure and thermodynamic properties of hydrated MMT. By comparing the relative thermodynamic energy of different systems, the steady-state water content of hydrated MMT can be determined. Furthermore, the mechanical properties of the steady-state system will be further simulated by the constant strain method. The constant strain method obtains the elastic constants via minimizing the energy of the system and deforming the supercell parameters in 12 directions in order. To improve the statistical accuracy, the simulations of mechanical properties are carried out based on the production stage trajectory files obtained above, that is, 20 frames of instant conformation are extracted from the corresponding trajectory files. The strain step is set as 10, and the maximum strain in each strain step is set to 1% according to the linear elastic assumption. The final mechanical parameters are obtained by the geometric average of the 20 frames.

2.3. Methods for Calculating the Elastic Constants

Acting on external forces, the hydrated MMT system will be in a state of stress. If the system is in equilibrium, the external forces must be exactly balanced by internal stress. In general, stress is a second rank tensor with nine components as in Equation (1):

$$\stackrel{\rightharpoonup }{\sigma }=\left[{\sigma }_{ij}\right]=\left[\begin{array}{ccc}{\sigma }_{xx}& {\sigma }_{xy}& {\sigma }_{xz}\\ {\sigma }_{xy}& {\sigma }_{yy}& {\sigma }_{yz}\\ {\sigma }_{xz}& {\sigma }_{yz}& {\sigma }_{zz}\end{array}\right]$$

(1)

In an atomistic calculation, the internal stress tensor can be obtained using the so-called virial expression as in Equation (2):

$$\stackrel{\rightharpoonup }{\sigma }=-\frac{1}{V_0}\left({\displaystyle \sum _{i=1}^N{m}_i({\stackrel{\rightharpoonup }{v}}_i{\stackrel{\rightharpoonup }{v}}_i^T)+{\displaystyle \sum _{i<j}{\stackrel{\rightharpoonup }{r}}_{ij}{\stackrel{\rightharpoonup }{F}}_{ij}^T}}\right)$$

(2)

where index i runs over all particles 1 through N; $m_i$, ${\stackrel{\rightharpoonup }{v}}_i$, and ${\stackrel{\rightharpoonup }{F}}_i$ denote the mass, velocity, and force acting on particle i, respectively; V₀ denotes the initial system volume.

The application of stress to the hydrated MMT system results in a change in the relative positions of particles within the system expressed quantitatively via the strain tensor as in Equation (3):

$$\stackrel{\rightharpoonup }{\epsilon }=\left[{\epsilon }_{ij}\right]=\left[\begin{array}{ccc}{\epsilon }_{xx}& {\epsilon }_{xy}& {\epsilon }_{xz}\\ {\epsilon }_{xy}& {\epsilon }_{yy}& {\epsilon }_{yz}\\ {\epsilon }_{xz}& {\epsilon }_{yz}& {\epsilon }_{zz}\end{array}\right]$$

(3)

For a parallelepiped (for example, a periodic simulation cell) characterized by the three column vectors $[{\stackrel{\rightharpoonup }{a}}_0 {\stackrel{\rightharpoonup }{b}}_0 {\stackrel{\rightharpoonup }{c}}_0]$ in some reference state, and by the vectors $[\stackrel{\rightharpoonup }{a} \stackrel{\rightharpoonup }{b} \stackrel{\rightharpoonup }{c}]$ in the deformed state, the strain tensor is given by Equation (4):

$$\stackrel{\rightharpoonup }{\epsilon }=\frac{1}{2}\left[{\left({\stackrel{\rightharpoonup }{h}}_0^T\right)}^{-1}\stackrel{\rightharpoonup }{G}{\stackrel{\rightharpoonup }{h}}_0^{-1}-1\right]$$

(4)

where ${\stackrel{\rightharpoonup }{h}}_0$ denotes the matrix formed from the three column vectors $[{\stackrel{\rightharpoonup }{a}}_0 {\stackrel{\rightharpoonup }{b}}_0 {\stackrel{\rightharpoonup }{c}}_0]$, $\stackrel{\rightharpoonup }{h}$ denotes the corresponding matrix forms from $[\stackrel{\rightharpoonup }{a} \stackrel{\rightharpoonup }{b} \stackrel{\rightharpoonup }{c}]$, $T$ denotes the matrix transpose, and $\stackrel{\rightharpoonup }{G}$ denotes the metric tensor ${\stackrel{\rightharpoonup }{h}}^T\stackrel{\rightharpoonup }{h}$.

Therefore, the elastic stiffness coefficients, relating the various components of stress and strain, are defined by Equation (5):

$$\stackrel{\rightharpoonup }{C}=\left[C_{ijkl}\right]={\frac{\partial {\sigma }_{ij}}{\partial {\epsilon }_{kl}}|}_{T,{\epsilon }_{kl}}={\frac{1}{V_0}\frac{{\partial }^{2}A}{\partial {\epsilon }_{ij}\partial {\epsilon }_{kl}}|}_{T,{\epsilon }_{ij},{\epsilon }_{kl}}$$

(5)

where $A$ denotes the Helmholtz free energy, and T denotes temperature.

For small deformations, the relationship between the stresses and strains may be expressed in terms of a generalized Hooke’s law:

$${\sigma }_{ij}=C_{ijkl}{\epsilon }_{kl}$$

(6)

and, alternatively:

$${\epsilon }_{ij}=S_{ijkl}{\sigma }_{kl}$$

(7)

where $C_{ijkl}$ and $S_{ijkl}$ denote stiffness tensor and compliance tensor, respectively.

Since both the stress and strain tensors are symmetric, it is often convenient to simplify these expressions by making use of Voigt vector notation. The stress is represented as in Equation (8):

$$\left[\begin{array}{ccc}{\sigma }_{xx}& {\sigma }_{xy}& {\sigma }_{xz}\\ {\sigma }_{xy}& {\sigma }_{yy}& {\sigma }_{yz}\\ {\sigma }_{xz}& {\sigma }_{yz}& {\sigma }_{zz}\end{array}\right]\Rightarrow \left[\begin{array}{ccc}{\sigma }_1& {\tau }_6& {\tau }_5\\ & {\sigma }_2& {\tau }_4\\ & & {\sigma }_3\end{array}\right]$$

(8)

The strain is represented as in Equation (9):

$$\left[\begin{array}{ccc}{\epsilon }_{xx}& {\epsilon }_{xy}& {\epsilon }_{xz}\\ {\epsilon }_{xy}& {\epsilon }_{yy}& {\epsilon }_{yz}\\ {\epsilon }_{xz}& {\epsilon }_{yz}& {\epsilon }_{zz}\end{array}\right]\Rightarrow \left[\begin{array}{ccc}{\epsilon }_1& {\gamma }_6& {\gamma }_5\\ & {\epsilon }_2& {\gamma }_4\\ & & {\epsilon }_3\end{array}\right]$$

(9)

The generalized Hooke’s law is thus often written as in Equation (10):

$$\left[\begin{array}{l}{\sigma }_1\\ {\sigma }_2\\ {\sigma }_3\\ {\tau }_4\\ {\tau }_5\\ {\tau }_6\end{array}\right]=\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& C_{14}& C_{15}& C_{16}\\ & C_{22}& C_{23}& C_{24}& C_{25}& C_{26}\\ & & C_{33}& C_{34}& C_{35}& C_{36}\\ & & & C_{44}& C_{45}& C_{46}\\ & & & & C_{55}& C_{56}\\ & & & & & C_{66}\end{array}\right]\left[\begin{array}{l}{\epsilon }_1\\ {\epsilon }_2\\ {\epsilon }_3\\ {\gamma }_4\\ {\gamma }_5\\ {\gamma }_6\end{array}\right]$$

(10)

alternatively:

$${\sigma }_i=C_{ij}{\epsilon }_j$$

(11)

where C_ij denotes the stiffness tensor and is a symmetric matrix of 6 × 6. The number of independent elastic constants for any system is 21. The mineral of MMT crystals belongs to the monoclinic system ($a\ne b\ne c, \alpha =\beta =90°\ne \gamma$) and has C2/M symmetry. It, then, only needs 13 independent parameters to describe its stiffness tensor. However, according to the study in [23], the MMT crystal can be considered approximately as orthotropic symmetry, thus further reducing the stiffness tensor parameters to 9 to characterize the sheet. These 9 independent elastic constants include 3 Young’s moduli (E₁, E₂, E₃), 3 shear moduli (G₂₃, G₃₁ and G₁₂), and 3 Poisson’s ratios (μ₁₂, μ₁₃ and μ₂₃). In terms of these parameters, the stiffness tensor C_ij can be further expressed as in Equation (12):

$$\begin{array}{l}\left[\begin{array}{cccccc}{C}_{11}& C_{12}& C_{13}& 0& 0& 0\\ C_{21}& C_{22}& C_{23}& 0& 0& 0\\ C_{31}& C_{32}& C_{33}& 0& 0& 0\\ 0& 0& 0& C_{44}& 0& 0\\ 0& 0& 0& 0& C_{55}& 0\\ 0& 0& 0& 0& 0& C_{66}\end{array}\right]=\\ \left[\begin{array}{cccccc}\frac{1-{\mu }_{23}{\mu }_{32}}{E_2{E}_3\Delta }& \frac{{\mu }_{21}+{\mu }_{31}{\mu }_{23}}{E_2{E}_3\Delta }& \frac{{\mu }_{31}+{\mu }_{21}{\mu }_{32}}{E_2{E}_3\Delta }& 0& 0& 0\\ \frac{{\mu }_{12}+{\mu }_{13}{\mu }_{32}}{E_3{E}_1\Delta }& \frac{1-{\mu }_{31}{\mu }_{13}}{E_3{E}_1\Delta }& \frac{{\mu }_{32}+{\mu }_{31}{\mu }_{12}}{E_3{E}_1\Delta }& 0& 0& 0\\ \frac{{\mu }_{13}+{\mu }_{12}{\mu }_{23}}{E_1{E}_2\Delta }& \frac{{\mu }_{23}+{\mu }_{13}{\mu }_{21}}{E_1{E}_2\Delta }& \frac{1-{\mu }_{12}{\mu }_{21}}{E_1{E}_2\Delta }& 0& 0& 0\\ 0& 0& 0& 2G_{23}& 0& 0\\ 0& 0& 0& 0& 2G_{31}& 0\\ 0& 0& 0& 0& 0& 2G_{12}\end{array}\right]\end{array}$$

(12)

where $\Delta =\frac{1-{\mu }_{12}{\mu }_{21}-{\mu }_{23}{\mu }_{32}-{\mu }_{31}{\mu }_{13}-2{\mu }_{12}{\mu }_{23}{\mu }_{31}}{E_1{E}_2{E}_3}$.

Therefore, the key to characterizing the elastic mechanical properties of the MMT–water-ion system at the nanoscale level is to obtain the parameters of the stiffness tensor or the compliance tensor.

3. Results and Discussion

3.1. Swelling and Thermodynamics Analysis

To obtain the water content corresponding to the thermodynamic steady-state in the hydration process of MMT, the strategy of simulating 23 different water content systems was adopted. The interlayer space of the MMT lamellae increases with the increase of water content, and this swelling process can be described by using the parameter of basal spacing. The basal spacing can be obtained from the view of water content by averaging the system volume during the production stage of simulation and defined as in Equation (13):

$$d_{001}=\frac{\langle V\rangle }{\langle a\rangle \langle b\rangle \sin \langle \alpha \rangle }$$

(13)

where <V> denotes the statistically averaged volume, and <a> <b> <α> are the statistically averaged parameters of the simulation supercell.

Figure 2 illustrates the variation of basal spacing of different hydrated montmorillonite systems as a function of water content from this study and available results from the literature. In the early stage of crystalline swelling, the basal spacing increases from step to step, while the step characteristic (plateau) disappears, and shows approximately linear variation in the following stage of osmotic swelling. The swelling curves of MMT with different compensation cations are different. The basal spacing of Cs-MMT is greater than Na- and Ca-MMT under the same water content. However, there is almost no difference between the latter two. This is mainly due to the fact that Cs⁺ ions are a heavy metal ion and their ion radius is about 1.7 times that of Na⁺ and Ca²⁺ ions, while the latter two’s ion radii are almost the same. The simulated results can well penetrate the distribution range of other available measurements.

It has been shown that the energy contribution to free energy dominates clay swelling and the entropy term only played a minor role. The local minima of the swelling energy curve correspond to the stable hydration state. Hence, the relative stabilities of different states can be determined by immersion energy and hydration energy. The immersion energy is defined as in Equation (14):

$$Q=\langle E(N)\rangle -\langle E(N^0)\rangle -\left(N-N^0\right)E_{bulk}$$

(14)

Here <E(N)> is the average potential energy of the hydrate with water content N, and <E(N⁰)> is the average energy of the referential hydration state (the state with the highest water content, 49%, was selected). E_bulk is the mean interaction potential of the bulk flexible SPC water; E_bulk = −41.5392 kJ/mol. The immersion energy Q is the energy released when MMT with water content N transforms into the referenced state with water content N⁰ by adsorbing water from bulk water.

The hydration energy is defined as in Equation (15)

$$\Delta E=\frac{\langle E(N)\rangle -\langle E(0)\rangle }{N}$$

(15)

Here <E(0)> is the average potential energy of the dry MMT. The hydration energy ΔE evaluates the energy change associated with water uptake by the dry clay.

The simulated hydration energy and immersion energy curves of MMT with different compensation cations are presented in Figure 3. For Na- and Cs-MMT the local minima occur when the interlayer water molecular number is 40 and 80 in the hydration energy curves, while Ca-MMT has only a global minimum at low water content. Correspondingly, in the immersion energy curves of Na- and Cs-MMT, local minima occur when the interlayer water molecular numbers are also 40 and 80. The global minima of Cs-MMT is found at n = 80, while Na-MMT global minima is located at n = 160, as shown by the enlarged illustration in Figure 3. Meanwhile, the variation of the Ca-MMT immersion energy curve decreases monotonically with the increase of water content, and the global minima occurs also at n = 160. Then, since the local minima of the swelling energy corresponds to the thermodynamic stable states, it can be deduced that the Na- and Cs-MMT will form the single- and double-layer stable hydrated states when the number of interlayer water molecules reach 40 and 80 respectively. The double-layer hydrate for Cs-MMT is in the most stable state. Thus, the Cs⁺ ions can inhibit swelling. The global minima of Na-MMT is located at n = 160. After forming the first and second hydration layer, Na-MMT will continue expansion to the osmotic swelling stage, namely that the Na⁺ ions promote swelling in contrast with Cs⁺. For Ca-MMT, its role in promoting expansion is even more pronounced. The intermediate state of the first or second hydration layer will not form during the swelling process, and the osmotic swelling stage will rapidly and directly develop. The nanoscale elastic properties of these stable state water content systems will be further simulated in the following sections.

3.2. Nanoscale Elastic Properties

To fully demonstrate the yield behavior of the MMT mineral, the maximum strain of each direction is expanded to 10% and the corresponding calculation step is increased to 100 steps. The stress–strain curves of the typical MMT system obtained from simulation are shown in Figure 4, where the interlayer compensation cation is Ca²⁺ ion, and the interlayer water number n is 80, corresponding to bilayer hydration.

The micromechanical properties of hydrated MMT show obvious anisotropy. The in-plane uniaxial tensile and compressive strength of MMT is approximately equal (for σ₁ and σ₂), while the tensile strength is slightly less than the compressive strength. When the tensile strain reaches about 7%, the material yields. Meanwhile, the compressive stress–strain relationship shows a good linear relationship within the maximum compressive strain (10%). The tensile and compressive strength (σ₃) of MMT orthogonal to the mineral plane is much lower than that in-plane strength (σ₁ and σ₂). This is mainly due to the appearance of the interlayer water in MMT lamellae, and the bonding strength between the MMT crystal and water molecule through non-bonding is much lower than that of covalent or ionic bonding in the mineral plane. In addition, the compression of MMT in the Z direction shows obvious non-linear characteristics and hardening under high pressure, which is different from the linear characteristics in the crystal plane.

The hydration molecular layer that mixed within MMT lamellae has a certain shear strength (τ4 and τ5) and also shows anisotropy along with different directions of the MMT crystal plane. The shear stiffness of the interlayer water along the Y direction (τ5) is slightly larger than that in the X direction (τ4), while its shear strength is slightly lower. The shear capability of MMT lamellae is very high (τ6) and shows good linear characteristics in the simulated strain ranges.

By comparing the quantitative values of in-plane strength, it was found that the in-plane compressive strength of hydrated MMT is the largest, followed by tensile strength and shear strength. Similarly, this variation also applies to the strengths that are perpendicular to the mineral plane.

All of the stress responses show good linearity in the range of 1% small strain, and this is illustrated in the lower right of each block in Figure 4. This also shows that the linear elastic constants calculated by the small strain of 1% are reliable. The stiffness tensor and linear elastic constants of typical MMT in the 1% strain range corresponding to Figure 4 are listed in

Table 1. Typical stiffness coefficients and elastic constants of hydrated Ca-MM (GPa).

${\mathit{C}}_{\mathit{i}\mathit{j}}$	1	2	3	4	5	6
1	189.6710 ± 1.0140	81.7433 ± 0.8172	10.1290 ± 0.4701	1.2430 ± 0.3777	−10.4597 ± 0.2771	0.4964 ± 0.4787
2		212.5464 ± 1.7814	8.4054 ± 0.4993	3.5123 ± 0.4520	−3.4307 ± 0.3618	1.1526 ± 1.0715
3			28.4335 ± 0.5057	0.3278 ± 0.2660	1.3468 ± 0.1654	0.1060 ± 0.2955
4				4.2630 ± 0.1881	−0.6221 ± 0.1371	−2.7480 ± 0.2758
5	Symmetric				5.4343 ± 0.1049	0.3864 ± 0.2236
6						45.4441 ± 0.6942
Young’s modulus	$E_1$	139.1605	$E_2$	174.2115	$E_3$	27.0289
Shear modulus	$G_{23}$	2.1315	$G_{31}$	2.7172	$G_{12}$	22.7221
Poisson’s ratios	${\mu }_{12}$	${\mu }_{21}$	${\mu }_{13}$	${\mu }_{31}$	${\mu }_{23}$	${\mu }_{32}$
	0.3463	0.4335	0.3436	0.0667	0.1199	0.0186
Bulk modulus	70.1340
Compressibility (1/TPa)	38.6319

.

The simulated stiffness tensor components in Table 1 are in good agreement with simulations and measurements in the literature. The in-plane components C₁₁ and C₂₂ of the MMT with water bilayer at 300 K are 189 GPa and 191 GPa, respectively, from [25], which are almost equal to the results of this study that C₁₁ = 189 GPa and C₂₂ = 212 GPa. Furthermore, there are differences in the in-plane Young’s moduli E₁ and E₂ which are derived from the stiffness tensor, and these differences are induced by the different topologies of trivalent cations in dioctahedral structure of MMT, as shown in Figure 5. However, the in-plane Young’s moduli (E₁ and E₂) are all much higher than the Young’s modulus perpendicular to the plane (E₃). On the contrary, the in-plane shear moduli (G₂₃ and G₃₁) are lower than the shear modulus perpendicular to the plane (G₁₂). The larger difference of Poisson’s ratio in different directions once again shows the anisotropy of the MMT structure. The boning strength in the mineral plane is much higher than that of non-bonding strength between mineral lamellae.

The evolution of elastic properties of hydrated MMT under different stable state water content will be further illustrated by analyzing the changes of some components in the stiffness tensor. To validate the study, the simulation results from [29] and UPV test results from [17] are adopted for comparison due to the high correlation with the work in this study. The authors in [29] gave the variation of stiffness components during the hydration process of Na-Wyoming MMT with a water content of 0–40%. The authors in [17] obtained the range of elastic stiffness constants of smectite minerals on the microscale based on UPV test results of clay composition in shale.

For compression of clay under high pressure, the most concerning component of the stiffness tensor is C₃₃, that is, the ability to resist deformation in the direction of perpendicular to the mineral plane. Figure 6 shows the variation of C₃₃ obtained upon hydration by [29], Na-, Cs-, Ca- MMT in this study under different stable state conditions, and the test range results of [17]. Since the stiffness tensor elements in this study are obtained by a statistical average of 20 frame trajectory files, the statistical standard deviation of data is also given in Figure 6.

Ebrahimi et al. (2012) revealed that C₃₃ decreased sharply from dry MMT to the formation of the first hydrated layer between lamellae, then increased to a local maximum with the formation of the stable hydrated layer, and finally maintained a constant until the basic spacing d > 15.6 Å [29]. The values of C₃₃ obtained in this study are in good agreement with the results of [29], apart from the first hydrated layer. The C₃₃ of Ca-MMT is slightly higher than that of Na-MMT, while the C₃₃ of Cs-MMT is the lowest among all the same hydrated states. Apart from the dry MMT, all of the simulated C₃₃ fall well into the given range of UPV measured by [17] (C₃₃ = 11.5~30.4 GPa).

The in-plane stiffness tensor components C₁₁, C₂₂, and C₁₂, which reflect the Poisson effect in the plane, are shown in Figure 7 with basal spacing. C₁₁, C₂₂, and C₁₂ all decrease nonlinearly with the increase of water content. These components are mainly controlled by the ionic or covalent bonding strength of mineral atoms and the geometric conformation of the hydrated MMT system. In addition, the simulation results of C₁₁ and C₁₂ are in good agreement with those of [29], but C₁₂ is slightly lower. The reason may lie in the difference between our simulation structure and the octahedral isomorphic substitutions of [29]. However, the measured range of UPV given by [17] (C₁₁ = 13.8~46.1 GPa and C₁₂ = 6.9~17.8 GPa) is much lower than the simulated results. This is because the MMT mineral lamellae are very thin, and the relative slip between lamellae easily occurs in the UPV test. However, the simulated results in [25] (C₁₁ = 189~231 GPa) and test results by Brillouin scattering on mica minerals in [26] (C₁₁ = 181 GPa and C₂₂ = 178 GPa) both show that the results of this study are credible.

For the shear of clay under high pressure, the most concerned is the shear resistance of limited interlayer water molecules in the hydrated MMT system, i.e., the in-plane shear stiffness tensor components C₄₄ and C₅₅. Their variations with basal spacing are shown in Figure 8; C₄₄ and C₅₅ fluctuate with the increase of water content which is like the variation of C₃₃. Therefore, the mechanism of shear strength of interlayer water is closely related to the formation of H-bonds and atom-free volume. The shear strength between dry MMT lamellae is much larger than that of hydrated MMT lamellae, that is to say, the addition of water plays a lubricating role in interlayer friction, which makes C₄₄ and C₅₅ decrease sharply. However, these limited interlayer water molecules do not completely show the zero shear strength characteristics of free water. As water content increases, the in-plane shear strength increases to a local maximum when the system approaches the thermodynamic stable state and minimum atom free volume occurs. In the range of simulated water content, C₄₄ and C₅₅ eventually tend to be constant. Furthermore, the results in Figure 8 shows that the components C₄₄ and C₅₅ are in good agreement with those obtained by [29], and the in-plane shear strength between different MMT systems is close except for the dry MMT system. Meanwhile, the simulated results also fall well into the UPV measurement range of [17] (C₄₄ = 2.9~8.9 GPa).

The shear of hydrated MMT along the direction perpendicular to the mineral plane is mainly controlled by the relative torsion of bonded atoms. The variations of C₆₆ with basal spacing are shown in Figure 9. With the increase of interlayer water content, C₆₆ decreases continuously, which is mainly due to the increase of interlayer water contribution in the weak shear zone. Considering that the ratio of width to height of a true MMT lamella is very large (about 100:1 to 1000:1), shear along the direction of orthogonal to the mineral plane is almost impossible, so the variations of C₆₆ will not be discussed here.

Young’s moduli of the different MMT systems under various hydration can be obtained by further conversion of the stiffness tensor, and the variation of Young’s moduli with basal spacing is shown in Figure 10. The in-plane Young’s moduli E₁ and E₂ are consistent with C₁₁ and C₁₂ that decrease nonlinearly with the increase of water content. For the Young’s modulus of E₃ which is perpendicular to the mineral plane, its value for the dry MMT system is 2–3 times that of the hydrated system for Na-MMT and Ca-MMT. E₃ of the hydrated MMT system also decreases nonlinearly with the increase of water content. Interestingly, for Cs-MMT, the Young’s modulus, E₃, of different hydrated MMT systems, even dry MMT systems, is almost the same. Preliminary analysis shows that this may be related to the larger radius of Cs⁺ ions. For the dry MMT system, the interlayer spacing is large, which induces a low compression strength. After hydration, the Cs⁺ ions always form an inner-sphere complex, which has little effect on the strength, thus making all Cs-MMT systems of E₃ almost the same within the range of simulated water content. However, this interesting phenomenon has not yet been reported in the literature, so the above corollary needs further study in the future.

From the above comparisons, molecular microscopic testing techniques, such as UPV, can only qualitatively analyze the strength of clay minerals and the test results are not accurate enough. Molecular dynamics simulation provides a good way to understand the properties of hydrated clay minerals on an atomic scale. In particular, unconstrained system techniques and thermodynamic stable-state water content analysis techniques are used to make the calculation results of molecular dynamics more reliable.

4. Conclusions

Molecular dynamic simulations were conducted to investigate the mechanical properties of montmorillonite with different compensation cations upon hydration with unconstrained system atoms, expanding further studies of nanoscale elastic properties. The conclusions are as follows:

(1)

The basal spacing of all MMT systems increases from step to step in the initial crystalline swelling stage. While this step characteristic (plateau) disappears, the basal spacing shows approximately linear variation in the subsequent osmotic swelling stage. However, the quantitative values of the basal spacing for different compensation cationic MMTs are different, which is related to the ion radius of cations.

(2)

The nanoscale elastic properties of hydrated MMT show obvious anisotropy. The tensile and compressive strength (σ₃) of MMT orthogonal to the mineral plane is much lower than the in-plane strengths (σ₁ and σ₂), and the hydration molecular layer mixed within MMT lamellae has certain shear strengths (τ₄ and τ₅). The order of strength values from large to small is the compressive strength, the tensile strength, and the shear strength, respectively.

(3)

The components of C₃₃, C₄₄, and C₅₅ all fluctuate with the increase of interlayer water content, while the variation of the in-plane stiffness tensor components C₁₁, C₂₂, and C₁₂, all decrease nonlinearly with the increase of water content. These components are mainly controlled by the bonding strength of mineral atoms and the geometric conformation of the hydrated MMT system. The C₆₆ component reflects the shear of hydrated MMT along the direction perpendicular to the mineral plane, which is mainly controlled by the relative torsion of bonded atoms. Due to the increase of interlayer water contribution in the weak shear zone, C₆₆ decreases continuously with the increase of interlayer water content. Young’s moduli all decrease nonlinearly with the increase of water content.

The simulation results reveal that the water molecules linked to the mineral surface will be arranged in an ordered structure (H-bond formation and thermodynamic dense stable state) and their shear strength are not zero. Preliminary discussion shows that this is related to the strength mechanism of clay at the microscale. Building a bottom-up procedure toward a constitutive law for soil mechanical still needs further research.