1. Introduction

Carbon nanoparticles, in particular graphene both in sheet-like and in cylindrically wrapped (carbon nanotube) form, have attracted tremendous scientific and technological interest, with applications ranging from nanoelectronics ^[1] to functional and structural composites ^[2]. In structural applications, a most attractive feature of graphene is the extreme in-sheet elastic stiffness and rupture strength of perfect sp² covalently bonded structures ^[3,4].

Significant research has been devoted to assessing the influence of defects on the mechanical properties of carbon nanoparticles. This is partly motivated by the fact that experimentally measured rupture strengths e.g. of carbon nanotubes may fall significantly below theoretical expectations ^[4] for the defect-free structures, and that they exhibit a very significant statistical scatter ^[5]. Theoretical investigations of failure of defected carbon nanotubes have mostly focused on molecular mechanics (MM) or molecular dynamics (MD) simulations, see e.g. ^[6,7], although a few ab-initio studies have also been reported (e.g. ^[8]). Besides investigations of the effect of single defects, such as vacancies, Stone-Wales defects, or cracks ^[6,7], several studies have addressed the effect of multiple, randomly distributed defects. Multiple defects may not only lead to a reduction in strength due to microcrack interactions, but may also give rise to a significant scatter in strength concomitant with statistical size effects as larger samples have an enhanced probability to contain weak local configurations. Both experiment ^[5] and simulation ^[9,10] indicate that the strength distribution of carbon nanotubes with multiple defects is well described by Weibull statistics and exhibits the typical size effects characteristic of weakest-link controlled failure. Similar investigations have been conducted for Graphene ^[11,12,13].

In the present paper we aim at identifying the statistical features of failure in graphene sheets containing multiple defects. To this end we conduct an investigation of sheets with randomly distributed vacancies, systematically varying both the sheet size and the vacancy concentration, and determining for each set of parameters the statistical distribution of failure stresses as obtained from an ensemble of simulations for different initial defect configurations. Our simulation methodology is detailed in Section 2.1, while the theoretical framework used for statistically analyzing the generated data is explained in Section 2.2. Results are presented in Section 3, and we close with a Discussion and Conclusions.

---

2. Materials and Methods

---

2.1. Molecular Simulation

In this study, we perform deformation simulations with the LAMMPS package, using both MD and MM methods. We stretch graphene sheets along the armchair direction until failure. To correctly represent the atomic-level mechanisms of graphene failure, we perform our simulations with a modified REBO potential proposed by Pastewka et al. ^[14]. Standard AIREBO potentials suffer from the well-known problem that the built-in cut-off function leads to an artificial stiffening of the bond response at large bond strains and thus leads to unrealistic deformation behavior and too-high failure stresses. The commonplace method to overcome this problem is to set the inner cut-off radius (the radius where the cut-off kicks in) equal to the outer cut-off (the radius where the interaction energy becomes zero), see e.g. ^[10,12]. This method has, however, the obvious disadvantage of rendering the potential non-analytic -- a feature which makes it impossible to perform meaningful molecular mechanics simulations based upon conjugate-gradient energy minimization and which introduces artefacts into the MD simulations. The potential of Pastewka et. al. overcomes this problem by using a local environment-dependent cut-off function which leads to smooth changes in forces and to a much better description of bond breaking and re-forming as compared with conventional truncated REBO potentials. To illustrate the advantages we note that, for vacancy defects as considered here, the environment-dependent cutoff of the Pastewka potential leads to an interaction between the atoms surrounding a vacancy and to spontaneous reconstruction of two of the three dangling bonds that are left once an atom is removed. In the standard AIREBO potential the atoms around a vacancy do not interact, which is clearly unphysical, and the extensive presence of dangling bonds in multiply defected graphene samples leads to interesting simulation artefacts, for an example see results reported in ^[12].

