Nonlinear fracture toughness measurement and crack propagation resistance of functionalized graphene multilayers

See allHide authors and affiliations

Science Advances  06 Apr 2018:
Vol. 4, no. 4, eaao7202
DOI: 10.1126/sciadv.aao7202


Despite promising applications of two-dimensional (2D) materials, one major concern is their propensity to fail in a brittle manner, which results in a low fracture toughness causing reliability issues in practical applications. We show that this limitation can be overcome by using functionalized graphene multilayers with fracture toughness (J integral) as high as ~39 J/m2, measured via a microelectromechanical systems–based in situ transmission electron microscopy technique coupled with nonlinear finite element fracture analysis. The measured fracture toughness of functionalized graphene multilayers is more than two times higher than graphene (~16 J/m2). A linear fracture analysis, similar to that previously applied to other 2D materials, was also conducted and found to be inaccurate due to the nonlinear nature of the stress-strain response of functionalized graphene multilayers. A crack arresting mechanism of functionalized graphene multilayers was experimentally observed and identified as the main contributing factor for the higher fracture toughness as compared to graphene. Molecular dynamics simulations revealed that the interactions among functionalized atoms in constituent layers and distinct fracture pathways in individual layers, due to a random distribution of functionalized carbon atoms in multilayers, restrict the growth of a preexisting crack. The results inspire potential strategies for overcoming the relatively low fracture toughness of 2D materials through chemical functionalization.


Two-dimensional (2D) materials [for example, graphene, graphene oxide (GO), boron nitride, and MoS2] have attracted great attention recently and have shown potential in a number of applications, including composites (13), sensors (46), energy storage devices (7, 8), and electronics (911). Mechanical failure is a major concern in many of these applications, particularly when materials are synthesized into large-area films. Although one commonality of this class of materials is that they have very high intrinsic in-plane strength (1215), the onset of failure is typically governed by local bond breakage at preexisting defects (13).

The capability of a material with preexisting cracks to resist fracture is defined as fracture toughness and is a critical mechanical property of interest because it determines the structural integrity and reliability of the device. Determining the fracture toughness of 2D materials, however, is quite challenging due to technical difficulties. For example, an attempt to measure the fracture toughness of MoSe2 was reported to be unsuccessful due to the inability to create artificial defects without damaging the film (16). Several experimental studies have investigated the fracture toughness of 2D materials, limited to graphene (17, 18), boronitrene (18), MoS2 (19), and MoSe2 (16). For graphene, the critical stress intensity factor has been measured to be 4 ± 0.6 MPa√m (bilayer) (17). However, owing to lacking information of the crack geometry, the intricacies of the crack tips were not taken into account. This led to the omission of crack blunting effects in these studies to which fracture toughness is known to be highly sensitive (20). To reconcile these effects, Wei et al. (18) applied finite element (FE) analysis to resolve the stresses at the crack tips and determined the critical stress intensity factors of multilayer graphene and boron nitride to be 12.0 ± 3.9 and 5.5 ± 0.7 MPa√m, respectively. A major issue, however, with this analysis is the assumption of linear elasticity for 2D materials that are instead known to exhibit nonlinear elastic (12) or plastic behavior (21), which can lead to the prediction of unrealistic fracture stresses at the crack tip (22).

Given that the fracture toughness of 2D materials is relatively low compared with conventional bulk materials (J/m2 versus kJ/m2) (23), because they typically fail in a brittle manner (1619), a key challenge is to engineer approaches to increase their fracture toughness. From past developments in bulk materials, it is well known that plastic deformation or crack arrest mechanisms can significantly enhance the fracture toughness by absorbing strain energy; however, such an attempt has yet to be demonstrated experimentally for 2D materials. Furthermore, at the nanoscale, defects play an increasingly important role in affecting fracture behavior. Recently, using molecular dynamics (MD) simulations, Wang et al. (19) proposed that by changing the defect density, the fracture in MoS2 can be modified from brittle to ductile and that the fracture toughness of defective MoS2 was predicted to exceed that of graphene. A practical alternative means of providing a crack arrest mechanism in 2D materials could come from functionalization of material surfaces with chemical groups that would pin the crack front and retard the crack propagation. Here, we investigated such a crack arresting effect and the resulting fracture toughness of GO, as a representative 2D material that contains functional groups, by using a microelectromechanical systems (MEMSs)–based in situ scanning transmission electron microscopy (STEM) technique. MD simulations were subsequently conducted to uncover details of the underlying crack arrest mechanism. By combining these investigations with nonlinear FE fracture analysis, informed through the material constitutive relation derived from atomistic simulation, the mechanisms resulting in increased fracture toughness of 2D films through chemical functionalization is demonstrated.


