Effects of hydrogen adsorption on the fracture properties of graphene

Effects of hydrogen adsorption on the fracture properties of graphene

Computational Materials Science 121 (2016) 151–158 Contents lists available at ScienceDirect Computational Materials Science journal homepage: www.e...

5MB Sizes 0 Downloads 36 Views

Computational Materials Science 121 (2016) 151–158

Contents lists available at ScienceDirect

Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci

Effects of hydrogen adsorption on the fracture properties of graphene Alireza Tabarraei ⇑, Xiaonan Wang, Dan Jia Department of Mechanical Engineering & Engineering Science, University of North Carolina at Charlotte, Charlotte, NC 28223, United States

a r t i c l e

i n f o

Article history: Received 25 February 2016 Received in revised form 16 April 2016 Accepted 25 April 2016

Keywords: Graphene Fracture Critical stress intensity factor Hydrogen embrittlement

a b s t r a c t Adsorption of environment molecules to graphene alters the atomic structure of graphene and impacts fracture properties of graphene. In this paper, we use molecular dynamics (MD) modeling to study the effects of hydrogen adsorption on the fracture mechanics of graphene. Both armchair and zigzag cracks under mode I and II fracture loadings are considered. Our molecular dynamics simulations predict that adsorption of hydrogen atoms to the crack tip or surface leads to a reduction in graphene toughness and can alter the crack propagation paths. Based on our MD simulations, the location of carbon atoms at which hydrogen atoms are attached is a main factor in the level of impact that hydrogen atoms have on the fracture properties of graphene. Ó 2016 Elsevier B.V. All rights reserved.

1. Introduction Fantastic mechanical [1], electrical [2,3] and physical [4] properties of graphene make it an ideal candidate for a wide spectrum of applications such as composite materials [5], nano-sensors and nano-devices [6], nanomedicine [7], and energy storage [8,9]. Understanding mechanical properties of graphene is essential toward reliable design and manufacturing of graphene-based nano-devices and nano-materials. Specially, it is important to predict and eventually prevent the fracture of graphene, which necessitates a thorough understanding of the fracture properties of graphene. Due to the difficulties in designing and performing experiments at nanoscale, experimental approach for determining the physical and mechanical properties of nanomaterials is very challenging. Atomistic modeling, on the other hand, can provide valuable insights into the properties of nanomaterials [10–15]. In particular, atomistic modeling techniques have proven themselves as invaluable tools in determining fracture properties of graphene and graphen-like two-dimensional materials [16–27]. Omeltchenko et al. [21] used the first generation of reactive empirical bond order (REBO) potential to study mode I crack propagation in graphene. Khare et al. [22] used a coupled quantum mechancis/molecular mechanics model to study the effect of defects and cracks on the mechanical properties of graphene. Xu et al. [23] used a coupled quantum/continuum mechanics approach to obtain the critical stress intensity factors of graphene under quasi-static loading. ⇑ Corresponding author. E-mail address: [email protected] (A. Tabarraei). http://dx.doi.org/10.1016/j.commatsci.2016.04.037 0927-0256/Ó 2016 Elsevier B.V. All rights reserved.

Zhang et al. [24] used molecular dynamics simulations to study fracture in graphene under mixed mode I and II loading. Dewapriya et al. [25,26] used atomistic simulations to study how nanocracks impact the fracture strength of graphene sheets. Tabarraei et al. [27] used density functional theory (DFT) to study the fracture mechanism of graphene nanoribbons with different chirality. Jung et al. [28] used molecular dynamics method to show that polycrystalline graphene releases up to 50% more fracture energy than pristine graphene. The previous studies on graphene fracture have focused on a graphene sheet with no adsorbed atoms of any kinds. In many applications, chemical functionalization [29] of graphene via attachment of atoms have been used as a mean for tuning the properties of graphene. Among different functional groups, hydrogen has attracted considerable attentions. Adsorption of hydrogen atoms can change the graphene atomic structure [30] which in turn can modify electronic and magnetic properties. For example, the adsorption of hydrogen atoms can induce a band gap around Fermi level which converts graphene from a semimetal to a semiconductor [31], or it can lead to the formation of a magnetic structure [30,32] in graphene. In addition to deliberate functionalization of graphene sheet with hydrogen, the presence of hydrogen in many environments such as water and air provides many opportunities for its introduction to graphene. The adsorption of hydrogen atoms alter the mechanical properties of graphene. Graphene sheets with adsorbed hydrogen atoms have a lower Young’s modulus, tensile strength and fracture strain than those of pristine graphene sheets [33]. Furthermore, it is shown that hydrogenation of graphene grain boundaries leads to substantial reduction in the strength