We consider sheets with different sizes ranging from 1200 graphene lattice sites (6.5 × 5.0 nm²) to 19176 lattice sites (25 × 20 nm²) contained in the xy plane. Periodic boundary conditions are imposed in both x and y directions. Point defects (vacancies) are created by randomly removing atoms with probability f from their lattice sites. Each atom is considered independently, and removed atoms are taken out of the system. Vacancy concentrations range from f = 0.25% to 8%. For each vacancy concentration and system size we simulate 50-100 samples with different, statistically independent defect configurations.

Deformation is imposed by stretching the system volume in y direction to impose a series of fixed strains (MM) or a fixed strain rate (MD). No cross-contraction is allowed. In MM at a given strain, the configuration is optimized using a conjugate-gradient (CG) method which yields a minimal potential energy configuration. In MD, we perform simulations in the NVT ensemble where the system volume is stretched in y direction at a fixed strain rate of 0.001 ps^-1. This value is low enough to prevent the results being affected by wave propagation. The system is thermostated by using a Nose-Hoover thermostat with a relaxation time equal to 0.1 ps. Our time step is 0.5 fs, and prior to deformation the system is equilibirated over a time of 20 ps.

---

2.2. Statistical Analysis

Our statistical analysis of simulation data builds upon the Duxbury-Leath-Beale (DLB) theory of load-driven breakdown in random media ^[15]. The theory is originally formulated for a network of fuses carrying a scalar load (a current), however, the results can be transferred to networks of springs or beams carrying vectorial or tensorial mechanical loads ^[17]. Here we summarize only the main ideas, following the presentation in ^[16], while for mathematical elaboration and further detail the reader is referred to the original paper ^[15].

Consider a 2D lattice of nodes connected by load-carrying links. The lattice is assumed of square shape with a load σ applied in y direction and periodic boundary conditions in x direction. The linear size of the lattice is L and there are N = (L/a₀)² lattice sites (thus, a₀ is of the order of the lattice constant).

A fraction f of nodes is randomly removed. A cluster of n such nodes aligned perpendicular to the loading direction is termed a microcrack of size αna₀ where α is a geometrical factor which depends on the lattice morphology. The probability that adjacent vacancies form a microcrack of length a or larger is given by

$P(a) \approx f^{\displaystyle\left[\frac{a}{\alpha a_0}\right]} = \exp\left[-\frac{a}{a_1}\right] \quad, \quad a_1(f) = - \frac{\alpha a_0}{\ln f}.$

(1)

Thus, a simple geometrical argument shows that the microcrack population has an exponential size distribution. Upon loading, a certain degree of subcritical crack growth and damage coalescence may take place prior to failure. As demonstrated by Manzato and co-workers, this may increase the characteristic crack length a₁ but does not change the exponential nature of the microcrack size distribution ^[16].

Now assume that elastic interactions between different microcracks can be neglected (this approximation becomes asymptotically exact in the limit f → 0). Under a stress σ, the stress concentration at the tip of a microcrack of length a can then be estimated as ^[15,16]

$\hat{\sigma} \propto \sigma \sqrt{\frac{a}{a_2}},$

(2)

where a₂ is again of the order of the lattice spacing. Crack propagation occurs once $\hat{\sigma} > \sigma_c$, hence if

$\sigma >{{\sigma }_{c}}\sqrt{\frac{{{a}_{2}}}{a}}, a>{{a}_{2}}\sqrt{\frac{{{\sigma }_{c}}}{\sigma }}.$

(3)

In this approximation, the problem of finding the failure stress is thus equivalent to identifying the size of the largest micro-crack in the sample. There are N=(L/a₀)² possible crack initiation sites, and the probability distribution of the largest crack can then be obtained using standard arguments of extreme value statistics: The probability of finding a largest microcrack of size > a in a system of size L derives from

$1- P_{\rm max}(a, L) = \left[1 - P(a)\right]^{(L/a_0)^2}.$

(4)

Now use Eq. 3 to express a in terms of the stress σ at which the largest microcrack becomes critical. This directly yields the survival probability distribution (the fraction of samples of size L that have not yet failed at stress σ):

