Molecular Dynamics Simulation of Nanocrack ... - ACS Publications

Article Cite This: J. Phys. Chem. C 2018, 122, 1351−1360

pubs.acs.org/JPCC

Molecular Dynamics Simulation of Nanocrack Propagation in Single-Layer MoS2 Nanosheets Hongwei Bao,† Yuhong Huang,‡ Zhi Yang,† Yunjin Sun,§ Yu Bai,∥ Yaping Miao,†,⊥ Paul K. Chu,⊥ Kewei Xu,†,# and Fei Ma*,†,⊥ †

State Key Laboratory for Mechanical Behavior of Materials, Xi’an Jiaotong University, Xi’an 710049, Shaanxi, China College of Physics and Information Technology, Shaanxi Normal University, Xi’an 710062, Shaanxi, China § Beijing Laboratory of Food Quality and Safety, Faculty of Food Science and Engineering, Beijing University of Agriculture, Beijing 102206, China ∥ Suzhou Research Institute, Xi’an Jiaotong University, Suzhou 215123, Jiangsu, China ⊥ Department of Physics and Materials Science, City University of Hong Kong, Tat Chee Avenue, Kowloon, Hong Kong, China # Department of Physics and Opt-electronic Engineering, Xi’an University of Arts and Science, Xi’an 710065, Shaanxi, China ‡

S Supporting Information *

ABSTRACT: Single-layer MoS2 (SLMoS2) nanosheets promise potential applications in flexible electronic and optoelectronic nanodevices for which the mechanical stability is crucial. However, the measured fracture strength is extremely dispersive, which might be due to the random crack configuration. In this work, molecular dynamics (MD) simulations are conducted to investigate the propagation of nanocracks in SLMoS2 nanosheets and the fracture mechanism at atomic scale, and the modified Griffith criterion developed by Yin et al. is adopted to fit the dependence of fracture stress on the initial crack length. Moreover, the fracture stress is highly dependent on the initial crack configuration, crack length, and crack angle. The energy release rate (GS) decreases with increasing initial crack length, crack angle, and temperature but is not sensitive to strain rate. The average propagation velocity of cracks (V̅ ) is substantially reduced with increasing initial crack length and crack angle but is almost independent of temperature and strain rate. The V̅ at lower GS is well predicted by linear elastodynamic theory but approaches 66% Rayleigh-wave speed at a higher GS of >5.78 J/m2. It is also found that fracture is preferred along the zigzag direction of SLMoS2 nanosheets. The results provide us a clear understanding on the dispersive data of measured fracture strength of SLMoS2 nanosheets. 270 ± 100 GPa was obtained.12 However, it is extremely challenging to study the deformation and fracture behaviors in 2D materials by experiments. Na et al.13 conducted in situ tensile tests on polycrystalline graphene grown on copper foil by using a scanning electron microscope. It was found that the cracks were initiated perpendicular to the loading direction throughout the graphene sheets, but wrinkles parallel to the loading direction emerged. Zhang et al.14 conducted an in situ tensile test of suspended graphene with a predefined central crack and found that brittle fracture occurred with substantially reduced fracture strength. Recently, Yang et al.15 developed a new nanomechanical device to quantitatively measure the mechanical properties of freestanding MoSe2 nanosheets through in situ tensile loading, and they found that MoSe2 was more brittle than graphene. The measured fracture strength varies from 2.2 to 9.9 GPa,

1. INTRODUCTION Since graphene was successfully fabricated, two-dimensional (2D) materials became one of the research focuses in the fields of material science, physics, and chemistry in recent years.1−3 In principle, any bulk materials with a laminated structure and bonded by van der Waals force can be exfoliated into graphenelike layers, including transition metal dichalcogenides (TMDs),4 transition metal oxides,5 boron nitride,6 black phosphorus,7 and even artificially produced honeycomb structures.8 There are more than 40 compounds in the family of TMDs, and their physical properties vary from insulating to superconducting.9 Single-layer MoS2 (SLMoS2), as a prototypical example of TMDs, has a direct band gap of 1.8 eV and promising potential in nextgeneration field-effect transistors, optoelectronics, energy harvesting devices, and flexible electronics nanodevices.10,11 The stability and reliability of these devices depend critically on the mechanical properties of SLMoS2. Atomic force microscopy (AFM) based nanoindentation has been commonly used to measure the effective elastic modulus of SLMoS2, and a value of © 2017 American Chemical Society

Received: October 11, 2017 Revised: December 21, 2017 Published: December 21, 2017 1351

DOI: 10.1021/acs.jpcc.7b10094 J. Phys. Chem. C 2018, 122, 1351−1360

Article

The Journal of Physical Chemistry C

Figure 1. (a) Typical geometrical configurations of the SLMoS2 nanosheet with edge or center cracks under mode I tensile loading. The insets show the atomic structures of AC-cen, AC-2S, AC-1Mo, and ZZ-asy cracks. (b) Top view and (c) side view of an atomic model of the SLMoS2 nanosheet with edge cracks under tearing loading (mode III).

been well understood at atomic scale. In this work, besides the above-mentioned research contents, the dependences of fracture strength on the temperature and strain rate have been systematically explored by MD simulations, and the results are compared with continuum theory. It is found that the fracture stress is highly dependent on the initial crack configuration, crack length, and crack angle, and the fracture is preferred along the zigzag direction of SLMoS2 nanosheets. More importantly, the relationship between energy release rate and crack velocity at the atomic scale is established.