152

A. Tabarraei et al. / Computational Materials Science 121 (2016) 151–158

and ductility of polycrystalline graphene sheets [34]. Recent experiments provide evidence that the attachment of hydrogen atoms to the crack tip or surface can leads to the environment assisted cracking of graphene sheets [35], as a consequence of which the graphene undergoes sub-critical crack growth. Attachment of hydrogen atoms to the carbon atoms located in the vicinity of crack tip or edges change the morphology of crack and subsequently affect the fracture properties of graphene. The impact of hydrogen atoms on the fracture properties of graphene sheets has not yet been fully investigated. In this paper, we use molecular dynamics simulations to probe how hydrogen atoms modify the crack configuration and impact the fracture mechanism of graphene sheets. 2. Simulation models and methodology Our molecular dynamics simulations are performed using LAMMPS (Large-Scale Atomic/Molecular Massively Paralle Simulator) package. The adaptive intermolecular reactive bond order (AIREBO) [36] potential is used to consider the interaction between atoms in the graphene sheet. AIREBO potential allows bond breaking and forming, hence it is suitable for modeling failure. The original AIREBO potential employs two cutoff distances with values of 1.7 and 2 to allow a smooth transition of cutoff function from one to zero. It has been observed that such a cutoff function generates non-physical strain hardening in the stress–strain behavior of carbon nanostructures which is due to the discontinuity in the second derivative of the cutoff function. As a remedy of this issue it is suggested to eliminate the transition zone; the cutoff function is 1.0 for distances less than R and be 0 otherwise. A value of 1.9–2.2 Å [24,34,37–39] has been used for the cutoff distance R. In this paper we use a value of 1.92 Å, as is suggested in Ref. [37] and has been used in the past to successfully capture the fracture properties of graphene [24,40]. We consider the crack propagation of zigzag and armchair cracks under mode I and II loading. For the purpose of computational efficiency, a small circular-shaped domain of radius 70 Å from the crack tip is used to model a semi-infinite crack in a graphene sheet. The size of the domain is chosen such that the whole domain falls within the K-dominant field. The initial zigzag and armchair cracks are generated by removing two or three layers of atoms, respectively. Using this geometry, our domain is composed of about 6000 carbon atoms. Mode I and II displacement field are applied as the boundary conditions to the outermost layers of atoms. The boundary atoms are shown in blue in Fig. 1. We have assumed that the far field behavior is well-represented by a linear, isotropic response because the crack-tip strain filed reduces with distance and graphene behavior is anisotropic only at large strain [41,42]. Atomic displacements on the boundary of the molecular dynamic (MD) zone are prescribed using the crack-tip asymptotic displacement field given by [43]

rffiffiffiffiffiffiffi rffiffiffiffiffiffiffi r h K II r h cos ½j  cos h þ sin ½j þ 2 þ cos h; 2p 2 2 2l 2p rffiffiffiffiffiffiffi rffiffiffiffiffiffiffi KI r h K II r h uy ¼ sin ½j  cos h  cos ½j  2 þ cos h; 2 2 2l 2p 2l 2p

ux ¼

KI 2l

ð1aÞ ð1bÞ

where r is the distance between a boundary atom and crack-tip, h is the polar angle depicted in Fig. 1, K I and K II are mode I and II stress m, intensity factors, respectively, l is the shear modulus, and j ¼ 3 1þm where m is the graphene Poisson’s ratio which is taken equal to 0.165. To apply a pure mode I or pure mode II displacement field, the values of K II or K I parameters are set equal to zero, respectively. The displacement field is gradually applied to the boundary atoms by incrementing the stress intensity factors in increments pffiffiffiffiffi of 0.01 MPa m. After each load increment, by keeping the

Fig. 1. An initial armchair crack in a graphene sheet. The K-field displacement field is applied to the boundary atoms (shown in blue) while the position of interior atoms (shown in green) are relaxed. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