In situ TEM fracture toughness measurements of multilayer GO

A transmission electron microscopy (TEM)–compatible monolithic MEMS device (24) was used to perform fracture toughness measurements of multilayer GO nanosheets under STEM imaging (Fig. 1). First, GO nanosheets were suspended over the MEMS actuation shuttles by drop casting a GO solution using a micropipette. Upon evaporation of the liquid, GO was left suspended over the MEMS shuttles. The GO films used for this study have an oxygen functionalization of 20% as measured via x-ray photoelectron spectroscopy (XPS), and their chemical structure has been reported in detail previously (25). The MEMS device was then mounted in a custom-made TEM holder and placed inside the TEM chamber (Hitachi HF-3300 TEM/STEM) for tensile testing. GO samples were examined before testing to make sure no significant flaw existed in the nanosheets. Before creating an artificial crack in each nanosheet, electron energy-loss spectroscopy (EELS) inelastic mean free path (IMFP) mapping was conducted for the measurement of sheet thickness (see Materials and Methods). An artificial crack was then introduced to the GO nanosheet via electron irradiation with a beam energy of 300 keV (Fig. 1, E to G). Note that an electron beam was used instead of focused ion beam (FIB) [previously used for precracking of graphene (17)] to eliminate any effects of Ga+ deposition and implantation. In addition, catastrophic failure has been reported to occur while using FIB to create artificial cracks in other 2D materials such as MoSe2 (16). The use of electron beam precracking ensured less unwanted damage in the regions adjacent to the desired area when creating artificial cracks. We have also developed a “hole-by-hole” irradiation procedure to create nanosized cracks with improved control over size and shape (see figs. S1 and S2). The ratio of lateral length of the crack-to-film width was maintained to be ~10% for all tested GO naonsheets.

Fig. 1 In situ TEM setup for tensile test of precracked GO nanosheet.

(A) MEMS-based in situ TEM/STEM/SEM (scanning electron microscopy) experimental setup. (B) MEMS device placed inside a custom-made TEM holder. (C) SEM image of the TEM-compatible monolithic MEMS device, as dash box in (B). (D) Higher magnification SEM image showing actuation shuttles with a ~2.5-μm gap [dash box in (C)]. (E) SEM image of pristine GO nanosheet suspended over the actuation shuttles in (D). (F) SEM image showing high-energy electron beam–irradiated GO nanosheet (34 nm; sample 3) in (E); a ~1.2-μm crack was introduced. (G) High-resolution SEM image of the artificial crack in (F). The entire experiment was conducted in STEM mode and Hitachi HF-3300 STEM/TEM, which has a secondary electron detector that allows capturing SEM images in STEM mode.

After electron beam irradiation, displacement-controlled uniaxial tensile tests were performed on four different precracked GO nanosheets. Stress-strain measurements (Fig. 2A) before the first subcritical propagation of the artificial crack were recorded following our previously reported method (14). The thicknesses of the four samples of GO nanosheets ranged from ~15 to ~132 nm. Detailed geometry of samples can be found in table S1. Snapshots of various stages of the GO nanosheet (34 nm thick; sample 3) during the tensile test are presented in Fig. 2B and movie S1. Under continued loading, the first subcritical crack, which was almost perpendicular to the tensile direction, propagated from the left tip of the crack to the left of the nanosheet but did not develop all the way to the end. Instead, the crack stopped halfway, and another crack began to initiate from the other crack tip, which also stopped halfway (movie S1) and subsequently propagated to the edge of the nanosheet. This behavior indicates an improved crack resistance compared with graphene (17). Furthermore, another crack front opened near the freshly formed crack tip on the left of the nanosheet, propagated in a direction not orthogonal to the loading direction toward the edge of the nanosheet causing complete fracture. Once strain was released, the gap in between the separated nanosheet was closed. Note that overlapping of films was observed after the load was completely released (Fig. 2B), suggesting that the plastic deformation within the nanosheet occurred during the tensile test. Similar fracture behavior was also observed in other GO nanosheets tested in this work (see figs. S1 and S2).

