Molecular dynamics simulations on deformation and fracture of bi-layer graphene with different stacking pattern under tension

Physics Letters A 380 (2016) 609–613

Contents lists available at ScienceDirect

Physics Letters A www.elsevier.com/locate/pla

Molecular dynamics simulations on deformation and fracture of bi-layer graphene with different stacking pattern under tension M.D. Jiao a , L. Wang a , C.Y. Wang b , Q. Zhang a,∗ , S.Y. Ye a , F.Y. Wang a a b

College of Mechanics and Materials, Hohai University, Nanjing 210098, China Faculty of Civil Engineering and Mechanics, Jiangsu University, Zhenjiang 212013, China

a r t i c l e

i n f o

Article history: Received 12 June 2015 Received in revised form 25 October 2015 Accepted 9 November 2015 Available online 27 November 2015 Communicated by R. Wu Keywords: Bi-layer graphene films Stacking pattern Molecular dynamics Tension Mechanical behavior

a b s t r a c t Based on AIREBO (Adaptive Intermolecular Reactive Empirical Bond Order) potential, molecular dynamics simulations (MDs) are performed to study the mechanical behavior of AB- and AA-stacked bilayer graphene films (BGFs) under tension. Stress–strain relationship is established and deformation mechanism is investigated via morphology analysis. It is found that AA-stacked BGFs show wavy folds, i.e. the structural instability, and the local structure of AB-stacked BGFs transforms into AA-stacked ones during free relaxation. The values of the Young’s modulus obtained for AA-stacked zigzag and armchair BGFs are 797.2 GPa and 727.4 GPa, and those of their AB-stacked counterparts are 646.7 GPa and 603.5 GPa, respectively. In comparison with single-layer graphene, low anisotropy is observed for BGFs, especially AB-stacked ones. During the tensile deformation, hexagonal cells at the edge of BGFs are found to transform into pentagonal rings and the number of such defects increases with the rise of tensile strain. © 2015 Elsevier B.V. All rights reserved.

1. Introduction Monolayer graphene is characterized by the honeycomb lattice structure of hybridized SP2 carbon atoms [1–3]. It exhibits superior property (e.g. super strength and extreme stiffness) that is promising for a wide range of applications in nanoelectronics, nanodevices and nanocomposites [4,5]. Graphene thus has attracted considerable attention from researchers worldwide. In particular, the fundamental mechanical property of graphene is a topic of great interest in the areas of both nanomechanics and nanomaterials. In 2008, nano-indentation was performed by Lee et al. [6] to study the mechanical and thermal properties of monolayer graphene, where the Young’s modulus E = 1.0 ± 0.1 TPa (associated with thickness of 0.335 nm), the breaking strength about 40 N/m, and the thermal conductivity 5000 W m−1 K−1 (room temperature) were reported. Through experiments on epoxy composites, Rafiee et al. [7] found the influence of graphene on mechanical property, e.g. elastic modulus, tensile strength, fracture toughness, and resistance to fatigue crack, of composite material is remarkable as compared with single- and multi-walled carbon nanotubes. Using density functional theory, Liu et al. [8] found that the Young’s modulus and Poisson’s ratio for monolayer graphene

*

Corresponding author. E-mail address: [email protected] (Q. Zhang).

http://dx.doi.org/10.1016/j.physleta.2015.11.018 0375-9601/© 2015 Elsevier B.V. All rights reserved.

