Clay minerals are layered aluminosilicates occurring in sediments, soils and soft rocks. Clays have various uses and applications in pottery, construction and the mining industry (Murray, Reference Murray2000). They have also been used in new technological fields over recent decades. For example, clay minerals are used as fillers in polymer nanocomposites (Chen et al., Reference Chen, Evans, Greenwell, Boulet, Coveney, Bowden and Whiting2008). These materials have excellent properties, such as high dimensional stability, reduced gas permeability and enhanced mechanical properties (Podsiadlo et al., Reference Podsiadlo, Kaushik, Arruda, Waas, Shim, Xu and Ramamoorthy2007; Ham et al., Reference Ham, Kim and Chang2013; Abdi-Khangah et al., Reference Abdi-Khangah, Barati and Zhang2018). Moreover, various molecules may be intercalated between clay layers, enlarging the volume of the clays and improving their traditional applications (Kakegawa & Ogawa, Reference Kakegawa and Ogawa2002; Rutkai et al., Reference Rutkai, Mako and Kristof2009). Recently, they have also been used as natural or artificial barriers in radioactive waste depositories due to their low permeability and long-term healing capacity (Mondol et al., Reference Mondol, Bjorlykke, Jahren and Hoeg2007; Wersin et al., Reference Wersin, Johnson and McKinley2007). Although these materials have numerous applications, little work has been carried out on their mechanical properties. This is mainly due to the fact that the small grain size makes it difficult to isolate an individual clay crystallite large enough to satisfy experimental requirements. In addition, their high affinity for water and fine-grained polycrystalline nature make it challenging to measure their physical and chemical properties. Therefore, to fully understand the complex structure of these materials and the critical mechanisms of intercalation, it has become necessary to use modern molecular simulation methods, especially because many of the layered materials are restricted to nano-sized morphologies and are less suitable for conventional experimental analysis.
The atoms in pure clay minerals are arranged regularly and mechanically more durable. However, defects are distributed randomly in clay minerals, which lead the physical and chemical properties of clays to exhibit great variability. For example, Nisar et al. (Reference Nisar, Århammar, Jämstorp and Ahuja2011) performed first-principles calculations for defective kaolinite in order to determine its electronic structure, atomic geometry and defect formation energies. The results implied that the presence of intrinsic defects will change the electronic structure of kaolinite and influence properties such as cation-retention capacity and drug release. He et al. (Reference He, Fang and Zhang2009) investigated the formation energetics and charge-doping properties of the extrinsic substitutional defects in kaolinite, which excluded the formation possibility of the defects CuAl and ZnAl in natural kaolinite. Li et al. (Reference Li, Li and Yang2015) studied comprehensively the adsorption behaviour of metal ions on defective kaolinite by using molecular dynamics (MD). They indicated that both the amount and stability of adsorption were enhanced by increasing the number of defects and their sizes. The effect of intrinsic defects, such as vacancies in kaolinite, on the mechanical properties of this mineral has, to our knowledge, not been investigated previously.
Molecular dynamics simulation is an effective tool to probe molecular-level structural, physical and mechanical properties, which are rather difficult to explore directly by macroscopic experiments. Density functional theory (DFT) calculations are also a powerful tool to investigate the physical and chemical properties of materials. However, a very limited number of atoms can be built up by using DFT calculations. In contrast, the MD method can simulate millions of atoms and extract deformation and failure details of a system at the atomic scale. Many researchers have used MD methods to investigate the micro-properties of clay minerals and have obtained satisfactory results (Cygan et al., Reference Cygan, Greathouse, Heinz and Kalinichev2009). For example, Teich-McGoldrick et al. (Reference Teich-McGoldrick, Greathouse, Jovécolón and Cygan2015) used MD and Monte Carlo simulations to study the effects of layer charge location, interlayer cation and temperature on the swelling of montmorillonite and beidellite clay minerals. Zhang et al. (Reference Zhang, Liu, Cheng, Li, Zeng and Frost2014) investigated the intercalation behaviour and structure of nanoconfined dodecylamine in kaolinite at the atomic scale. Many studies have also focused on the investigation of the mechanical properties of clay minerals by using MD simulations. For example, Heinz et al. (Reference Heinz, Vaia and Farmer2006) computed cleavage energies for different clay minerals, and these results showed agreement with experimental results. Hantal et al. (Reference Hantal, Brochard, Laubie, Ebrahimi, Pellenq, Ulm and Coasne2014) used the CLAYFF (Cygan et al., Reference Cygan, Liang and Kalinichev2004) and REAXFF (van Duin et al., Reference van Duin, Dasgupta, Lorant and Goddard2001) force fields to study the elastic and failure properties of illite at the atomic scale. Teich-McGoldrick et al. (Reference Teich-McGoldrick, Greathouse and Cygan2012) used the CLAYFF force field to obtain the mechanical properties of muscovite and observed better agreement with experiment results than earlier DFT calculations. Few of these simulations have involved kaolinite, one of the most abundant clay minerals in the Earth's crust, although Benazzouz and Zaoui (Reference Benazzouz and Zaoui2012) determined the elastic constant tensors of kaolinite under pressure based on MD methods. However, there is still a lack of comprehensive research on the microstructural details of kaolinite under deformation. A fundamental question is whether or not the defects would lead to changes in the mechanical properties of kaolinite.
In this study, the effects of defects and temperature on the mechanical properties of kaolinite were studied using classical MD simulations. Kaolinites with four types of defect were constructed. The uniaxial tensile stress–strain responses of different defect types, locations and densities were obtained. In addition, the basic mechanical properties, such as Young's modulus and maximum stress, for different defect densities were derived and are discussed in detail. Finally, the temperature effects of kaolinites with different defect types were investigated.
Simulation details
Crystal structures
Kaolinite is a 1:1 clay mineral that is abundant in the Earth's crust. The kaolinite structure consists of a siloxane tetrahedral sheet and a gibbsite-type octahedral sheet (Fig. 1). In the Si tetrahedral sheet (Fig. 1a), each Si atom is connected to three other Si atoms by sharing oxygen atoms, forming the siloxane surface. This arrangement creates hexagonal cavities on the surface. In the Al octahedral sheet (Fig. 1b), each Al atom is coordinated by two vertex oxygen atoms of the siloxane surface and shares four hydroxyls with neighbouring Al atoms. These two sheets are bonded through hydrogen bonds (Fig. 1c).
Structure of kaolinite viewed from (a) Si tetrahedral sheets, (b) Al octahedral sheets and (c) side directions.
Defects in kaolinite and montmorillonite have been investigated experimentally, and different types of defect have been observed to be associated with silanol nests (Yamagishi et al., Reference Yamagishi, Namba and Yashima1991; Fickel et al., Reference Fickel, Shough, Doren and Lobo2010; Bayahia et al., Reference Bayahia, Kozhevnikova and Kozhevnikov2013). Silanol defects have also been detected in other clay minerals (Libowitzky & Beran, Reference Libowitzky and Beran1995; Khraisheh et al., Reference Khraisheh, Al-Ghouti, Allen and Ahmad2005; Bobon et al., Reference Bobon, Christy, Kluvanec and Illasova2011), and it has been proposed that a significant fraction of carbon resources reside in the defects of minerals in the Earth's mantle. Here, MD simulations are employed to simulate the mechanical properties of kaolinite with different silanol defects (Fig. 2). The defect sizes follow the order Si1 < Si2 < Si3 < Si4, with removal of 1, 4, 6 and 10 lattice Si atoms, respectively. The dangling O atoms of defect sites were saturated with H atoms to ensure electroneutrality. In this study, models with 512 unit cells (16 × 8 × 8 along OA, OB and OC directions) were built and then different defect types were created on the kaolinite models by deleting Si atoms and saturating with O atoms.
Enlarged microstructures of various defect types (Si1, Si2, Si3 and Si4) created at the tetrahedral SiO4 surface. Oh is from the hydroxyl of silanol and Oa (bridging O) is also bonded to the silanol Si atom, while Ob (bridging O) has no relation to the defect sites.
Method
Selecting an appropriate interatomic potential is a key task in MD simulations for describing accurately the interactions between atoms. The CLAYFF force field (Cygan et al., Reference Cygan, Liang and Kalinichev2004) is highly efficient at modelling the crystal structures of kaolinite and montmorillonite (Larentzos et al., Reference Larentzos, Greathouse and Cygan2007; Cygan et al., Reference Cygan, Romanov and Myshakin2012). In this study, the CLAYFF force field was used to investigate the mechanical behaviour of kaolinite with defects.
The total potential energy includes contributions from Coulombic interactions, van der Waals interactions and any bonded interactions:
$$\eqalign{& E_{{\rm total}} = E_{{\rm Coul}} + E_{{\rm VDW}} + E_{{\rm bond} \hbox {-} {\rm strech}} \cr &= \displaystyle{{e^2} \over {4\pi \varepsilon _0}}\sum\limits_{i\ne j} {\displaystyle{{q_iq_j} \over {r_{ij}}}} + \sum\limits_{i\ne j} {D_{0,ij}} \left[\left(\displaystyle{{R_{0,ij}} \over {r_{ij}}}\right)^{12}-2\left(\displaystyle{{R_{0,ij}} \over {r_{ij}}}\right)^6\right] + k_1(r_{ij}-r_0)^2} $$where r ij is the distance between atoms i and j. The first term is Coulombic energy, where e is the charge of an electron, ε0 is the dielectric permittivity of a vacuum (8.85419 × 10–12 F/m) and q i and q j are the atomic partial charges, which can be derived from quantum mechanics calculations. The second term is van der Waals energy, D 0,ij and R 0,ij, which can be derived as:
$$R_{0,ij} = \displaystyle{1 \over 2}(R_{0,i} + R_{0,j}) $$
$$D_{0.ij} = \sqrt {D_{0,i}D_{0,j}} $$The third term represents the bonded stretch energy (O–H), where k 1 = 554.1349 kcal/mol Å2 is a force constant and r 0 = 1 Å is the equilibrium bond length. The force field parameters of saturated Oa–H pairs are the same as the hydroxyl groups (O–H) in the CLAYFF force field.
All of the simulations were performed using open-source parallel code LAMMPS (Plimpton, Reference Plimpton1995). The OA, OB and OC coordination axes correspond to the crystallographic orientations [1 0 0], [0 1 0] and [0 0 1], respectively. The periodic boundary condition is used for all directions to remove possible surface effects. The NVT (with a constant particle number–volume–temperature at 300 K) and NPT (with a constant particle number–pressure–temperature at 300 K and 1.0 atmosphere) ensembles were implemented here with a Nosé–Hoover thermostat and barostat.
The mechanical properties of regular and defective kaolinite were investigated using MD simulations. First, the regular kaolinite was built by reduplicating the unit cell in OA, OB and OC directions. Then, defective kaolinite models with different defect densities were generated by removing Si atoms in the SiO4 surface based on the regular kaolinite model. Finally, tensile tests were carried out by continuously increasing deformation along their corresponding axes; three uniaxial tension tests were performed on an equilibration configuration. The stress was determined by summing up the internal forces at individual particle points and dividing by the initial system area. To obtain a freestanding configuration, a given system was relaxed for 50 ps through NVT (with a constant particle number–volume–temperature at 300 K) and NPT (with a constant particle number–pressure–temperature at 300 K and 1.0 atmosphere) ensembles sequentially, and it reached an equilibrium state. Previous studies (Rottler & Robbins, Reference Rottler and Robbins2003; Lyulin et al., Reference Lyulin, Li, Mulder, Vorselaars and Michel2006; Mahajan & Basu, Reference Mahajan and Basu2010) have indicated that MD is a useful tool for investigating qualitatively the relationship between micro-scale phenomena and macro-scale mechanical properties and, to some extent, to determine quantitatively these terms. However, there is significant discrepancy between the simulated and experimental strain rates. Due to the computational costs and the time required to produce results that are qualitatively similar to experimental results, MD simulations typically use strain rates that are about six to ten orders of magnitude larger than the highest strain rates commonly used in experiments (Sahputra & Echtermeyer, Reference Sahputra and Echtermeyer2013). We tested the effect of strain rate on mechanical behaviour using three different strain rates (5 × 108, 5 × 109 and 5 × 1010 s–1). Negligible change was observed in the results as a function of strain rate (Fig. S1 in Supplementary Information). In order to save computational time and to obtain sufficient configurations under different strains, one cell direction was continuously deformed by a strain rate of 5 × 109 s–1. The timestep is set to be 1.0 fs and the motion of atoms was integrated by using a Verlet leapfrog algorithm.
Results and discussion
Effect of defect location
Kaolinite is a 1:1 layered clay mineral with two types of defect in each layer. In type I, the locations of the defects are the same in each layer, and in type II, the defects are distributed randomly in each layer (Fig. 3). The stress–strain curves of Si1 kaolinite with different locations of defects are shown in Fig. 4, where point A refers to the initial system without strain and points B, C and D are the strained states at the maximum stress, under inflection and after failure, respectively. Three types of defect location are selected randomly for type II kaolinite. The stress–strain curves of defect location for type II kaolinite are nearly coincident (Fig. 4). In addition, the tensile strength and Young's modulus of kaolinite are almost the same for both defect types, suggesting that the location of defects does not affect the mechanical properties of kaolinite. In addition, the stress increases linearly with increasing strain in the elastic stage before decreasing sharply at the maximum stress, which is due to the fact that the fracture occurred in each layer at the same time. The microstructure change in kaolinite was tracked at different defect locations during tension (Fig. 5). The crack was initiated at the defect sites and propagated very quickly once crack nucleation appeared, which led to the stress drop. Moreover, for the OA and OB directions, the crack growth trajectory of type I kaolinite is linear since the defect locations are the same in each layer. However, the cracks propagated randomly in each layer for type II kaolinite. The deformation processes of types I and II defective kaolinites are the same when loading along the OC direction as only hydrogen bonds between layers are broken.
Schematic diagram of kaolinite with defect locations of (a) type I, and (b) type II. The crosses within circles represent defect sites.
Stress–strain curves of Si1 kaolinite along the (a) OA, (b) OB and (c) OC directions.
The corresponding microstructures of kaolinite with type I and II defect locations at four specific points in the (a) OA, (b) OB and (c) OC directions.
Effect of defect density
Si1 defective kaolinite models with different defect densities were built. Here, the defect density is defined as the ratio of the number of defective nests generated to the total number of unit nests in the kaolinite structure. The Young's modulus and tensile strength of regular kaolinite and those of kaolinite with different defect densities of Si vacancies at room temperature are shown in Table 1. The Young's moduli and maximum stress decreased with increasing defect density, which indicates that more vacancy defects degrade the stiffness of materials. Sato et al. (Reference Sato, Ono, Johnston and Yamagishi2005) calculated the Young's modulus of regular kaolinite at three directions by using DFT calculations, with the corresponding values of Young's moduli in the OA, OB and OC directions being 92, 170 and 31 GPa, respectively, which are in accord with the MD simulation results of this study. Furthermore, in comparison with kaolinite without defects, the Young's moduli of defective kaolinite decrease more sharply in the OB direction, suggesting that the Si vacancies affect severely the mechanical properties of kaolinite in the OB direction. This is because the interaction force decreased most significantly in the OB direction (Fig. 6). However, the Young's moduli in the OC direction hardly changed with the variation of defect density, because the OH···O hydrogen bonds, which are the dominant interactions between layers and Si vacancies, have little effect in the OC direction.
Young's moduli and maximum stresses of kaolinite with different defect densities in the OA, OB, and OC directions.
Interaction force in the OA, OB and OC directions.
Effect of defect type
The atom density of defective kaolinite is defined as:
$$\rho {\rm}={\rm } \displaystyle{m \over V} $$where m is the mass of atoms and V is the volume of the system. The defect density of 50% was selected to investigate the effect of the defect type of kaolinite. The numbers of Si atoms removed in Si1, Si2, Si3 and Si4 systems are 256, 1024, 384 and 768, respectively. The atom density and Young's modulus of kaolinite with different defect types in the OA and OB directions both follow the order Si2 < Si4 < Si3 < Si1 < Si0 (Table 2). Here, Si0 represented the system with no defects. The normalized Young's moduli are shown in Fig. 7 with increasing numbers of removed atoms. Here, the moduli were normalized by using the Young's modulus of the Si0 system. The results indicate that the more atoms that were removed, the weaker the mechanical properties of kaolinite. This is because removing more Si atoms leads to lower atom density, which weakens interactions between atoms and decreases the Young's modulus and strength of kaolinite.
Atom densities and Young's moduli of kaolinite with different defect types.
Normalized Young's moduli of systems with different numbers of removed atoms.
Effect of temperature
Four temperatures (300, 400, 500 and 600 K) and a defect density of 50% were selected to investigate the temperature effect of kaolinite with different types of defect. The normalized Young's moduli in the OA, OB and OC directions are shown in Fig. 8. Here, the normalized moduli of the three directions were calculated by using their respective values at 300 K. The moduli decrease with increasing temperature in all of the directions. This is due to the fact that the volume of kaolinite expanded with increasing temperature, which weakened the interactions between atoms (Fig. 9a). Moreover, the greater the volume expansion, the faster the modulus decrease, and thus the temperature sensitivity of defective kaolinite follows the order Si2 > Si4 > Si3 > Si1, which has the same trend as that for atom density. In addition, the modulus in the OC direction decreased the most. As is shown in Table 2, Si vacancies affected most significantly the mechanical properties in the OA and OB directions. The OC direction is dominated by hydrogen bonds, which indicate that basal d spacing is more sensitive to temperature compared with the Si defects.
Normalized Young's moduli of kaolinite at different temperatures in the (a) OA, (b) OB and (c) OC directions.
(a) Normalized volumes and (b) d spacings of kaolinite at different temperatures.
Conclusions
The mechanical properties of kaolinite with Si vacancies were investigated using MD simulations. Uniaxial tension tests of defective kaolinite models under various conditions, including defect type, location and density effects, were carried out, and the stress–strain curves in three crystal directions were obtained. The stress–strain relationships of kaolinite with different defect locations in each layer are almost the same, and the defect sites are the weakest points, thereby initiating and propagating the cracks. The Young's modulus and maximum stress exhibit significant decreases as defect density or temperature increase. The decrease with temperature seen in the mechanical properties might be attributed to the thermal expansion of defective kaolinite in the OA, OB and OC directions. The results obtained will be helpful in understanding related defective minerals.
Supplementary material
To view supplementary material for this article, please visit https://doi.org/10.1180/clm.2019.22
Acknowledgements
This work was supported by the National Key Research and Development Program (2016YFC0600901) and the Chinese Postdoctoral Foundation (2017M611043).