boundary atoms fixed, first the energy of the system is minimized, and then the system is equilibrated for 6000 time step at a temperature of 300 K using Nosé-Hoover thermostat [44]. The velocity-Verlet scheme with a time step of 1 fs is utilized for the purpose of time integration. The stress intensity factor at which the initial crack advances or a new crack branches off the original crack is considered as the critical stress intensity factor of the domain. To check the validity of our MD modeling approach and for the purpose of comparison, we first find mode I and II critical stress intensity factors of armchair and zigzag cracks when no hydrogen atoms is adsorbed. Our MD simulations predict a mode I critical pffiffiffiffiffi stress intensity factor ðK cI Þ of 3.22 MPa m for armchair cracks pffiffiffiffiffi and 2.88 MPa m for zigzag crack, whereas mode II critical stress intensity factors ðK cII Þ of armchair and zigzag cracks are pffiffiffiffiffi pffiffiffiffiffi 3.95 MPa m and 3.65 MPa m, respectively. These values are in good agreement with the values reported previously [23,24,45]. The crack propagation path of armchair and zigzag cracks under mode I and II loading are shown in Fig. 2. These results show that zigzag and armchair cracks propagate along a zigzag path under both mode I and II loading. This is consistent with the previous results [23,24] and confirms that in graphene the surface energy of zigzag edges is smaller than that of armchair edges. Moreover, our MD simulations predict that under mode II loading the top surface of the crack which is under a compressive load will undergo a buckling (out-of-plane) deformation, a phenomenon which has been observed in thin cracked plates under tensile and shear loading [46–48] as well. The buckling of a graphene sheet with an armchair crack under mode II loading is shown in Fig. 3. The impact of the out-of-plane deformation of graphene on fracture toughness of graphene is studied in [49] and it is shown that the out-of-plane deformation of the graphene can be used as a topological defect to increase the fracture toughness of graphene. To investigate the impact of the domain size on the results, we have modeled the crack propagation paths and obtained the crack intensity factors using domains with 60,000 atoms. The results

A. Tabarraei et al. / Computational Materials Science 121 (2016) 151–158

153

Fig. 2. Mode I and II crack propagation path of armchair and zigzag cracks and zoom-in of their crack tips. The broken bond at the crack tip is shown in green. (a) Zigzag crack under mode I; (b) armchair crack under mode I; (c) zigzag crack under mode II; and (d) armchair crack under mode II. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 3. Buckling of a crack top surface under mode II loading.

obtained from this domain are not significantly different from the previous results, indicating the domain size did not impact the results obtained from the smaller domain. Therefore, to reduce the computational costs we will use the small domain in the following modelings. 3. Results and discussion The influence of hydrogen adsorption on the fracture properties of graphene is studied by adding hydrogen atoms to the carbon atoms located at the crack tip or crack edges, as are shown in Fig. 4. Depending on the position of the carbon atom at which a hydrogen is attached, the relaxed position of the hydrogen atom can be in or out of the graphene plane. If a hydrogen atom is adsorbed to an edge carbon atom with a dangling bond (a carbon

atom with only two C–C covalent bonds), the relaxed position of the hydrogen atom will be in the same plane as graphene; whereas the relaxed position of a hydrogen atom attached to a bulk carbon atom is out of the graphene plane. The bulk and edge carbon atoms are shown respectively in red and green in Fig. 4. In this paper we study the effect of both in-plane and out-ofplane hydrogen atoms on the fracture properties of zigzag and armchair cracks. The impact of the number of adsorbed hydrogen atoms are examined by considering models with different number of hydrogen atoms. For mode I fracture loading, we consider two cases: hydrogen atoms are attached to one surface or both crack surfaces. Since under mode II loading the domain deformation is not symmetric, three cases are considered for mode II: hydrogen atoms attached to the top surface, bottom surface or both surfaces of the crack.

154

A. Tabarraei et al. / Computational Materials Science 121 (2016) 151–158

(a)

(b)

Fig. 4. The tip structure of (a) zigzag and (b) armchair cracks.