are 1.05 TPa and 0.186, respectively. Based on orthogonal tightbinding method and MDs, Zhao et al. [9] studied the effects of size and chirality on the elastic property of graphene, where the Young’s modulus of zigzag sheets was found to be larger than that of armchair ones. Similarly, size-dependent mechanical properties were studied to predict the non-linear behavior of graphene nanoribbons [10]. Results show that both linear and non-linear properties are strongly dependent on the structure as well as on the size of the graphene strip tested. By MD simulations, Zhang and Gu [11] investigated the effects of layer number, temperature and isotope on mechanical property of graphene which is found to be sensitive to temperature but insensitive to the layer numbers. Up to now, most studies are focused on monolayer graphene [9,12–19], whereas multilayer graphene films (MGFs) have not received enough attention in spite of the fact that MGFs also hold great promise for various applications in nanotechnology [7,20,21]. Using molecular structural mechanics methods, Hosseini Kordkheili et al. [22] calculated the mechanical properties of double-layer graphene, where non-linear beam and truss elements were introduced to modeling different bonds. In particular, the stacking pattern and the interlayer van der Waals (vdW) interaction may lead to different mechanical property of MGFs from those of monolayer counterpart. Here it is noted that the bi-layer graphene film (BGF) is the simplest case and thus provide an excellent platform for the study of MGF mechanics. BGFs exhibit some prominent attributes, e.g. their optical [20] and electrical properties [21], which enable

610

M.D. Jiao et al. / Physics Letters A 380 (2016) 609–613

Fig. 1. Three typical stacking patterns of BGFs: (a) Mis-oriented; (b) AA-stacked; (c) AB-stacked. Yellow and red atoms represent top-level and underlying graphene respectively. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

their use in field effect transistors. Furthermore, their double-layer structures can substantially reduce the influence of noise in the transistors by overcoming Hooge’s rule which states that as transistors become smaller, the tiny electron charges inside them will threaten to overwhelm a desired signal with unwanted interference [23]. In this paper, MDs are performed to study the free relaxation and deformation mechanism of both AA- and AB-stacked BGFs subject to uniaxial tension. The whole paper is organized as follows: In section 2, geometrical models and simulation details used in the present study are briefly introduced. Results and discussions are included in section 3 where the stress–strain evolutionary relationship of BGFs is achieved, showing the tensile deformation process for BGFs of different stacking patterns. Finally, conclusions are summarized in section 4. 2. Geometrical models and simulation details Two single graphene layers can be stacked in an arbitrary pattern as shown in Fig. 1(a). Such stacking order is referred to as mis-oriented or twisted pattern which is extremely unstable due to its high free energy. Thus, BGFs are usually found in the following two stable stacking modes, i.e. AA- and AB-stacked structures shown in Figs. 1(b) and 1(c), respectively. In the AA-stacked structure shown in Fig. 1(b), the vertical projections of the corresponding carbon atoms on the two layers completely coincide with each other. While for the AB-stacked structure shown in Fig. 1(c), half carbon atoms of the top layer sit directly on the top of the bottom layer atoms while the other half sit at the centers of the hexagonal rings in the bottom layer. Both AA- and AB-stacked structures are stable due to their relatively low energy and thus are the focus of some research [23]. In nanomechanics, MDs is considered as a standard technique that can be efficiently used to quantify the mechanical property and capture essential characteristics of nanostructures and materials. The basic principle of MD is to work out the force acting on individual atoms through the atom–atom interaction potential, and calculate their trajectory, acceleration and velocity by solving Newtonian mechanics equations. Here a certain number of molecules (atoms) under initial conditions are considered in the calculation which will be done with selected time step and boundary conditions. Subsequently, the statistical average of the obtained results will be outputted, from which the required macroscopic physical and mechanical quantities can be obtained. Also, it is a key issue to select a proper atomic potential function which determines the simulation accuracy for the calculation of MDs. AIREBO (adaptive intermolecular reactive empirical bond order) potential function [24] is mainly for carbons atoms and used to describe the atomic interaction within the same layer as follows:



E=

1  2

i

LJ

E iRB j + Eij +

i = j

 

 TORS E ki jl

(1)

k= j l=i , j ,k

where the second and third items in bracket represent a term for long-range interaction and a twisting term relying on dihedral angle, respectively. The interlayer cohesion maintained by vdW force is characterized by 12-6 Lennard-Jones potential [25]:

 12 V (r ) = 4ε

σ r

 6  −

σ r

(2)