dependent on the samples, that is, the experimental data are extremely dispersive. Randomly distributed defects, grain boundaries (GBs), and cracks in single-layer sheets might be the reasons;16−19 however, it is still difficult to quantitatively measure and evaluate the influences on the deformation and fracture behaviors from the experimental aspect.20−22 In comparison, molecular dynamics (MD) simulations have been widely adopted to directly predict the mechanical properties as well as the atomic evolution during mechanical loading.23,24 Based on MD simulations, Jiang et al.,25 Xiong et al.,26 and Zhao et al.27 found that the elastic modulus of SLMoS2 is in the range of 110−220 GPa, which is close to the measured value of 270 ± 100 GPa. Yin et al.28 found that the Griffith criterion is still valid in describing the crack propagation in graphene and proposed a geometrical factor to modify the Griffith criterion. Liu et al.29 verified the modified Griffith criterion in single-layer black phosphorus (SLBP) and found that the geometrical factor α is almost a constant. Budarapu et al.30 studied the initiation and growth of cracks in graphene sheets and discussed the physical mechanism, focusing on the potential energy, stress, bond length, and reorientation at the tip of cracks. Abadi et al.31 studied the temperature and grain size dependences of mechanical properties of polycrystalline hexagonal boronnitride (h-BN) nanosheets with precracks; the influences of the crack size on the crack propagation speed and orientation were also discussed. In the work of Xiong et al.,32 the uniaxial tensile loading of precracked SLMoS2 sheets with different crack tips, locations, lengths, and angles was simulated. They found that the fracture behaviors of monolayer MoS2 sheets are significantly dependent on the crack tip configuration, but little on the crack location. Wang et al.33 mainly focused on the fracture toughness and crack propagation path in SLMoS2 sheets in a mixed loading mode, and obtained the fracture toughness in the range of 1−1.8 MPa·m1/2. However, the dynamic characteristics for the crack propagation in the natural brittle 2D materials have not

2. SIMULATION MODEL AND METHOD In SLMoS2, the Mo monolayer is sandwiched by two S monolayers, and the Mo and S atoms are covalently bonded with each other. Based on a lattice constant of 2H-MoS2 (a = 3.16 Å), the initial model of the SLMoS2 sheet is established.10,11 The cracks with different sizes and configurations are predefined by deleting some atoms. Figure 1a shows the model of 50 nm × 50 nm SLMoS2 nanosheet, and the atomic configurations at crack tips of given edges are displayed in the insets. Type 1 (AC-cen) is a crack at the center of the SLMoS2 nanosheet with armchair edges, of which two atomic configurations at the tips are involved: in one case, two S atoms are included, and in the other case, one Mo atom is included. Types 2 and 3 (AC-2S or AC-1Mo) are cracks of armchair edges at the boundaries of the SLMoS2 nanosheet in which two S atoms and one Mo atom are involved at the crack tips, respectively. Type 4 (ZZ-asy) is a crack of asymmetrical zigzag edges at the boundaries of SLMoS2 nanosheet. Tensile loading is applied perpendicular to the cracks, that is, along the y axis, which is denoted as mode I. The initial length of the cracks changes from 3 to 20 nm. In addition, some other models including angle θ (0°, 30°, 60°) between the crack and tensile loading direction are also simulated for comparison. In order to study the in-plane fracture during mechanical exfoliation,34,35 taking AC-2S and ZZ-asy cracks as examples, 1352


Article

The Journal of Physical Chemistry C the tearing loading (mode III) is also simulated, as Figure 1b,c shows. LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) is used to conduct the MD simulation.36 A new Stillinger−Webber (SW) potential developed by Jiang et al.25,37 is adopted to describe the interatomic interactions, bond breaking, and bond reforming in SLMoS2. An additional cutoff for the three-body interaction is embedded, and bond stretching and angle rotation are described as V2 = A e[ρ / r − rmax](B /r 4 − 1)

(1)

V3 = A e[ρ1 /(r12 − rmax12) + ρ2 /(r13− rmax13)](cos θ − cos θ0)2

(2)

The simulated elastic properties, deformation behavior, thermal conductivity, and phonon spectrum based on the SW potential show good agreement with the DFT and experimental results.38,39 The initial SLMoS2 nanosheets are relaxed in the NPT (number of atoms N, pressure P, and temperature T are constant) ensemble at T = 10 K and P = 0 Pa for 2 ns, and the nanosheets are thermally equilibrated using the Nosé−Hoover thermostat. The standard Newton equations of atomic motion are integrated using the velocity Verlet algorithm with a time step of 1 fs, and the uniaxial tensile loading is applied by uniform displacement of lattice in the nanosheets. The periodic boundary conditions (PBCs) are applied in all directions, and a vacuum 30 Å in thickness is added along all directions to avoid the interactions of atoms near the opposite edges. The two edges parallel to the loading direction are free to relax, but the other two edges perpendicular to the loading direction are clamped. A strain rate of 108 S−1 and a temperature of 10 K are applied. For comparison, the temperature in the range of 10−400 K and the strain rate in the range of 5 × 107−109 S−1 are also considered. For the tearing loading, as schematically shown in Figure 1b, in the 50 nm × 50 nm SLMoS2 nanosheet, the region A with 0.4 nm in width is fixed, but the regions B and C are dragged along the opposite out-of-plane directions at the velocities of 0.5 nm and −0.5 nm·ps−1, respectively, as indicated by the arrows in Figure 1c. A vacuum region with 20 nm in width is added in the out-of-plane direction to avoid the nonphysical atomic interaction near the opposite edge, and the temperature is set as 10 K. The atomic stress is calculated using the virial definition40 σij =

⎛ 1 ⎜1 V ⎝⎜ 2

N

∑ α=1

⎛ Δx αβ Δx αβ i j U ′ − ∑ ⎜⎜γ αβ αβ γ ⎝ β≠α N

⎞⎞ α β ⎟⎟ m x x ∑ α i̇ j̇ ⎟⎟ ⎠⎠ α=1 N