3.1. In-plane hydrogen atoms The critical stress intensity factors of armchair and zigzag cracks as a function of the number of adsorbed hydrogen atoms are shown in Fig. 5. These plots show that the attachment of hydrogen atoms to the crack tip and surface leads to a reduction in mode I and II critical stress intensity factors of both armchair and zigzag cracks. The attachment of a hydrogen atoms to only one surface of the cracks, leads to a reduction of the critical stress intensity factor pffiffiffiffiffi pffiffiffiffiffi ðK cri m to 2.71 MPa m and I Þ of the zigzag crack from 2.88 MPa pffiffiffiffiffi pffiffiffiffiffi for armchair cracks from 3.22 MPa m to 2.84 MPa m. For the case when hydrogen atoms are adsorbed to both surfaces, the mode I critical stress intensity factor of armchair cracks is the same as when hydrogen atoms are adsorbed to one surface. However, for zigzag cracks the critical stress intensity factor is further reduced pffiffiffiffiffi to 2.49 MPa m when hydrogen atoms are attached to both surfaces. Plots of Fig. 5 show that adsorption of hydrogen atoms to the top or bottom surfaces have different effects on mode II critical stress intensity factors ðK cri II1 Þ. While hydrogen attachment to the top surface (compressive surface) does not have an important impact on mode II critical stress intensity factor, hydrogen attachment to the bottom (tensile) surface significantly reduce mode II critical stress intensity factor. Due to the attachment of hydrogen atoms to the bottom surface, K cri II1 of zigzag cracks reduces from pffiffiffiffiffi pffiffiffiffiffi 3.70 MPa m to 3.03 MPa m, while K cri II1 of armchair cracks pffiffiffiffiffi pffiffiffiffi ffi reduces from 4.0 MPa m to 3.68 MPa m. For both armchair and zigzag cracks, when hydrogen atoms are adsorbed to both surfaces, K cri II1 is almost the same as when hydrogens are attached to only the bottom surface, confirming that the attachment of hydrogen atoms to the top surface does not impact K cri II1 . The plots of Fig. 5

Number of H atoms on each crack edge (a)

also show that for both armchair and zigzag cracks under mode I or II loading, only the first few hydrogen atoms which are close to the crack tip impact the critical stress intensity factors and the adsorption of hydrogen atoms to the carbon atoms located always from the tip does not affect the critical stress intensity factors. The mode I and II crack propagation paths of zigzag and armchair cracks with six hydrogen atom adsorbed to each crack face are shown in Fig. 6. Comparison of Fig. 6 with Fig. 2 indicates that in-plane hydrogen atoms do not alter the crack propagation paths of armchair or zigzag cracks. This is a results of the fact that inplane hydrogen attachment to carbon atoms does not change the morphology of the crack tip and surface, i.e. bond stretching and rotation does not occur, hence crack growth path does not change.

3.2. Out-of-plane hydrogen atoms Opposed to in-plane hydrogen atoms, out-of-plane hydrogen atoms significantly alter the morphology of the crack tip and crack edges. This is demonstrated in Fig. 7 for a zigzag crack with eight hydrogen atoms attached to its top edge. This figure shows that relaxed position of carbon atoms at which a hydrogen is adsorbed is out of the graphene plane. For example, for the case shown in Fig. 7, the out-of-plane displacement of the carbon atom located at the crack tip is 1:30 Å. This amount of out-of-plane deformation induces large stretching and rotation of C–C bonds located in the vicinity of crack tip and edges. Our MD simulations predict a similar change of morphology for armchair cracks as well. The critical stress intensity factors of zigzag and armchair cracks with out-of-plane hydrogen atoms are presented in Fig. 8. Comparing the plots of Fig. 8 with those of Fig. 5 indicate that the adsorption of out-of-plane hydrogen atoms have a more sever impact on the critical stress intensity factors than the adsorption of in-plane hydrogen atoms.

Number of H atoms on each crack edge (b)

Fig. 5. Critical stress intensity factors of zigzag and armchair cracks with adsorbed in-plane hydrogen atoms. (a) Zigzag cracks and (b) armchair cracks.

155

A. Tabarraei et al. / Computational Materials Science 121 (2016) 151–158

Fig. 6. Propagation path of zigzag and armchair cracks with six adsorbed hydrogen atoms on each face under mode I and II loading. A zoom in of the crack-tip in its relaxed form is shown for each case on which the first broken bond is shown in green. Carbon and hydrogen atoms are shown in black and purple respectively. (a) Zigzag crack under mode I, (b) zigzag crack under mode II, (c) armchair crack under mode I, and (d) armchair crack under mode II fracture loading. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 7. Crack tip and crack surface morphology of a zigzag crack with eight out-of-plane hydrogen atoms adsorbed to its top surface. Atoms are colored based on their z-coordinate. (a) Top view of the crack, (b) zoom-in of the crack tip, (c) side view of the top crack edge, and (d) zoom-in of the side view in the vicinity of crack. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Number of H atoms on each crack edge