where ε = 2.968 meV and σ = 0.3407 nm. In the present simulations, a uniformly distributed tensile force is applied on one side of rectangular BGFs in the perpendicular direction and the atoms on the opposite side are fixed. The geometrical model is square with length of 7.0 nm and the total number of atoms is 3944. Periodic boundary conditions are imposed on the two sides along horizontal direction, as shown in Fig. 2. Nose–Hoover method [26,27] is applied for isothermal adjustment at the temperature of 0.001 K. The low temperature is considered in order to avoid complex changes of mechanical property caused by atomic thermal activation. During the simulation process, unconstrained relaxation is first carried out on initial configuration to keep the system in an equilibrium state with the lowest energy. Then uniform tensile deformation is conducted on the BGFs. Here a time step 1.0 fs is selected and an increment of tensile strain 0.001 is imposed, followed by 1000 steps of relaxation for 1.0 ps. The calculated strain rate is of 1.0 × 109 S−1 . The stretching and relaxation process will be repeated to deform the BGFs. All simulations in this paper were conducted with the LAMMPS code developed by Sandia National Laboratory in US [28]. Young’s modulus, strain and stress are traditionally defined in the framework of continuum mechanics. In the present work, these concepts will be extended to discrete BGFs consisting of two oneatom thick sheets. The equivalent thickness of graphene needs to be defined in an appropriate manner. The thickness of monolayer graphene is assumed to be 0.335 nm which is the interlayer spacing of multilayer graphene. Following this way, here we define the equivalent thickness of BGFs as d = 0.670 nm. 3. Results and discussion To validate the MDs program described in section 2, we first studied the tensile behavior of a monolayer graphene and obtained its Young’s modulus and tensile strength. Assuming the thickness of single-layer graphene is 0.335 nm and based on the stress– strain curves achieved for the sample graphene sheet, the obtained value of Young’s modulus is 897.0 GPa, which is in agreement with available theoretical calculations of 0.8–1.1 TPa [6,29,30]. Tensile strength and fracture strain obtained for zigzag graphene are

M.D. Jiao et al. / Physics Letters A 380 (2016) 609–613

611

Fig. 2. Boundary condition and force application: (a) Zigzag AA-stacked BGFs; (b) Armchair AA-stacked BGFs.

Fig. 3. BGFs after free relaxation: (a) Top view of AA-stacked BGFs; (b) Top view of AB-stacked BGFs; (c) Front view of AA-stacked BGFs; (d) Front view of AB-stacked BGFs.

142.7 GPa and 0.35 respectively, which are also close to the results 130–180 GPa and 0.3428 reported in [6,31]. The good agreements achieved confirm that the simulation details in the present study are reliable in studying uniaxial stretching of graphene sheets. 3.1. Morphology analysis of BGFs after free relaxation In the process of free relaxation, small corrugated wrinkles are found at the edge of AA-stacked BGFs. The wavelength of the wrinkles and their transverse deformation then grow gradually as time elapses. Wrinkles thus become more and more pronounced, as shown in Fig. 3(c). The BGFs with rippled edges finally reach a stable equilibrium state associated with the lowest energy. Such intrinsic wrinkles of AA-stacked BGFs observed in the present study are consistent with the TEM observations by Meyer et al. [2,32] and Ishigami et al. [4], and the Monte Carlo simulation of Fasolino et al. [33]. However, in the existing experiments, wrinkles occur not only at the edge of the BGFs but also in their internal area. Such discrepancy between present MDs and previous experiments can possibly be attributed to the different boundary constraints imposed on the BGFs. Indeed, during the free relaxation the four

