Nanoscale Tensile, Shear, and Failure Properties of Layered Silicates

J. Phys. Chem. C 2010, 114, 1763–1772

1763

Nanoscale Tensile, Shear, and Failure Properties of Layered Silicates as a Function of Cation Density and Stress Gregory D. Zartman,† Hua Liu,† Brahim Akdim,‡ Ruth Pachter,‡ and Hendrik Heinz*,† Department of Polymer Engineering, UniVersity of Akron, Akron, Ohio 44325 and Nanostructured and Biological Materials Branch, Materials and Manufacturing Directorate, Wright-Patterson AFB, Ohio 45433 ReceiVed: July 23, 2009; ReVised Manuscript ReceiVed: December 2, 2009

Mechanical properties of layered silicates on the nanometer scale have been associated with large uncertainty. We attempt to clarify the linear elastic properties including tensile moduli, shear moduli, and potential failure mechanisms for the minerals pyrophyllite, montmorillonite, and mica in the order of increasing cation exchange capacity (CEC) under a broad range of stress using electronic structure calculations, semiempirical classical molecular dynamics simulation, and the comparison to available macroscopic experimental data. In-plane tensile moduli (xx and yy) are ∼160 GPa independent of CEC and stress, whereas perpendicular tensile moduli (zz) range from 5 to 60 GPa as a function of CEC at low stress (0.01 to 1 GPa) and approach in-plane values at high stress. In-plane shear moduli (xy) are ∼70 GPa independent of CEC and the shear strength increases from ∼1 to ∼3 GPa with increasing cation density. Shear moduli parallel to the layers (xz and yz) are between 2 and 20 GPa as a function of CEC, with a shear strength of 0.2 to 1 GPa beyond which the layers exhibit lateral shear flow. Tensile zz moduli, shear moduli, and shear strength in the xz and yz direction reach a local minimum at a cation density of 0.3 relative to mica. The simulation suggests sliding of the layers, in-plane kinks, and cation intrusion into the layers as potential failure mechanisms equal to amorphization on the macroscale. The anisotropy and stress-dependence of the mechanical properties is determined by the presence of rigid layers and flexible interlayer spaces of variable cation density. Current classical simulation models tend to overestimate in-plane moduli (xx, yy, xy) in a systematic way relative to electronic structure (DFT) and experimental results. 1. Introduction Aluminosilicates form an important fraction of the earth’s surface and are widely used in polymer-clay nanocomposites for automotive and aerospace parts, biocompatible composite materials, and packaging, as well as in drilling liquids and cosmetics.1–6 Dispersion of clay minerals and organically modified clay minerals7 in synthetic and biological polymer matrices2–6 can lead to improved mechanical properties such as higher modulus and yield strength. Gas diffusion rates in thin polymer films for packaging applications can also be reduced by several orders of magnitude through the inclusion of a few volume percent of exfoliated clay platelets of high aspect ratio. Moreover, layered silicates can be employed to modify the action of detergents and cosmetics on fabrics and human skin.8 Easy access to even surfaces on the scale of nanometers to micrometers as well as the widely tunable range of the cation density renders clay minerals such as pyrophyllite, montmorillonite, and mica as interesting substrates for the assembly of complex interfaces (Figure 1).9–13 The minerals consist of well-defined layers separated by interlayer spaces, which may contain alkali, earth alkali, and other cations. Each layer consists of two outer tetrahedral silica sheets covalently connected to a central octahedral aluminum oxide hydroxide sheet, thus often called 2:1 clay minerals. The presence of Si f Al and Al f Mg substitution leads to negative charges within the layers, which are compensated by potassium * To whom correspondence should be addressed. E-mail: hendrik.heinz@ uakron.edu. † University of Akron. ‡ Wright-Patterson AFB.

interlayer ions in mica and by sodium as well as other interlayer cations in montmorillonite. Therefore, mica contains predominantly SiO2 f AlO2- · · · K+ substitutions in the outer tetrahedral silica sheets and montmorillonite contains predominantly AlO(OH) f MgO(OH)- · · · Na+ substitutions in the central octahedral aluminum oxide hydroxide sheet. Pyrophyllite contains no such defect sites and thus no alkali cations in the interlayer space. The cation density in the interlayer and the cation exchange capacity (CEC) is dictated by the concentration of defect sites7,14 and their spatial distribution has been characterized by solid-state NMR methods (Figure 1).15,16 In the exfoliated state, the thickness of individual layers of the minerals varies between 0.919 nm (pyrophyllite) and ∼1.3 nm (montmorillonite and mica with superficial cations on both sides).17,18 Equilibrium cell parameters and the equilibrium gallery spacing in experiment agree with those in force-fieldbased simulation models17 within 1% (Table 1) and with those in DFT simulation within 3%. Mechanical properties of clay minerals such as responses to tensile, shear, and bending strain as well as failure mechanisms at the nanometer scale are important to understand the complex mechanical behavior of polymer nanocomposite materials in automotive, aerospace, biomedical, and thin film applications. The rigidity of soft polymeric host materials increases through the presence of harder, well-dispersed layered silicates, particularly in the case of high aspect ratios in excess of 102. In spite of intense research efforts on clay-hybrid materials over the past decades and the availability of mechanical data for selected clay minerals on the macroscale,25–37 nanoscale mechanical properties remained rather speculative.37 The interpreta-

