Fracture behavior dependent on crack-tip shapes in nanoscale crack-defect monolayer boron nitride sheets

ABSTRACT Nanoscale defects, including cracks, circular holes, and the triangular-shaped defects, often occur in the growth of boron nitride nanosheets (BNNS). In this study, the fracture behavior of chiral BNNS with different crack-tip shapes and the interactions of nanoscale crack-defects are studied using molecular dynamics (MD) simulations and finite element (FE) analysis. Both MD and FE results indicate that the fracture strength of BNNS with two crack tips (t = 2) is significantly higher than that with one crack tip (t = 1), in which the difference in zigzag (ZZ) direction is more obvious than that in armchair (AC) direction, mainly due to the fact that the change of bond angles near the crack tips is more substantial in the ZZ direction than those in the AC direction. Our results show that the fracture strength of BNNS strongly depends on crack-tip shapes, chiral angles, the defect-to-crack tip spacing and deflection angles. Checking against the current MD simulations and FE analysis shows the present results are reasonable. This study should be of great importance for enhancing the fracture performance of BNNS by modulating their crack-tip shapes and the interactions of nanoscale crack-defects. Graphical abstract


Introduction
In recent years, two-dimensional (2D) nanomaterials, such as graphene [1,2], hexagonal boron nitride (h-BN) [3][4][5], have attracted extensive attention due to their excellent mechanical properties and application prospects in the field of electronic components and nanomaterials [6][7][8][9]. Although the fracture strength of these 2D nanomaterials has been extensively investigated [10][11][12][13], the genetic defects obtained from the preparation of boron nitride nanosheets (BNNS) during the chemical vapor deposition (CVD) [14,15], such as cracks, holes and triangular-shaped defects, greatly reduce the mechanical strength and properties of the corresponding materials, where the mechanism behind the association between defects and fracture strength is not clear enough.
In comparison to other types of defects, cracks are usually more harmful to the mechanical properties of the materials. Therefore, it is always necessary to fully understand the mechanical properties of BNNS with triangular-shaped defects and cracks before putting them into practical applications. Liu et al. [16] reported the epitaxial growth of a single-crystal monolayer h-BN on the low symmetric Cu (110) surface, which is obtained by annealing of industrial copper foil. Xiong et al. [17] researched the impact of different defects on monolayer h-BN fracture behavior through molecular dynamics (MD) simulations. Zhang et al. [18,19] studied the branching and oscillation phenomena of BN, with the crack propagation velocity variation. Wei et al [20]. studied the fracture strength variation of pristine graphene with the density functional theory (DFT) and energy minimization calculations, and it was found that the fracture strength was monotonous with the variation of chiral angle. Meng et al. [21,22] researched the shielding effect of monolayer graphene under loading when 5-7 defects were found near the crack, which provided a basis for crack propagation protection. Yao et al. [23] investigated the interactions of nanoscale crack-hole in chiral graphene nanosheets (GNS) using finite element (FE) analysis and MD simulations.
Despite the fact that the mechanical properties of BNNS have been reported by experiments and MD simulations [24][25][26], the effect of different crack-tip shapes on the chiral BNNS fracture behavior has been rarely investigated. In this paper, both FE analysis and MD simulations are used to study the fracture behavior of BNNS at multiple scales [27]. To compare with the MD results, a FE model based on Tersoff potentials [28] is established to obtain the fracture behaviors of chiral BNNS.
In this study, the effect of crack-tip shapes on the mechanical behavior of BNNS is mainly studied, which has not been reported in the previous research. The triangularshaped defects are commonly existed in BNNS, so it is of great importance for studying the interactions between crack and triangular-shaped defects, which also has not been mentioned in the previous articles. It is found that the mechanical properties of BNNS can be properly improved by modulating the crack-tip shapes and the relative position of the crack and defects. Checking against MD and FE results shows that the present results are reasonable. This study should be of great importance for improving the fracture performance of BNNS by modulating their crack-tip shapes and interactions of nanoscale crack-defects [29,30]. Figure 1 shows the structures of BNNS with one crack and four symmetric defects. The lengths in BNNS are around 170 Å along the armchair (AC) direction and 150 Å along the zigzag (ZZ) direction. The crack length is 20 Å and the circular defects diameter are 7.6 Å. The triangular-shaped defects are formed by deleting four atoms.

