- Email: [email protected]
A molecular dynamics study of nanofracture in monolayer boron nitride
Materials Science & Engineering A 641 (2015) 225–230
Contents lists available at ScienceDirect
Materials Science & Engineering A journal homepage: www.elsevier.com/locate/msea
A molecular dynamics study of nanofracture in monolayer boron nitride Alireza Tabarraei n, Xiaonan Wang Department of Mechanical Engineering & Engineering Science, University of North Carolina at Charlotte, Charlotte, NC 28223, United States
art ic l e i nf o
a b s t r a c t
Article history: Received 14 March 2015 Received in revised form 29 May 2015 Accepted 2 June 2015 Available online 17 June 2015
In this paper, we use molecular dynamics (MD) modeling to study the fracture properties of monolayer hexagonal boron nitride (h-BN) under mixed mode I and II loading. We investigate the impact of crack edge chirality, crack tip configuration and loading phase angle on the crack propagation path and critical stress intensity factors. The MD results predict that under all the loading phase angles cracks prefer to propagate along a zigzag direction and the critical stress intensity factors of zigzag cracks are higher than those of armchair cracks. Under mixed mode loading, the h-BN sheets can undergo out-of-plane deformations due to the buckling induced by compressive stresses. The out-of-plane deformations can be significant when mode II loading is dominant. An excessive amount of out-of-plane deformation can induce buckling cracks. Depending on the loading phase angle and crack configurations, buckling cracks can nucleate before or after the propagation of the original cracks. & 2015 Elsevier B.V. All rights reserved.
Keywords: Two dimensional materials Monolayer hexagonal boron nitride Fracture mechanics Molecular dynamics
1. Introduction The current interest to graphene and the need to introduce an electronic bandgap to the gapless pristine graphene has brought intense interests to other graphene-like two-dimensional materials such as hexagonal boron nitride [1–3]. Hexagonal boron nitride (hBN) has a honeycomb atomic structure in which boron (B) and nitrogen (N) atoms occupy alternating sites. Despite similar morphology with graphene, the binary atomic structure of hexagonal boron nitride leads to mixed ionic-covalent atomic bonding. Such atomic bonding which is different from the covalent sp2 bonding of graphene gives h-BN distinct physical and mechanical properties than graphene. Opposed to graphene which displays a zero bandgap and is semimetallic [4], h-BN shows a band gap and is an insulator with a wide band gap of 5.6 eV [5]. Such differences in the physical properties of h-BN and graphene, has motivated the synthesize of hybrid graphene–boron nitride sheets with properties which are complementary to both graphene and boron nitride [3,6]. Besides being complementary to graphene, the remarkable physical and mechanical properties h-BN such as low dielectric constant [7], high temperature stability [8], high thermal conductivity [9] and high strength [10,11] make it appealing for a wide spectrum of applications in its own right. The wide spectrum of applications of h-BN expose it to different thermal, electrical and n
Corresponding author. Fax: þ 1 704 687 8345. E-mail address: [email protected] (A. Tabarraei).
http://dx.doi.org/10.1016/j.msea.2015.06.012 0921-5093/& 2015 Elsevier B.V. All rights reserved.
vibrational loadings. For a reliable usage of h-BN in such sensitive devices it is essential to ensure that h-BN can sustain such loadings. This necessitates a fundamental understanding regarding the mechanical and fracture properties of h-BN. Opposed to graphene whose mechanical and fracture has been studied in the past [12–14], the studies on the fracture properties of graphene-like two-dimensional materials is scarce. Considering that graphene-like two-dimensional materials such as h-BN have a more complex atomic structure than graphene, their fracture properties can be quite different from graphene. Such differences necessitate a separate study of the fracture of h-BN. Ideally, the fracture and mechanical properties of h-BN should be characterized experimentally, e.g., by uniaxial test. However, designing and conducting test at nanoscales is very complicated; to date no uniaxial test on two-dimensional materials has been reported. Computational studies such as molecular dynamics, on the other hand, can provide valuable insights regarding the behavior of twodimensional materials. In this paper, we use molecular dynamics simulations using a Tersoff potential [15] to study the fracture properties of single layer boron nitride sheets under mixed mode I and II loading.
2. Tersoff potential We use molecular dynamics simulations to study the toughness and crack propagation path of monolayer boron nitride under
226
A. Tabarraei, X. Wang / Materials Science & Engineering A 641 (2015) 225–230
mixed mode I and II loading. The molecular dynamics simulations are conducted using LAMMPS package [16,17]. The interatomic interaction between boron and nitrogen atoms in the BN sheets are prescribed using a Tersoff potential [15]. The Tersoff potential can be written as
Vij = fC (rij )⎡⎣fR (rij ) + bij f A (rij )⎤⎦,
(1)
where fC is a cut-off function defined as
⎧1, r≤R−D ⎪ ⎪ ⎛π r − R⎞ ⎟, R − D < r < R + D fC (r ) = ⎨ 0.5 − 0.5 sin⎜ ⎝2 D ⎠ ⎪ ⎪ ⎩ 0, r≥R+D
ux =
(2)
The function fR is a two body term and represents a repulsive pair potential whereas fA is a three body term and represents the attractive pair potential due to the atomic bondings. The functions fR (rij ) and fA (rij ) are given by
fR (rij ) = A exp(λ1rij )
(3a)
f A (rij ) = − B exp(λ2rij ).
(3b)
Function bij represents a measure of the bond order which depends on the local coordination of the neighbors of atoms i and the angle θijk between atoms i, j and k. Function bij is defined by
(
bij = 1 + β nζijn
ζij =
−1/2n
)
,
(4a)
⎛ fC (rij )g (θijk ) exp⎜⎜λ33 rij − rik ⎝ k ≠ i, j
∑
g (θ ) = 1 +
(
3
⎞
) ⎟⎟, ⎠
c2 c2 − 2 . d2 d + (cos θ − h)2
is shown in Fig. 1. The circular domain is chosen large enough that its boundary falls in the K-dominant zone. The equilibrium configuration of the cracked domain is obtained by first applying the crack tip asymptotic field to all the atoms in the domain and then while the boundary atoms are kept fixed the position of the interior atoms is relaxed. The crack tip asymptotic displacement fields for a linear isotropic material under combined mode I and II loading are given by [23]
(4b)
(4c)
In the above equations, the summation is taken over all neighbors j and k of atoms i which are located within a cutoff distance of Rþ D. Several sets of Tersoff potential parameters for interaction between boron and nitrogen atoms have been developed in the past [9,18–20]. Sekkal et al. [18] modified the Tersoff set of parameters of carbon [15] to describe the B and N interactions in cubic boron nitride (c-BN) systems. Verma et al. [19] modified this set of parameters to describe the interactions of B and N in hexagonal BN (h-BN) systems. More recently, adjusted set of Tersoff potential parameters for boron nitride are developed by fitting the obtained interatomic forces, bond lengths, cohesive energy and phonon dispersion curves to the experimental and ab initio modeling data [9,20,21]. Furthermore, Tersoff potential parameters have been extended to describe the interaction among carbon, boron and nitrogen [20,22]. Such parameters are used to study the thermal transport and mechanical properties of more complex materials such as hybrid graphene-boron nitride [20] or graphitic carbon nitride [22]. In this paper, we use the potential parameters and their corresponding values for BN provided by Çağın et al. [9,20]. This potential parameter set has been successfully used in accurate reproduction of structural, mechanical, vibrational and thermal transport properties of hexagonal boron nitride nanostructures [9] and hybrid graphene-boron nitride nanoribbons [20].
3. Molecular dynamics simulations Our MD model is a circular domain cut around the crack tip as
uy =
⎡ θ⎛ 2 θ⎞ ⎢⎣KI cos ⎝⎜κ − 1 + 2 sin ⎠⎟ 2 2 ⎤ θ ⎛⎜ 2 θ⎞ ⎟ +KII sin κ + 1 + 2 cos ⎥ 2⎝ 2 ⎠⎦
(5a)
⎞ θ⎛ 1+v r ⎡ 2 θ⎟ ⎢KI sin ⎜⎝κ + 1 − 2 cos E 2π ⎣ 2 2⎠
(5b)
1+v r E 2π
−KII cos
θ ⎛⎜ θ ⎞⎤ κ − 1 − 2 sin2 ⎟⎥ 2⎝ 2 ⎠⎦
(5c)
where r and θ are the polar coordinates shown in Fig. 1, ux and uy are the displacement components in the x and y directions, E¼ 925 GPa [24] is Young's modulus of pristine monolayer h-BN, ν ¼ 0.23 [24] is its Poisson's ratio and κ is the Kolosov constant which is equal to (3 − v ) /(1 + v ) for plane stress condition. KI and KII are the mode I and II stress intensity factors; respectively. The effective stress intensity factor is Keff =
KI2 + KII2 , and the loading
−1
phase angle defined as ϕ = tan (KI/KII) describes the ratio of the mode I and II loading. Based on this definition a loading phase angle of zero degree corresponds to a pure mode I loading and a loading phase angle of 90° corresponds to a pure mode II loading. The initial zigzag and armchair cracks are generated by eliminating respectively four and three rows of atoms as are shown in Fig. 2. We investigate the effect of crack edge chirality on the fracture properties by considering both armchair (AC) and zigzag (ZZ) cracks. As are shown in Fig. 2a, the edges of zigzag cracks are not alike; in one of the edges, boron atoms occupy the outermost layer while in the other edge nitrogen atoms are located at the outermost layer. When the loading phase angle is not zero (i.e. loading is not pure mode I), the deformation of the top and bottom edges of the crack is not symmetric, e.g. under pure mode II loading one edge is under compression and the other edge is under tension. To take into account the crack edge unsymmetry the loading phase angle is varied from 90° to 90°. This allows us to consider the effect of tension or compression in either of the edges on the critical stress intensity factor, Kcr .
Fig. 1. Molecular dynamic domain. The boundary atoms are shown in cyan. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
A. Tabarraei, X. Wang / Materials Science & Engineering A 641 (2015) 225–230
227
Fig. 2. Schematic of the five crack types: (a) ZZ crack, (b) a blunt and c) sharp AC crack with a B atom at their tip, (d) a blunt and (e) sharp AC crack with a N atom at their tip.
The edges of cracks with armchair chirality are symmetric, however armchair cracks can have different tip configurations. The impact of tip configuration on the critical stress intensity factors and crack propagation paths are studied by considering the four different crack tips shown in Fig. 2b–e. Depending on the location of the crack, a boron (B) or a nitrogen atom (N) can occupy the crack tip. Furthermore, the crack tip can have a blunt shape (Fig. 2b and d) or can have a sharp configuration as are shown in Fig. 2c and e. We refer to blunt and sharp cracks with a boron atom at the tip as AC-BB and AC-BS and to the blunt and sharp cracks with a nitrogen atom at the tip as AC-NB and AC-NS; respectively. Our goal is to model the crack propagation under quasi-static loading. For this purpose, the loading is increased in increments of ΔKeff = 0.01 MPa m . After each loading increment the position of boundary atoms are kept fixed while the positions of interior atoms are relaxed using the conjugate gradient method. The temperature of interior domain is then increased to 300 K using the velocity-rescaled Berendsen thermostat. After relaxing for 3 ps using a microcanonical (NVE) ensemble a Nosé–Hoover thermostat is used to maintain the temperature at 300 K for 60 ps in a canonical (NVT) ensemble. The velocity-Verlet algorithm with a time step of 1 fs is used for the purpose of time integration of atoms trajectory. The critical stress intensity factor reported in this paper corresponds to the stress intensity factor at which the first bond at the crack tip or crack edge breaks. A bond is considered broken if the average distance of the atoms of the bond over the last 10 ps of simulation exceeds the potential cutoff distance, i.e, if the average bond length becomes larger than R þD.
4. Results and discussion The crack propagation path of zigzag and armchair cracks are shown in Fig. 3. These figures show that loading phase angle, crack edge chirality and crack tip configuration affect the crack growth path. An important observation is that regardless of the loading phase angle or crack geometry, cracks tend to advance along a zigzag direction. A similar observation has been made previously for crack propagation in graphene sheets [12,13]. The propagation paths shown in Fig. 3 indicate that as the cracks grow, they can kink. In that case the crack path might include segments of armchair direction. To clearly separate the armchair and zigzag segments of the propagation paths, in Fig. 3 the zigzag atoms are shown in green and armchair atoms are shown in red. These figures show that the length of armchair segments are much smaller than the zigzag segments of crack path. The only condition in which crack propagates along a self-similar path is a zigzag crack under pure mode I loading (loading phase equal to zero). Armchair cracks kink even under mode I loading to propagate in a zigzag direction. Moreover, Fig. 3 shows that in addition to crack edge chirality and loading phase angle, crack tip configuration also impacts the crack propagation path. Out-of-plane deformation occurs in the boron nitride sheet when the loading phase angle is not zero as is shown in Fig. 4. The out of plane deformation is due to the buckling of BN sheet under compressive loading. Since BN sheet is very thin, compressive
stress can generate local buckling in the sheet. The out-of-plane deformation releases some energy through bending and can postpone the crack propagation. Such a phenomenon has been observed previously in cracked thin plates under tensile or shear loading [25–28] and is predicted to occur in some other two-dimensional materials such as MoS2 [29]. The maximum out-ofplane deformation as a function of the loading phase angle is shown in Fig. 5. The plots of this figure indicate that the out-ofplane deformation is negligible when the loading phase angle is close to zero and increases as the loading phase angle approaches 90°, i.e. out-of-plane deformation is considerable when mode II loading is dominant. Excessive out-of-plane deformation can induce buckling cracks (BC) in the sheet. Since increase in the loading phase angle leads to increase in the out-of-plane deformation, it is expected that buckling cracks initiate for the loading phase angles that are close to 90°. A buckling crack is shown in Fig. 4. As shown, buckling cracks start at the crack edge near the crack tip, propagate in a direction orthogonal to the initial crack and end on the same crack edge. Hence, due to the buckling crack, a segment of the boron nitride sheet is separated from the rest of the sheet. eff The effective critical stress intensity factor (Kcr ) associated with the propagation of the original crack (CP) and nucleation of buckling cracks (BC) are shown in Figs. 6–8 for zigzag and armchair cracks. The plots of Fig. 6 show that the mode I critical stress intensity factor for a zigzag crack is about 5.56 MPa m . This can be used to find the surface energy of a zigzag surface by using the Griffith relation for brittle materials which relates the critical stress intensity factor to the surface energy density γ by
Kcr =
2γ E .
(6)
Using Eq. (6), the surface energy density in a zigzag direction is 18.06 J/m2 which is about 50% higher than the zigzag surface energy of 12 J/m2 for graphene [30]. The critical stress intensity factors are larger when the loading phase angle is less than zero, i.e. the edge with nitrogen atoms at its outermost layer are in compression. The pure mode II (ϕ ¼90° or 90°) critical stress intensity factor for a zigzag crack is about 6.60 MPa m . This indicates that crack propagation under shearing mode demands more energy compared to the energy required for crack propagation under opening load. This is partly due to the out-of-plane deformation generated when mode II is dominant. Some of the energy of the applied loading is used by the elastic energy of buckling bending, hence more external energy should be applied to the system before a bond at the crack tip breaks. The dashed curve in Fig. 6 corresponds to the stress intensity factor at which a buckling crack nucleates. These curves show that when loading phase angle is less than 45°, no buckling crack is generated. This is in consistency with the out-of-plane deformation curves shown in Fig. 5, i.e. the amount of out-of-plane deformation is negligible when the loading phase angle is less than 45°. Moreover, the solid and dashed curves of Fig. 6 intersect at a phase angle of 52° indicating that the buckling cracks nucleate before the propagation of the original crack only when the loading phase angle is larger than 52°. The effective critical stress intensity factors for armchair cracks
228
A. Tabarraei, X. Wang / Materials Science & Engineering A 641 (2015) 225–230
Fig. 3. Crack propagation paths in a BN sheet. (a) ZZ crack, (b) blunt AC crack with a B atom at the tip, (c) sharp AC crack with a B atom at the tip, (d) blunt AC crack with a N atom at the tip and (e) sharp AC crack with a N atom at the tip.
Fig. 4. A buckling crack under pure mode II loading in boron nitride sheet.
with boron and nitrogen atoms at their tips are shown in Figs. 7 and 8; respectively. Comparing the plots of these figures with the plots of Fig. 6 shows that in general the critical stress
intensity factors of zigzag cracks are higher than those of armchair cracks. The plots of these figures indicate that crack tip configuration can affect the critical stress intensity factor. Mode I stress intensity factor is almost the same for sharp and blunt configurations. The mode I critical stress intensity factor of an armchair crack is 5.35 MPa m and 5.15 MPa m for cracks with a boron and nitrogen atom at their tips, respectively. The mode II critical stress intensity factor is higher than mode I and depending on the tip configuration varies from 6.2 MPa m to 7.2 MPa m . For both cases of cracks with a boron or nitrogen atom at the tip, when the loading phase angle is less than 30°, the critical stress intensity eff of cracks with a factor of cracks with a blunt tip is larger than Kcr sharp tip and when the loading phase angle is larger than 30° cracks with a sharp tip configuration have a higher critical stress intensity factor. The crack tip configuration also impacts the formation of buckling cracks. Buckling cracks are not generated for armchair cracks with a sharp crack tip. Buckling cracks initiate only if the crack tip is blunt and the loading phase angle is larger than 75°.
A. Tabarraei, X. Wang / Materials Science & Engineering A 641 (2015) 225–230
Fig. 5. Maximum out-of-plane deformation of cracked boron nitride sheet as a function of loading phase angle.
229
Fig. 8. Critical stress intensity as a function of loading phase angle for an armchair crack with a nitrogen atom at its tip.
5. Conclusion
Fig. 6. Critical stress intensity as a function of loading phase angle for a zigzag crack.
In summary, we have studied the fracture properties of monolayer boron nitride sheet using molecular dynamics method. The results show that the crack edge chirality, crack tip configuration and loading phase angle affect the critical stress intensity factor and crack propagation path of the boron nitride sheet. In all the cases studied in this paper, the critical stress intensity factor corresponding to pure mode II loading is higher than that of pure mode I loading. Furthermore, all the cracks propagate along a zigzag direction, although some cracks kink during growth in which case a short path of crack growth can be along the armchair direction. Besides the propagation of the main crack, excessive out-of-plane deformation of BN sheet under the mixed mode loading can lead to the formation of buckling cracks. Buckling cracks are more important in sheets with an initial zigzag crack since in such sheets buckling cracks can form when the loading phase angle is larger than 45°. In BN sheets containing an armchair crack, buckling cracks can only initiate when the crack tip configuration is blunt and the loading phase angle is larger than 75°.
Acknowledgment The authors would like to thank the University of North Carolina at Charlotte for providing the fund for partial support of this research.
References
Fig. 7. Critical stress intensity as a function of loading phase angle for an armchair crack with a boron atom at its tip.
Comparing the dashed curves of Figs. 7 and 8 with those of Fig. 6 indicate that the formation of buckling cracks is more a severe issue for sheets with initial zigzag cracks in them.
[1] G. Giovannett, P.A. Khomyakov, G. Brocks, P.J. Kelly, J. van den Brink, Phys. Rev. B 76 (2007) 073103. [2] C. Dean, A. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. Shepard, et al., Nat. Nanotechnol. 5 (2010) 722. [3] L. Ci, L. Song, C. Jin, D. Jariwala, D. Wu, Y. Li, A. Srivastava, Z. Wang, K. Storr, L. Balicas, et al., Nat. Mater. 9 (2010) 430. [4] B. Partoens, F. Peeters, Phys. Rev. B 74 (2006) 075404. [5] Y. Kubota, K. Watanabe, O. Tsuda, T. Taniguchi, Science 317 (2007) 932. [6] Z. Liu, L. Ma, G. Shi, W. Zhou, Y. Gong, S. Le, X. Yang, J. Zhang, J. Yu, K. Hackenberg, A. Babakhani, J. Idrobo, R. Vajtai, J. Lou, P. Ajayan, Nat. Nanotechnol. 8 (2013) 119. [7] L. Britnell, R. Gorbachev, R. Jalil, B. Belle, F. Schedin, M. Katsnelson, L. Eaves, S. Morozov, A. Mayorov, N. Peres, A. Castro Neto, J. Leist, A. Geim, L. Ponomarenko, K. Novoselov, Nano Lett. 12 (2012) 1707. [8] J. Eichler, C. Lesniak, J. Eur. Ceram. Soc. 28 (2008) 1105. [9] C. Sevik, A. Kinaci, J.B. Haskins, T. Çağın, Phys. Rev. B 84 (2011) 085409. [10] Q. Peng, A. Zamiri, W. Ji, S. De, Acta Mech. 223 (2012) 2591. [11] T. Han, Y. Luo, C. Wang, J. Phys. D: Appl. Phys. 47 (2014) 025303.
230
A. Tabarraei, X. Wang / Materials Science & Engineering A 641 (2015) 225–230
[12] M. Xu, A. Tabarraei, J. Paci, T. Belytschko, Int. J. Fract. 173 (2012) 163. [13] B. Zhang, L. Mei, H. Xiao, Appl. Phys. Lett. 101 (2012) 121915. [14] M.A. Rafiee, J. Rafiee, I. Srivastava, Z. Wang, H. Song, Z.-Z. Yu, N. Koratkar, Small 6 (2010) 179. [15] J. Tersoff, Phys. Rev. B 37 (1988) 6991. [16] S. Plimpton, J. Comput. Phys. 117 (1995) 1. [17] 〈http://lammps.sandia.gov〉. [18] W. Sekkal, B. Bouhafs, H. Aourag, M. Certier, J. Phys.: Condens. Matter 10 (1998) 4975. [19] V. Verma, V. Jindal, K. Dharamvir, Nanotechnology 18 (2007) 435711. [20] A. Kınacı, J.B. Haskins, C. Sevik, T. Çağın, Phys. Rev. B 86 (2012) 115410. [21] M.-L. Liao, Y.-C. Wang, S.-P. Ju, T.-W. Lien, L.-F. Huang, J. Appl. Phys. 110 (2011) 054310.
[22] B. Mortazav, G. Cuniberti, T. Rabczuk, Comput. Mater. Sci. 99 (2015) 285. [23] T.L. Anderson, Fracture Mechanics: Fundamentals and Applications, CRC Press, Boca Raton, FL, 2005. [24] Q. Peng, W. Ji, S. De, Comput. Mater. Sci. 56 (2012) 11. [25] K. Markstrom, B. Storakers, Int. J. Solid Mech. 16 (1980) 217. [26] G. Sih, Y. Lee, Theoret. Appl. Fract. Mech. 6 (1986) 129. [27] R. Brighenti, Thin-walled Struct. 43 (2005) 209. [28] R. Brighenti, Int. J. Mech. Mater. Des. 6 (2010) 73. [29] X. Wang, A. Tabarraei, D. Spearot, Nanotechnology 26 (2015) 175703. [30] H. Yin, H. Qi, F. Fan, T. Zhu, B. Wang, Y. Wei, Nano Lett. 15 (2015) 1918–1924.
Contents lists available at ScienceDirect
Materials Science & Engineering A journal homepage: www.elsevier.com/locate/msea
A molecular dynamics study of nanofracture in monolayer boron nitride Alireza Tabarraei n, Xiaonan Wang Department of Mechanical Engineering & Engineering Science, University of North Carolina at Charlotte, Charlotte, NC 28223, United States
art ic l e i nf o
a b s t r a c t
Article history: Received 14 March 2015 Received in revised form 29 May 2015 Accepted 2 June 2015 Available online 17 June 2015
In this paper, we use molecular dynamics (MD) modeling to study the fracture properties of monolayer hexagonal boron nitride (h-BN) under mixed mode I and II loading. We investigate the impact of crack edge chirality, crack tip configuration and loading phase angle on the crack propagation path and critical stress intensity factors. The MD results predict that under all the loading phase angles cracks prefer to propagate along a zigzag direction and the critical stress intensity factors of zigzag cracks are higher than those of armchair cracks. Under mixed mode loading, the h-BN sheets can undergo out-of-plane deformations due to the buckling induced by compressive stresses. The out-of-plane deformations can be significant when mode II loading is dominant. An excessive amount of out-of-plane deformation can induce buckling cracks. Depending on the loading phase angle and crack configurations, buckling cracks can nucleate before or after the propagation of the original cracks. & 2015 Elsevier B.V. All rights reserved.
Keywords: Two dimensional materials Monolayer hexagonal boron nitride Fracture mechanics Molecular dynamics
1. Introduction The current interest to graphene and the need to introduce an electronic bandgap to the gapless pristine graphene has brought intense interests to other graphene-like two-dimensional materials such as hexagonal boron nitride [1–3]. Hexagonal boron nitride (hBN) has a honeycomb atomic structure in which boron (B) and nitrogen (N) atoms occupy alternating sites. Despite similar morphology with graphene, the binary atomic structure of hexagonal boron nitride leads to mixed ionic-covalent atomic bonding. Such atomic bonding which is different from the covalent sp2 bonding of graphene gives h-BN distinct physical and mechanical properties than graphene. Opposed to graphene which displays a zero bandgap and is semimetallic [4], h-BN shows a band gap and is an insulator with a wide band gap of 5.6 eV [5]. Such differences in the physical properties of h-BN and graphene, has motivated the synthesize of hybrid graphene–boron nitride sheets with properties which are complementary to both graphene and boron nitride [3,6]. Besides being complementary to graphene, the remarkable physical and mechanical properties h-BN such as low dielectric constant [7], high temperature stability [8], high thermal conductivity [9] and high strength [10,11] make it appealing for a wide spectrum of applications in its own right. The wide spectrum of applications of h-BN expose it to different thermal, electrical and n
Corresponding author. Fax: þ 1 704 687 8345. E-mail address: [email protected] (A. Tabarraei).
http://dx.doi.org/10.1016/j.msea.2015.06.012 0921-5093/& 2015 Elsevier B.V. All rights reserved.
vibrational loadings. For a reliable usage of h-BN in such sensitive devices it is essential to ensure that h-BN can sustain such loadings. This necessitates a fundamental understanding regarding the mechanical and fracture properties of h-BN. Opposed to graphene whose mechanical and fracture has been studied in the past [12–14], the studies on the fracture properties of graphene-like two-dimensional materials is scarce. Considering that graphene-like two-dimensional materials such as h-BN have a more complex atomic structure than graphene, their fracture properties can be quite different from graphene. Such differences necessitate a separate study of the fracture of h-BN. Ideally, the fracture and mechanical properties of h-BN should be characterized experimentally, e.g., by uniaxial test. However, designing and conducting test at nanoscales is very complicated; to date no uniaxial test on two-dimensional materials has been reported. Computational studies such as molecular dynamics, on the other hand, can provide valuable insights regarding the behavior of twodimensional materials. In this paper, we use molecular dynamics simulations using a Tersoff potential [15] to study the fracture properties of single layer boron nitride sheets under mixed mode I and II loading.
2. Tersoff potential We use molecular dynamics simulations to study the toughness and crack propagation path of monolayer boron nitride under
226
A. Tabarraei, X. Wang / Materials Science & Engineering A 641 (2015) 225–230
mixed mode I and II loading. The molecular dynamics simulations are conducted using LAMMPS package [16,17]. The interatomic interaction between boron and nitrogen atoms in the BN sheets are prescribed using a Tersoff potential [15]. The Tersoff potential can be written as
Vij = fC (rij )⎡⎣fR (rij ) + bij f A (rij )⎤⎦,
(1)
where fC is a cut-off function defined as
⎧1, r≤R−D ⎪ ⎪ ⎛π r − R⎞ ⎟, R − D < r < R + D fC (r ) = ⎨ 0.5 − 0.5 sin⎜ ⎝2 D ⎠ ⎪ ⎪ ⎩ 0, r≥R+D
ux =
(2)
The function fR is a two body term and represents a repulsive pair potential whereas fA is a three body term and represents the attractive pair potential due to the atomic bondings. The functions fR (rij ) and fA (rij ) are given by
fR (rij ) = A exp(λ1rij )
(3a)
f A (rij ) = − B exp(λ2rij ).
(3b)
Function bij represents a measure of the bond order which depends on the local coordination of the neighbors of atoms i and the angle θijk between atoms i, j and k. Function bij is defined by
(
bij = 1 + β nζijn
ζij =
−1/2n
)
,
(4a)
⎛ fC (rij )g (θijk ) exp⎜⎜λ33 rij − rik ⎝ k ≠ i, j
∑
g (θ ) = 1 +
(
3
⎞
) ⎟⎟, ⎠
c2 c2 − 2 . d2 d + (cos θ − h)2
is shown in Fig. 1. The circular domain is chosen large enough that its boundary falls in the K-dominant zone. The equilibrium configuration of the cracked domain is obtained by first applying the crack tip asymptotic field to all the atoms in the domain and then while the boundary atoms are kept fixed the position of the interior atoms is relaxed. The crack tip asymptotic displacement fields for a linear isotropic material under combined mode I and II loading are given by [23]
(4b)
(4c)
In the above equations, the summation is taken over all neighbors j and k of atoms i which are located within a cutoff distance of Rþ D. Several sets of Tersoff potential parameters for interaction between boron and nitrogen atoms have been developed in the past [9,18–20]. Sekkal et al. [18] modified the Tersoff set of parameters of carbon [15] to describe the B and N interactions in cubic boron nitride (c-BN) systems. Verma et al. [19] modified this set of parameters to describe the interactions of B and N in hexagonal BN (h-BN) systems. More recently, adjusted set of Tersoff potential parameters for boron nitride are developed by fitting the obtained interatomic forces, bond lengths, cohesive energy and phonon dispersion curves to the experimental and ab initio modeling data [9,20,21]. Furthermore, Tersoff potential parameters have been extended to describe the interaction among carbon, boron and nitrogen [20,22]. Such parameters are used to study the thermal transport and mechanical properties of more complex materials such as hybrid graphene-boron nitride [20] or graphitic carbon nitride [22]. In this paper, we use the potential parameters and their corresponding values for BN provided by Çağın et al. [9,20]. This potential parameter set has been successfully used in accurate reproduction of structural, mechanical, vibrational and thermal transport properties of hexagonal boron nitride nanostructures [9] and hybrid graphene-boron nitride nanoribbons [20].
3. Molecular dynamics simulations Our MD model is a circular domain cut around the crack tip as
uy =
⎡ θ⎛ 2 θ⎞ ⎢⎣KI cos ⎝⎜κ − 1 + 2 sin ⎠⎟ 2 2 ⎤ θ ⎛⎜ 2 θ⎞ ⎟ +KII sin κ + 1 + 2 cos ⎥ 2⎝ 2 ⎠⎦
(5a)
⎞ θ⎛ 1+v r ⎡ 2 θ⎟ ⎢KI sin ⎜⎝κ + 1 − 2 cos E 2π ⎣ 2 2⎠
(5b)
1+v r E 2π
−KII cos
θ ⎛⎜ θ ⎞⎤ κ − 1 − 2 sin2 ⎟⎥ 2⎝ 2 ⎠⎦
(5c)
where r and θ are the polar coordinates shown in Fig. 1, ux and uy are the displacement components in the x and y directions, E¼ 925 GPa [24] is Young's modulus of pristine monolayer h-BN, ν ¼ 0.23 [24] is its Poisson's ratio and κ is the Kolosov constant which is equal to (3 − v ) /(1 + v ) for plane stress condition. KI and KII are the mode I and II stress intensity factors; respectively. The effective stress intensity factor is Keff =
KI2 + KII2 , and the loading
−1
phase angle defined as ϕ = tan (KI/KII) describes the ratio of the mode I and II loading. Based on this definition a loading phase angle of zero degree corresponds to a pure mode I loading and a loading phase angle of 90° corresponds to a pure mode II loading. The initial zigzag and armchair cracks are generated by eliminating respectively four and three rows of atoms as are shown in Fig. 2. We investigate the effect of crack edge chirality on the fracture properties by considering both armchair (AC) and zigzag (ZZ) cracks. As are shown in Fig. 2a, the edges of zigzag cracks are not alike; in one of the edges, boron atoms occupy the outermost layer while in the other edge nitrogen atoms are located at the outermost layer. When the loading phase angle is not zero (i.e. loading is not pure mode I), the deformation of the top and bottom edges of the crack is not symmetric, e.g. under pure mode II loading one edge is under compression and the other edge is under tension. To take into account the crack edge unsymmetry the loading phase angle is varied from 90° to 90°. This allows us to consider the effect of tension or compression in either of the edges on the critical stress intensity factor, Kcr .
Fig. 1. Molecular dynamic domain. The boundary atoms are shown in cyan. (For interpretation of the references to color in this figure caption, the reader is referred to the web version of this paper.)
A. Tabarraei, X. Wang / Materials Science & Engineering A 641 (2015) 225–230
227
Fig. 2. Schematic of the five crack types: (a) ZZ crack, (b) a blunt and c) sharp AC crack with a B atom at their tip, (d) a blunt and (e) sharp AC crack with a N atom at their tip.
The edges of cracks with armchair chirality are symmetric, however armchair cracks can have different tip configurations. The impact of tip configuration on the critical stress intensity factors and crack propagation paths are studied by considering the four different crack tips shown in Fig. 2b–e. Depending on the location of the crack, a boron (B) or a nitrogen atom (N) can occupy the crack tip. Furthermore, the crack tip can have a blunt shape (Fig. 2b and d) or can have a sharp configuration as are shown in Fig. 2c and e. We refer to blunt and sharp cracks with a boron atom at the tip as AC-BB and AC-BS and to the blunt and sharp cracks with a nitrogen atom at the tip as AC-NB and AC-NS; respectively. Our goal is to model the crack propagation under quasi-static loading. For this purpose, the loading is increased in increments of ΔKeff = 0.01 MPa m . After each loading increment the position of boundary atoms are kept fixed while the positions of interior atoms are relaxed using the conjugate gradient method. The temperature of interior domain is then increased to 300 K using the velocity-rescaled Berendsen thermostat. After relaxing for 3 ps using a microcanonical (NVE) ensemble a Nosé–Hoover thermostat is used to maintain the temperature at 300 K for 60 ps in a canonical (NVT) ensemble. The velocity-Verlet algorithm with a time step of 1 fs is used for the purpose of time integration of atoms trajectory. The critical stress intensity factor reported in this paper corresponds to the stress intensity factor at which the first bond at the crack tip or crack edge breaks. A bond is considered broken if the average distance of the atoms of the bond over the last 10 ps of simulation exceeds the potential cutoff distance, i.e, if the average bond length becomes larger than R þD.
4. Results and discussion The crack propagation path of zigzag and armchair cracks are shown in Fig. 3. These figures show that loading phase angle, crack edge chirality and crack tip configuration affect the crack growth path. An important observation is that regardless of the loading phase angle or crack geometry, cracks tend to advance along a zigzag direction. A similar observation has been made previously for crack propagation in graphene sheets [12,13]. The propagation paths shown in Fig. 3 indicate that as the cracks grow, they can kink. In that case the crack path might include segments of armchair direction. To clearly separate the armchair and zigzag segments of the propagation paths, in Fig. 3 the zigzag atoms are shown in green and armchair atoms are shown in red. These figures show that the length of armchair segments are much smaller than the zigzag segments of crack path. The only condition in which crack propagates along a self-similar path is a zigzag crack under pure mode I loading (loading phase equal to zero). Armchair cracks kink even under mode I loading to propagate in a zigzag direction. Moreover, Fig. 3 shows that in addition to crack edge chirality and loading phase angle, crack tip configuration also impacts the crack propagation path. Out-of-plane deformation occurs in the boron nitride sheet when the loading phase angle is not zero as is shown in Fig. 4. The out of plane deformation is due to the buckling of BN sheet under compressive loading. Since BN sheet is very thin, compressive
stress can generate local buckling in the sheet. The out-of-plane deformation releases some energy through bending and can postpone the crack propagation. Such a phenomenon has been observed previously in cracked thin plates under tensile or shear loading [25–28] and is predicted to occur in some other two-dimensional materials such as MoS2 [29]. The maximum out-ofplane deformation as a function of the loading phase angle is shown in Fig. 5. The plots of this figure indicate that the out-ofplane deformation is negligible when the loading phase angle is close to zero and increases as the loading phase angle approaches 90°, i.e. out-of-plane deformation is considerable when mode II loading is dominant. Excessive out-of-plane deformation can induce buckling cracks (BC) in the sheet. Since increase in the loading phase angle leads to increase in the out-of-plane deformation, it is expected that buckling cracks initiate for the loading phase angles that are close to 90°. A buckling crack is shown in Fig. 4. As shown, buckling cracks start at the crack edge near the crack tip, propagate in a direction orthogonal to the initial crack and end on the same crack edge. Hence, due to the buckling crack, a segment of the boron nitride sheet is separated from the rest of the sheet. eff The effective critical stress intensity factor (Kcr ) associated with the propagation of the original crack (CP) and nucleation of buckling cracks (BC) are shown in Figs. 6–8 for zigzag and armchair cracks. The plots of Fig. 6 show that the mode I critical stress intensity factor for a zigzag crack is about 5.56 MPa m . This can be used to find the surface energy of a zigzag surface by using the Griffith relation for brittle materials which relates the critical stress intensity factor to the surface energy density γ by
Kcr =
2γ E .
(6)
Using Eq. (6), the surface energy density in a zigzag direction is 18.06 J/m2 which is about 50% higher than the zigzag surface energy of 12 J/m2 for graphene [30]. The critical stress intensity factors are larger when the loading phase angle is less than zero, i.e. the edge with nitrogen atoms at its outermost layer are in compression. The pure mode II (ϕ ¼90° or 90°) critical stress intensity factor for a zigzag crack is about 6.60 MPa m . This indicates that crack propagation under shearing mode demands more energy compared to the energy required for crack propagation under opening load. This is partly due to the out-of-plane deformation generated when mode II is dominant. Some of the energy of the applied loading is used by the elastic energy of buckling bending, hence more external energy should be applied to the system before a bond at the crack tip breaks. The dashed curve in Fig. 6 corresponds to the stress intensity factor at which a buckling crack nucleates. These curves show that when loading phase angle is less than 45°, no buckling crack is generated. This is in consistency with the out-of-plane deformation curves shown in Fig. 5, i.e. the amount of out-of-plane deformation is negligible when the loading phase angle is less than 45°. Moreover, the solid and dashed curves of Fig. 6 intersect at a phase angle of 52° indicating that the buckling cracks nucleate before the propagation of the original crack only when the loading phase angle is larger than 52°. The effective critical stress intensity factors for armchair cracks
228
A. Tabarraei, X. Wang / Materials Science & Engineering A 641 (2015) 225–230
Fig. 3. Crack propagation paths in a BN sheet. (a) ZZ crack, (b) blunt AC crack with a B atom at the tip, (c) sharp AC crack with a B atom at the tip, (d) blunt AC crack with a N atom at the tip and (e) sharp AC crack with a N atom at the tip.
Fig. 4. A buckling crack under pure mode II loading in boron nitride sheet.
with boron and nitrogen atoms at their tips are shown in Figs. 7 and 8; respectively. Comparing the plots of these figures with the plots of Fig. 6 shows that in general the critical stress
intensity factors of zigzag cracks are higher than those of armchair cracks. The plots of these figures indicate that crack tip configuration can affect the critical stress intensity factor. Mode I stress intensity factor is almost the same for sharp and blunt configurations. The mode I critical stress intensity factor of an armchair crack is 5.35 MPa m and 5.15 MPa m for cracks with a boron and nitrogen atom at their tips, respectively. The mode II critical stress intensity factor is higher than mode I and depending on the tip configuration varies from 6.2 MPa m to 7.2 MPa m . For both cases of cracks with a boron or nitrogen atom at the tip, when the loading phase angle is less than 30°, the critical stress intensity eff of cracks with a factor of cracks with a blunt tip is larger than Kcr sharp tip and when the loading phase angle is larger than 30° cracks with a sharp tip configuration have a higher critical stress intensity factor. The crack tip configuration also impacts the formation of buckling cracks. Buckling cracks are not generated for armchair cracks with a sharp crack tip. Buckling cracks initiate only if the crack tip is blunt and the loading phase angle is larger than 75°.
A. Tabarraei, X. Wang / Materials Science & Engineering A 641 (2015) 225–230
Fig. 5. Maximum out-of-plane deformation of cracked boron nitride sheet as a function of loading phase angle.
229
Fig. 8. Critical stress intensity as a function of loading phase angle for an armchair crack with a nitrogen atom at its tip.
5. Conclusion
Fig. 6. Critical stress intensity as a function of loading phase angle for a zigzag crack.
In summary, we have studied the fracture properties of monolayer boron nitride sheet using molecular dynamics method. The results show that the crack edge chirality, crack tip configuration and loading phase angle affect the critical stress intensity factor and crack propagation path of the boron nitride sheet. In all the cases studied in this paper, the critical stress intensity factor corresponding to pure mode II loading is higher than that of pure mode I loading. Furthermore, all the cracks propagate along a zigzag direction, although some cracks kink during growth in which case a short path of crack growth can be along the armchair direction. Besides the propagation of the main crack, excessive out-of-plane deformation of BN sheet under the mixed mode loading can lead to the formation of buckling cracks. Buckling cracks are more important in sheets with an initial zigzag crack since in such sheets buckling cracks can form when the loading phase angle is larger than 45°. In BN sheets containing an armchair crack, buckling cracks can only initiate when the crack tip configuration is blunt and the loading phase angle is larger than 75°.
Acknowledgment The authors would like to thank the University of North Carolina at Charlotte for providing the fund for partial support of this research.
References
Fig. 7. Critical stress intensity as a function of loading phase angle for an armchair crack with a boron atom at its tip.
Comparing the dashed curves of Figs. 7 and 8 with those of Fig. 6 indicate that the formation of buckling cracks is more a severe issue for sheets with initial zigzag cracks in them.
[1] G. Giovannett, P.A. Khomyakov, G. Brocks, P.J. Kelly, J. van den Brink, Phys. Rev. B 76 (2007) 073103. [2] C. Dean, A. Young, I. Meric, C. Lee, L. Wang, S. Sorgenfrei, K. Watanabe, T. Taniguchi, P. Kim, K. Shepard, et al., Nat. Nanotechnol. 5 (2010) 722. [3] L. Ci, L. Song, C. Jin, D. Jariwala, D. Wu, Y. Li, A. Srivastava, Z. Wang, K. Storr, L. Balicas, et al., Nat. Mater. 9 (2010) 430. [4] B. Partoens, F. Peeters, Phys. Rev. B 74 (2006) 075404. [5] Y. Kubota, K. Watanabe, O. Tsuda, T. Taniguchi, Science 317 (2007) 932. [6] Z. Liu, L. Ma, G. Shi, W. Zhou, Y. Gong, S. Le, X. Yang, J. Zhang, J. Yu, K. Hackenberg, A. Babakhani, J. Idrobo, R. Vajtai, J. Lou, P. Ajayan, Nat. Nanotechnol. 8 (2013) 119. [7] L. Britnell, R. Gorbachev, R. Jalil, B. Belle, F. Schedin, M. Katsnelson, L. Eaves, S. Morozov, A. Mayorov, N. Peres, A. Castro Neto, J. Leist, A. Geim, L. Ponomarenko, K. Novoselov, Nano Lett. 12 (2012) 1707. [8] J. Eichler, C. Lesniak, J. Eur. Ceram. Soc. 28 (2008) 1105. [9] C. Sevik, A. Kinaci, J.B. Haskins, T. Çağın, Phys. Rev. B 84 (2011) 085409. [10] Q. Peng, A. Zamiri, W. Ji, S. De, Acta Mech. 223 (2012) 2591. [11] T. Han, Y. Luo, C. Wang, J. Phys. D: Appl. Phys. 47 (2014) 025303.
230
A. Tabarraei, X. Wang / Materials Science & Engineering A 641 (2015) 225–230
[12] M. Xu, A. Tabarraei, J. Paci, T. Belytschko, Int. J. Fract. 173 (2012) 163. [13] B. Zhang, L. Mei, H. Xiao, Appl. Phys. Lett. 101 (2012) 121915. [14] M.A. Rafiee, J. Rafiee, I. Srivastava, Z. Wang, H. Song, Z.-Z. Yu, N. Koratkar, Small 6 (2010) 179. [15] J. Tersoff, Phys. Rev. B 37 (1988) 6991. [16] S. Plimpton, J. Comput. Phys. 117 (1995) 1. [17] 〈http://lammps.sandia.gov〉. [18] W. Sekkal, B. Bouhafs, H. Aourag, M. Certier, J. Phys.: Condens. Matter 10 (1998) 4975. [19] V. Verma, V. Jindal, K. Dharamvir, Nanotechnology 18 (2007) 435711. [20] A. Kınacı, J.B. Haskins, C. Sevik, T. Çağın, Phys. Rev. B 86 (2012) 115410. [21] M.-L. Liao, Y.-C. Wang, S.-P. Ju, T.-W. Lien, L.-F. Huang, J. Appl. Phys. 110 (2011) 054310.
[22] B. Mortazav, G. Cuniberti, T. Rabczuk, Comput. Mater. Sci. 99 (2015) 285. [23] T.L. Anderson, Fracture Mechanics: Fundamentals and Applications, CRC Press, Boca Raton, FL, 2005. [24] Q. Peng, W. Ji, S. De, Comput. Mater. Sci. 56 (2012) 11. [25] K. Markstrom, B. Storakers, Int. J. Solid Mech. 16 (1980) 217. [26] G. Sih, Y. Lee, Theoret. Appl. Fract. Mech. 6 (1986) 129. [27] R. Brighenti, Thin-walled Struct. 43 (2005) 209. [28] R. Brighenti, Int. J. Mech. Mater. Des. 6 (2010) 73. [29] X. Wang, A. Tabarraei, D. Spearot, Nanotechnology 26 (2015) 175703. [30] H. Yin, H. Qi, F. Fan, T. Zhu, B. Wang, Y. Wei, Nano Lett. 15 (2015) 1918–1924.