Fig. 2 Fracture behavior of GO nanosheet under tension.

(A) Stress-strain response before subcritical crack propagation (captured during in situ TEM uniaxial tensile testing). Data were recorded for all the four tested samples of GO nanosheets with different thickness. Inset: Electron irradiated precracks of four different samples. (B) Representative snapshots of propagation of the artificial crack (34-nm-thick sample). The last image is a high-magnification image showing film overlap after load release.

Linear fracture analysis of multilayer GO

The fracture toughness of 2D materials such as graphene (17), h-BN (18), MoS2 (19), and MoSe2 (16) has recently been studied using the postulates of linear elastic fracture mechanics via critical stress intensity factor (KIC) and critical energy release rate (GIC), which are suitable measures for isotropic linear elastic brittle materials. Here, the same approach was applied first to calculate KIC and GIC of the GO nanosheets (see the Supplementary Materials) for purposes of comparison and demonstration. KIC and GIC of all the four tested GO nanosheets are summarized in Table 1. KIC for GO nanosheets was found to be similar to that of graphene (~4 MPa√m) (17); however, because GO nanosheets have a significantly lower Young’s modulus (E) as compared with graphene, the GO nanosheets were estimated to have significantly higher GIC values (GIC = KIC2/E, 76 to 262 J/m2), than those previously reported for graphene (15.9 J/m2) (17). Accordingly, this suggests that GO nanosheets are able to absorb a higher amount of strain energy than graphene before critical fracture. Note that critical strain energy release rates are reported as a range due to the uncertainty associated with the magnitudes of E as function of the thickness of the GO nanosheets. A lower bound of GIC is estimated by considering the modulus of monolayer GO (26), whereas the upper bound is determined using the modulus of GO nanosheets (14).

Table 1 Summary of fracture toughness measurements.
View this table:

To determine the critical stresses in the vicinity of the crack tips during fracture initiation and to take into account crack tip blunting effects, we performed FE-based fracture analysis (samples 1 and 2) using ANSYS 16. An image processing algorithm was used to create 3D cracked models (see fig. S3) of GO that were identical in shape and size to samples 1 and 2 (Fig. 3, A and B). MD simulations were performed to obtain the stress-strain relations of pristine GO nanosheets with up to nine layers, without preexisting cracks (see figs. S4 and S5A). The simulated stress-strain response of defect-free monolayer GO and GO nanosheets under uniaxial tensile loading shows two distinct regimes of pre-ultimate tensile strength (UTS) (see fig. S5A), including a linear Hookean response followed by a nonlinear stress-strain response with softening at large strains. Although this observation is not directly evident from the experimental stress-strain data due to the limitation in force resolution of our MEMS devices, the overlap of fractured films after complete load release and the crack arrest behavior (Fig. 2B) indicate the occurrence of plastic deformation. These results are also consistent with previous studies that suggested nonlinear behavior for hydroxyl-dominated monolayer GO using density functional theory calculations (13).

Fig. 3 Fracture toughness analysis of GO nanosheet.

(A) SEM image of a precracked GO nanosheet before loading (sample 1) for which J integral calculations were performed. (B) FE analysis stress contour (GPa) in the loading direction of sample 1 corresponding to the maximum experimentally measured far-field stress immediately before the propagation of the crack. (C) High magnification of the stress distribution near the crack tip in (B). The region colored orange experienced stress in the nonlinear stress-strain regime. (D) Stress as a function of distance from the crack tip. A nonlinear zone from the crack tip to 5.6 nm is shaded. (E) Variation of JIC along the crack front as a function of the number of contours considered for the domain integral method.