Finite element method
The FE analysis in this paper is carried out by the ABAQUS/Explicit code and the material properties of the beam model are defined by VUMAT [31]. The Timoshenko beams model is suitable to describe the variation of bond length and B-N bond fracture behavior before fracture. The geometry of BNNS used in FE analysis is the same as that after energy minimization in MD simulations. Since slight spatial wrinkles in BNNS can be found after energy minimization in MD simulations, the three-dimensional (3D) beam elements of B31 are used to study the fracture behavior of BNNS.
In this study, both FE analysis and MD simulations are used to research the fracture behaviors of BNNS with the Tersoff potential, which can effectively characterize the large deformation and nonlinear behaviors of B-N bonds. In order to obtain the stress-strain relationship of B-N bonds conveniently, the potential energy U ij [29,32] between the ith B atom and the jth N atom based on Tersoff is expressed as The functions V A and V R , which represent the attractions and repulsions, respectively, are pair-additive interactions. r ij is the distance between the ith B atom and the jth N atom. b ij is a function of the bond angle and f c is the cutoff function. θ ijk is the angle between atoms i, j, and k. All other parameters are presented in Table 1. According to Eq. (1), the interatomic force F can be written as The stress and strain of B-N bond equivalent to the beam model can be expressed as where d 0 and S are, respectively, the initial cross-sectional diameter and area of B-N bonds, and r 0 is the initial bond length. The above equations and parameters are used to analyze the force and geometry balance of BN beam model [22]. A series of specific parameters for B31 are listed in Table  2. Due to the quite complex calculation of Eq. (9), it is difficult to describe the stress-strain relation of B-N bonds in VUMAT. The nonlinear constitutive relation of the B-N bond before broken can be derived as  [ 1] h = -0.89000 In Figure 2a, the accuracy of Eq. (12) is verified by the comparison between Eq. (9) and Eq. (12) for a B-N bond. The relative error is less than 5% before fracture, indicating that the Eq. (9) and Eq. (12) results are considerably consistent. In this study, the constitutive relation in FE analysis is selected as Eq. (12). The elongation of the bonds and dihedral angle are also considered.

Molecular dynamics simulations
In MD simulations, the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [33] is used to research the BNNS mechanical behavior and the Tersoff potential is employed to describe the B-N bond rupture [34]. The cutoff distance of the B-N bond is 1.778 Å [35], as shown in Figure 2a. Free boundary conditions are applied along the x and y directions, while the periodic boundary condition is applied along the z direction and the thickness of vacuum layer is around 20 Å along the z direction. BNNS is stretched under a constant displacement after energy minimization. Firstly, the structure of BNNS is optimized by the energy minimization based on conjugate-gradient algorithm, where energy and force tolerances are both 10 − [8]. After initial structural optimization, BNNS is stretched by the deformationcontrol method and the time step is chosen as 1 fs [36]. Afterward, the structure is optimized after each displacement increment (0.01 Å) and the optimized structure is taken as the initial geometry for the next calculations. Figure 2b shows the stress-strain curves of BNNS under uniaxial tension with FE methods and MD simulations based on Tersoff potential. The stress-strain curves between two methods are distinct in the AC and ZZ directions. The major reason is that the B-N bond is considered as Timoshenko beam model and the fracture condition is determined by the ultimate stress in FE method. When the beam stress reaches the ultimate stress, the beam length is less than 1.778 Å (the cutoff distance of MD simulations). Therefore, the ultimate stresses in MD simulations are always higher than those of FE analysis in Figure 2b.