10.1021/jp907012w  2010 American Chemical Society Published on Web 01/07/2010

1764

J. Phys. Chem. C, Vol. 114, No. 4, 2010

Zartman et al.

Figure 1. (a) Models of the (water-free) surface of the aluminosilicate layers as a function of the average number of cations per nm2. Cation distributions are patchy, consistent with 29Si NMR data on the statistical distribution of SiO2 f AlO2- · · · K+ and AlO(OH) f MgO(OH)- · · · Na+ defects (refs 15 and 16). The surface area shown equals ∼2.60 × ∼2.70 nm2. (b) Side view of the corresponding 3D supercells of pyrophyllite (5 × 3 × 2), two montmorillonites of different CEC, and mica (5 × 3 × 1) under atmospheric pressure. The (ideal) cation exchange capacity (CEC in meq/100 g) is indicated. These model cells were employed in classical molecular dynamics simulation of the compressive elastic, shear, and failure properties.

TABLE 1: Equilibrium Cell Parameters in NPT Classical Molecular Dynamics Simulation and in Experiment at 298 K, as well as the Corresponding Equilibrium Basal Plane Spacing h between Individual Layers of the Nonhydrated Minerals cell dimension pyrophyllite montmorillonite CEC 91 meq/100 g montmorillonite CEC 143 meq/100 g mica

sim exptla sim exptlb sim exptlb sim exptlc

5 5 5 5 5 5 5 5

× × × × × × × ×

3 3 3 3 3 3 3 3

× × × × × × × ×

2 2 2 2 2 2 1 1

a (nm)

b (nm)

c (nm)

89.67

95-100

90.10

h (nm)

2.583 2.580 2.578 2.6 2.582 2.6 2.587 2.596

2.691 2.690 2.681 2.7 2.682 2.7 2.687 2.705

1.858 1.869 1.916 ∼2.0 1.856 ∼2.0 2.002 2.005

89.67 91.18 90.82 90 90.25 90 89.44 90

101.28 100.46 94.54 95-100 93.78 95-100 95.79 95.74

90.10 89.64 89.98 90 89.89 90 89.93 90

0.911 0.919 0.955d ∼1.0 0.926d ∼1.0 0.996d 0.997

a Ref 19. b Refs 20–23. The exact compositions of the samples is not known and data are therefore approximate. c Ref 24. d The thickness of individual, isolated layers is higher (∼1.3 nm) due to the presence of cations on both sides.

tion of available data for layered silicates suggests Young’s moduli in a range of 5 to 400 GPa,36–39 including prior simulation results reporting in-plane moduli as high as 20038 and 400 GPa.39 Verifiable experimental data on the macroscale include the 13 elastic constants of muscovite,26,28–31 phlogopite,32 and phengite micas33,34 at low stress, and the observation of increasing bulk moduli of muscovite mica,31 phengite,34 and pyrophyllite25,27,35 under increasing pressure. Elastic properties of smectite clays of intermediate CEC such as montmorillonites were difficult to obtain due to smaller grain size and ease of hydration;36,37 nevertheless, tentative bulk moduli between 6 and 12 GPa were reported. Amorphization of macroscopic mica samples was observed under 27 GPa triaxial pressure and above.30 The linear mechanical properties thus leave a number of open questions, and we attempt a clarification of the tensile properties, shear properties, and potential failure mechanisms as a function of the cation exchange capacity (CEC) and stress over a wide range. We rely on density functional theory (DFT) at the full

electronic structure level for small unit cells of pyrophyllite and mica, on classical molecular dynamics simulation with carefully derived force field models for all systems (Figure 1) as a function of applied stress,14,17,40 and on the interpretation of experimental data in the context of the simulation results. We will not consider bending deformations of single and multiple aluminosilicate layers in this contribution. A survey of the most important results for the linear mechanical properties at low stress is given in Table 2, which indicates a close correspondence between macroscale and nanoscale properties. Among the results, systematic overestimates of in-plane moduli (xx, yy, xy) in the classical force field model become apparent; however, the analysis of quantum-mechanical models, classical models, and experiment allows the formulation of a correction factor and the determination of the linear mechanical behavior for all layered silicate minerals in quantitative precision. The outline of the article is therefore as follows. We describe the molecular models, DFT methods, classical methods, and their limitations in section 2. Elastic constants, Young’s moduli,