sides of the BGFs are free, while in the experiments the tested BGFs have to be placed on a substrate where the vdW interaction can firmly hold the BGFs and thus change the stress state in their internal area. The relaxation process of AB-stacked BGFs is substantially different from their AA-stacked counterparts as no wrinkle is observed as shown in Fig. 3(d). Interestingly, it is noticed that in some areas, stocking order transforms from AB-stacked to the stacking pattern quite close to AA-stacking. 3.2. Tensile behavior of BGFs In the tensile deformation of BGFs, the stress–strain relationships are established, where important material property, e.g. fracture strain and elastic modulus, can be obtained. Since AA- and AB-staked BGFs show similar stress–strain curves, only the curves of AB-stacking pattern BGFs are plotted in Fig. 4, where two different atomic structures are considered, i.e. zigzag and armchair. As shown in Fig. 4, the stress–strain curve of BGFs can be divided into four stages, i.e. linear elastic deformation, yielding, hardening and fracture of BGFs, which are analogous to the tensile behaviors of low-carbon steel bar. For AB-stacked BGFs, the

612

M.D. Jiao et al. / Physics Letters A 380 (2016) 609–613

that the anisotropy of multilayer graphene depends sensitively on the number of the graphene layers and their stacking patterns. In Table 1, the values of Young’s modulus obtained for zigzag and armchair AB-stacked BGFs are 646.7 GPa and 603.5 GPa, while the corresponding values of AA-stacked BGFs are 797.2 GPa and 727.4 GPa. These values are lower than 897 GPa obtained in the present study for monolayer graphene. From here it follows that the Young’s modulus of graphene decreases with the increase of the layer number. This trend in the present MDs is found to be in agreement with experimental data [36], where the elastic modulus is 891.0 GPa for monolayer graphene and it reduces to 393.0 GPa for BGFs. The lower experimental value of BGFs is possibly due to the defects (e.g. SP3 hybrid orbitals) formed during the fabrication of BGFs [6]. It is noted that the dependency of the Young’s modulus of BGFs on the number of layers is qualitatively similar to the trend found in multi-walled carbon nanotubes [37]. It is believed that the interlayer vdW interaction and the SP3 hybrid orbitals formed between adjacent layers contribute to the decreasing elastic modulus of graphene with increasing number of layers.

Fig. 4. Stress–strain curves of AB-stacked BGFs. Table 1 Comparison on Young’s modulus of single and bi-layer graphene with different stacking patterns. Graphene type

Young’s modulus (GPa) Zigzag

Armchair

Single-layer [31] Single-layer AA-stacked AB-stacked

1160.0 897.0 797.2 646.7

1050.0 808.6 727.4 603.5

Differences (%)

10.5 11.0 9 .6 7 .2

Table 2 Comparison on ultimate stress of single and bi-layer graphene with different stacking patterns. Graphene type

Ultimate stress (GPa) Zigzag

Armchair

Single-layer [31] Single-layer AA-stacked AB-stacked

180 142.7 149.5 200.3

210 170.6 173.6 230.5

Differences (%)

16.7 19.6 16.1 13.3

fracture strength of armchair BGF is 230.5 GPa and zigzag counterpart 200.3 GPa. The corresponding values obtained for AA-stacked BGFs in the present study are 173.6 GPa and 149.5 GPa, respectively. Thus, the strength of BGFs is found to be higher than that of monolayer graphene, which is due to the coupling and cohesion between two constituent films via vdW interaction. This theory is also used to explain the higher strength of multi-walled carbon nanotubes relative to that of single-walled carbon nanotubes [34]. It is also noted in Fig. 4 that the fracture strain associated with zigzag structure is 0.541, which is greater than 0.35–0.38 reported for monolayer graphene in this study and Ref. [30]. Furthermore, the anisotropy is also observed for BGFs in the present study. Similar observation was also reported for monolayer graphene previously [35]. In Tables 1 and 2, comparison has been made among the graphene in terms of their anisotropic property. From Tables 1 and 2, it can be seen that, as compared with monolayer graphene, low anisotropy is observed for BGFs, especially for AB-stacked ones. Here the vdW interaction between the two layers should be responsible for the low anisotropy of the BGFs. Specifically, the top view of AB-stacked BGFs (Fig. 1(c)) show that the atoms are more densely located in the system than their AA-staked counterparts, which leads to even stronger interlayer vdW interaction and thus further decreases the anisotropy of the AB-stacked BGFs. Based on the above results, it may be concluded