(a)

Number of H atoms on each crack edge

(b)

Fig. 8. Critical stress intensity factors of (a) zigzag and (b) armchair cracks with adsorbed out-of-plane hydrogen atoms.

156

A. Tabarraei et al. / Computational Materials Science 121 (2016) 151–158

Fig. 9. Crack propagation path of zigzag cracks with out-of-plane adsorbed hydrogen atoms under mode I loading.

Fig. 10. Crack propagation path of armchair cracks with out-of-plane adsorbed hydrogen atoms under mode I loading.

If out-of-plane hydrogen atoms are adsorbed to only one surface of cracks, the mode I critical stress intensity factors of zigzag pffiffiffiffiffi pffiffiffiffiffi and armchair cracks reduces to 2.25 MPa m and 1.60 MPa m, respectively. Attachment of hydrogen atoms to both crack surfaces does not further reduce K cri I1 of armchair cracks, while it reduces the pffiffiffiffiffi cri K I1 of zigzag cracks further to 1.50 MPa m. Similar to what is observed for in-plane hydrogen atoms, adsorption of out-of-plane hydrogen atoms to the top (compressive) face of cracks, does not alter their mode II critical stress intensity factors K cri II1 . However, if hydrogen atoms are adsorbed to the bottom faces of the cracks, their K cri II1 is significantly reduced. For such a case, mode II critical stress intensity factors of zigzag and pffiffiffiffiffi pffiffiffiffiffi armchair cracks are 2.45 MPa m and 2.50 MPa m, respectively. The mode II critical stress intensity factors of armchair and zigzag cracks does not further alter if hydrogen atoms are attached to both surfaces which support the fact that hydrogen attachment to the top surface which is in compression under mode II does not impact the value of K cri II1 . The crack propagation path of zigzag and armchair cracks under mode I and II loadings are shown in Figs. 9–11. The mode I and

mode II propagation path of zigzag cracks shown in Fig. 9 show that despite change in the morphology of cracks induced by hydrogen atoms, the propagation path of zigzag cracks does not alter. The bond stretching and rotation at the crack tip and edges do not lead to any extra bond failure at the tip or edges of the crack and growth path is the same as is observed for zigzag cracks with no hydrogen atoms. Out-of-plane hydrogen atoms have a more sever impact on the propagation paths of armchair cracks. As is shown in Fig. 10, under mode I loading multiple C–C bonds located on the crack edges fail before the initial crack propagates. The failure of edge bonds generates a non-smooth face for the crack, however in all the cases we have studied the initial crack eventually propagates and failure of edge bonds does not lead to crack branching. When hydrogen atoms are adsorbed to only one crack surface, the first bond that fails under mode I loading is the bond between the crack tip and the surface with adsorbed hydrogen atoms. This bond is shown in green in Fig. 10. When hydrogen atoms are attached to both surfaces, a chain of atoms form at the crack tip; indicating that some of the bonds in the crack path do not break and crack propagates by breaking the bonds behind the crack tip.

A. Tabarraei et al. / Computational Materials Science 121 (2016) 151–158

157

Fig. 11. Crack propagation path of armchair cracks under mode II loading. Out-of-plane hydrogen atoms are bonded to the compressive surface.

Similar to mode I loading, under mode II fracture loading also multiple C–C bonds located on the edges of armchair cracks break. This is true even for the case when hydrogen atoms are attached to the compressive surface, although in this case crack propagation path does not change. When hydrogen atoms are attached to the bottom surface or both surfaces, the failure of the crack edge C–C bonds completely alter the crack propagation path of armchair cracks. In these cases, the initial crack does not propagate any more, but a new crack perpendicular to the original crack nucleates and develop. Therefore, out-of-plane hydrogen atoms can significantly alter the propagation path of armchair cracks if they are attached to the tensile surface of the crack. The change in the crack growth path of graphene in the presence of chemical additives such as oxygen has been reported previously [50]. Our modelings also show that the adsorption of hydrogen atoms can modify the crack propagation path. Such results suggest that chemomechanical conditions can be used as a tool to tune and regulate the fracture path of graphene. 4. Conclusions In this paper, we have employed molecular dynamics simulations to study the effects of hydrogen adsorption on the fracture properties of graphene sheets. Both armchair and zigzag cracks under mode I and II fracture loading were considered. The MD