Properties of Layered Silicates

J. Phys. Chem. C, Vol. 114, No. 4, 2010 1765

TABLE 2: Summary of Most Important Results: Young’s Moduli E, Bulk Moduli K, Shear Moduli G, and Failure Strength of Layered Silicates as a Function of CEC at Low Stress in GPaa Pyrophyllite (CEC 0) exptl Ex Ey Ez K Gxy Gxz Gyz

NA NA ∼49f 37h,I (1) NA 11h,l

comp

Montmorillonite (CEC 91) est

m

400 [3], 170 400 [10], 150m 38 (1), 70 (NA)m 34 (1), 39 (NA)m 190 [3] 3-6 [0.2]

160 [1] 160 [4] 38 (1) 37 (1) 71 [1] 3-6 [0.2]

exptl

comp

NA NA 6g (0.01) NA NA 4-6f

400 [8] 400 [10] 32 (1) 29 (1) 190 [3] 2-2.5 [0.25]

est 160 [3] 160 [4] 32 (1) 29 (1) 71 [1] 2-2.5 [0.25]

Montmorillonite (CEC 143) exptl NA NA ∼60f NA NA 4-6f

comp 400 [10] 400 [50] 60 (1) 43 (1) 190 [7] 3.6-5 [0.25]

est 160 [3] 160 [20] 60 (1) 43 (1) 71 [3] 3.6-5 [0.25]

Mica (CEC 251) exptl b,c,d,e

160 160b,c,d,e 59c,d (1) 59j,k (2) 71c,d 17c,d 15c,d

comp

est m

400 [15], 160 400 [50], 150m 60 (1), 70 (NA)m 54 (2) 190 [8] 20 [1] 14 [1]

160 [6] 160 [20] 59 (1) 59 (2) 71 [3] 17 [1] 15 [1]

a Experimental, computational, and estimated values are shown. Estimated values are derived from experiment, DFT, and classical MD simulation (including a correction factor 0.4 for in-plane elastic properties) for highest consistency. For stress-dependent values, the reference stress is shown in round brackets. Failure strengths are shown in square brackets where known. b Ref 26. c Ref 28. d Ref 29. e Ref 32. f Ref 37. g Ref 36. h Ref 27. I Ref 35. j Ref 30. k Ref 31. l Ref 25. m DFT (all other computational results force field-based).

and Poisson ratios derived from electronic structure calculations follow in section 3. The analysis of compressive tensile, shear, and failure properties as a function of the CEC and a broad range of stress by classical molecular dynamics simulation is presented in section 4. The article concludes with a summary in section 5. 2. Molecular Models and Simulation Details 2.1. Molecular Models. Models of the four layered silicates pyrophyllite, montmorillonite of CEC 91 and 143 meq/100 g, and mica cover the full range of cation density on the surface and were derived from X ray data (Figure 1)17 using the Materials Studio graphical interface.41 For the DFT calculation of elastic constants and in-plane elastic moduli, the smallest possible (1 × 1 × 1) unit cells of pyrophyllite and of mica were chosen for an affordable computational expense (Table 1). The mica unit cell contained four SiO2 f AlO2- · · · K+ defects in regular order. Montmorillonite models were not prepared for DFT calculations because adjustment of the cation density (CEC) requires larger models. For the classical molecular dynamics (MD) simulation of stress-strain relations, larger (5 × 3 × 2) and (5 × 3 × 1) periodic supercells of approximately 2.6 nm × 2.7 nm × 2 nm size were constructed (Table 1). Carefully justified atomic charges in agreement with an extended Born model were assigned to represent the polarity of the models,14 and the distribution of the SiO2 f AlO2- · · · K+ and AlO(OH) f MgO(OH)- · · · Na+ defects was chosen consistent with solidstate 29Si NMR data.15,16 2.2. DFT Simulation. The calculation of elastic constants and Young’s moduli at the electronic structure level was performed on models of the unit cell of pyrophyllite and mica using DFT. We employed the LDA approximation, PAW pseudopotentials, an energy cutoff at 600 eV, and plane waves (4 × 2 × 2 k points for the pyrophyllite unit cell of 40 atoms and 4 × 2 × 1 k points for the mica unit cell of 84 atoms) as implemented in the Vienna Ab-initio Simulation Package (VASP).42,43 First, 130 steps of geometry optimization of the respective unit cell were performed using the conjugate gradient method. The final energy (convergence 10-6 eV) and cell parameters were recorded. Second, a sequence of strains was applied to the optimized unit cell, represented by a strain tensor 5 S

5 S )

(

1 + δxx δxy δxz δxy 1 + δyy δyz δxz δyz 1 + δzz

)