In the FE simulations, strain-controlled quasi-static load was imposed normal to the crack face until the far-field stress reached the experimentally measured critical stress of 5.4 GPa for sample 1 and 5.2 GPa for sample 2 (Fig. 2A). For sample 1, a linear stress-strain analysis using a Young’s Modulus value of 532 GPa, obtained from our MD simulated stress-strain response, and a Poisson’s ratio of 0.165 (13, 14) predicted an unrealistically large stress of 96 GPa in the elements ahead of the crack tip. This large stress is significantly higher than the UTS of multilayer GO (~46 GPa) obtained using MD simulations. The contours of the principle stress component parallel to the loading direction obtained from the linear analysis are shown in fig. S6A. These results suggest that an FE-based linear analysis can produce unrealistically high stress fields near the crack tip for materials with nonlinear mechanical properties such as GO. A similarly unrealistic prediction was previously reported when a linear fracture mechanics analysis was applied to bilayer graphene; wherein a linear FE analysis predicted a stress of 360 GPa at the crack tip (22), which is much higher than the in-plane strength of graphene. In agreement with our findings, the authors attributed the unrealistic high stress to the linear material property-based FE analysis they used and suggested that a nonlinear material constitutive law should be used.

Nonlinear fracture analysis (J integral) of multilayer GO

We performed a nonlinear stress analysis using the MD simulation–obtained constitutive law of a nine-layer GO nanosheet. It is well known that the mechanical properties of GO nanosheets are sensitive to the sample thickness (14, 21), as the elastic modulus and strength deteriorate steeply with increasing number of layers. The thickness of sample 1 (14.6 nm) is closest to that of nine-layer GO (~5.5 nm). In addition, the uniaxial stress-strain relationship obtained from MD simulations predicted similar mechanical response for multilayer GO samples. Figure 3 (B and C) shows the contours of principle stress at the experimental fracture strength of sample 1 obtained from the nonlinear analysis. It can be seen that when the model of the GO nanosheet was loaded to the experimentally measured far-field stress of 5.4 GPa (Fig. 2A), the stresses in the elements around the crack tip reached approximately 50 GPa, which is in good agreement with the estimates of UTS obtained from our MD simulations. The nonlinear fracture analysis also enabled us to capture the nonlinear zone at the crack tip. As shown in Fig. 3D, a 5.6-nm nonlinear zone (shaded) can be observed (also shown by orange color in Fig. 3C) as the stress-strain response of multilayer GO nanosheets deviates from linearity above a stress level of 26.3 GPa (fig. S5B).