Figure 2. (a) Stress−strain curves of SLMoS2 nanosheet with AC-cen, AC-2S, AC-1Mo, and ZZ-asy cracks of different initial length. (b) Fracture stress as a function of initial crack length. (c) Fracture stress as a function of crack angle; the inset shows the fracture stress as a function of initial crack length for given crack angles.

(3)

in which V is the total volume of the SLMoS2 nanosheets with a thickness of 0.61 nm, N is the total number of atoms, ẋαi is the ith α β component of the velocity of atom α, Δxαβ j = ẋj − ẋj , mα is the αβ mass of atom α, γ is the distance between atoms α and β, and U′ is the potential energy function. The atomic structures at different times are visualized by the OVITIO package.41

51.7%, and from 6.65 to 2.99 GPa by 55.03% for AC-cen, AC-2S, AC-1Mo, and ZZ-asy typed cracks, respectively. According to the modified Griffith criterion,28 the fracture stress is given as

3. RESULTS AND DISCUSSION Figure 2a plots the stress−strain curves of SLMoS2 nanosheets with cracks of AC-cen, AC-2S, AC-1Mo, and ZZ-asy types. After elastic deformation, the stress decreases sharply at a critical strain, characteristic of brittle fracture. The fracture stress is plotted in Figure 2b. As the initial length of cracks is increased from 3 to 20 nm, the fracture stress decreases from 8.23 to 3.3 GPa by 59.9%, from 6.72 to 3.18 GPa by 52.6%, from 7.26 to 3.5 GPa by

σf =

α F(ϕ)

2E Γ πa

(4) 12,25−27

in which E is Young’s modulus of 270 GPa, a is the initial length of cracks, and Γ is the apparent fracture resistance of the crack plane (it is the surface energy for bulk materials and the edge energy for 2D materials). Based on the DFT calculations, the values of 1.94 and 1.84 J/m2 are chosen for armchair and 1353


Article

The Journal of Physical Chemistry C zigzag edges, respectively,42 α is geometrical correction, and F(ϕ) is a geometrical factor given by

Griffith criterion remains valid even in a large crack length range from hundreds nanometers (experimental) to several nanometers (MD simulations). However, commonly, the cracks of different sizes are randomly distributed in CVD-grown 2D samples, and they will interact with each other, which makes the fracture complicated. In the experiments, it is difficult to directly visualize the atomic structure at the crack tip. Therefore, it is mandatory to develop experimental technologies and large-scale MD simulations to discover the complicated fracture mechanism in the future. In fracture mechanics, the energy release rate (GS) is usually adopted to evaluate the brittleness of material.44 Accordingly, during the crack propagation, the energy stored at the crack tip is partly converted into edge energy along the cracks and partly dissipated as atomic motion. The external loading promotes the crack extension, corresponding to the release of mechanical energy per unit crack, GS,45

⎛ πϕ ⎞ W F(ϕ) = (1 − 0.025ϕ2 + 0.06ϕ4) sec⎜ ⎟ , and ϕ = ⎝ 2 ⎠ a (5)

in which W is the width of the crack. The solid blue lines in Figure 2b show the fitted results, and they match well with the MD results. The geometrical corrections, α, are around 1.225 and 1.343, which are much lower than that of graphene28 but much larger than that of SLBP.29 The modified Griffith criterion can be used to predict the fracture stress of SLMoS2 nanosheets with a crack length even down to 3 nm. Taking the AC-cen crack as an example, the crack propagation in SLMoS2 nanosheets with different crack angles and initial crack length is investigated. Pure tensile loading takes place at the crack angles θ = 0° and 90°, while mixed loading occurs at the crack angles θ = 30° and 60°. The crack is parallel to the tensile loading direction at the crack angle of 0°, but perpendicular to the tensile loading direction at the crack angle of 90°. As shown in Figure 2c, the fracture stress decreases by 76.5−353.6% for the SLMoS2 nanosheets, dependent on the crack length, as the crack angle is increased from 0° to 90°. As evidenced from the inset in Figure 2c, the modified Griffith criterion is still valid if tensile loaded along a direction of 30°, 60°, and 90° with respect to the cracks. Apparently, the fracture strength is highly dependent on the initial crack length, crack configuration, and crack angle, so that the measured fracture strength of CVD-grown SLMoS2 will be extremely dispersive. The available results of experimental tests and MD simulations on graphene nanosheets,14,28,43 SLBP nanosheets,29 and SL/BLMoSe215 are summarized in Figure 3.

GS = −

W − Wa +Δa dW = a 2t da 2t Δa

(6)

in which t is the effective thickness taking the value of 0.61 nm, and Wa and Wa+Δa are the total external work required for an initial crack length a and a + Δa, respectively, and can be calculated from the stress−strain curves. Figure S1 shows the crack length as a function of strain. As shown in Figure 4a, the GS decreases from 6.17 to 2.06 J/m2 by 67.7%, from 3.36 to 1.02 J/m2 by 69.6%, from 3.45 to 1.71 J/m2 by 50.4%, and from 2.71 to 1.16 J/m2 by 57.2% for AC-cen, AC-2S, AC-1Mo, and ZZ-asy cracks, respectively, as the initial crack length is increased from 3 to 20 nm. The GS of center cracks is always larger than that of boundary cracks by 200%. The average GS of SLMoS2 is much smaller than that of graphene (15.9 J/m2),14 indicating that SLMoS2 is more brittle than graphene. Taking the center cracks as examples, Figure 4b,c shows that the GS decreases with increasing crack angle and elevating temperature substantially. As shown in Figure 4d, the GS is not sensitive to strain rate, and the strain rate effect is negligible. In fact, GS is not only a function of the loading force required to create fresh fracture surface but also a function of the crack itself, and GS is given by46 GS =

σ 2W (1 − ν 2) 2E