(1)

b′,c which resulted in a series of strained unit cells (a b′,b b′)T ) T 5 b,c bi′,y bi′,z bi′)T S(a b,b b) with updated coordinates of all atoms i (x T 5 ) S(x bi,y bi,z bi) . For example, to evaluate the elastic constant C11, δxx was varied in a sequence of 10 strains (excluding zero) as -0.05, -0.04, ..., +0.05 while all other δij ) 0 (δxx ) 0 known from geometry optimization). To evaluate the elastic constant C22, δyy was varied in the same sequence while all other δij ) 0. To evaluate the elastic constant C12, δxx ) δyy were varied at the same time while all other δij ) 0. In general, to obtain an elastic constant Cijkl, the strains δij ) δkl were varied in the order -0.05, -0.04, ..., +0.05, resulting in several sets of 10 strained unit cells. Third, for each of these strained cells, 60 steps of geometry optimization were performed to allow the atoms to find their new equilibrium positions while the cell parameters (and thus the adjusted cell volume) remained fixed. The equilibrium energy of these strained unit cells was then recorded. Fourth, the computation of the elastic constants was carried out using the second derivatives of the energy for each set of strained cells as a function of strain. The energy is lowest for 5 S ) 0, and in each series of strain values the energy increases with increasing nonzero positive and negative strain. To obtain the elastic constants, a Taylor expansion of the energy E at the energy at zero strain E0 was considered assuming small strains:

E(ε11, ε12, ..., ε33) ) E0 +

1 1!

∂E εij + ∑ ∂ε ij ij

1 2!

E(ε11, ε12, ..., ε33) ) E0 +

V0 1!

2

∑ ∂ε∂ij∂εE kl εijεkl + ...

(2)

ij,kl

∑ σijεij + ij

V0 2!

∑ Cijklεijεkl + ...

(3)

ij,kl

In eqs 2 and 3, the first term is the equilibrium energy E0 at zero strain, the second term represents the contribution of bij (using σij ) -Fij/Aij, dxij ) x0,ijεij, a restoring force - b Fijdx V0 ) Aijx0,ij), and the third term characterizes the convex line shape represented by the elastic constants as a derivative of the stress with respect to deformation. The stress σij and the

1766

J. Phys. Chem. C, Vol. 114, No. 4, 2010

Zartman et al.

TABLE 3: Performance of the Phyllosilicate Force Field in PCFF17 and of CLAYFF49 for the Deformation of Micas in Comparison to Experimental Compression Data under Isotropic Stress;31,34 a Comparison of the Strains εxx, εyy, εzz, and of the Bulk Modulus K is Shown εxx

εyy

εzz

K (GPa)

-0.0258 (2) -0.030 -0.023

57 56 66

a

exptla sim - PFF-PCFFb sim - CLAYFFc

muscovite mica, 2 GPa -0.0046 (5) -0.0049 (5) -0.0028 -0.0031 -0.0036 -0.0041

exptld sim - PFF-PCFFb sim - CLAYFFc

phengite mica, 7.44 GPad -0.014 (2) -0.0186 (13) -0.056 (2) -0.0087 -0.0097 -0.0775 -0.0133 -0.0119 -0.057

a

86 79 98

Ref 31. b Ref 17. c Ref 49. d Ref 34.

second order adiabatic elastic constants Cijkl can also be written as:

σij )

1 ∂E 1 ∂2E and Cijkl ) V0 ∂εij V0 ∂εij∂εkl

(4)

At zero strain, the restoring force (the second term in eqs 2 and 3) is zero so that the nonzero second derivative (third term in eqs 2 and 3) can be directly analyzed. Computation of the second derivative of the energy at 5 S ) 0 then yields the elastic constants. For example, C11 was obtained from ∂2E/∂ε2xx using the series of δxx strain, C22 from ∂2E/∂ε2yy using the series of δyy strain, and C12 from the derivative ∂2E/∂ε2xx and the corresponding sum of elastic constants C11 + 2C12 + C22 using the series of δxx ) δyy strain. Moreover, Young’s moduli Exx, Eyy, and Ezz were calcu5 of 3 × 3 elastic constants and the lated using a reduced set C 5(εxx,εyy,εzz)T. For example, to obtain Exx equation (σxx,σyy,σzz)T ) C

Exx )

σxx εyy εzz ) C11 + C12 + C13 εxx εxx εxx

(5)