3.3. Atomic structure changes In this section, further study is conducted on both the structural changes and fracture process of BGFs under tension. To this end, snapshots in Fig. 5 show the structural changes of AA- and AB-stacked BGFs as the tensile strain increases. It can be seen from Fig. 5 that, in the tensile deformation of AA-stacked BGFs, elongation of hexagonal carbon atom rings is observed at small axial strain, i.e. linear elastic stage. When axial strain exceeds the elastic limit, the hexagonal carbon rings at the edges turn out to be irregular (right panel in Fig. 5(b)). With further increase of the axial strain to a certain value, some of the carbon–carbon bonds are found to break, showing the starting point of fracture process. Differently, in the same process, ABstacked BGFs first changes to similar AA-stacked pattern in three places. Along with the development of stretching, the area that firstly converts into AA-stacked is the place where tensile failure occurs. For AA-stacked BGFs, with the rise of tensile strain, folds along the stretching direction are gradually flattened while those are perpendicular to the stretching direction do not appear similar phenomenon. With the increasing of strain and the changing of the distance between layers, the folding phenomenon is more obvious, this can be attributed to the disappearance of folds in one direction which releases some strain energy and enhance the folds in orthogonal direction [38]. The stretching process also shows that folds similar to corrugated structure also appear in AB-stacked BGFs, but the fold degree is lower than that of AA-stacked counterparts, as shown in Fig. 6. 4. Conclusions MDs are performed to study the uniaxial tensile behavior of both AA- and AB-stacked BGFs, from which the stress–strain evolving curve is achieved and the deformation and fracture mechanisms are investigated via morphology analysis. Simulations results show that: AA-stacked BGFs do not have perfect planar structure after free relaxation and the surface is not completely flat. Instead, folding phenomenon occurs. For AB-stacked structure, stacking pattern transforms into the one similar to AA-stacking and folding phenomenon is not observed. MDs results show that the values of Young’s modulus for AAstacked armchair and zigzag BGFs are 797.2 GPa and 727.4 GPa respectively. While those of AB-stacked armchair and zigzag counterparts are 646.7 GPa and 603.5 GPa respectively. In comparison with single-layer graphene, low anisotropy is observed for

M.D. Jiao et al. / Physics Letters A 380 (2016) 609–613

613

Fig. 5. Structural changes under tension: (a) Linear elastic deformation of AA-stacked BGFs; (b) Hardening of AA-stacked BGFs; (c) Linear elastic deformation of AB-stacked BGFs; (d) Hardening of AB-stacked BGFs and amplified local area.

Fig. 6. Front view of different stacked BGFs under tension: (a) AA-stacked BGFs; (b) AB-stacked BGFs.

BGFs, especially AB-stacked ones. During the stretching process, the hexagonal cells at edges of AA-stacked BGFs firstly transform into pentagon cellular. Under the action of tensile load, some local areas transform from AB-stacked to AA-stacked, with analogous AA-stacked being the initial position of destruction. Acknowledgements This work is supported by the National Natural Science Foundation of China (Nos. 11102058, 11472098 and 11372099) and the Program for New Century Excellent Talents in University of China (NCET-0773). References [1] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, et al., Science 306 (2004) 666. [2] J.C. Meyer, A.K. Geim, M.I. Katsnelson, K.S. Novoselov, T.J. Booth, S. Roth, Nature 446 (2007) 60. [3] A.K. Geim, K.S. Novoselov, Nat. Mater. 6 (2007) 183. [4] M. Ishigami, J.H. Chen, W.G. Cullen, M.S. Fuhrer, E.D. Williams, Nano Lett. 7 (2007) 1643. [5] F. Schedin, A.K. Geim, E.W. Hill, P. Blake, M.I. Katsnelson, K.S. Novoselov, Nat. Mater. 6 (2007) 652. [6] C. Lee, X. Wei, J.W. Kysar, J. Hone, Science 321 (2008) 385. [7] M.A. Rafiee, J. Rafiee, Z. Wang, H. Song, Z.Z. Yu, N. Koratkar, ACS Nano 3 (2009) 3884. [8] F. Liu, P.B. Ming, J. Li, Phys. Rev. B 76 (2007) 064120. [9] H. Zhao, K. Min, N.R. Aluru, Nano Lett. 9 (2009) 3012. [10] S.K. Georgantzinos, G.I. Giannopoulos, D.E. Katsareas, P.A. Kakavas, N.K. Anifantis, Comput. Mater. Sci. 50 (2011) 2057. [11] Y.Y. Zhang, Y.T. Gu, Comput. Mater. Sci. 71 (2013) 197.