(7)

in which σ is the critical driving force, that is, the fracture stress, W is the crack width, the elastic modulus E is 270 GPa, and the Poisson’s ratio ν is 0.267. Commonly, the fracture stress decreases with increasing initial crack length, crack angle, and temperature, and GS decreases accordingly. It was reported that the strain rate has a little influence on the fracture stress at low temperature.47 Therefore, GS is not sensitive to strain rate, especially at low temperature (10 K), as reported recently on SLBP in the strain rate range of 5 × 107−109 S−1.29 According to the linear elastodynamic theory of fracture, the crack propagation velocity should match the energy release rate. The average velocity V̅ can be calculated as V̅ = d(a − a0)/dt. As shown in Figure S1, the crack length increases with the strain linearly. As shown in Figure 5, the value of V̅ in SLMoS2 is in the range of 0.55−3.13 km/s, which is smaller than that of graphene (7.13−8.82 km/s).48 According to Yoffe’s model, the propagation speed of cracks under mode I loading cannot exceed the Rayleigh wave speed, that is, the velocity of surface acoustic waves.49 Based on the linear elastodynamic theory, the Rayleigh

Figure 3. Fracture stress of different 2D materials as a function of initial crack length. The results of graphene nanosheets (Zhang et al.,14 Yin et al.,28 and Liu et al.43), SLBP nanosheets (Liu et al.29), and SL/BLMoSe2 nanosheets (Yang et al.15) reported previously are also included for comparison.

Although the initial crack length, temperature, and strain rate in experiments and MD simulations are different, both illustrate the applicability of the Griffith criterion in describing the brittle fracture of 2D materials.14 In the work of Yang et al.,15 the crack length in SL/BLMoSe2 is simply estimated in the range of 77.5−3.6 nm by the Griffith criterion in which only the dominant fracture-producing crack is considered. In such a case, the initial crack length estimated by the Griffith criterion is well matched with the MD results, as shown in Figure 3. Furthermore, the

wave velocity, CR, can be calculated as C R ≈ 0.9265 1354

E 2ρ(1 + ν)

in


Article


Figure 4. Average crack velocities as a function of (a) initial crack length, (b) crack angle, (c) temperature, and (d) strain rate.

Figure 5. Energy release rate as a function of (a) initial crack length, (b) crack angle, (c) temperature, and (d) strain rate.

which the elastic modulus E of 270 ± 100 GPa, the Poisson’s ratio ν of 0.267, and the density ρ of 4.80 are adopted.50,51 The calculated CR of the SLMoS2 nanosheet is about 5.11 km/s. In this work, the Yoffe’s model is still valid, and V̅ < 0.66CR. A model based on continuum fracture mechanics is given to describe the crack propagation velocity in SLMoS2 under tensile loading:52

a ⎞ ⎛ V0 = 0.66κC R ⎜1 − 0 ⎟ ⎝ a⎠

(8)

in which V0 is the instantaneous velocity, κ is a constant, CR is the Rayleigh wave velocity, a0 is the initial crack length, and a is the current crack length. Figure S2 shows the calculated velocity with CR of 5.11 km/s and κ of 1. Apparently, the crack propagation 1355


Article


the crack velocity at higher energy release rate evaluated from MD simulations substantially deviates what is predicted by the continuum theory. Similar results were found in the other brittle materials, such as poly(methyl methacrylate) (PMMA)56 and single-crystal silicon.57 Experimentally, crack propagation can be measured by means of high-speed photography or optical methods of photoelasticity and caustics.58,59 In principle, crack velocity can be measured according to the crack propagation at a given time, and the energy release rate be evaluated from the contours and the stress field at the crack tips. Based on this, the crack propagation behaviors in the samples of millimeters in size have been widely explored.58,59 However, it is still intractable to measure the energy release rate and crack velocity of nanocracks in 2D materials. Hence, MD simulation is an effective approach to provide an insight into the crack propagation at nanoscale. Taking SLMoS2 with an initial crack length of 5 nm as an example, Figure 7 shows the atomic configuration evolution and crack propagation path for AC-cen, AC-2S, AC-1Mo, and ZZ-asy type cracks at 10 K in which tensile loading is applied perpendicular to the cracks, that is, in mode I. The distribution of stress per atom along the y axis, σy, is displayed, and the atomic configurations near the crack tip are plotted in the insets. As the tensile strain is increased, stress concentration takes place around the crack tip, and the crack tip becomes blunted; as a result, the S−Mo bonds near the crack tip break first, and then the crack propagates far away. Remarkably, crack kinks are evidenced at the early stage of crack propagation for the armchair-edged cracks (AC-cen, AC-2S, and AC-1Mo) but not for the zigzag-edged cracks (ZZ-asy). More interestingly, it can be seen that the newly formed cracks will deviate from initial armchair orientation to the zigzag one, forming a zigzaged crack path with only a few nanometers, as shown in Figure 7a−c. The angles between the different parts of kinked cracks are close to 120°, as shown in Figure S4. So the kinked cracks are mainly composed of zigzag edges, as observed in experiments.60,61 Fracture-induced edges in SLMoS2 nanosheets with ZZ-asy cracks are atomically flat. The difference in morphology depends on the atomic configurations around the crack tips. For a given initial crack length, the fracture strength of SLMoS2 nanosheet with a center crack is much higher than that with a crack at the edges, and the fracture strength of nanosheets with an armchair-edged crack (AC-cen, AC-2S, and AC-1Mo) is higher than that with zigzag-edged crack (ZZ-asy). This can be attributed to the distinct bond arrangement at the crack tip. If σi is used to denote the breaking strength of an individual S−Mo bond, the stress to break the S−Mo bond at the crack tip, σa will be σi σa = sin(90° − φ) (11)