the set of 3 linear equations for a nonzero stress σxx in xx direction yields the two Poisson ratios νyx ) -εyy/εxx, νzx ) -εzz/εxx, and Exx in good approximation (neglecting tensile-shear coupling described by C14, C15, C16, which are close to zero). 2.3. Semiempirical Classical Simulation. Classical molecular dynamics simulation was carried out using the Phyllosilicate Force Field (PFF)17 embedded in the Polymer Consistent Force Field (PFF-PCFF). In contrast to other models such as CLAYFF with deviations up to 5% in cell parameters and up to 500% in surface energies in comparison to experiment,17 deviations in cell parameters using PFF-PCFF are 100 GPa), the interlayer space is compressed, cations (magnified for visibility) travel close to the silicate surface (mica, mont 143), and may enter the layers at low CEC (mont 91). In pyrophyllite, the void interlayer space is compressed. Under extreme triaxial tensile stress, similar observations were made, although cations did not enter into the layers for any CEC. (c) Under terminal xy shear stress, mica layers bend and fail (50 ps shown at 9 GPa). In contrast, all other systems show the onset of horizontal sliding similar to (a). (d) Above terminal yz and xz shear stress, shear flow in the lateral direction is observed in all minerals (60 ps shown for pyrophyllite under 0.5 GPa yz shear stress).

is observed for the clay minerals of higher CEC with interlayer space of highest cation density. At higher stress (3-15 GPa), pyrophyllite leads in bulk modulus due to the compression of the thinnest interlayer space and lack of cations, and mica is lowest due to the association of cations with negative charge defects (SiO2 f AlO2- · · · K+) directly located on the surface. At stresses above 20 GPa, mica displays the highest modulus and at extreme stress above 40 GPa, bulk moduli eventually become proportional to CEC in a narrow range (Figure 3). Experimental bulk moduli for mica at low stress were found to be 49-61 GPa depending on method and actual stress,28,30,31,33,34,37 including 61 GPa31 and 57 GPa31,34

at 2 GPa stress with an increment of 7-9 GPa per GPa additional stress. The simulation agrees quantitatively with 48 at 1 GPa stress and with 56 at 2 GPa stress. For montmorillonite, K values between 6 and 50 GPa are reported in experiment36,37 but the stress is not known. The simulation falls within this range with 43 GPa (143 meq/100 g) and 29 GPa (91 meq/100 g) at 1 GPa stress. Even lower values are possible around 20 GPa), however, true bulk moduli may decrease down to ∼40% of the computed value in Figure 3 due to the systematic overestimate of in-plane (xx and yy) stiffness (section 2.3). For pyrophyllite, a change in cell angle β from 102° to 112° was observed at ∼13 GPa stress (Figure 3 and Table 1); the modified angle remained stable up to 65 GPa. 4.3. Uniaxial Stress and Stress-Induced Failure. The application of uniaxial tensile stress leads to an expansion of the cell in the two directions not under compression and failure of the minerals under lower stress compared to triaxial compression can be observed (parts a and b of Figure 4). We find a lower resistance to failure in the xx and yy directions compared to the zz direction (Figure 5) where the cations may closely approach the surface or enter into the layer next to negative charge defects. In-plane elastic moduli (xx and yy) are lower than the slope under triaxial compression (Figure 2) and independent of CEC with a numerical value of 400 GPa. This value, after correction for overestimates of in-plane moduli and strengths (section 2.3), amounts to ∼160 GPa as on macroscopic mica26,28,29 and phlogopite32 samples and in agreement with DFT results (Table 4). Interesting qualitative trends include the similarity of the xx and yy modulus of all four clay minerals regardless of stress. Failure of the minerals was seen above 3 (xx) GPa/10 GPa (yy) for pyrophyllite, 7 GPa (xx)/10 GPa (yy) for montmorillonite of CEC 91, 10 GPa (xx)/50 GPa (yy) for

montmorillonite of CEC 143, and above 15 GPa (xx)/50 GPa (yy) for mica. The anticipated failure strength corresponds to 0.4 times these computed values (summary in Table 2). In contrast, Young’s modulus in the zz direction is highly dependent on stress and CEC. The deformation in the zz direction is significantly nonlinear, for example, the computed modulus for CEC 91 (part c of Figure 5) is ∼5 GPa near 10 MPa stress, ∼20 GPa near 100 MPa stress, and ∼32 GPa near 1 GPa stress due to compression of the alkali-containing interlayer. In the same curve, an irregularity εzz ) 0.20 is associated with a sudden approach of interlayer cations to the Si-surface level plane (CEC 91 meq/100 g), prior to entering the layers (part c of Figure 5). As a function of CEC, we observe values in the order 32 GPa (CEC 91), 38 GPa (CEC 0), 60 GPa (CEC 143), 60 GPa (CEC 251) at ∼1 GPa. The precise determination of experimental values at low stress is challenging except for mica; the closest reported data for the same minerals are ∼6 GPa (atomic force acoustic microscopy for a similar montmorillonite at ∼10 MPa),36 49 GPa,37 51-69 GPa (comparable smectites),37 and 59-61 GPa28,29 near 1 GPa. These trends are closely reflected in the simulation and reminiscent of bulk moduli as mostly the interlayer space in the zz direction is involved in the initial stages of compression. The smallest modulus of 5 to 32 GPa is associated with the wide, least-occupied interlayer space at low cation density in montmorillonite of CEC 91. In comparison, pyrophyllite possesses a smaller interlayer without cations and thus exhibits a higher zz modulus. 4.4. Shear Stress and Shear-Induced Failure. Responses to shear were analyzed in the range 0.0001 to 10 GPa in the xy,