A nonlinear description of fracture toughness was subsequently obtained by calculating the J integral approach introduced by Rice (27) and Hutchinson (28). It is defined as a path-independent integral that can be imagined around the crack tip and can be viewed also as an energy release rate parameter. For a 3D crack problem, assuming a planar crack surface at any point on the crack front, J is defined locally, and it varies along the crack front. In a 3D form, J integral for points along the crack front is given byEmbedded Image(1)where W is the strain energy density, nk is the unit normal vector to the integration path in the outward direction, ti is Cauchy stress tensor, ui is the displacement vector u of node i, ui,k is the derivative of ui versus the kth element, s is the isoprametric coordinate, σij is the stress tensor, δ is the Kronecker δ, γ is the 2D integration path surrounding the crack tip, and A is the domain bounded by γ. Therefore, in a 3D formulation, the magnitude of the vector J(η) is given by the sum of a line integral term (first term on the right-hand side of Eq. 1) and a surface integral term (second term on the right-hand side of Eq. 1). Displacements and stresses are calculated directly by FE analysis and therefore, the line integral term is readily calculated. On the other hand, the calculation of the surface term requires second derivatives of displacements. The surface term can be insignificant if J is constant with respect to the crack front coordinate; however, the magnitude of this term can be large if there are strong gradients in stresses and strains, that is, at the specimen surface (29). To avoid the influences of sharp curvatures, we consider the value of J inside the crack front (see fig. S8, A and B). The magnitude of JIC (JIC represents the J integral corresponding to crack initiation) for the crack tip nodes away from the surface (that is, inside the bulk), as shown in Fig. 3E, is evaluated as 39 J/m2 for crack tip T1 and 34 J/m2 for crack tip T2 in sample 1. This magnitude of JIC is significantly higher than graphene, which has a GIC of 15.9 J/m2 (JIC = GIC for a linear material) (30). Note that the 3D stress-strain response of GO is dependent on the assumption of thickness, as 3D stress values are obtained by dividing 2D stresses by the thickness of the sample. For example, Meng et al. (31) assumed a thickness of 0.75 nm for monolayer GO, whereas our MD simulations predict the thickness to be 0.62 nm. Therefore, when directly comparing the magnitude of J, the thickness assumption must be considered. Furthermore, note that there is a dramatic difference in fracture toughness values obtained by linear versus nonlinear analysis (80 to 92 J/m2 versus 34 to 39 J/m2; Table 1), which additionally justifies the need for nonlinear fracture analysis for GO nanosheets. Magnitudes of J integral can be found in figs. S8 and S9). Note that the simplistic MD model used for evaluating material stress-strain behavior was pristine in nature and devoid of crystallographic defects (point defects, dislocations, and grain boundaries) and irregularities which were found to exist in the experimental samples, such as wrinkles and folds. In addition, the experimental samples could have contained grain boundaries which are typically found in 2D nanosheets. Previous research has shown that these defects or irregularities affect both modulus and strength (14). In addition, the thickness of the MD model was smaller than both samples 1 and 2. These factors could have resulted in a systematic error in the J integral estimates in coupling atomistic simulations and continuum-based FE simulations.

Atomic origins of crack arresting mechanism

The observed higher fracture toughness of multilayered GO nanosheets as compared to that of graphene stems from the difference in their fracture mechanisms. MD simulations were performed to further investigate the fracture behavior of precracked multilayered GO and multilayer graphene nanosheets. A central crack (2.5 nm × 0.5 nm) was introduced in the MD models of one-, three-, four-, and six-layered GO (fig. S10). The crack face was taken parallel to the zigzag direction. As loading was applied, the width of the cracks in all the samples gradually increased, resulting in crack blunting. Under continued loading, mechanical fracture of GO nanosheets was observed to initiate from the crack tip and propagate through the basal plane (Fig. 4, A to G), conforming to a Mode 1 fracture similar to our experimental finding. As shown in Fig. 4J, the precracked GO samples exhibit a nonlinear stress-strain response under uniaxial tension. Fractured monolayer GO has the largest critical stress and undergoes sudden brittle fracture. The critical fracture strength of multilayer GO samples is lower than that of the monolayer, and the multilayer GO samples exhibit progressive fracture during crack propagation.

Fig. 4 Atomistic origin of enhanced fracture toughness of GO nanosheet.

(A) The initial configuration of a layer in a precracked four-layered GO sample and the schematic of the constraints. (B to G) Snapshots of the distinct stages of fracture path through one of the layers in a precracked four-layered GO sample during uniaxial loading. Atoms in the precracked layer are colored on the basis of their shear strain magnitudes. For clarity, the functional groups are hidden in these images. Here, atomic strain is defined as the von Mises shear strain invariant of the atomic Green-Lagrangian strain tensor, which may be derived directly from the definition of the local deformation gradient. (H) Uniaxial stress-strain response of the layer and labels in stress-strain curve refer to MD snapshot panels in this figure. (I) Fracture pathways of constitutive layers in a four-layer GO sample. (J) MD simulations: Uniaxial stress-strain response of precracked monolayer and multilayer GO nanosheets subjected to uniaxial tensile loading. (K) Energy absorbed per unit volume by GO nanosheets of varying thicknesses during crack propagation. It is calculated as the area under the stress-strain curve in (J) from the point of fracture initiation to the 50% of their peak stress, where the sheet is considered to have failed completely.