$P_L(\sigma) = \left(1 - \exp\left[-\frac{a_2}{a_1}\left(\frac{\sigma_c}{\sigma}\right)^2\right] \right)^{(L/a_0)^2} .$

(5)

In the limit of large system sizes this expression is well approximated by the DLB distribution

${{D}_{L}}(\sigma )=\exp \left[-{{\left( \frac{L}{{{a}_{0}}} \right)}^{2}}\exp \left[-{{\left( \frac{{\tilde{\sigma }}}{\sigma } \right)}^{2}} \right] \right]$

(6)

where $\tilde{\sigma} = \sigma_{\rm c}a_0/a_1$. This distribution has several remarkable features. In the limit of very large system sizes, the distribution flows not towards the Weibull distribution but towards the Gumbel distribution of extremal statistics ^[16]. The mean failure stress is of the form

$\langle \sigma \rangle =\tilde{\sigma }\frac{C}{\sqrt{\ln {{(L/{{a}_{0}})}^{2}}}}={{\sigma }_{\text{c}}}\frac{{C}'\ln (1/f)}{\sqrt{\ln {{(L/{{a}_{0}})}^{2}}}}.$

(7)

Here C and C' are numerical factors which depend on the lattice morphology, and σ_c is the critical stress for breaking bonds at the crack tip. We thus expect according to the DLB argument a logarithmic size dependence and also a logarithmic dependence of the failure stress on vacancy concentration f.

---

3. Results

---

3.1. Stress-Strain Curves and Failure Patterns

Figure 1 shows stress-strain curves simulated by MD and MM for pristine as well as for defected graphene samples. While the failure stresses of the defect-free sheets are almost identical, we see that for defected graphene MD predicts a lower failure stress than MM. At the same time it is evident that the MD stress-strain curve has a much more brittle appearance: failure is indicated by an abrupt load drop and there is little straining in the regime past the peak stress. In MM simulations, by contrast, we see a much more gradual decrease in stress concomitant with significant post-peak strain accumulation, hence a more ductile response.

To understand the origin of these differences we investigate the concomitant evolution of the crack patterns. Figure 2 shows emergent crack patterns during a simulated deformation experiment on a graphene sheet containing N= 19176 lattice sites with f=0.5% vacancies. The physical dimension of the sheet is 25 × 20 nm². The top series represents a MM simulation, the lower series a MD simulation performed at a temperature of T=0.1 K. Clearly the simulated failure patterns are very different: In MM, we see evidence of multiple cracking and significant interactions between microcracks. After the first crack starts to expand, the concomitant growing stress concentration activates further microcracks while sometimes, the first crack gets arrested as the reduction of the overall stress decreases the stress intensity at the crack tip. Thus system failure is not always associated with expansion of the first activated crack, and the simultaneous growth of multiple cracks leads to a quasi-ductile response where a significant amount of strain can accumulate before the sheet finally ruptures. The thus emerging multi-crack configurations represent metastable energy minima, and the quasi-ductile behavior in a MM simulation arises as the system passes through a sequence of such minima. In MD, by contrast, failure is governed by the activation and propagation of a single crack which takes the system directly to its global energy minimum, i.e., a broken sample.

These differences are at first glance astonishing given that, at T=0.1 K, thermal fluctuations are negligible as compared to typical bond energies in Graphene. However, the low-temperature MD simulation is not simply tantamount to an energy relaxation: bond breakage releases kinetic energy which cannot be instantaneously dissipated. In the MD simulation, the finite time needed for kinetic energy to dissipate is reflected by the finite response time of the MD thermostat. The presence of excess kinetic energy at the crack tip is evidenced in Figure 3, and this local excess kinetic energy prevents crack arrest and leads to a macroscopically brittle response.

We can thus conclude that, in MD simulation, dynamic effects associated with inertia may reduce the complexity of the energy landscape as such effects can help to overcome small local energy barriers and favor single-crack propagation over multiple-crack activation. Given that such dynamic effects are surely physically meaningful, this leads us to regard quasi-static simulation methods such as MM, which implicitly assume instantaneous energy dissipation, with caution.

---

3.2. Statistical Analysis

We now proceed to statistically analyze the failure stresses derived from our MM and MD simulations. We first investigate the dependency of the mean failure stress on damage fraction f and system size L. According to Eq. (7), the mean failure stress should be proportional to the logarithm of 1/f, and the inverse square of the mean failure stress should be a linear function of the logarithm of system size L² (or, equivalently, of the number of atoms N). These relations are indeed fulfilled as shown in Figure 4.

Comparing the results of MM and MD simulations, we find that for both types of simulations the f dependence of the mean failure strength is well described by the predictions of DLB statistics. The only difference between the MM and the MD results is that, in fitting Eq. (7) to the data, we must assume for MD a critical stress value σ_c that is about 8% lower than for MM (note that all other parameters refer to system size and lattice geometry and must thus be identical). This finding is interesting for two reasons: (i) The reduced stress for crack propagation cannot be attributed to assistance of thermal equilibrium fluctuations: MD simulations were at a temperature of 0.1K where such fluctuations are negligible. Rather we are dealing with a dynamic effect: Because of the finite time required for the system to thermalize, the energy released at the crack tip remains to some extent available for crack propagation and prevents the system from settling into metastable energy minima. (ii) The very significant differences in failure patterns between MD and MM do not affect the applicability of the DLB argument, which is based upon consideration of isolated cracks - this despite the fact that in MM multiple cracking and crack interactions are clearly in evidence from the damage patterns. It thus appears that the presence of multiple metastable minima in the energy landscape can be represented in the DLT picture -- which considers a simple energy landscape with a single relevant saddle point -- by the simple means of introducing a slightly increased effective propagation threshold.

We now turn to the statistics of failure thresholds and its size dependence. Figure 5 shows survival probability distributions P_N(σ) determined for a vacancy concentration of 2% and system sizes ranging from 1200 to 19176 atoms (sheet sizes ranging from 6.5 × 5 to 25 × 20 nm²). With increasing system size the distributions shift to lower stresses, as already expected from the behavior of the mean values shown in the inset of Figure 4. To assess whether the behavior of the distributions is consistent with DLB statistics we note that, according to Eq. (6), a plot of ln(-ln(D_L)/N) vs. 1/σ² (note that (L/a₀)² = N) should collapse distributions that pertain to different L (or N) onto a common straight line. Performing such a plot for P_N (inset of Figure 5) shows that the distributions indeed collapse, within the limitations of the available statistics, onto a common straight line, thus corroborating the applicability of the DLB distribution to our simulation data.

---

4. Discussion and Conclusions

Our investigation shows for the case of rupture of strongly defected Graphene sheets that simple lattice models of failure in disordered materials, such as the scalar fuse model considered by DLB and others ^[15,16,17], can provide important insights even into much more complex material systems. The energy landscape that controls microstructure evolution in strongly defected Graphene under increasing strain is characterized by multiple metastable minima indicating barriers to crack propagation, multiple microcracking, and microcrack interactions. Despite this complexity, the main statistical characteristics of failure can be envisaged within the DLB paradigm which focuses on a single critical configuration associated with geometrically adjacent defects forming a critical crack. Interestingly, the comparison of MD and MM simulations shows that dynamic effects in MD -- and most likely in reality -- may help the system to overcome intermediate energy barriers, thus effectively reducing the complexity of the energy landscape and making the problem of failure in randomly defected materials amenable to simple statistical arguments at least in the regime of not-too-strong disorder. An important physical question which we have not addressed here concerns the influence of temperature on the fracture process. Investigating this question will require an extension of the -- in essence athermal -- DLB framework to include effects of thermal fluctuations on fracture nucleation. This remains an important task for future investigations.

---

Acknowledgments

The authors acknowledge support of ZISC in carrying out the molecular dynamics simulations.

---

Conflict of Interest

All authors declare no conflict of interest in this paper.

---