velocity strongly depends on the initial crack length, that is, the shorter the initial crack, the larger the crack velocity (Figure 5a). If there is an inclusive angle between the crack and the loading direction, both tensile and shearing loading modes will be involved, resulting in larger stress concentration and higher crack propagation velocity (Figure 5b). However, according to eq 8, the crack velocity is independent of the temperature and strain rate (Figure 5c,d). Generally speaking, the cracks propagate along with the stress wave propagation (SWP) under mode I loading. As illustrated in Figure S3, the stress concentration around the crack tip decays in an approximately semicircular wavefront. This is similar to laser-induced surface acoustic waves in single MoS2 flake and thin films.53,54 For an approximately straight crack extension, the dynamic fracture energy Γ can be expressed as55 ⎛ V̅ ⎞ Γ = GS⎜1 − ⎟ CR ⎠ ⎝

(9)

and for 2D materials, Γ is a constant and equal to the edge energy for the simplest case. Based on the DFT calculations, the value of 1.89 J/m2 is chosen.46 Equation 9 can be rewritten as V̅ = 5.11 −

9.89 GS

(10)

Accordingly, the average velocity can be fitted as a function of energy release rate; the result is plotted as a black solid line in Figure 6 in which CR (5.11 km/s) and 0.66CR (3.37 km/s) are

Figure 6. Average crack velocity as a function of energy release rate. The symbols show the MD results, and the dark solid line shows the fitted curve.

highlighted by the red and blue solid lines, respectively. The black solid line and the blue solid line intersect with each other at the critical GS value of 5.78 J/m2. If GS is smaller than that value, the average velocity is well predicted by the continuum theory. However, if GS is larger than that value, the crack velocity from MD simulation is much lower than what is predicted by the continuum theory but is close to 0.66 CR (3.37 km/s), as predicted by Yoffe’s model. For GS < 5.78 J/m2, fracture results in nearly perfect cleavage with mirror-like fresh crack edges, that is, the energy released at the crack tip is completely dissipated for crack propagation. However, for GS > 5.78 J/m2, fracture edges become misty and hackle, that is, the energy released at the crack tip is dissipated not only to create new cracks but also to roughen the atomic structure. However, the changes in atomic configuration cannot be described by the continuum theory, and thus,

in which φ is the angle between the loading direction and the S−Mo bond direction, as showed in Figure S5. For SLMoS2 nanosheets with armchair-edged cracks (AC-cen, AC-2S, and AC-1Mo crack), φ = 30° and σa = 1.15σi, and for SLMoS2 nanosheets with zigzag-edged cracks (ZZ-asy crack), φ = 0° and σa = σi. Hence, the stress required to break the ZZ-asy crack is smaller than that used to break the AC-cen, AC-2S, or AC-1Mo cracks. In addition, the S−Mo bonds at AC-cen, AC-2S, and AC-1Mo crack tips are oriented with respect to the loading direction by 30°, while the S−Mo bonds at ZZ-asy crack tip are oriented parallel to the loading direction. Hence, the S−Mo bonds at AC-cen, AC-2S, and AC-1Mo crack tips should be rotated and stretched before breaking. Furthermore, there are eight S−Mo 1356


Article


Figure 7. Atomic configuration evolution during crack extension for (a) AC-cen, (b) AC-2S, (c) AC-1Mo, and (d) ZZ-asy cracks with an initial crack length of 5 nm. The closeups around the crack tip are plotted in the insets, and the atoms are colored according to the stress per atom along the tensile loading direction.

bonds at an AC-cen crack, but four at AC-2S and AC-1Mo cracks. So, larger energy is needed to break S−Mo bonds at AC-cen cracks, resulting in the highest fracture strength. For comparison, the detailed atomic configurations at crack tips of graphene and SLBP are also given in Figure S5. The ideal hexagonal lattice in graphene leads to the similar results, but the in-plane wrinkling in SLBP induces large anisotropy in fracture strength.28,29 The results are consistent with the experimental observations. For instance, Guo et al.60 found that it was preferentially cleaved along the zigzag orientation on monolayer MoS2, and Wang et al.61 observed the cracks with zigzag edges by using aberrationcorrected TEM. Figure 8a−c shows the typical atomic configuration and crack propagation path of SLMoS2 nanosheets tensile loaded along a direction of 60°, 30°, and 0°, respectively, with respect to the cracks. Armchair edges are mainly distributed at the cracks in the cases with a crack angle of 30° and 90°, but zigzag edges in the cases with a crack angle of 0° and 60°. The stress concentration around the crack tip promotes the breaking of S−Mo bonds, and the cracks extend along zigzag direction at early stage; afterward, the newly formed crack edges become misty and hackle. So, fracture along the zigzag direction is preferred in SLMoS2 nanosheets with cracks under tensile loading. In order to study the in-plane fracture during mechanical exfoliation, taking AC-2S and ZZ-asy cracks as examples, the tearing loading (mode III) is also simulated, and the results are shown in Figure 9. Figure 9a shows the potential energy change per atom (ΔPe) in the SLMoS2 nanosheet with ZZ-asy cracks. After elastic deformation, ΔPe increases stepwise owing to the propagation and kink of cracks. Then ΔPe decreases sharply at a tearing time of 2.42 ns,

Figure 8. Atomic configuration evolution during crack extension for a crack angle of θ = (a) 60°, (b) 30°, and (c) 0°. The closeups around the crack tip are shown in the insets, and the atoms are colored according to the stress per atom along the tensile loading direction.

corresponding to the cleavage of the SLMoS2 nanosheet. It is difficult to measure the crack length because of the out-of-plane displacement during tearing, and thus, it is impossible to calculate 1357


Article


cracks are always parallel to the zigzag direction, as observed in mechanical exfoliation.