The significant difference between multilayer and monolayer GO in the stress-strain response after the onset of crack initiation was further studied to understand the role of interlayer interactions provided by the functional groups. Figure 4K shows the area under the stress-strain curves of monolayer and multilayer GO samples from the point of initial fracture to a level that is 50% of their critical stress magnitude (we considered a 50% loss in tensile strength as equivalent to the final fracture of the nanosheet). This area represents the energy absorbed by GO nanosheets per unit volume during fracture propagation before catastrophic failure. There are two components to this energy: the energy required to create additional surfaces as the fracture propagates and the energy spent to deform the entire structure. It is evident that the higher the number of GO layers, the higher the energy required to propagate cracks to fracture. In particular, multilayers of GO require significantly higher energy as compared to monolayers of GO. For example, the magnitude of energy absorbed by the six-layer GO is four times higher than that of the monolayer GO.

To explore the mechanism responsible for this higher amount of strain energy in multilayer GO, we analyzed the stress-strain response of individual constituent layers in the precracked multilayer GO samples (the results for the four-layered case is shown in fig. S11). For all of the multilayered GO samples (that is, bilayer, three layer, four layer, and six layer), stress at fracture initiation and strain to failure are observed to be different among the individual GO layers. To understand the origin of this difference among individual layers, we analyzed the atomistic configurations of each of the constituent layers individually from the onset of loading to ultimate fracture. Distinct stages of crack growth in one of the constituent layers in a four-layer GO is shown in Fig. 4 (A to G). It can be seen that fracture initiates from one of the crack tips, nearly perpendicular to the loading direction, and before this crack becomes critical, another crack front opens from a secondary crack tip and continues to grow. Tiny peaks and corresponding stress drops in the tensile responses correspond to these subcritical cracking events, before these microcracks grow unstably or combine to form a critical crack. For example, the stress-strain response corresponding to the atomistic configurations in Fig. 4 (A to G) is shown in Fig. 4H. Ultimate fracture in this case occurs when both crack fronts grow completely and meet the edge of the nanosheet. This development of crack growth is in very good agreement with the experimental observations using STEM (see Fig. 4). Despite that the GO nanosheets used for fracture experiments contained wrinkles and defects, although the MD models used did not have any of such complexities, MD simulations can still capture the crack arrest behavior during propagation. These results demonstrate that the plasticity in multilayer GO originates from the complex interactions among the functional atoms sandwiched between carbon basal planes. The plasticity observed in this class of materials differs from a classical description of plastic deformation governed by the motion of dislocations, twins, and grain boundaries. However, combination and variation of these defects may potentially affect the fracture toughness measurement reported here because it was previously shown by molecular simulations that the size of grain boundaries and the presence of triple junctions can affect the fracture toughness of graphene (3234). In those studies, toughness was predicted to be inversely affected by the presence of grain boundaries in relation to the crack tip. Comprehensive studies of the effect of different types of defects and grain boundaries on the fracture toughness of functionalized graphene multilayers are worthy of future investigations. Furthermore, it is interesting to note that the fracture pathways of each individual layer were found to be significantly different (Fig. 4I). Depending on the distribution of the random functionalization groups, once initiated, the crack grows because of the combined effect of global loading condition and localized strain fields, tracing the carbon atoms attached to these functionalized carbon atoms. Furthermore, during crack growth, the strain energy is released through the neighboring hydrogen bonds in the functional groups in the adjacent carbon layers (14). Previous studies have shown that the interactions between functionalized atoms (attractive forces) in GO lead to movements similar to stick slip–like motion (35). As a result, these atomic interactions between individual atomic layers result in the crack arresting mechanism unique to GO. For example, a similar crack-resisting behavior would not be expected for graphene because there is no functionalization but carbon atoms connected only by sp2 bonds. For example, previous experimental fracture studies of graphene do not exhibit this crack arresting behavior (17). In addition, the tensile response of multilayer graphene with a precrack studied herein (fig. S12) reveals that the fracture pathways of individual atomic layers within the multilayer are identical in contrast to that shown for GO.