Properties of Layered Silicates

J. Phys. Chem. C, Vol. 114, No. 4, 2010 1771

Figure 6. Shear stress-strain relationship for (a) pyrophyllite (CEC ) 0 meq/100 g), (b) montmorillonite (CEC ) 91 meq/100 g), (c) montmorillonite (CEC ) 143 meq/100 g), and (d) mica (CEC ) 251 meq/100 g). Failure can be associated with sliding of the layers rather than fracture (Figure 4). Note that the true in-plane stiffness (xy) only amounts to 0.4 times the computed value (limitations in section 2.3).

xz, and yz planes (Figure 6). All minerals are about an order of magnitude more resistant to shear in the xy plane compared to the xz and yz planes. The xy shear modulus was computed as ∼190 GPa independent of shear stress and of CEC of the mineral. The experimental xy shear modulus from macroscopic samples (mica) amounts to 71 GPa,28,29 which agrees with the computed value after correction (×0.4) of 76 GPa (section 2.3). The computed xy shear strength is 3, 3, 7, 8 GPa in the order of increasing CEC (Figure 6) and failure occurs by severe deformation (amorphization) (part c of Figure 4) or lateral sliding of the layers in the simulation (part d of Figure 4). The anticipated, true xy failure strength corresponds to 0.4 times the computed values (Table 1), that is, 1, 1, 3, 3 GPa in the order of increasing CEC. Approximate shear moduli in the xz and yz planes parallel to the aluminosilicate layers in the computation are 2 to 2.5 GPa for montmorillonite (CEC 91), 3 to 6 GPa for pyrophyllite, 3.6 to 5 GPa for montmorillonite (CEC 143), and 14 to 20 GPa for mica (Figure 6). Experimental data on macroscopic pyrophyllite samples are ∼11 GPa,25,27 4 to 6 GPa are suggested for montmorillonite,36 and 17 to 20 GPa28/13 to 15 GPa29 for mica. Computed data agree with experiment due to the reproduction of surface and interface properties in the model.17,40 Nevertheless, a lower shear modulus is computed for pyrophyllite in comparison to experiment, possibly related to imperfect crystallinity and orientation in the sample. The onset of sustained shear flow between the aluminosilicate layers was observed at stresses of 0.2 GPa, 0.2-0.3 GPa, 0.2-0.3 GPa, and 1 GPa in the order

of increasing CEC (part d of Figure 4) in the computation. Therefore, the comparatively low shear strength increases from pyrophyllite without interlayer cations over montmorillonites to a maximum for mica where all superficial cavities are filled with alkali cations (Figure 1). 5. Conclusions We presented a systematic investigation of the tensile and shear properties of mica-type silicates over the range of possible cation densities and a wide range of stress at the molecular scale. We quantified all major nanomechanical properties such as Young’s moduli in-plane and perpendicular, shear moduli inplane and perpendicular, as well as the strength and possible failure mechanisms that have been previously associated with uncertainties exceeding 1 order of magnitude or unknown (Table 2 and Figure 4). The analysis has been enabled by the use of hierarchical simulation techniques at the electronic structure (DFT) level and at the classical semiempirical (MD) level in combination with available experimental data. DFT calculations lead to accurate estimates of the full set of elastic constants for pyrophyllite and mica at zero strain, whereas the classical force field-based simulation describes dynamical mechanical properties for significantly larger model structures as a function of interlayer composition and stress. In the forcefield-based simulation, however, corrections for a systematic overestimate of the in-plane stiffness were necessary by comparison to experimental data and electronic structure results.

1772

J. Phys. Chem. C, Vol. 114, No. 4, 2010