simulations predict that the adsorption of hydrogen atoms to the crack tip or crack surfaces reduce the critical stress intensity factors of graphene sheets. Based on our MD simulation results, the impact of hydrogen atoms on critical stress intensity factors depends on the position of the atom sites at which hydrogen atoms are adsorbed. For example, adsorption of hydrogen atoms to a crack surface which is under compressive stress in a mode II loading does not have any significant impact on mode II critical stress intensity factors, whereas attachment of hydrogen atoms to the surface which is under tensile stress significantly reduces the mode II critical stress intensity factors. The effects of hydrogen atoms on crack growth path is not the same in all the cases. The in-plane hydrogen atoms do not affect the crack surface or tip morphology and do not impact the crack propagation path. On the other hand, the adsorption of out-ofplane hydrogen atoms leads to out-of-plane deformation of carbon atoms, hence the crack tip and edges become distorted. As a result of this, out-of-plane hydrogen atoms reduce the critical stress intensity factors of graphene sheets more significantly. Out-ofplane hydrogen atoms do not alter the crack propagation paths of zigzag cracks, however they remarkably alter the propagation path of armchair cracks under mode II loading. For example, in some instances due to the presence of out-of-plane hydrogen atoms the initial armchair crack does not propagate under mode II loading; but a new crack perpendicular to the initial crack

158

A. Tabarraei et al. / Computational Materials Science 121 (2016) 151–158

branches and propagates in a path perpendicular to the initial crack. In general, based on our vast number of MD simulations we can conclude that the influence of out-of-plane hydrogen atoms on armchair cracks under mode II loading is the most significant among all the cases studied in this paper. Acknowledgment This research is supported by the NC Space Grants. The authors also acknowledge the support from the University of North Carolina at Charlotte. References [1] C. Lee, X. Wei, J.W. Kysar, J. Hone, Science 5887 (2008) 385–388. [2] K. Bolotin, K. Sikes, Z. Jiang, M. Klima, G. Fudenberg, J. Hone, P. Kim, H. Stormer, Solid State Commun. 146 (2008) 351–355. [3] E. Hwang, S. Adam, S.D. Sarma, Phys. Rev. Lett. 98 (2007) 186806. [4] A.C. Neto, F. Guinea, N. Peres, K. Novoselov, A. Geim, Rev. Mod. Phys. 81 (2009) 109. [5] S. Stankovich, D.A. Dikin, G.H.B. Dommett, K.M. Kohlhaas, E.J. Zimney, E.A. Stach, R.D. Piner, S.T. Nguyen, R.S. Ruoff, Nature 442 (2006) 282–285. [6] S. Wu, Q. He, C. Tan, Y. Wang, H. Zhang, Small 9 (2013) 1160–1172. [7] H. Mao, S. Laurent, W. Chen, O. Akhavan, M. Imani, A. Ashkarran, M. Mahmoudi, Chem. Rev. 113 (2013) 3407–3424. [8] M. Pumera, Energy Environ. Sci. 4 (2011) 668–674. [9] X. Yang, C. Cheng, Y. Wang, L. Qiu, D. Li, Science 341 (2013) 534–537. [10] A. Tabarraei, X. Wang, S. Shadalou, Mechanical properties of monolayer molybdenum disulfide, in: ASME 2014 International Mechanical Engineering Congress and Exposition, American Society of Mechanical Engineers, 2014. pp. V009T12A070–V009T12A070. [11] J. Hu, X. Ruan, Y.P. Chen, Nano Lett. 9 (7) (2009) 2730–2735. [12] A. Bagri, S.-P. Kim, R.S. Ruoff, V.B. Shenoy, Nano Lett. 11 (9) (2011) 3917–3921. [13] B. Mortazavi, S. Ahzi, Carbon 63 (2013) 460–470. [14] A. Tabarraei, Comput. Mater. Sci. 108 (2015) 66–71. [15] A. Tabarraei, X. Wang, Appl. Phys. Lett. 108 (2016) 181904. [16] A. Tabarraei, X. Wang, Mater. Sci. Eng.: A 641 (2015) 225–230. [17] B. Zhang, L. Mei, H. Xiao, Appl. Phys. Lett. 101 (12) (2012) 121915. [18] M.-Q. Le, R.C. Batra, Comput. Mater. Sci. 69 (2013) 381–388.