Fracture toughness of BNNS with crack
Based on the classical fracture mechanics, the fracture toughness K I under mode I loading can be obtained from the stress distribution at the crack tip. The fracture toughness K I can be defined as [37,38] Where Y is a dimensionless parameter (Y ≈ 1), σ m is the fracture stress and a is half of the crack length. As shown in Figure 3a, there are four kinds of cracks with different shapes in MD models. Figure 3b-c shows the distribution of K I with the variation of crack length in armchair and zigzag BNNS with FE analysis and MD simulations. The fracture toughness K I of FE analysis is consistent with that of MD simulations. As the result of the different cracktip shapes, the fracture toughness along the AC and the ZZ direction shows an obvious distinction. When the crack has two tips (t = 2), the fracture toughness in AC and ZZ directions is larger than that of one tip (t = 1), and the difference in ZZ direction is more evident. As shown in Figure 3d, the K I of GNS with four different cracks is obtained by MD simulations. Different from BNNS, under the condition of a certain crack length, the K I of graphene in AC direction is greater than that in ZZ direction, which are in agreement with previous research [39,40]. It can be seen that GNS has the same properties as BNNS in Figure 3b-c. Figure 3e shows the relative difference of K I between two kinds of tip-shape cracks in different chiral BNNS and GNS. The variation of the fracture toughness in the ZZ direction is more pronounced than that in the AC direction at the two crack tips whether it is GNS or BNNS. However, it is also worth noting that the difference of K I in BNNS is greater than that in GNS.
The mechanism of the phenomenon in Figure 3 is further revealed based on the variation of bond length and bond angle. Figure 4 shows the variation of bond length and bond angle in BNNS with crack length of 2a = 15 Å. Due to the geometric symmetry of the atomic structure, the bond length and bond angle which play a major role in the model are only selected. The variation of the crack with two kinds of tip shapes (t = 1 and t = 2) in AC and ZZ directions are different. In order to intuitively analyze the variation of bond length and bond angle, a relative change rate of bond length k L is introduced as: where the ΣΔL 1 , ΣΔL 2 are the sum of ΔL that the crack with one tip and two tips at the moment of fracture, respectively. The ΣΔθ 1 , ΣΔθ 2 are the sum of Δθ that the crack with one tip and two tips at the moment of fracture, respectively. According to the Eq. (14)(15), the k L in AC direction and ZZ direction are both about 35% in the case of two shapes crack tips. However, the k θ with the two kinds of cracks in AC direction is about 11.8%, while the k θ in ZZ direction is about 41.2%. These results indicate that the variation of bond angle in BNNS dominates the fracture difference in Figure 3 Figure 5 indicates fracture stress distribution of BNNS with different crack widths in the AC and ZZ directions with two kinds of crack tips. The fracture stress increases slightly with increasing crack width, indicating that the crack width has little effect on its fracture behavior. However, the fracture stresses of two tips are greater than those of one tip, and the difference along the ZZ direction are more obvious. Therefore, the tip shape still has a significant effect on the fracture behavior for different width cracks.

