- Email: [email protected]
Thickness and chirality effects on tensile behavior of few-layer graphene by molecular dynamics simulations
Computational Materials Science 53 (2012) 298–302
Contents lists available at SciVerse ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Thickness and chirality effects on tensile behavior of few-layer graphene by molecular dynamics simulations Bohayra Mortazavi a,b, Yves Rémond a, Said Ahzi a,⇑, Valérie Toniazzo b a b
Institut de Mécanique des Fluides et des Solides, University of Strasbourg/CNRS, 2 Rue Boussingault, 67000 Strasbourg, France Centre de Recherche Public Henri Tudor, Department of Advanced Materials and Structures, 66 Rue de Luxembourg, BP 144, L-4002 Esch/Alzette, Luxembourg
a r t i c l e
i n f o
Article history: Received 27 June 2011 Accepted 15 August 2011 Available online 22 October 2011 Keywords: Molecular dynamics Few-layer graphene Chirality Uniaxial tension Brittle
a b s t r a c t The mechanical response of few-layer graphene (FLG), consisting of 2–7 atomic planes and bulk graphite is investigated by means of molecular dynamics simulations. By performing uniaxial tension tests at room temperature, the effects of number of atomic planes and chirality angle on the stress–strain response and deformation behavior of FLG were studied using the Tersoff potential. It was observed that by increasing of the FLG number of layers, the increase of bonding strength between neighboring layers reduce the elastic modulus and ultimate strength. It was found that, while the chirality angle of FLG showed a significant effect on the elastic modulus and ultimate tensile strength of two and three graphene layers, it turns to be less significant when the numbers of layers are more than four. Finally, by plotting the deformation behavior, it was concluded that FLGs present brittle performance. Ó 2011 Elsevier B.V. All rights reserved.
1. Introduction A two dimensional structure of carbon atoms, graphene [1] is an ongoing research interest in low dimensional materials due to its exceptionally high crystal and electronic quality [2–4]. As a member of carbon family, graphene materials present extraordinary thermal [5], mechanical [6] and electrical [7] properties. In this way, the graphene is considered as an excellent candidate for building block in nano-electro mechanical systems (NEMS) [8,9], nanoelectronics [10,11] and nanocomposites [12,13]. Thus, the fundamental understanding of graphene properties is the critical issue in its future applications. One of the key properties of graphene is its mechanical properties that have considerable effects on its application for the NEMS devices and reinforcement part in polymer nanocomposites. Among the limited available experimental studies, by using atomic force microscope nanoindentation experiment, Frank et al. [14] reported the elastic modulus of stacks of graphene sheets as 0.5 TPa. Using the same experimental technique, Lee et al. [6] acquired the elastic modulus of 1 TPa and breaking strength of 0.13 TPa for monolayer and free of defect graphene sheets. Since the experimental characterization techniques at nanoscale are considerably complex, today the numerical simulations are gaining conspicuous interest as an alternative for experimental studies [15]. Molecular dynamics (MD) is one of the most popular techniques that have been used for atomic-scale
⇑ Corresponding author. Tel.: +33 (0) 3 68 85 29 52; fax: +33 (0) 3 68 85 29 36. E-mail address: [email protected] (S. Ahzi). 0927-0256/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2011.08.018
characterization of material properties. MD simulation computes the motion of individual molecules in the models of solid, liquid, or gas. The key idea here is the motion [16], then physical and chemical properties can be evaluated based on the statistical mechanics theories. There have been numerous MD studies on the mechanical properties of bulk and nanoribbons graphene. Bu et al. [17] and Ni et al. [18] studied the mechanical behavior of graphene nanoribbons and sheets respectively by MD simulations. Tsai and Tu [19] calculated the mechanical properties of single graphene layer and bulk graphite by MD simulation. One of the alternatives of MD simulations in evaluation of mechanical properties of carbon nanotubes and sheets is structural mechanics modeling. In the structural mechanics, the atoms are considered as balls and their interaction force fields are represented by beams or spring like finite elements. Recently, Georgantzinos et al. [20] studied the size-dependent non-linear mechanical properties of graphene nanoribbons. Furthermore, Sakhaee-Pour [21] and Scarpa et al. [22] used structural mechanics modeling techniques for the evaluation of elastic properties of the single layer graphene sheets. Georgantzinos et al. [23] also studied the single and multilayer graphene structures mechanical properties by structural mechanics modeling. The main privilege of structural mechanics is its certainty and computational efficiency in the evaluation of elastic constants. On the other hand MD seems to be a more versatile technique for assessment of large deformations, failure mechanism and temperature effects. In the current study, uniaxial MD tension tests were carried out in order to investigate the directional and thickness effect on deformation and mechanical response of FLG sheets using Tersoff potential at room temperatures.
B. Mortazavi et al. / Computational Materials Science 53 (2012) 298–302
2. Molecular dynamics modeling LAMMPS (Large-Scale Atomic/Molecular Massively Parallel Simulator) [24], developed at Sandia National Laboratories, is used for molecular dynamics simulations. The carbon–carbon bonded interaction is modeled by Tersoff potential [25,26]. In case of nonbonded van der Waals interaction between individual carbon atoms, there is a variety of used potentials in the literature. Commonly, the nonbonded interactions are expressed as Lennard– Jones (LJ) or Morse potentials. In the current study, the LJ potential was used which is expressed by:
12 6 r r /ðrÞ ¼ 4e r r
ð1Þ
where r refers to the distance between carbon atoms at different atomic layers. e is the depth of the potential well and r is the equilibrium distance at which the inter-particle potential is zero. A parameterization frequently employed for LJ interactions in the modeling of multi-layer carbon structures is e = 2.4 m eV and r = 3.4 Å which is based on the compressibility and graphene layer spacing [27] which also satisfy the density of the graphite (2.2 g/ cm3). As a common assumption for the LJ potentials, interactions were ignored at distances higher than 2.5r. As a result of the van der Waals interactions between carbons atoms in different layers, in the plane properties of FLG, are expected to be dependent of the number of layers. The MD simulation model for uniaxial tensile test of FLG with the chirality angle of 10° is illustrated in Fig. 1. The chirality angle of 0° and 30° are commonly called armchair and zigzag directions, respectively. A FLG is constructed by putting graphene sheets up together by a spacing distance of 3.4 Å. The periodic boundary conditions were applied in x and y directions, thus these simulations are representative of an infinite FLG sheet and not the FLG nanoribbons. The properties of bulk graphite were also evaluated by considering eight layers of graphene sheets while the periodic boundary condition was also applied in the graphene sheets normal direction (z direction). Due to high computational cost of MD simulations, the tensile tests are usually performed at high strain rates, in the order of 109 s1. The limitation in approaching lower strain rates is mainly due to the small interatomic distances and high velocities of atoms. In this way, to accurately count for the deformations of FLG at considerably high applying strain rates, small time increment of 0.5 fs (=0.5 1015 s) was used. Prior to applying the loading condition, the FLG structure was left to anisotropically relax to zero pressure in x and y directions at a temperature of 300 K for 30,000 time increments using constant pressure–temperature (i.e. NPT
299
ensemble) by means of Nosé–Hoover barostat and thermostat method. The loading condition was applied by extending the periodic simulation box size in the x direction under constant strain rates. In order to satisfy the uniform uniaxial stress conditions, the periodic boundary condition was applied in the loading direction and the simulation box in the y direction was altered to reach zero stress in this direction using NPT Nosé–Hoover method. Virial stresses in the x direction (uniaxial stress) were calculated at each strain level. The engineering strain at each time step was calculated by multiplying of total time of loading simulation by the applied strain rate. Using the Hooke’s law, the slope of the initial linear part of uniaxial stress–strain curve is equivalent to the elastic modulus of FLG along the loading direction. 3. Results and discussions Since the applied strain rates are considerably high, the simulations were performed at different strain rate to verify the level of strain rate dependency of the obtained properties. In Fig. 2, the acquired engineering stress–strain diagrams of three layers zigzag graphene and five layers armchair graphene at different strain rates (SR) are plotted. At low strain levels (up to 0.05), a linear relation between stress and strain values is obtained. This linear response is commonly assumed as the linear elasticity part, in which the corresponding slope is equivalent to the elastic modulus. Based on the Fig. 2 observations, the elastic modulus of the armchair and zigzag FLG are reasonably independent of the strain rate values. At higher strain level, the stress–strain response presents a nonlinear behavior which continues up to the maximum tensile strength, the point in which rupture occurs in the first atomic plane. In case of zigzag FLG, the nonlinear part is also found to be acceptably independent of strain rate where the rupture occurs at approximately the same stress and strain levels. On the other hand, while the nonlinear part of armchair FLG stress–strain curves also present acceptable level of strain rate independency up to the tensile strength. However, the strains at rupture decreases as the strain rate decreases. Furthermore, by decreasing of the loading strain rate in both armchair and zigzag loading directions, the rupture is accompanied with sharper decline in the stress–strain diagram which consequently means more brittle behavior. Based on the obtained results, the stain rate of 1 109 s1 is selected for the further MD simulations in this study, because of its computational efficiency as well as predicting a stress–strain trend, closer to the lower strain rate of 1 108 s1. In Fig. 3a, the engineering stress–strain response of zigzag FLG containing different number of layers of 2, 3, 6 and bulk graphite
Fig. 1. MD model of uniaxial tension test of FLG with chirality angle of 10° in respect to armchair direction. The periodic boundary conditions were applied in the x and y direction. The loading condition was applies by the extending the periodic simulation box in the x direction.
300
B. Mortazavi et al. / Computational Materials Science 53 (2012) 298–302
Fig. 2. Uniaxial engineering stress–strain relation of three layers zigzag graphene (dashed line) and five layers armchair graphene (continuous line) at different loading strain rates (SR).
Fig. 4. (a) Uniaxial engineering stress–strain relation of three layers FLG as function of different loading chirality angle (ChA) directions and (b) the effect of loading direction on the elastic modulus for FLGs with different number of layers and bulk graphite.
Fig. 3. (a) Uniaxial engineering stress–strain curve of zigzag FLG for different number of layers and bulk graphite and (b) the dependence of FLG number of layers on the elastic modulus for armchair and zigzag loading directions.
(containing eight layers in which periodic boundaries are also applied along the thickness) is presented. As an overall trend, by increasing of the number of atomic layers the elastic modulus and ultimate tensile strength decrease (see Fig. 3a and b). Ni
et al. [18] reported the elastic modulus of around 1.1 TPa for a single layer graphene sheet. In our case, due to the high level of uncertainty for the results of a single layer graphene, we are not reporting any value for this case. However, if we take into account the value of 1.1 TPa for the single layer case, we can conclude that the trend of the elastic modulus variation as function of the number of atomic layers of grapheme (see Fig. 3b) is similar to the trend reported by Liew et al. [28] for multi-walled carbon nanotubes (CNTs). These authors [28] found that the double wall CNTs present higher elastic modulus and ultimate strength than the single layer CNTs, and by increasing of CNTs number of walls, both of these mechanical properties decline. As it can be observed from Fig. 3a in case of zigzag FLGs, the ultimate tensile strength, where the ruptures occur in the first atomic layer, is obtained at approximately the same strain values for all zigzag FLG cases. Since the subsequent ruptures at adjacent layers may not occur at the same time, series of stress fall and raise can be observed after the ultimate tensile strength point. In Fig. 3b the effect of FLG number of layers on the elastic modulus for zigzag and armchair loading directions are shown. In the current study the elastic modulus of bulk graphite was found to be around 0.98 TPa, which is in fine agreement with the experimental ranges of 1.02 ± 0.03 TPa [29]. Based on the obtained results, it can be concluded that as the number of FLG layers reach the number of six, the obtained elastic properties become acceptably independent of number of layers, where the elastic modulus is close to the ranges of experimental bulk properties. This matter can be interesting in the modeling of polymer
B. Mortazavi et al. / Computational Materials Science 53 (2012) 298–302
301
Fig. 5. The brittle fracture process at different engineering strain level of (a–c) two layers zigzag graphene and (d–f) three layers graphene with chirality angle of around 10°. The initial void shapes in (a) and (c) insets and the fracture line in (b) and (f) insets.
nanocomposites, reinforced by multi-layer graphene. Thus, as an accurate assumption in the modeling of the graphene with higher number of layers than six, we can use the graphite bulk properties. It is worth mentioning that, in-plane thermal conductivity of FLG was experimentally proven to become independent of the number of layers as the number of layers is more than four [30]. Fig. 4a illustrates the effects of loading direction chirality angle (ChA) on the engineering stress–strain response of three layers FLG. As an overall trend, both elastic modulus and ultimate strength of FLG increase by increasing of the chirality angle from 0° (Armchair) to 30° (Zigzag). As presented in Fig. 4b, the chirality dependency of mechanical properties of FLG is more obvious for two and three layers. For the case of FLG with higher number of layers than four, the elastic modulus dependency on chirality reduces considerably. It was also found that in case of bulk graphite the elastic modulus is independent of loading direction. However it was found that the ultimate tensile strength of armchair is still less than zigzag direction by around 15%.
In Fig. 5, the deformation behavior of two layers zigzag graphene (a–c) and three layers graphene with chirality angle of around 10° (d–f) is shown. It was found, that up to the high strain levels (about 30% strain) the specimen extends uniformly along the loading direction. At a point, very close to the ultimate strength, several voids are formed along one of the graphene sheets. These formed voids have the shapes illustrated in Fig. 5a,c insets. The formations of these voids take place during the short time, which results in fast subsequent failure in the first atomic sheet while the remaining sheets are undamaged. The beginning of voids formation and subsequent failure occur at considerably close strain levels. This observation is suggesting that the deformation mechanism of FLG present a high level of void sensitivity; which accordingly means brittle failure behavior. Furthermore, as the rupture takes place in each graphene sheet, the sheets tend to retrieve their original structure and contract rapidly (see Fig. 5d,e – red sheet). However as shown in Fig. 5c inset, the initially formed voids are stabilized in the structure. This damage/failure mechanism extends from one
302
B. Mortazavi et al. / Computational Materials Science 53 (2012) 298–302
sheet to the other and the final rupture of the FLG occurs at the failure of the last sheet. Moreover, as presented in Fig. 5b,f insets, the cut edge of the graphene sheets are often along the zigzag direction which reveals that the debonding between the carbon atoms happened along the armchair direction. This observation is in agreement with the fact that the ultimate tensile strength of FLG and mono layer graphene sheets [18] along the armchair directions falls lower than that in the zigzag direction.
4. Conclusion Molecular dynamics simulations of uniaxial tension were performed in order to investigate the stress–strain and failure behavior of few-layer graphene. The simulations were carried out at different strain rates ranging between 1010 and 108 s1 and it was observed that both linear and nonlinear part of the stress strain response of zigzag FLG are convincingly strain rate independent and only the strain at failures of armchair FLG decreases by decreasing the loading rate. The modeling technique was validated by comparison of the obtained elastic modulus of the simulation and the bulk graphite experiments. The simulated results show that the properties are sensitive to the number of layers and loading direction as well. It was observed that by increasing the number of graphene layers the elastic modulus and ultimate tensile strength dropped. However, the obtained elastic modulus for FLG containing six atomic planes is close to the bulk graphite values. By performing the simulations at various loading directions, it was found that the directional sensitivity of the elastic modulus is considerable only for the FLG with two and three layers. We also found that the ultimate tensile strength increases when the loading direction chirality angle goes from 0° to 30°. Finally by plotting the deformation process of the FLG, it was concluded that the FLG presented a brittle behavior and the debonding between carbon atoms occurred, frequently, along the armchair direction. Acknowledgment The authors acknowledge the financial support of the FNR of Luxembourg via the AFR Grants (PHD-09-016). The first author
greatly acknowledges Michael Essa from University of Strasbourg for his valuable helps. References [1] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov, Science 306 (2004) 666. [2] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, M.I. Katsnelson, I.V. Grigorieva, S.V. Dubonos, A.A. Firsov, Nature 438 (2005) 197. [3] Y.B. Zhang, Y.W. Tan, H.L. Stormer, P. Kim, Nature 438 (2005) 201. [4] A.K. Geim, K.S. Novoselov, Nat. Mater. 6 (2007) 183. [5] A.A. Balandin, S. Ghosh, W.Z. Bao, I. Calizo, D. Teweldebrhan, F. Miao, et al., Nano. Lett. 8 (2008) 902. [6] C. Lee, X. Wei, J.W. Kysar, J. Hone, Science 321 (2008) 385. [7] J.R. Williams, L. DiCarlo, C.M. Marcus, Science 317 (2007) 638. [8] J.S. Bunch, A.M. VanderZande, S.S. Verbridge, I.W. Frank, D.M. Tanenbaum, J.M. Parpia, H.G. Craighead, P.L. McEuen, Science 315 (2007) 490. [9] B. Standley, W.Z. Bao, H. Zhang, J. Bruck, C.N. Lau, M. Bockrath, Nano. Lett. 8 (2008) 3345. [10] D.A. Areshkin, C.T. White, Nano. Lett. 7 (2007) 3253. [11] R.M. Westervelt, Science 320 (2008) 324. [12] 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. [13] T. Ramanathan, A.A. Abdala, S. Stankovich, D.A. Dikin, M. Herrera-Alonso, R.D. Piner, D.H. Adamson, H.C. Schniepp, X. Chen, R.S. Ruoff, S.T. Nguyen, I.A. Aksay, R.K. Prud’homme, L.C. Brinson, Nat. Nanotechnol. 3 (2008) 327. [14] I.W. Frank, D.M. Tanenbaum, A.M. Van der Zande, P.L. McEuen, J. Vac. Sci. Technol. B 25 (2007) 2558. [15] B Mortazavi, A.A. Khatibi, C. Politis, J. Comput. Theor. Nanosci. 6 (2009) 644. [16] J.M. Haile, Molecular Dynamic Simulation, John Wiley & Sons, New York, 1992. [17] H. Bu, Y. Chen, M. Zou, H. Yia, K. Bi, Z. Ni, Phys. Lett. A 373 (2009) 3359. [18] Z. Ni, H. Bu, M. Zou, H. Yi, K. Bi, Y. Chen, Physica B 405 (2010) 1301. [19] J.L. Tsai, J.F. Tu, Mater. Des. 31 (2010) 194. [20] S.K. Georgantzinos, G.I. Giannopoulos, D.E. Katsareas, P.A. Kakavas, N.K. Anifantis, Comput. Mater. Sci. 50 (2011) 2057. [21] A. Sakhaee-Pour, Solid State Commun. 149 (2009) 91. [22] F. Scarpa, S. Adhikari, A. Srikantha Phani, Nanotechnology 20 (2009) 065709. [23] S.K. Georgantzinos, G.I. Giannopoulos, N.K. Anifantis, Mater. Des. 31 (2010) 4646. [24] S. Plimpton, J. Comput. Phys. 117 (1995) 1. [25] J. Tersoff, Phys. Rev. B 37 (1988) 6991. [26] J. Tersoff, Phys. Rev. Lett. 61 (1988) 2879. [27] L.A. Girifalco, M. Hodak, R.S. Lee, Phys. Rev. B 62 (2000) 13104. [28] K.M. Liew, X.Q. He, C.H. Wong, Acta Mater. 52 (2004) 2521. [29] O.L. Blakslee, D.G. Proctor, E.J. Seldin, G.B. Spence, T. Weng, J. Appl. Phys. 41 (1970) 3373. [30] S. Ghosh, W. Bao, D.L. Nika, S. Subrina, E.P. Pokatilov, C.N. Lau, A.A. Balandin, Nat. Mater. 9 (2010) 555.
Contents lists available at SciVerse ScienceDirect
Computational Materials Science journal homepage: www.elsevier.com/locate/commatsci
Thickness and chirality effects on tensile behavior of few-layer graphene by molecular dynamics simulations Bohayra Mortazavi a,b, Yves Rémond a, Said Ahzi a,⇑, Valérie Toniazzo b a b
Institut de Mécanique des Fluides et des Solides, University of Strasbourg/CNRS, 2 Rue Boussingault, 67000 Strasbourg, France Centre de Recherche Public Henri Tudor, Department of Advanced Materials and Structures, 66 Rue de Luxembourg, BP 144, L-4002 Esch/Alzette, Luxembourg
a r t i c l e
i n f o
Article history: Received 27 June 2011 Accepted 15 August 2011 Available online 22 October 2011 Keywords: Molecular dynamics Few-layer graphene Chirality Uniaxial tension Brittle
a b s t r a c t The mechanical response of few-layer graphene (FLG), consisting of 2–7 atomic planes and bulk graphite is investigated by means of molecular dynamics simulations. By performing uniaxial tension tests at room temperature, the effects of number of atomic planes and chirality angle on the stress–strain response and deformation behavior of FLG were studied using the Tersoff potential. It was observed that by increasing of the FLG number of layers, the increase of bonding strength between neighboring layers reduce the elastic modulus and ultimate strength. It was found that, while the chirality angle of FLG showed a significant effect on the elastic modulus and ultimate tensile strength of two and three graphene layers, it turns to be less significant when the numbers of layers are more than four. Finally, by plotting the deformation behavior, it was concluded that FLGs present brittle performance. Ó 2011 Elsevier B.V. All rights reserved.
1. Introduction A two dimensional structure of carbon atoms, graphene [1] is an ongoing research interest in low dimensional materials due to its exceptionally high crystal and electronic quality [2–4]. As a member of carbon family, graphene materials present extraordinary thermal [5], mechanical [6] and electrical [7] properties. In this way, the graphene is considered as an excellent candidate for building block in nano-electro mechanical systems (NEMS) [8,9], nanoelectronics [10,11] and nanocomposites [12,13]. Thus, the fundamental understanding of graphene properties is the critical issue in its future applications. One of the key properties of graphene is its mechanical properties that have considerable effects on its application for the NEMS devices and reinforcement part in polymer nanocomposites. Among the limited available experimental studies, by using atomic force microscope nanoindentation experiment, Frank et al. [14] reported the elastic modulus of stacks of graphene sheets as 0.5 TPa. Using the same experimental technique, Lee et al. [6] acquired the elastic modulus of 1 TPa and breaking strength of 0.13 TPa for monolayer and free of defect graphene sheets. Since the experimental characterization techniques at nanoscale are considerably complex, today the numerical simulations are gaining conspicuous interest as an alternative for experimental studies [15]. Molecular dynamics (MD) is one of the most popular techniques that have been used for atomic-scale
⇑ Corresponding author. Tel.: +33 (0) 3 68 85 29 52; fax: +33 (0) 3 68 85 29 36. E-mail address: [email protected] (S. Ahzi). 0927-0256/$ - see front matter Ó 2011 Elsevier B.V. All rights reserved. doi:10.1016/j.commatsci.2011.08.018
characterization of material properties. MD simulation computes the motion of individual molecules in the models of solid, liquid, or gas. The key idea here is the motion [16], then physical and chemical properties can be evaluated based on the statistical mechanics theories. There have been numerous MD studies on the mechanical properties of bulk and nanoribbons graphene. Bu et al. [17] and Ni et al. [18] studied the mechanical behavior of graphene nanoribbons and sheets respectively by MD simulations. Tsai and Tu [19] calculated the mechanical properties of single graphene layer and bulk graphite by MD simulation. One of the alternatives of MD simulations in evaluation of mechanical properties of carbon nanotubes and sheets is structural mechanics modeling. In the structural mechanics, the atoms are considered as balls and their interaction force fields are represented by beams or spring like finite elements. Recently, Georgantzinos et al. [20] studied the size-dependent non-linear mechanical properties of graphene nanoribbons. Furthermore, Sakhaee-Pour [21] and Scarpa et al. [22] used structural mechanics modeling techniques for the evaluation of elastic properties of the single layer graphene sheets. Georgantzinos et al. [23] also studied the single and multilayer graphene structures mechanical properties by structural mechanics modeling. The main privilege of structural mechanics is its certainty and computational efficiency in the evaluation of elastic constants. On the other hand MD seems to be a more versatile technique for assessment of large deformations, failure mechanism and temperature effects. In the current study, uniaxial MD tension tests were carried out in order to investigate the directional and thickness effect on deformation and mechanical response of FLG sheets using Tersoff potential at room temperatures.
B. Mortazavi et al. / Computational Materials Science 53 (2012) 298–302
2. Molecular dynamics modeling LAMMPS (Large-Scale Atomic/Molecular Massively Parallel Simulator) [24], developed at Sandia National Laboratories, is used for molecular dynamics simulations. The carbon–carbon bonded interaction is modeled by Tersoff potential [25,26]. In case of nonbonded van der Waals interaction between individual carbon atoms, there is a variety of used potentials in the literature. Commonly, the nonbonded interactions are expressed as Lennard– Jones (LJ) or Morse potentials. In the current study, the LJ potential was used which is expressed by:
12 6 r r /ðrÞ ¼ 4e r r
ð1Þ
where r refers to the distance between carbon atoms at different atomic layers. e is the depth of the potential well and r is the equilibrium distance at which the inter-particle potential is zero. A parameterization frequently employed for LJ interactions in the modeling of multi-layer carbon structures is e = 2.4 m eV and r = 3.4 Å which is based on the compressibility and graphene layer spacing [27] which also satisfy the density of the graphite (2.2 g/ cm3). As a common assumption for the LJ potentials, interactions were ignored at distances higher than 2.5r. As a result of the van der Waals interactions between carbons atoms in different layers, in the plane properties of FLG, are expected to be dependent of the number of layers. The MD simulation model for uniaxial tensile test of FLG with the chirality angle of 10° is illustrated in Fig. 1. The chirality angle of 0° and 30° are commonly called armchair and zigzag directions, respectively. A FLG is constructed by putting graphene sheets up together by a spacing distance of 3.4 Å. The periodic boundary conditions were applied in x and y directions, thus these simulations are representative of an infinite FLG sheet and not the FLG nanoribbons. The properties of bulk graphite were also evaluated by considering eight layers of graphene sheets while the periodic boundary condition was also applied in the graphene sheets normal direction (z direction). Due to high computational cost of MD simulations, the tensile tests are usually performed at high strain rates, in the order of 109 s1. The limitation in approaching lower strain rates is mainly due to the small interatomic distances and high velocities of atoms. In this way, to accurately count for the deformations of FLG at considerably high applying strain rates, small time increment of 0.5 fs (=0.5 1015 s) was used. Prior to applying the loading condition, the FLG structure was left to anisotropically relax to zero pressure in x and y directions at a temperature of 300 K for 30,000 time increments using constant pressure–temperature (i.e. NPT
299
ensemble) by means of Nosé–Hoover barostat and thermostat method. The loading condition was applied by extending the periodic simulation box size in the x direction under constant strain rates. In order to satisfy the uniform uniaxial stress conditions, the periodic boundary condition was applied in the loading direction and the simulation box in the y direction was altered to reach zero stress in this direction using NPT Nosé–Hoover method. Virial stresses in the x direction (uniaxial stress) were calculated at each strain level. The engineering strain at each time step was calculated by multiplying of total time of loading simulation by the applied strain rate. Using the Hooke’s law, the slope of the initial linear part of uniaxial stress–strain curve is equivalent to the elastic modulus of FLG along the loading direction. 3. Results and discussions Since the applied strain rates are considerably high, the simulations were performed at different strain rate to verify the level of strain rate dependency of the obtained properties. In Fig. 2, the acquired engineering stress–strain diagrams of three layers zigzag graphene and five layers armchair graphene at different strain rates (SR) are plotted. At low strain levels (up to 0.05), a linear relation between stress and strain values is obtained. This linear response is commonly assumed as the linear elasticity part, in which the corresponding slope is equivalent to the elastic modulus. Based on the Fig. 2 observations, the elastic modulus of the armchair and zigzag FLG are reasonably independent of the strain rate values. At higher strain level, the stress–strain response presents a nonlinear behavior which continues up to the maximum tensile strength, the point in which rupture occurs in the first atomic plane. In case of zigzag FLG, the nonlinear part is also found to be acceptably independent of strain rate where the rupture occurs at approximately the same stress and strain levels. On the other hand, while the nonlinear part of armchair FLG stress–strain curves also present acceptable level of strain rate independency up to the tensile strength. However, the strains at rupture decreases as the strain rate decreases. Furthermore, by decreasing of the loading strain rate in both armchair and zigzag loading directions, the rupture is accompanied with sharper decline in the stress–strain diagram which consequently means more brittle behavior. Based on the obtained results, the stain rate of 1 109 s1 is selected for the further MD simulations in this study, because of its computational efficiency as well as predicting a stress–strain trend, closer to the lower strain rate of 1 108 s1. In Fig. 3a, the engineering stress–strain response of zigzag FLG containing different number of layers of 2, 3, 6 and bulk graphite
Fig. 1. MD model of uniaxial tension test of FLG with chirality angle of 10° in respect to armchair direction. The periodic boundary conditions were applied in the x and y direction. The loading condition was applies by the extending the periodic simulation box in the x direction.
300
B. Mortazavi et al. / Computational Materials Science 53 (2012) 298–302
Fig. 2. Uniaxial engineering stress–strain relation of three layers zigzag graphene (dashed line) and five layers armchair graphene (continuous line) at different loading strain rates (SR).
Fig. 4. (a) Uniaxial engineering stress–strain relation of three layers FLG as function of different loading chirality angle (ChA) directions and (b) the effect of loading direction on the elastic modulus for FLGs with different number of layers and bulk graphite.
Fig. 3. (a) Uniaxial engineering stress–strain curve of zigzag FLG for different number of layers and bulk graphite and (b) the dependence of FLG number of layers on the elastic modulus for armchair and zigzag loading directions.
(containing eight layers in which periodic boundaries are also applied along the thickness) is presented. As an overall trend, by increasing of the number of atomic layers the elastic modulus and ultimate tensile strength decrease (see Fig. 3a and b). Ni
et al. [18] reported the elastic modulus of around 1.1 TPa for a single layer graphene sheet. In our case, due to the high level of uncertainty for the results of a single layer graphene, we are not reporting any value for this case. However, if we take into account the value of 1.1 TPa for the single layer case, we can conclude that the trend of the elastic modulus variation as function of the number of atomic layers of grapheme (see Fig. 3b) is similar to the trend reported by Liew et al. [28] for multi-walled carbon nanotubes (CNTs). These authors [28] found that the double wall CNTs present higher elastic modulus and ultimate strength than the single layer CNTs, and by increasing of CNTs number of walls, both of these mechanical properties decline. As it can be observed from Fig. 3a in case of zigzag FLGs, the ultimate tensile strength, where the ruptures occur in the first atomic layer, is obtained at approximately the same strain values for all zigzag FLG cases. Since the subsequent ruptures at adjacent layers may not occur at the same time, series of stress fall and raise can be observed after the ultimate tensile strength point. In Fig. 3b the effect of FLG number of layers on the elastic modulus for zigzag and armchair loading directions are shown. In the current study the elastic modulus of bulk graphite was found to be around 0.98 TPa, which is in fine agreement with the experimental ranges of 1.02 ± 0.03 TPa [29]. Based on the obtained results, it can be concluded that as the number of FLG layers reach the number of six, the obtained elastic properties become acceptably independent of number of layers, where the elastic modulus is close to the ranges of experimental bulk properties. This matter can be interesting in the modeling of polymer
B. Mortazavi et al. / Computational Materials Science 53 (2012) 298–302
301
Fig. 5. The brittle fracture process at different engineering strain level of (a–c) two layers zigzag graphene and (d–f) three layers graphene with chirality angle of around 10°. The initial void shapes in (a) and (c) insets and the fracture line in (b) and (f) insets.
nanocomposites, reinforced by multi-layer graphene. Thus, as an accurate assumption in the modeling of the graphene with higher number of layers than six, we can use the graphite bulk properties. It is worth mentioning that, in-plane thermal conductivity of FLG was experimentally proven to become independent of the number of layers as the number of layers is more than four [30]. Fig. 4a illustrates the effects of loading direction chirality angle (ChA) on the engineering stress–strain response of three layers FLG. As an overall trend, both elastic modulus and ultimate strength of FLG increase by increasing of the chirality angle from 0° (Armchair) to 30° (Zigzag). As presented in Fig. 4b, the chirality dependency of mechanical properties of FLG is more obvious for two and three layers. For the case of FLG with higher number of layers than four, the elastic modulus dependency on chirality reduces considerably. It was also found that in case of bulk graphite the elastic modulus is independent of loading direction. However it was found that the ultimate tensile strength of armchair is still less than zigzag direction by around 15%.
In Fig. 5, the deformation behavior of two layers zigzag graphene (a–c) and three layers graphene with chirality angle of around 10° (d–f) is shown. It was found, that up to the high strain levels (about 30% strain) the specimen extends uniformly along the loading direction. At a point, very close to the ultimate strength, several voids are formed along one of the graphene sheets. These formed voids have the shapes illustrated in Fig. 5a,c insets. The formations of these voids take place during the short time, which results in fast subsequent failure in the first atomic sheet while the remaining sheets are undamaged. The beginning of voids formation and subsequent failure occur at considerably close strain levels. This observation is suggesting that the deformation mechanism of FLG present a high level of void sensitivity; which accordingly means brittle failure behavior. Furthermore, as the rupture takes place in each graphene sheet, the sheets tend to retrieve their original structure and contract rapidly (see Fig. 5d,e – red sheet). However as shown in Fig. 5c inset, the initially formed voids are stabilized in the structure. This damage/failure mechanism extends from one
302
B. Mortazavi et al. / Computational Materials Science 53 (2012) 298–302
sheet to the other and the final rupture of the FLG occurs at the failure of the last sheet. Moreover, as presented in Fig. 5b,f insets, the cut edge of the graphene sheets are often along the zigzag direction which reveals that the debonding between the carbon atoms happened along the armchair direction. This observation is in agreement with the fact that the ultimate tensile strength of FLG and mono layer graphene sheets [18] along the armchair directions falls lower than that in the zigzag direction.
4. Conclusion Molecular dynamics simulations of uniaxial tension were performed in order to investigate the stress–strain and failure behavior of few-layer graphene. The simulations were carried out at different strain rates ranging between 1010 and 108 s1 and it was observed that both linear and nonlinear part of the stress strain response of zigzag FLG are convincingly strain rate independent and only the strain at failures of armchair FLG decreases by decreasing the loading rate. The modeling technique was validated by comparison of the obtained elastic modulus of the simulation and the bulk graphite experiments. The simulated results show that the properties are sensitive to the number of layers and loading direction as well. It was observed that by increasing the number of graphene layers the elastic modulus and ultimate tensile strength dropped. However, the obtained elastic modulus for FLG containing six atomic planes is close to the bulk graphite values. By performing the simulations at various loading directions, it was found that the directional sensitivity of the elastic modulus is considerable only for the FLG with two and three layers. We also found that the ultimate tensile strength increases when the loading direction chirality angle goes from 0° to 30°. Finally by plotting the deformation process of the FLG, it was concluded that the FLG presented a brittle behavior and the debonding between carbon atoms occurred, frequently, along the armchair direction. Acknowledgment The authors acknowledge the financial support of the FNR of Luxembourg via the AFR Grants (PHD-09-016). The first author
greatly acknowledges Michael Essa from University of Strasbourg for his valuable helps. References [1] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, Y. Zhang, S.V. Dubonos, I.V. Grigorieva, A.A. Firsov, Science 306 (2004) 666. [2] K.S. Novoselov, A.K. Geim, S.V. Morozov, D. Jiang, M.I. Katsnelson, I.V. Grigorieva, S.V. Dubonos, A.A. Firsov, Nature 438 (2005) 197. [3] Y.B. Zhang, Y.W. Tan, H.L. Stormer, P. Kim, Nature 438 (2005) 201. [4] A.K. Geim, K.S. Novoselov, Nat. Mater. 6 (2007) 183. [5] A.A. Balandin, S. Ghosh, W.Z. Bao, I. Calizo, D. Teweldebrhan, F. Miao, et al., Nano. Lett. 8 (2008) 902. [6] C. Lee, X. Wei, J.W. Kysar, J. Hone, Science 321 (2008) 385. [7] J.R. Williams, L. DiCarlo, C.M. Marcus, Science 317 (2007) 638. [8] J.S. Bunch, A.M. VanderZande, S.S. Verbridge, I.W. Frank, D.M. Tanenbaum, J.M. Parpia, H.G. Craighead, P.L. McEuen, Science 315 (2007) 490. [9] B. Standley, W.Z. Bao, H. Zhang, J. Bruck, C.N. Lau, M. Bockrath, Nano. Lett. 8 (2008) 3345. [10] D.A. Areshkin, C.T. White, Nano. Lett. 7 (2007) 3253. [11] R.M. Westervelt, Science 320 (2008) 324. [12] 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. [13] T. Ramanathan, A.A. Abdala, S. Stankovich, D.A. Dikin, M. Herrera-Alonso, R.D. Piner, D.H. Adamson, H.C. Schniepp, X. Chen, R.S. Ruoff, S.T. Nguyen, I.A. Aksay, R.K. Prud’homme, L.C. Brinson, Nat. Nanotechnol. 3 (2008) 327. [14] I.W. Frank, D.M. Tanenbaum, A.M. Van der Zande, P.L. McEuen, J. Vac. Sci. Technol. B 25 (2007) 2558. [15] B Mortazavi, A.A. Khatibi, C. Politis, J. Comput. Theor. Nanosci. 6 (2009) 644. [16] J.M. Haile, Molecular Dynamic Simulation, John Wiley & Sons, New York, 1992. [17] H. Bu, Y. Chen, M. Zou, H. Yia, K. Bi, Z. Ni, Phys. Lett. A 373 (2009) 3359. [18] Z. Ni, H. Bu, M. Zou, H. Yi, K. Bi, Y. Chen, Physica B 405 (2010) 1301. [19] J.L. Tsai, J.F. Tu, Mater. Des. 31 (2010) 194. [20] S.K. Georgantzinos, G.I. Giannopoulos, D.E. Katsareas, P.A. Kakavas, N.K. Anifantis, Comput. Mater. Sci. 50 (2011) 2057. [21] A. Sakhaee-Pour, Solid State Commun. 149 (2009) 91. [22] F. Scarpa, S. Adhikari, A. Srikantha Phani, Nanotechnology 20 (2009) 065709. [23] S.K. Georgantzinos, G.I. Giannopoulos, N.K. Anifantis, Mater. Des. 31 (2010) 4646. [24] S. Plimpton, J. Comput. Phys. 117 (1995) 1. [25] J. Tersoff, Phys. Rev. B 37 (1988) 6991. [26] J. Tersoff, Phys. Rev. Lett. 61 (1988) 2879. [27] L.A. Girifalco, M. Hodak, R.S. Lee, Phys. Rev. B 62 (2000) 13104. [28] K.M. Liew, X.Q. He, C.H. Wong, Acta Mater. 52 (2004) 2521. [29] O.L. Blakslee, D.G. Proctor, E.J. Seldin, G.B. Spence, T. Weng, J. Appl. Phys. 41 (1970) 3373. [30] S. Ghosh, W. Bao, D.L. Nika, S. Subrina, E.P. Pokatilov, C.N. Lau, A.A. Balandin, Nat. Mater. 9 (2010) 555.