Layered silicates of any cation density exhibit a large mechanical anisotropy, which is characterized by in-plane moduli of ∼160 GPa independent of CEC and stress and by smaller perpendicular moduli between 30 and 60 GPa upon initial compression. At high stress >10 GPa, however, the perpendicular moduli converge to a similar value as the in-plane modulus. Novel observations are also potential failure mechanisms such as kinks in the aluminosilicate layers, cation intrusion, and shear flow. In particular, lateral shear flow between the aluminosilicate layers can be a likely failure mechanism in aggregates of layered silicates before the internal, chemically bonded structure is modified. The onset of shear flow between the clay mineral layers requires a lateral shear stress between 0.2 and 1.0 GPa (in the xz or yz direction) with increasing CEC and then leads to a continuous horizontal offset of the individual clay mineral layers mediated by the interlayer cations or an empty interlayer space (in case of pyrophyllite). A nearly quantitative correlation between the macroscale and nanoscale mechanical properties is observed with a trend toward slightly increased tensile strength on the smaller scale. The insight into mechanical data of clay minerals aids in understanding their interfaces in polymeric and biocompatible composites, thin films, and natural environments. The results may also be utilized to parametrize micromechanical models of composites at coarse-grain and finite element levels. Future challenges include the 3D visualization of local stress tensors,58 the understanding of bending energies and bending stability of layered minerals.59 DFT calculations may be applied to increasingly larger systems of related materials as they demonstrate consistency with experiment. Acknowledgment. We are grateful for support by the University of Akron, the Air Force Research Laboratory, and the allocation of computational resources at the Ohio Supercomputing Center. G.D.Z. acknowledges a Goodyear fellowship and the NSF-REU program. References and Notes (1) Hydrous Phyllosilicates. In ReViews in Mineralogy; Bailey, S. W., ed.; Mineralogical Society of America: Chelsea, MI, 1988; Vol. 19. (2) Ishida, H.; Campbell, S.; Blackwell, J. Chem. Mater. 2000, 12, 1260–1267. (3) Yalcin, B.; Cakmak, M. Polymer 2004, 45, 6623–6638. (4) Osman, M. A.; Mittal, V.; Morbidelli, M.; Suter, U. W. Macromolecules 2004, 37, 7250–7257. (5) Drummy, L. F.; Koerner, H.; Farmer, K.; Tan, A.; Farmer, B. L.; Vaia, R. A. J. Phys. Chem. B 2005, 109, 17868–17878. (6) Osman, M. A.; Rupp, J. E. P.; Suter, U. W. Polymer 2005, 46, 1653–1660. (7) Heinz, H.; Vaia, R. A.; Farmer, B. L. Langmuir 2008, 24, 3727– 3733. (8) Wang, C. C.; Juang, L. C.; Lee, C. K.; Lee, J. F.; Huang, F. C. J. Colloid Interface Sci. 2004, 273, 80–86. (9) (a) Ruiz-Hitzky, E. Chem. Record 2003, 3, 88–100. (b) Salles, F.; Devautour-Vinot, S.; Bildstein, O.; Jullien, M.; Maurin, G.; Giuntini, J. C.; Douillard, J. M.; Van Damme, H. J. Phys. Chem. C 2007, 112, 14001– 14009. (10) Roiter, Y.; Ornatska, M.; Rammohan, A. R.; Balakrishnan, J.; Heine, D. R.; Minko, S. Nano Lett. 2008, 8, 941–944. (11) Stepanenko, V.; Stocker, M.; Muller, P.; Buchner, M.; Wurthner, F. J. Mater. Chem. 2009, 19, 6816–6826. (12) Podsiadlo, P.; Michel, M.; Critchley, K.; Srivastava, S.; Qin, M.; Lee, J. W.; Verploegen, E.; Hart, A. J.; Qi, Y.; Kotov, N. A. Angew. Chem., Int. Ed. 2009, 48, 7073–7077. (13) Sun, X. P.; Ko, S. H.; Zhang, C. A.; Ribbe, A. E.; Mao, C. D. J. Am. Chem. Soc. 2009, 131, 13248–3733. (14) Heinz, H.; Suter, U. W. J. Phys. Chem. B 2004, 108, 18341–18352. (15) Heinz, H.; Suter, U. W. Angew. Chem., Int. Ed. 2004, 43, 2239– 2243. (16) Herrero, C. P.; Sanz, J. J. Phys. Chem. Solids 1991, 52, 1129– 1135.