Interactions of nanoscale crack-defects in BNNS
In this study, the stress field at the singular crack tip without defects is described with the stress intensity factor (SIF) K I (BN) [41]. Similarly, K I (BN-h) can describe the stress field that the crack tip with defects, which can be written as where ΔK I is the toughness variation due to the presence of defects. Figure 6 shows the normalized stress intensity factor K I (BN-h) /K I (BN) for the BNNS with different crack-tip shapes in AC and ZZ directions using FE analysis and MD simulations when circular defects with deflection angle θ = 0° are contained. The results of the two methods are in good agreement. It can be seen from the Figure 6 that the normalized critical K I (BN-h) /K I (BN) decrease with increasing distance between the circular defects and the crack. In the two cases of AC direction cracks, the K I (BN-h) /K I (BN) diversity between two kinds of crack-tip is not obvious. As shown in Figure 6d, the crack in ZZ direction contains two tips and the normalized critical K I (BN-h) /K I (BN) is larger than the other three ones. The initial fracture occurs near the circular defect based on the stress nephogram (B atoms and N atoms are colored in the light of their von Mises stress). This situation in Figure 6d can be regarded as the coexistence of one crack with two tips and two small cracks. The fracture toughness of small cracks with one tip is smaller than that of large crack with two tips, especially in ZZ direction. As a result, the fracture occurs from the circular defect in the first place during the tensile process. Therefore, the fracture strength of the structure with circular defects in Figure 6d is greatly different from that without circular defects, and thus leads to a larger K I (BN-h) /K I (BN) .
Compared with the circular defects, the triangular-shaped defects with losing four atoms are more common in the preparation of BNNS due to the influence of substrate and preparation process. The interaction between the triangular-shaped defects and the crack with different tip shapes is investigated with FE analysis and MD simulations (see Figure 7). The defect-to-crack tip spacing r is 11 Å and the deflection angles θ (θ is the angle between the defect-to-crack tip direction and the horizontal direction) are 10°, 30°, 45°, 60° and 90°, respectively. The normalized critical K I (BN-h) /K I (BN) gradually decreases under uniaxial tension with increasing θ. In Figure 7a, the crack in the AC direction interacts with the triangular-shaped defects and the K I (BN-h) /K I (BN) of one tip crack is  greater than that of two tips crack. It indicates that the triangular-shaped defects have a more significant effect on the crack with one tip. Whereas the effect of triangular-shaped defects on the two kinds of crack-tip has no significant difference in ZZ direction (see in Figure 7b). Before θ = 90°, the values of K I (BN-h) /K I (BN) are always bigger than 1, which demonstrates that the existence of triangular-shaped defects will accelerate the fracture. However, for θ = 90°, the values of K I (BN-h) /K I (BN) are smaller than 1, that is to say, triangular-shaped defects can restrain the crack growth in the tensile process.

Interactions of nanoscale crack-defect in chiral BNNS
The chiral BNNS with nanoscale crack-defects interactions is investigated with FE analysis and MD simulations in Figure 8, where chiral angle α are 0°, 7°, 12°, 17°, 22° and 30°. In Figure 8a, the definition of BN chiral angle can be intuitively understood. Figure 8b shows the variation of fracture stress with chiral angle for complete BNNS and BNNS with horizontal crack. The fracture stress decreases with increasing chiral angle. When the chiral angle α = 30°, there is an upward trend in complete monolayer chiral BNNS. However, this phenomenon does not exist in chiral BNNS with horizontal cracks. Based on the Figure 4, the reasons for this case is mainly attributed to the variation of the bond angle, which plays a vital role under uniaxial tension when the chiral angle α = 30°. Four symmetrically distributed triangular-shaped defects are located near the horizontal crack of the chiral BNNS, and the deflection angles θ are 10°, 45°, and 90°, respectively (see Figure 8c). The values of K I (BN-h) /K I (BN) present an arched trend with the variation of chiral angle under uniaxial tension. Therefore, if the defects distributed near the crack tip in chiral BNNS, they have a marked impact on the fracture toughness when the chiral angle is around 15°. For the same chiral angle, the smaller the deflection angle θ, the larger the normalized K I (BN-h) /K I (BN) . As shown in Figure 8d, the surface energies of both BNNS and GNS generally decrease with increasing chiral angle. It is predicted that the tensile strength breaking armchair planes is bigger than that tearing apart zigzag planes, which is consistent with the stress-strain curves in Figure 8b.

Conclusion
In this study, fracture properties of the chiral BNNS with different crack-tip shapes and the nanoscale crack-defect interactions have been studied with MD simulations and FE analysis. In the FE method, the B-N bond has been equivalent to a nonlinear Timoshenko beam. The fracture toughness of BNNS at two tips is always higher than that with one tip, and the two kinds of crack-tip shapes in the ZZ direction have more obvious difference on fracture toughness. The main reason is that the bond angle have more variations in two tips cracks. The influence of bond angle variation is also a nonnegligible factor in BN fracture process. Furthermore, in view of the interaction between defects and the main crack, the fracture strength of BNNS depends on deflection angle θ, the defect-to-crack tip spacing r and the chiral angle α. Our work is of great help for improving the fracture strength of chiral BNNS with cracks and understanding the interaction between defects and crack, and then providing the basis on the fracture engineering of BNNS.