[12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38]

P.A. Thrower, Chem. Phys. 5 (1969) 217. A.J. Stone, D.J. Wales, Chem. Phys. Lett. 128 (1986) 501. Q.B. Zheng, Y. Geng, S.J. Wang, Z.G. Li, J.K. Kim, Carbon 48 (2010) 4315. M. Neek-Amal, F.M. Peeters, Appl. Phys. Lett. 97 (2010) 153118. M. Neek-Amal, F.M. Peeters, Phys. Rev. B 82 (2010) 085432. W.H. Duan, C.M. Wang, Nanotechnology 20 (2009) 075702. W.H. Duan, C.M. Wang, W.X. Tang, J. Appl. Phys. 107 (2010) 074303. S.S. Gupta, R.C. Batra, J. Comput. Theor. Nanosci. 7 (2010) 2151. D.S.L. Abergel, V.I. Falko, Phys. Rev. B 75 (2007) 155430. T. Ohta, A. Bostwick, T. Seyller, K. Horn, E. Rotenberg, Science 313 (2006) 951. S.A. Hosseini Kordkheili, H. Moshrefzadeh-Sani, Comput. Mater. Sci. 69 (2013) 335. Y.M. Lin, P. Avouris, Nano Lett. 8 (2008) 2119. D.W. Brenner, O.A. Shenderova, J.A. Harrison, S.J. Stuart, B. Ni, S.B. Sinnott, J. Phys. Condens. Matter 14 (2002) 783. R. Saito, R. Matsuo, T. Kimura, G. Dresslhaus, M.S. Dresslhaus, Chem. Phys. Lett. 348 (2001) 187. S. Nosé, Mol. Phys. 52 (1984) 255. W.G. Hoover, Phys. Rev. A 31 (1984) 1695. S.J. Plimpton, J. Comput. Phys. 117 (1995) 1. C. Lee, X. Wei, Q. Li, R. Carpick, J.W. Kysar, Phys. Status Solidi B 246 (2009) 2562. I.W. Frank, D.M. Tanenbaum, A.M. van der Zande, P.L. McEuen, J. Vac. Sci. Technol. B 25 (2007) 2558. Z.H. Ni, H. Bu, M. Zou, H. Yi, K. Bi, Y.F. Chen, Science 405 (2010) 1301. J.C. Meyer, A.K. Geim, M.I. Katsnelson, K.S. Novoselov, D. Obergfelle, S. Rothe, C. Girit, A. Zettl, Solid State Commun. 143 (2007) 101. A. Fasolino, J.H. Los, M.I. Katsnelson, Nat. Mater. 6 (2007) 858. E.M. Byme, M.A. McCarthy, Z. Xia, W.A. Curtin, Phys. Rev. Lett. 49 (2009) 04502. Y. Gao, P. Hao, Physica E 41 (2009) 1561. Y.P. Zhang, C.X. Pan, Diam. Relat. Mater. 24 (2012) 1. B. Peng, M. Locascio, P. Zapol, S.Y. Li, S.L. mielke, G.C. Schatz, H.D. Espinosa, Nat. Nanotechnol. 3 (2008) 626. T. Zhang, X.Y. Li, H.J. Gao, J. Mech. Phys. Solids 67 (2014) 2.