Zartman et al. (17) Heinz, H.; Koerner, H.; Anderson, K. L.; Vaia, R. A.; Farmer, B. L. Chem. Mater. 2005, 17, 5658–5669. (18) Heinz, H.; Vaia, R. A.; Krishnamoorti, R.; Farmer, B. L. Chem. Mater. 2007, 19, 58–69. (19) Rothbauer, R. Neues Jahrb. Mineral., Monatsh. 1971, 143–154. (20) Brown, G. The X-ray Identification and Crystal Structures of Clay Minerals; Mineralogical Society: London, 1961. (21) ReViews in Mineralogy; Bailey, S. W., Ed.; Mineralogical Society of America: Chelsea, MI, 1988; Vol. 19 See also http://www.webmineral. com. (22) Tsipurski, S. I.; Drits, V. A. Clay Mineral. 1984, 19, 177–193, Mg positions not correct. (23) The exact crystal structure of montmorillonite depends on the nature of the cations (Na+, K+, Ca2+), charge density (in ref 22 K0.58), and presence of crystal water. Mainly the parameters c (roughly 9.9-10.3 Å) and β (roughly 95-100°) are affected. (24) Lee, J. H.; Guggenheim, S. Am. Mineral. 1981, 66, 350–357. (25) Giardini, A. A. J. Bas. Eng. 1967, 89, 554–560. (26) Simmons, G.; Wang, H. Single Crystal Elastic Constants and Calculated Aggregate Properties: A Handbook, 2nd ed.; MIT Press: Cambridge, MA, 1971. (27) Sachse, W.; Ruoff, A. L. J. Appl. Phys. 1975, 46, 3725–3730. (28) Vaughan, M. T.; Guggenheim, S. J. Geophys. Res. 1986, 91, 4657– 4664. (29) McNeil, L. E.; Grimsditch, M. J. Phys.: Condens. Matter 1993, 5, 1681–1690. (30) Faust, J.; Knitte, E. J. Geophys. Res. 1994, 99, 19785–19792. (31) Catti, M.; Ferraris, G.; Hull, S.; Pavese, A. Eur. J. Mineral. 1994, 6, 171–178. (32) Habelitz, S.; Carl, G.; Russel, C.; Theil, S.; Gerth, U.; Schnapp, J. D.; Jordanov, A.; Knake, H. J. Non-Cryst. Solids 1997, 220, 291–298. (33) Pavese, A.; Ferraris, G.; Pischedda, V.; Mezouar, M. Phys. Chem. Miner. 1999, 26, 460–467. (34) Smyth, J. R.; Jacobsen, S. D.; Swope, R. J.; Angel, R. J.; Arlt, T.; Domanik, K.; Holloway, J. R. Eur. J. Mineral. 2000, 12, 955–963. (35) Pawley, A. R.; Clark, S. M.; Chinnery, N. J. Am. Mineral. 2002, 87, 1172–1182. (36) Vanorio, T.; Prasad, M.; Nur, A. Geophys. J. Int. 2003, 155, 319– 326. (37) Chen, B.; Evans, J. R. G. Scripta Mater. 2006, 54, 1581–1585. (38) Collins, D. R.; Catlow, C. R. A. Am. Mineral. 1992, 77, 1172– 1181. (39) Manevitch, O. L.; Rutledge, G. C. J. Phys. Chem. 2004, 108, 1428– 1435. (40) Heinz, H.; Vaia, R. A.; Farmer, B. L. J. Chem. Phys. 2006, 124, 224713. (41) Materials Studio and DiscoVer program; Accelrys, Inc. 2006. (42) Kresse, G.; Hafner, J. Phys. ReV. B 1993, 47, 558–561. (43) Kresse, G.; Joubert, D. Phys. ReV. B 1999, 59, 1758–1775. (44) Cygan, R. T.; Greathouse, J. A.; Heinz, H.; Kalinichev, A. J. Mater. Chem. 2009, 19, 2470–2481. (45) Heinz, H.; Vaia, R. A.; Koerner, H.; Farmer, B. L. Chem. Mater. 2008, 20, 6444–6456. (46) The same high in-plane Young’s modulus is also found when scaling equilibrium bond lengths by a factor of 1.05 (ref 17) to counteract Coulomb attraction is not employed (equal to a factor of 1.00). (47) Frankland, S. J. V.; Harik, V. M.; Odegard, G. M.; Brenner, D. W.; Gates, T. S. J. Compos. Sci. Tech. 2003, 63, 1655–1661. (48) Heinz, H. Mol. Sim. 2007, 33, 748–758. (49) Cygan, R. T.; Liang, J. J.; Kalinichev, A. G. J. Phys. Chem. B 2004, 108, 1255–1266. (50) Hill, J. R.; Sauer, J. J. Phys. Chem. 1995, 99, 9536–9550. (51) Teppen, B. J.; Rasmussen, K.; Bertsch, P. M.; Miller, D. M.; Schafer, L. J. Phys. Chem. B 1997, 101, 1579–1587. (52) Greathouse, J. A.; Refson, K.; Sposito, G. J. Am. Chem. Soc. 2000, 122, 11459–11464. (53) Pospisil, M.; Capkova, P.; Merinska, D.; Malac, Z.; Simonik, J. J. Colloid Interface Sci. 2001, 236, 127–131. (54) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182–7190. (55) CASTEP program; Accelrys, Inc., 2006. (56) Craig, M. S.; Warren, M. C.; Dove, M. T.; Gale, J. D.; SanchezPortal, D.; Ordejon, P.; Soler, J. M.; Artacho, E. Phys. Chem. Miner. 2004, 31, 12–21. (57) Below 1 GPa, the change in cell volume is very small and falls within the error in mean volume measurements. Numerical errors in computed bulk moduli are