[19] X. Wang, A. Tabarraei, D.E. Spearot, Nanotechnology 26 (17) (2015) 175703. [20] A. Tabarraei, X. Wang, A. Sadeghirad, J. Song, Finite Elem. Anal. Des. 92 (2014) 36–49. [21] A. Omeltchenko, J. Yu, R. Kalia, P. Vashishta, Phys. Rev. Lett. 78 (1997) 2148. [22] R. Khare, S. Mielke, J. Paci, S. Zhang, R. Ballarini, G. Schatz, T. Belytschko, Phys. Rev. B 75 (2007) 075412. [23] M. Xu, A. Tabarraei, J. Paci, J. Oswald, T. Belytschko, Int. J. Fract. 173 (2012) 163–173. [24] B. Zhang, L. Mei, H. Xiao, Appl. Phys. Lett. 101 (2012) 121915. [25] M. Dewapriya, A.S. Phani, R. Rajapakse, Modell. Simul. Mater. Sci. Eng. 21 (2013) 065017. [26] M. Dewapriya, R. Rajapakse, A. Phani, Int. J. Fract. 187 (2014) 199–212. [27] A. Tabarraei, S. Shadalou, J.-H. Song, Comput. Mater. Sci. 96 (2015) 10–19. [28] G. Jung, Z. Qin, M.J. Buehler, Extreme Mech. Lett. 2 (2015) 52–59. [29] V. Georgakilas, M. Otyepka, A. Bourlinos, V. Chandra, N. Kim, K.C. Kemp, P. Hobza, R. Zboril, K. Kim, Chem. Rev. 112 (2012) 6156–6214. [30] D. Boukhvalov, M. Katsnelson, A. Lichtenstein, Phys. Rev. B 77 (2008) 035427. [31] R. Balog, B. Jørgensen, L. Nilsson, M. Andersen, E. Rienks, M. Bianchi, M. Fanetti, E. Lægsgaard, A. Baraldi, S. Lizzit, Z. Sljivancanin, F. Besenbacher, B. Hammer, T. Pedersen, P. Hofmann, L. Hornekær, Nat. Mater. 9 (2010) 315. [32] S. Casolo, O. Løevvik, R. Martinazzo, G. Tantardini, J. Chem. Phys. 130 (2009) 054704. [33] Q. Pei, Y. Zhang, V. Shenoy, Carbon 48 (2010) 898–904. [34] N. Li, Z.S. abd Q. Pei, Y. Zhang, J. Phys. Chem. C 118 (2014) 13769–13774. [35] Y. Hwangbo, C. Lee, S. Kim, J. Kim, K. Kim, B. Jang, H. Lee, S. Lee, S. Kim, J. Ahn, S. Lee, Sci. Rep. 4 (2014) 4439. [36] S. Stuart, A. Tutein, J. Harrison, J. Chem. Phys. 112 (2000) 073405. [37] R. Grantab, V. Shenoy, R. Ruoff, Science 330 (2010) 946–948. [38] A. Cao, J. Qu, Appl. Phys. Lett. 102 (2013) 071902. [39] H. Zhao, N. Aluru, J. Appl. Phys. 108 (2010) 064321. [40] B. Zhang, G. Yang, H. Xu, Physica B 434 (2014) 145–148. [41] F. Liu, P. Ming, Phys. Rev. B 7 (2002) 064120. [42] M. Xu, J. Paci, J. Oswald, T. Belytschko, Int. J. Solids Struct. 49 (2012) 25822589. [43] T. Anderson, Fracture Mechanics: Fundamentals and Applications, Taylor & Francis Group, 2004. [44] S. Nosé, J. Chem. Phys. 81 (1984) 511. [45] S. Zhang, T. Zhu, T. Belytschko, Phys. Rev. B 76 (2007) 094114. [46] G. Sih, Y. Lee, Theoret. Appl. Fract. Mech. 6 (1986) 129–138. [47] K. Markstrom, B. Storakers, Int. J. Solid Mech. 16 (1980) 217–229. [48] R. Brighenti, Thin-Wall. Struct. 43 (2005) 209–224. [49] T. Zhang, X. Li, H. Gao, Extreme Mech. Lett. 1 (2014) 3–8. [50] X. Huang, H. Yang, A.C. van Duin, K.J. Hsia, S. Zhang, Chemomechanics control of tearing paths in graphene, Phys. Rev. B 85 (19) (2012) 195453.