4. CONCLUSION In summary, MD simulations are conducted to investigate the nanocrack propagation in the SLMoS2 nanosheet with different crack types. The brittle fracture of the SLMoS2 nanosheet is verified. The temperature, strain rate, crack type, and loading mode dependent behaviors are explored in detail. It is found that the fracture strength decreases with increasing initial crack length, which could be well described by the modified Griffith criterion. The GS decreases with increasing initial crack length, crack angle, and temperature but was not sensitive to strain rate. The V̅ is substantially reduced with increasing initial crack length and crack angle but is almost independent of temperature and strain rate. The V̅ at lower GS is well predicted by linear elastodynamic theory, but the V̅ approaches 66% Rayleigh-wave speed at a higher GS of >5.78 J/m2. The cracks prefer to fracture along the zigzag direction of SLMoS2 nanosheets, as predicted by the static mechanics. The results provide us a clear understanding on the dispersive data of measured fracture stress and the distinctive edge profiles of SLMoS2 nanosheets.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.7b10094. Lattice parameters, thickness, elastic modulus, Poisson’s ratio and edge energy of SLMoS2 and SLMoSe2; crack length as a function of the applied strain on SLMoS2 nanosheet; crack velocity as a function of current crack length; approximate semicircular stress wave around crack tip fronts of SLMoS2 nanosheet; newly formed zigzag edges in SLMoS2 nanosheet; detailed atomic structure at crack tip for SLMoS2, graphene, and SLBP (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Hongwei Bao: 0000-0001-7086-1971 Fei Ma: 0000-0002-3911-7121 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was jointly supported by National Natural Science Foundation of China (Grant No. 51771144, 51471130, 51501012), Natural Science Foundation of Shaanxi Province (No. 2017JZ015), the Fund of the State Key Laboratory of Solidification Processing in NWPU (SKLSP201708), and Fundamental Research Funds for the Central Universities.

Figure 9. (a) Potential energy change per atom (ΔPe) in SLMoS2 nanosheet with ZZ-asy crack as a function of tearing time. (b,c) Atomic configurations of SLMoS2 nanosheets with AC-2S and ZZ-asy cracks, respectively.

■

the energy release rate and crack velocity. As shown in the inset of Figure 9a, the maximum ΔPe is not sensitive to the initial crack length during tearing loading, owing to the random deviation from the initial crack direction into zigzag one at the crack tip. Figure 9b,c shows the atomic configurations of the SLMoS2 nanosheets with different initial crack length for AC-2S and ZZ-asy cracks, respectively. Interestingly, the newly formed

REFERENCES

(1) Wu, Z.; Fang, B.; Wang, Z.; Wang, C.; Liu, Z.; Liu, F.; Wang, W.; Alfantazi, A.; Wang, D.; Wilkinson, D. P. MoS2 Nanosheets: A Designed Structure With High Active Site Density for The Hydrogen Evolution Reaction. ACS Catal. 2013, 3, 2101−2107. (2) Fontana, M.; Deppe, T.; Boyd, A. K.; Rinzan, M.; Liu, A. Y.; Paranjape, M.; Barbara, P. Electron-hole Transport and Photovoltaic Effect in Gated MoS2 Schottky Junctions. Sci. Rep. 2013, 3, 1634. 1358


Article


(26) Zhao, J. H.; Kou, L. Z.; Jiang, J. W.; Rabczuk, T. Tension-induced Phase Transition of Single Layer Molybdenum Disulphide (MoS2) at Low Temperatures. Nanotechnology 2014, 25, 295701. (27) Xiong, S.; Cao, G. X. Molecular Dynamics Simulations of Mechanical Properties of Monolayer MoS2. Nanotechnology 2015, 26, 185705. (28) Yin, H. Q.; Qi, H. J.; Fan, F. F.; Zhu, T.; Wang, B. L.; Wei, Y. J. Griffith Criterion for Brittle Fracture in Graphene. Nano Lett. 2015, 15, 1918−1924. (29) Liu, N.; Hong, J. W.; Pidaparti, R.; Wang, X. Q. Fracture Patterns and The Energy Release Rate of Phosphorene. Nanoscale 2016, 8, 5728. (30) Budarapu, P. R.; Javvaji, B.; Sutrakar, V. K.; Mahapatra, D. R.; Zi, G.; Rabczuk, T. Crack Propagation in Graphene. J. Appl. Phys. 2015, 118, 064307. (31) Abadi, R.; Uma, R. P.; Izadifar, M.; Rabczuk, T. Investigation of Crack Propagation and Existing Notch on The Mechanical Response of Polycrystalline Hexagonal Boron-nitride Nanosheets. Comput. Mater. Sci. 2017, 131, 86−99. (32) Xiong, Q. L.; Li, Z. H.; Tian, X. G. Fracture Behaviors of Precracked Monolayer Molybdenum Disulfide: A Molecular Dynamics Study. Beilstein J. Nanotechnol. 2016, 7, 1411−1420. (33) Wang, X. N.; Tabarraei, A.; Spearot, D. E. Fracture Mechanics of Monolayer Molybdenum Disulfide. Nanotechnology 2015, 26, 175703. (34) Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Electric Field Effect in Atomically Thin Carbon Films. Science 2004, 306, 666−669. (35) Tang, D. M.; Kvashnin, D. G.; Najmaei, S.; Bando, Y.; Kimoto, K.; Koskinen, P.; Ajayan, P. M.; Yakobson, B. I.; Sorokin, P. B.; Lou, J.; Golberg, D. Nanomechanical Cleavage of Molybdenum Disulphide Atomic Layers. Nat. Commun. 2014, 5, 3631. (36) Lammps. http://lammps.sandia.gov/index.html. (37) Jiang, J. W.; Park, H. S.; Rabczuk, T. Molecular Dynamics Simulations of Single-layer Molybdenum Disulphide (MoS2): StillingerWeber Parametrization, Mechanical properties, and Thermal Conductivity. J. Appl. Phys. 2013, 114, 064307. (38) Jiang, J. W. Phonon Bandgap Engineering of Strained Monolayer MoS2. Nanoscale 2014, 6, 8326. (39) Singh, S. K.; Neek-Amal, M.; Costamagna, S.; Peeters, F. M. Rippling, Buckling, and Melting of Single-and Multilayer MoS2. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 014101. (40) Subramaniyan, A. K.; Sun, C. T. Continuum Interpretation of Virial Stress in Molecular Simulation. Int. J. Solids Struct. 2008, 45, 4340−4346. (41) Stukowski, A. Visualization and Analysis of Atomistic Simulation Data with OVITO-The Open Visualization Tool. Modell. Simul. Mater. Sci. Eng. 2010, 18, 015012. (42) Qi, Z. N.; Cao, P. H.; Park, H. S. Density Functional Theory Calculation of Edge Stresses in Monolayer MoS2. J. Appl. Phys. 2013, 114, 163508. (43) Liu, F.; Tang, Q. H.; Wang, T. C. Intrinsic Notch Effect Leads to Breakdown of Griffith Criterion in Graphene. Small 2017, 13, 1700028. (44) Griffith, A. A. The Phenomena of Rupture and Flow in Solids. Philos. Trans. R. Soc., A 1921, 221, 163−198. (45) Zhang, Z.; Wang, X. Q.; Lee, J. D. An Atomistic Methodology of Energy Release Rate for Graphene at Nanoscale. J. Appl. Phys. 2014, 115, 114314. (46) Sharon, E.; Gross, S. P.; Fineberg, J. Energy Dissipation in Dynamic Fracture. Phys. Rev. Lett. 1996, 76, 2117−2120. (47) Chen, M. Q.; Quek, S. S.; Sha, Z. D.; Chiu, C. H.; Pei, Q. X.; Zhang, Y. W. Effects of Grain size, Temperature and Strain Rate on The Mechanical Properties of Polycrystalline Graphene-A Molecular Dynamics Study. Carbon 2015, 85, 135−146. (48) Zhang, B.; Yang, G.; Xu, H. M. Instability of Supersonic Crack in Graphene. Phys. B 2014, 434, 145−148. (49) Abraham, F. F.; Gao, H. J. How Fast Can Cracks Propagate? Phys. Rev. Lett. 2000, 84, 3113−3116. (50) Buehler, M. J.; Gao, H. J. Dynamical Fracture Instabilities Due to Local Hyperelasticity at Crack Tips. Nature 2006, 439, 307−310.

(3) Shanmugam, M.; Durcan, C. A.; Yu, B. Layered Semiconductor Molybdenum Disulfide Nanomembrane Based Schottky-barrier Solar Cells. Nanoscale 2012, 4, 7399−7405. (4) Xu, M.; Liang, T.; Shi, M.; Chen, H. Graphene-like Twodimensional Materials. Chem. Rev. 2013, 113, 3766−3798. (5) Ataca, C.; Şahin, H.; Ciraci, S. Single-layer MX2 Transition-metal Oxides and Dichalcogenides In A Honeycomb-like Structure. J. Phys. Chem. C 2012, 116, 8983−8999. (6) Li, L. F.; Wang, Y. L.; Xie, S. Y.; Li, X. B.; Wang, Y. Q.; Wu, R. T.; Sun, H. B.; Zhang, S. B.; Gao, H. J. Two-dimensional Transition Metal Honeycomb Realized: Hf on Ir (111). Nano Lett. 2013, 13, 4671−4674. (7) Jin, C.; Lin, F.; Suenaga, K.; Iijima, S. Fabrication of A Freestanding Boron Nitride Single Layer and Its Defect Assignments. Phys. Rev. Lett. 2009, 102, 195505. (8) Brent, J. R.; Savjani, N.; Lewis, E. A.; Haigh, S. J.; Lewis, D. J.; O’Brien, P. Production of Few-layer Phosphorene by Liquid Exfoliation of Black Phosphorus. Chem. Commun. 2014, 50, 13338−13341. (9) Meng, L.; Wang, Y. L.; Zhang, L. Z.; Du, S. X.; Wu, R. T.; Li, L. F.; Zhang, Y.; Li, G.; Zhou, H. T.; Hofer, W. A. Buckled Silicene Formation on Ir(111). Nano Lett. 2013, 613, 685−690. (10) Chhowalla, M.; Shin, H. S.; Eda, G.; Li, L. J.; Loh, K. P.; Zhang, H. The Chemistry of Two-dimensional Layered Transition Metal Dichalcogenide Nanosheets. Nat. Chem. 2013, 5, 263−275. (11) Mak, K. F.; Lee, C.; Hone, J.; Shan, J.; Heinz, T. F. Seeing Manybody Effects in Single- and Few-layer Graphene: Observation of 2Dimensional Saddle Point Excitons. Phys. Rev. Lett. 2010, 105, 136805. (12) Bertolazzi, S.; Brivio, J.; Kis, A. Stretching and Breaking of Ultrathin MoS2. ACS Nano 2011, 5, 9703−9709. (13) Na, S. R.; Wang, X. H.; Piner, R. D.; Huang, R.; Willson, C. G.; Liechti, K. M. Cracking of Polycrystalline Graphene on Copper under Tension. ACS Nano 2016, 10, 9616−9625. (14) Zhang, P.; Ma, L. L.; Fan, F. F.; Zeng, Z.; Peng, C.; Loya, P. E.; Liu, Z.; Gong, Y. J.; Zhang, J. N.; Zhang, X. X.; et al. Fracture Toughness of Graphene. Nat. Commun. 2014, 5, 3728. (15) Yang, Y. C.; Li, X.; Wen, M. R.; Hacopian, E.; Chen, W. B.; Gong, Y. J.; Zhang, J.; Li, B.; Zhou, W.; Ajayan, P. M.; et al. Brittle Fracture of 2D MoSe2. Adv. Mater. 2017, 29, 1604201. (16) Sundaram, R. S.; Engel, M.; Lombardo, A.; Krupke, R.; Ferrari, A. C.; Avouris, P.; Steiner, M. Electroluminescence in Single Layer MoS2. Nano Lett. 2013, 13, 1416−1421. (17) Jeong, H. Y.; Lee, S. Y.; Ly, T. H.; Han, G. H.; Kim, H.; Nam, H.; Jiong, Z.; Shin, B. G.; Yun, S. J.; Kim, J.; et al. Visualizing Point Defects in Transition-metal Dichalcogenides Using Optical Microscopy. ACS Nano 2016, 10, 770−777. (18) Hao, S.; Yang, B. C.; Gao, Y. L. Fracture-induced Nanoscrolls from CVD-grown Monolayer Molybdenum Disulfide. Phys. Status Solidi RRL 2016, 10, 549−553. (19) Komsa, H. P.; Krasheninnikov, A. V. Native Defects in Bulk and Monolayer MoS2 from First Principles. Phys. Rev. B: Condens. Matter Mater. Phys. 2015, 91, 125304. (20) Zhang, B.; Mei, L.; Xiao, H. F. Nanofracture in Graphene under Complex Mechanical Stresses. Appl. Phys. Lett. 2012, 101, 121915. (21) Zhang, T.; Li, X. Y.; Gao, H. J. Fracture of Graphene: A Review. Int. J. Fract. 2015, 196, 1−31. (22) Wang, S. S.; Li, H. S.; Sawada, H.; Allen, C. S.; Kirkland, A. I.; Grossman, J. C.; Warner, J. H. Atomic Structure and Formation Mechanism of Sub-nanometer Pores in 2D Monolayer MoS2. Nanoscale 2017, 9, 6417−6426. (23) Ma, F.; Sun, Y. J.; Ma, D. Y.; Xu, K. W.; Chu, P. K. Reversible Phase Transformation in Graphene Nano-ribbons: Lattice Shearing Based Mechanism. Acta Mater. 2011, 59, 6783−6789. (24) Bao, H. W.; Huang, Y. H.; Yang, Z.; Miao, Y. P.; Chu, P. K.; Xu, K. W.; Ma, F. Tensile Loading Induced Phase Transition and Rippling in Single-layer MoS2. Appl. Surf. Sci. 2017, 404, 180−187. (25) Jiang, J. W. Parametrization of Stillinger-Weber Potential Based on Valence Force Field Model: Application to Single-layer MoS2 and Black Phosphorus. Nanotechnology 2015, 26, 315706. 1359


Article

The Journal of Physical Chemistry C (51) Cooper, R. C.; Lee, C.; Marianetti, C. A.; Wei, X. D.; Hone, J.; Kysar, J. W. Nonlinear Elastic Behavior of Two-Dimensional Molybdenum Disulfide. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 87, 035423. (52) Anderson, T. L. Fracture Mechanics: Fundamentals and Applications; CRC Press: Boca Raton, FL, 2005. (53) Mezil, S.; Chonan, K.; Otsuka, P. H.; Tomoda, M.; Matsuda, O.; Lee, S. H.; Wright, O. B. Extraordinary Transmission of Gigahertz Surface Acoustic Waves. Sci. Rep. 2016, 6, 33380. (54) McKenna, A. J.; Eliason, J. K.; Flannigan, D. J. Spatiotemporal Evolution of Coherent Elastic Strain Waves in A Single MoS2 Flake. Nano Lett. 2017, 17, 3952−3958. (55) Cramer, T.; Wanner, A.; Gumbsch, P. Energy Dissipation and Path Instabilities in Dynamic Fracture of Silicon Single Crystals. Phys. Rev. Lett. 2000, 85, 788−791. (56) Scheibert, J.; Guerra; Célarié, C.; Dalmas, D.; Bonamy, D. BrittleQuasibrittle Transition in Dynamic Fracture: An Energetic Signature. Phys. Rev. Lett. 2010, 104, 045501. (57) Swadener, J. G.; Baskes, M. I.; Nastasi, M. Molecular Dynamics Simulation of Brittle Fracture in Silicon. Phys. Rev. Lett. 2002, 89, 085503. (58) Rosakis, A. J.; Samudrala, O.; Coker, D. Cracks Faster Than The Shear Wave Speed. Science 1999, 284, 1337−1440. (59) Livne, A.; Bouchbinder, E.; Svetlizky, I.; Fineberg, J. The Near-tip Fields of Fast Cracks. Science 2010, 327, 1359−1363. (60) Guo, Y.; Liu, C. R.; Yin, Q. F.; Wei, C. R.; Lin, S. H.; Hoffman, T. B.; Zhao, Y. D.; Edgar, J. H.; Chen, Q.; Lau, S. P.; et al. Distinctive Inplane Cleavage Behaviors of Two-dimensional Layered Materials. ACS Nano 2016, 10, 8980−8988. (61) Wang, S. S.; Qin, Z.; Jung, G. S.; Martin-Martinez, F. J.; Zhang, K.; Buehler, M. J.; Warner, J. H. Atomically Sharp Crack Tips in Monolayer MoS2 and Their Enhanced Toughness by Vacancy Defects. ACS Nano 2016, 10, 9831−9839.

1360


Molecular Dynamics Simulation of Nanocrack ... - ACS Publications

Recommend Documents