In summary, a nonlinear fracture toughness analysis (J integral) was applied to measure the fracture toughness of a multilayered GO to be ~39 J/m2, which is more than two times higher than those previously reported for graphene. The nonlinear fracture analysis was shown to be a more accurate approach for 2D materials that exhibit nonlinear stress-strain behavior. In addition, a unique fracture arresting behavior of multilayered GO was experimentally observed. Atomistic simulation was performed to identify the mechanism responsible for this behavior and the resulting high fracture toughness. It was found that interactions among functionalized carbon atoms, as well as distinct fracture pathways in individual layers, result in localized strain fields, which require more strain energy and inhibit crack growth. This fracture arresting behavior is unique to 2D films with functional groups that bridge multilayers and is not predicted for multilayer graphene. The results inspire potential strategies for overcoming the relatively low fracture toughness of 2D materials through chemical functionalization.

While the current manuscript was in the review process, an article discussing the fracture toughness of multilayer graphene was published (36). It is interesting that in the study, multilayer graphene was reported to have a higher fracture toughness than its monolayer counterpart due to weak interactions between adjacent layers; and fracture asynchronization was observed.


In situ TEM on-chip thickness measurement of GO nanosheet

IMFP mapping of GO nanosheets was captured using EELS with beam energy of 300 keV in TEM mode. By taking the average number of IMFP across a line profile on the map, the average IMFP of each GO nanosheet was calculated (24). The IMFP (λ) of GO was estimated using equation (37) Embedded Image, where E0 was the incidental beam energy (300 keV), β was the collection semi-angle (35 mrad), Em was the mean energy loss of a material (14.2 for arc-evaporated carbon film) (37), and F was a relativistic factor that was calculated to be 0.514 using equation Embedded Image. IMFP (λ) was calculated to be 158 nm. Upon multiplying the average IMFP by λ, the thickness of each set was estimated.


MD simulation. We have performed MD simulations to study fracture pathways in multilayer GO samples with preexisting cracks. In addition, MD simulations have also been used to investigate the material behavior of multilayer GO under uniaxial tension along the armchair direction. The GO nanosheets were 22% functionalized in the basal plane with an oxygen-to-carbon (O/C) ratio of approximately 1:4, a composition determined from XPS measurements. In the multilayer samples, different layers were stacked in an ABAB atomic configuration to be consistent with experimental measurements of stacking in bulk GO (38). The Lerf-Klinowski algorithm has been used to generate monolayer and multilayer GO. Periodic boundary conditions were applied along the in-plane directions, and in the thickness direction, a vacuum of 20 Å was used on either side. The temperature of the GO sample was raised to 300 K and subsequently thermally equilibrated in an NPT ensemble using a simulation time step of 0.25 fs. Here, N stands for number of atoms, P represents pressure, and T is temperature, and inside the ensemble, N, P, and T, were maintained constant. A Nose-Hoover thermostat was used for temperature control with a damping parameter of 100 fs. MD tensile simulations were performed using LAMMPS (39), and interatomic interactions were accounted for using the ReaxFF (40) potential, which has been used for predicting the mechanical properties of GO with a wide variety of composition and coverage (1, 25, 41).

Accuracy of interatomic potential. The ReaxFF potential used here was a part of a series of variable charge bond-order potentials in which the total energy of the system is described by bonding including, Coloumbic, overcoordination, and van der Waals energies. ReaxFF potentials have previously been used in a number of experimental coupled computational works for studying the mechanical properties of GO systems (21, 25). The estimates of various mechanical properties of GO at 300 K reported were in excellent agreement with the ground-state density functional theory based estimates reported by Meng et al. (31). For example, our MD simulations based on ReaxFF potential predicted a Young’s modulus (E) and UTS of 330 and 29.7 N/m, respectively, whereas Meng et al. (31) reported an E of 300 N/m and a UTS of 30 N/m. Moreover, a recent paper has also reported a successful use of similar atomistic-coupled FE simulations informed by experimental fracture stresses to calculate JIC of bilayer graphene samples (42).

Size-dependent constitutive law. The stress-strain response obtained from MD, which formed as an input to macroscopic FE model, was also size-dependent (especially dependent upon the number of GO layers); MD could not currently access large sample sizes present in experiments.

Effect of strain rate. Previous computational studies have pointed out that the mechanical properties (strain to failure, strength) of graphene were sensitive to the imposed strain rate (4345). Here, the MD simulations were performed using a strain rate of 109 s, which was significantly larger than the strain rates used in our experiments. It is possible that significantly higher strain rates in simulations may lead to overprediction of strength, while plasticity effects may get suppressed because sufficient time was not provided in simulations for nucleation and growth of failure processes. This strain rate sensitivity may result in a systematic error in the FE-based stress analysis and consequently in the reported value of J. However, estimation of the stress-strain response of GO nanosheets at strain rates comparable to those used in experiments using molecular simulations was computationally too expensive to carry out, and the consequent error is thus presently unknown.

Multiscale error propagation. Note that the MD model used for evaluating material stress-strain behavior considered that the material was pristine in nature and devoid of crystallographic defects (point defects, dislocations, and grain boundaries) and irregularities that were known to exist in the experimental samples, such as wrinkles and folds. In addition, the experimental samples likely contained grain boundaries that were typically found in 2D nanosheets. Previous research has shown that these defects or irregularities can affect both modulus and strength (14). In addition, the thickness of the MD model used was smaller than both samples 1 and 2. Generally, topological and grain boundary defects were expected to decrease stiffness and strength properties of GO, and thus, our calculations may have overestimated JIC.


Supplementary material for this article is available at

section S1. Hole-by-Hole crack creation via electron beam and additional experimental results

section S2. Applying Griffith theory of brittle fracture to calculate fracture toughness

section S3. Summary of sample geometry

section S4. FE-based stress analysis and calculation of J

section S5. Uniaxial tensile simulations of pristine multilayered GO

section S6. Linear and nonlinear stress analysis

section S7. J integral calculations

section S8. MD simulations of fracture in multilayered GO

fig. S1. Crack creation via electron beam and tensile testing of sample 4.

fig. S2. Crack creation via electron beam and tensile testing of sample 1.

fig. S3. FE models and mesh.

fig. S4. MD simulation model.

fig. S5. Nonlinear mechanical properties of GO nanosheets from MD simulations.

fig. S6. Stress analysis ahead of the crack tip in sample 1 using FE simulations.

fig. S7. Stress analysis ahead of the crack tip in sample 2 using FE simulations.

fig. S8. Crackfront nodes and corresponding JIC values for sample 1.

fig. S9. JIC for sample 2.

fig. S10. Atomic configurations of cracked AB-stacked GO sheets.

fig. S11. Stress-strain response for four-layer GO sheet.

fig. S12. Fracture in four-layered graphene.

table S1. Summary of sample geometry.

movie S1. In situ TEM tensile test of functionalized graphene multilayers.

References (4655)

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank M. Daly and Y. Zuo for helpful discussions on modeling. Funding: We acknowledge financial support from the Ontario Ministry of Research and Innovation Early Researcher Award, the Erwin Edward Hart Endowed Professorship, the Canada Research Chairs Program, the Natural Sciences and Engineering Research Council of Canada, Mitacs, and the Canada Foundation for Innovation. TEM analysis was carried out at Ontario Center for the Characterization of Advanced Materials. MD calculations were performed at the SciNet and Calculquebec consortia. SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute Canada, the Government of Ontario, Ontario Research Fund-Research Excellence, and the University of Toronto. Author contributions: C.C., S.M., Y.S., C.V.S., and T.F. conceived the idea and designed the project. Y.S. and T.F. directed C.C. to conduct sample preparation, in situ experiments, and analyze the experimental data. C.V.S. directed S.M. to perform MD and FE simulations. D.D.P. directed J.Y.H. to conduct TEM imaging and TEM analysis. C.C. and S.M. wrote the manuscript under guidance from Y.S., C.V.S., and T.F. All authors discussed the results and edited the manuscript. C.C. and S.M. contributed equally to this work. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article