Atomistically derived cohesive zone model of intergranular fracture in polycrystalline graphene

Since its discovery by Novoselov and Geim,¹ graphene, a two-dimensional allotrope of carbon, has generated extensive interest owing to its extraordinary electric and mechanical properties, and the variety of potential applications which they may allow. Pristine graphene is the strongest 2D material ever measured with a Young's modulus of 348 N/m and intrinsic strength of 39.5 N/m.² Upon normalizing by the distance between graphite basal planes (0.335 nm), the Young's modulus is equivalent to ≈1 TPa and the intrinsic strength to about 100 GPa. These properties are, however, sensitive to the presence of defects and, in particular, one-dimensional defects such as grain boundaries (GBs). The latter are inherent to scalable methods of graphene production such as chemical vapor deposition (CVD).³ The specimens produced in this way are polycrystalline in nature but retain the two-dimensional character of graphene. Their size may be millimetric to metric,⁴ and the grain size may vary from very small (0.5 to 1 μm) to large (tens to hundreds of micrometers). They would be more suitable for a variety of applications than pristine graphene obtained by mechanical exfoliation, provided their mechanical properties remain attractive.

The mechanical properties of defect-free graphene have been investigated experimentally by Lee et al.,⁵ who performed indentation experiments of free-standing circular monatomically thin membranes of graphene using the diamond tip of an atomic force microscope (AFM). Later, Wei et al.² established from density functional theory (DFT) calculations a continuum constitutive law for pristine graphene which is modeled as an anisotropic nonlinear elastic material. This constitutive law has been validated referring to the aforementioned experiments by the development of a multiscale model combining the continuum description provided by the DFT calculations and a finite element method (FEM) simulation of the indentation of the circular graphene membrane.⁶ 

In 2013, Lee et al.⁷ performed a set of indentation experiments to investigate the mechanical properties of polycrystalline graphene produced by CVD with small and large grain sizes. The monatomically thin graphene samples were transferred onto a silicon substrate patterned with an array of circular wells (with diameters of 1 and 1.5 μm). The free-standing membranes were indented with a commercial nanoindenter up to rupture, see Fig. 1 adapted from Ref. 7 and notice in Fig. 1(c) the false-colored graphene grains separated by GBs as revealed by dark field transmission electron microscopy (DF-TEM).

These experiments provide two important results for our own study. First, indentation experiments were performed on large grain CVD graphene for which the grain size is significantly larger than the well diameter that yield free-standing graphene specimens without GBs. Thus, the specimens did not contain grain boundaries, but they contained uncharacterized zero-dimensional point (i.e., atomic) defects such as atomic vacancies or substitutional atoms. Statistical analysis of the results reveals no significant difference in both the elastic stiffness and the breaking strength between pristine exfoliated and CVD produced large grain graphene. This suggests that the point defect density may be neglected when modeling polycrystalline graphene with sufficiently low point defect density. Second, indentation experiments of small grain CVD graphene, for which the free-standing specimens contain GBs, show a statistically significant reduction in strength on the order of a few percent. Moreover, when GBs are directly indented, the weakening effect may reach 15% and TEM images suggest that cracks are likely to initiate at the grain boundaries. Thus, the GBs are identified as constituting the main defects necessary to model mechanical failure in CVD graphene.

As a way to model crack initiation and propagation at the GB, a cohesive zone model (CZM) may be introduced in the FEM calculations of polycrystalline graphene deformation. The first step of this multiscale approach is to establish the characteristics of the CZM and this will be done at the atomistic scale, using molecular dynamics (MD) to simulate crack propagation along the GB of graphene bicrystals and, thus, derive a traction-separation law (TSL) characterizing the properties of the interface. The modeling of the mechanical properties of the GBs of graphene through the development of a characteristic TSL constitutes the main objective of this paper.

The CZM^8,9 is a classical concept in fracture mechanics which has found its implementation in FEM simulations to account for crack initiation and propagation.¹⁰ Intergranular fracture may be predicted by embedding cohesive surface elements along grain boundaries that incorporate a traction-separation law characterizing the interface properties.

Various functional forms of the TSL (bilinear, trapezoidal, exponential, and polynomial) have been proposed for both ductile and brittle materials.^8,9,11–13 They account for the fracture process at a macroscopic scale and are often phenomenological laws derived from macroscopic experiments. In the case of polycrystalline graphene, the approach has to be somewhat different since direct experimentation is not available. However, owing to the 2D nature of the material, the grain boundaries are one dimensional. They can be imaged by high resolution TEM and can be idealized in a fairly accurate way by periodic patterns of aligned defects in the honeycomb crystal lattice that characterizes graphene. The TSL can then be established by modeling, at the atomic scale, the crack propagation along grain boundaries. Molecular dynamics simulations that account for the aggregate behavior of hundreds of thousands of atoms are a powerful tool to establish a TSL for graphene.

Yamakov et al.¹⁴ have proposed a methodology to derive from MD simulations a CZM for intergranular fracture processes in aluminum and to incorporate it in continuum simulations. We shall use the main elements of that methodology with the objective of obtaining a quantitative TSL for intergranular fracture in graphene. Our need for a CZM that may be readily and consistently implemented in the FEM leads us to complement existing work^14–16 with, first, an energetic validation of the thermomechanical process accounted for by the TSL and, second, an analysis of the scale effect on the parameters of the TSL with regards to the FEM mesh size. To achieve these goals, the fracture process that we simulate by MD is the displacement controlled fracture of a monatomically thin bicrystal of graphene in the form of a double cantilever beam (DCB).

The paper is organized as follows. We review, in Sec. II, fundamental concepts for describing pristine and polycrystalline graphene and some important results on their mechanical properties. In Sec. III, we present the methodology to derive a cohesive zone model from molecular dynamics and implement it via simulation of controlled crack propagation along a GB. The results of our simulations on a few representative GBs are presented and analyzed in terms of energy balance and mesh size effect in Sec. IV. We conclude by discussing in Sec. V the range of validity of this original approach and its potential extensions.

Graphene is the newest experimentally accessible allotrope of carbon. A sheet of graphene is similar to a tiling of benzene where the hydrogen is replaced by carbon atoms to form neighboring hexagons. The atoms are arranged in a 2D regular honeycomb lattice due to their sp² hybridization. This lattice is not strictly a Bravais lattice since two neighboring sites are not equivalent. It may be viewed as a hexagonal Bravais lattice with a two-atom basis. Fig. 2(a) shows the two vectors (⁠$a1$ and $a2$⁠) that constitute along with an out of plane vector, the commonly used basis of the 3D hexagonal Bravais lattice. The distance between nearest neighbor carbon atoms is d = 0.142 nm, which is the average of single and double covalent bond distance for C. Hence, the lattice spacing that corresponds to the norm of basis vectors $a1$ and $a2$ is $a0=3d=0.246 nm$⁠. A third vector $a3$⁠, in the plane of (⁠$a1, a2$⁠), is usually introduced for the definition of the Miller–Bravais notation of crystallographic planes and directions.¹⁷ 

Two orthogonal directions within the crystal lattice referred to as the zigzag and armchair directions can be expressed using the Miller–Bravais notation as $[112¯0]$ and $[11¯00]$⁠, which are parallel to the x-direction and y-direction, respectively, in Fig. 2(a).

A graphene polycrystal is an assemblage of single crystals separated by grain boundaries that are 1D line defects. We are interested in the simplest of such structures: bicrystals, which are two crystalline domains linked by a GB. In order to characterize a bicrystal, two parameters are necessary. For instance, one may use θ_L and θ_R (⁠$0≤θL,θR≤30°$⁠) defined as the angles between the unit normal vector to the GB and a particular crystallographic direction of the left and right crystals, respectively (see Fig. 3). The misorientation angle θ between the two grains is expressed as: $θ=θL+θR$ if $θL+θR≤π/6$ and $θ=π/3−(θL+θR)$ if $θL+θR>π/6$⁠. The equality $θL=θR$ defines the symmetric GBs.

In 2D materials, a GB is a 1D chain of edge dislocations. Yazyev and Louie¹⁸ have shown that the atomic structure of a dislocation in graphene can be considered as a pair of positive and negative disclinations which consist of a five-ring and seven-ring atomic core, respectively. Hence, graphene GBs contain pentagonal and heptagonal carbon rings.

GBs observed with high-resolution transmission electron microscopy (HR-TEM)^3,19,20 appear to be composed of mixed regions of periodic and aperiodic sequences of dislocations resulting in sinuous GBs.²⁰ On the other hand, theoretical studies of graphene GBs^18,21,22 deal with the idealized periodic structures. It is important that the periodic structures consist of the periodic subsequences observed in TEM images of GBs. In particular, when the idealized periodic GB retains the sinuous character of a real GB, as is the case of our two model GBs, it has been shown that the idealized model accounts very well for the mechanical properties of the original GB.^20,23

Alternately, a GB may be characterized by the components (n_L, m_L) and (n_R, m_R) of the two periodic translation vectors $dL$ and $dR$ of the two grains expressed in the respective $(u,v)$ basis of the underlying hexagonal lattices (Figs. 2(b) and 2(c)). These vectors should match each other along the GB line to constitute the repeating vector of the GB. The GB is thus denoted as $(nL,mL)|(nR,mR)$⁠. This two-vector nomenclature, previously used by Refs. 20, 21, 23, and 24, encompasses the two angles that characterize a grain boundary and make apparent the symmetric or asymmetric character of the GB; in addition, one may distinguish GBs for which the matching between the two crystal lattices is exact or not. In fact, the length L of each periodic translation vector is computed from its components (n, m) and the lattice spacing a₀: $L=a0n2+nm+m2$⁠. For all symmetric GBs (⁠$(nL,mL)=(nR,mR)$⁠) and some asymmetric GBs (e.g., $(5,3)|(7,0)$ in Ref. 21), these periodic vectors respect the commensurability condition in which they have the same norm. For others asymmetric GBs, there is a mismatch between the two vectors, as for instance $(7,0)|(4,4)$ in Fig. 2(c) for which $LL=49a0$ and $LR=48a0$⁠. This small mismatch, responsible for a higher GB energy, is accommodated by local lattice distortions and the resulting repeating vector lies in between the two original vectors.

In our simulations, we shall consider two GBs. The symmetric tilt GB $(3,1)|(3,1)$ and the asymmetric tilt GB $(7,0)|(4,4)$ as presented in Figs. 2(b) and 2(c).

Graphene exhibits a nonlinear, anisotropic elastic behavior. From density functional theory (DFT) calculation, Wei et al.² derived a continuum constitutive relationship suitable for incorporation into the finite element method. This relationship results from a Taylor expansion of the elastic strain energy in strain truncated after the fifth-order term (Fig. 4). It provides a continuum description of graphene, valid for finite and arbitrary in-plane deformation, which constitutes the basis of a multiscale model of graphene.

For modeling polycrystalline graphene, as grown by CVD for instance, in addition to the constitutive behavior of the bulk, the mechanical properties of the GBs need to be characterized and several studies have dealt with this issue, either by numerical simulations using molecular dynamics (MD) or by experiments using nanoindentation.

When one is interested in the failure mechanics of graphene GBs, two quantities are of interest. The cohesive strength, defined as the maximum stress that a GB can sustain and the fracture toughness, which is relevant to predict fracture propagation. Some authors have estimated the cohesive strength of GBs by MD simulations^25–27 looking for dependence with respect to the misorientation angle of the grain boundary. However, recent studies²⁶ suggest that the detailed arrangement of defects at the grain boundary and the orientation of the GB line are also determinant factors for the strength. Experimental nanoindentation studies^7,28 are inconclusive regarding the misorientation angle dependence of the strength of GBs.

The fracture toughness of graphene (single crystalline and polycrystalline) has been the object of several studies, both theoretical and experimental. Theoretical studies are based on MD and coupled quantum/molecular mechanics simulations. As summarized in the review paper of Zhang et al.,²⁹ various studies have predicted fracture toughness values ranging from 2 to 4 eV/Å.

On the experimental side, Zhang et al.³⁰ recently measured the mode I fracture toughness of polycrystalline graphene by performing pioneering tensile loading experiments in pre-cracked sheets of polycrystalline graphene. By using a microelectromechanical (MEMS) device, they were able to test several bilayer membranes of polycrystalline graphene in which cracks were initially introduced by focused ion beam (FIB) cutting with initial crack lengths ranging from tens of nanometers to 1 μm. Their results show that in the range of crack length studied, the Griffith criterion of fracture holds in the sense that the product of the critical stress with the square root of the crack length is constant. They measured a fracture toughness of 3.3 eV/Å. The typical size of the grains that constituted the polycrystalline graphene ranged from hundreds of nanometers to a few microns. Thus, the fracture toughness they measured corresponds to that within a graphene crystallite rather than the fracture toughness of the constituting GBs. In contrast, in our multiscale approach, we are interested in the specific properties of particular GBs.

The approach that we propose in this work, deriving a cohesive zone model of the GBs in graphene, encompasses these two notions of strength and fracture toughness. Moreover, it provides a description of the GB that may be incorporated in finite element models of polycrystalline graphene.

In the present section, we introduce the concepts and the methods that we use to build a model that characterizes the mechanical properties of GBs in graphene. This model is a cohesive zone model of the GBs, for which our ultimate goal is the implementation in a FEM which simulates the indentation experiments of polycrystalline graphene. In the general case, a CZM should be mixed-mode in order to be able to take into account fracture processes that are a combination of different fracture modes (the opening, sliding, and tearing modes), see, for instance, Park et al.³¹ In addition, the influence of stress triaxiality on the CZM may be important as discussed by Siegmund and Brocks.³² While the CZM cannot originally account for stress triaxiality, some authors (see, for instance, Remmers et al.³³) have proposed an extension of the CZM that incorporates that effect.

In the present study, our aim is to develop a single-mode CZM that is adapted to conditions typical of the indentation experiments of graphene films. In their FEM simulations of indentation, Wei and Kysar⁶ have shown that the highest stresses are concentrated under the indenter tip and that the corresponding stress state there is equibiaxial. Also, near the indenter tip, i.e., where a GB is likely to fail since stresses are high, the stress state is approximately equibiaxial tensile and the in-plane shear stress is very small as compared to the normal stresses. Therefore, the operative mode of fracture is the opening separation (mode I fracture) and we can neglect the mode II contributions. Thus, we aim herein at building a mode I CZM valid in a range of stress states close to equibiaxial. To quantify the notion of closeness to an equibiaxial stress state, we note that our simulations suggest that the cohesive zone parameters may change by 20% between a uniaxial tension and an equibiaxial tension. Therefore, when considering biaxial stress states for which the ratio of the principal stresses is less than two to one, we can expect that the CZM will be valid within 10%.

In the first part of this section, we introduce the concepts that allow to build the CZM. We then develop the principles of the numerical tool that we use to derive our CZM, i.e., the molecular dynamics, and review the state of the art in the MD to CZM scale bridging. Finally, we describe the simulation that we use for deriving our CZM as well as the averaging procedure to transition from the atomistic scale to the continuum.

The concept of cohesive zone was first introduced by Barenblatt⁹ for brittle fracture and Dugdale⁸ for ductile fracture. We focus on the Barenblatt approach, well suited for the brittle fracture of graphene. It consists in considering a region near the tip of the crack where the two opposites sides are subjected to cohesive tractions as depicted in Fig. 5(a). The cohesive traction, for small crack openings, takes its origin from the interactions at the atomic scale, and allows to overcome the problem, present in the linear elastic fracture mechanics approach, of stress singularity at the tip of the crack. The region over which cohesive tractions are exerted is called the cohesive zone and its length, the cohesive zone length is often assumed to be small with respect to the crack length. Moreover, the distribution of tractions in the cohesive zone is described by the so-called traction-separation law (TSL); see Fig. 5(b) for a typical TSL. The TSL relates the local opening displacement to the cohesive traction along the cohesive zone and characterizes the interface. The TSL is expressed as a function of the separation that increases up to a maximum traction t_m with corresponding separation δ_m beyond which failure occurs irreversibly and the force decreases to zero, which corresponds to the critical separation δ_c for which the crack is fully opened. The TSL is characterized by its shape and the values of t_m and δ_c. An important quantity is the area under the traction–separation curve which is called the work of separation, w_sep, and corresponds to the energy absorbed by the fracture process per unit length of crack growth. Depending on the shape of the TSL, two or more parameters may be appropriate for its identification. If one chooses a bilinear law, as we shall explain later, the two main quantities to extract from these experimental curves are t_m and δ_c.

A common way to implement a CZM in an FEM model of fracture consists in deriving the parameters of the TSL from experiments. In the context of graphene GBs, the nanometric scale at which fracture is studied suggests to extract the TSL's parameters by performing MD simulations.

MD simulations consist in solving Newton's equations of motion at the level of the atoms that compose the material. It is a phenomenological method of material modeling where the interactions between the atoms are seen from a classical mechanics perspective. In MD, time is typically discretized at the scale of the femtosecond (fs) or 10^–15 s, which corresponds to the smallest time scale that we need to resolve, i.e., atomic vibrations.

The concepts of statistical physics allow then to relate the trajectories of the atoms and the interatomic forces to macroscopic quantities such as the strain, the stress, the potential, and kinetic energies. The primary outputs of MD are, at each time step, the position and the velocity of each atom in the system. From these quantities one may derive, given the potential of interaction between the atoms, the force exerted on each atom; the potential energy of the system U_pot that results from the sum of the potential energy of all the atomic interactions; the kinetic energy of the system U_kin which is the sum of the kinetic energy over all the atoms and which is directly related to the temperature. For less straightforward macroscopic quantities such as the stress, there exist different methods to derive them from the primary MD outputs. In Sec. III E, we discuss the Virial definition of stress. For an exhaustive discussion on the different atomic-based definitions of stress, see for instance Refs. 34–36.

The atomic interactions are accounted for by a phenomenological interatomic potential that embeds all the physics of the material. Therefore, the accuracy of the results derived from MD simulations is completely dependent on the validity of the chosen interatomic potential.

The most widely used interatomic potentials for MD simulations of graphene are the second generation reactive empirical bond order (REBO) potential³⁷ and its variant, the adaptive intermolecular REBO (AIREBO).³⁸ These potentials are designed for carbon and hydrocarbon molecules, and are based on the Abell–Tersoff bond order formalism. REBO and AIREBO are pair potentials where only nearest neighbors interactions are taken into account while many-body effects are introduced through a bond order function that parametrizes the properties of each bond with regards to its environment. The bond order term accounts for the modification in atomic hybridization when bonds are breaking and forming. However, it has been noted that the process of fracture in graphene, as modeled by REBO or AIREBO, is significantly different from the predictions of the density functional theory.^39,40 This problem has been addressed by changing the cut-off functions that turn off the atomic interactions beyond the nearest neighbors.^25–27,30

Alternately, other groups have developed improved versions of the REBO potential, especially designed for the study of bond-breaking processes in carbon-based materials.^40,41 As noted by Perriot et al.,⁴¹ the use of REBO or AIREBO potential, as well as the modified cut-off function versions, has not been properly validated and fails to describe the bond breaking phenomena. The problem, encountered in large deformations, in the REBO potential comes from the contradiction between the need to increase the cut-off distance to fully describe the nearest-neighbor interactions and to exclude the second-nearest-neighbor interactions originally not accounted for in the REBO potential.

The improvement proposed by Pastewka et al.⁴⁰ and Perriot et al.⁴¹ to the REBO potential consists in incorporating an environment-dependent screening function that allows to increase the cut-off distance while only keeping the nearest-neighbor interactions. When the local environment is such that second- and further-nearest-neighbor interactions are detected, a function screens these interactions. These screened environment-dependent potentials have been compared to the predictions of DFT computations, which are the benchmark results, in the range of large deformations and for the processes of fracture, and have proven to yield very good predictions with regards to energy and force profiles as a function of the interatomic distance. The development and validation of these potentials are of fundamental interest for the MD simulations of fracture processes in graphene.

In the present work, we have used the screened environment-dependent potential of Pastewka et al.,⁴⁰ called REBO2+S, as implemented in the atomistica package of the MD software LAMMPS.

Gall et al.⁴² and Spearot et al.⁴³ were pioneering in using MD to investigate the cohesive interface constitutive relations for a bi-material interface and a GB interface, respectively. They simulated tensile tests perpendicular to the interface on small samples of typical size 4 to 8 nm. By defining a measure of the opening and the tension in the sample, they derived interface separation relations. However, as noticed by Yamakov et al.,¹⁴ the small size of the sample and the boundary conditions are such that these authors tend to model interface adhesion rather than crack propagation.

Instead of interface adhesion, the CZM has originally been developed to describe the system of cohesive forces in the terminal region of a crack. Yamakov et al.¹⁴ extended the previous work by switching from MD simulations of homogeneous interface decohesion to simulations of crack propagation. The simulated cell, of typical size 100 nm, represents a bicrystal with a central crack along the GB. This bicrystal is loaded in order to trigger crack propagation along the GB. Yamakov et al.¹⁴ proposed a methodology that takes local measures of the stress and the opening displacement along the propagating crack to compute the TSL. The main elements of this methodology have been used again by Zhou et al.¹⁵ and Krull and Yuan¹⁶ who have investigated the TSL of interfaces in other contexts.

Herein, we adapt the methodology mentioned above to the case of mode I intergranular fracture in graphene and complement it by a quantitative analysis of the TSL for the purpose of incorporation in a FEM. For our simulations, we devise boundary conditions that allow an energetic verification of the derived TSL, and we discuss the mesh size effect for the implementation of the TSL in FEM simulations.

Following Yamakov et al.,¹⁴ we derive the TSL characterizing GBs in graphene by defining local measures of stress and strain in a MD simulation of crack propagation. We use the MD software LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator)⁴⁴ with the REBO2+S potential.⁴⁰ Unlike previous studies,^14–16 we choose the sample, the boundary conditions, and the loading conditions in order to control the rate of crack propagation. The above cited studies are based on simulations where, because of an elastodynamic instability, the crack propagates dynamically. In contrast, the controlled propagation allows us to perform an energetic verification of the TSL.

We simulate a finite size specimen that has the shape of a double cantilever beam (30 nm in width, 120 nm in length, one atomic layer of thickness) made of a bicrystal of graphene as seen in Fig. 6. The specimen contains, depending on the GB that it models, around 140 000 atoms of carbon. The motion of atoms is constrained in the (x, y) plane through reflective walls located on each side of the sheet so that, while the local motion of the atoms remains 3D, the graphene sheet cannot ripple or fold out of plane. This constraint does not significantly affect either the qualitative behavior of the fracture or the quantitative values of the observables. The edges at y = ±h are free. On the right end of the beam, atoms within a 0.5 nm wide strip are fixed. On the left end of the beam, and independently for the two grains, the atoms belonging to a 0.5 nm wide strip are bound together in order to form two rigid bodies. These rigid bodies constitute useful subsets of the structure where one can apply kinematic boundary conditions. Indeed, the displacement along the x direction and the rotation around the z axis of these two rigid bodies are free while the displacement of their center of mass along the y direction is prescribed by applying uniform equal and opposite velocities of magnitude v = 4 nm/ns. Thus, vertical prescribed displacements are $q(t)=±vt$⁠. With such boundary conditions, a controlled fracture process is simulated. The resulting local stress state at the tip of the crack is quasi-equibiaxial, thus reproducing the stress state under and near the indenter in the indentation experiment that we aim at modeling.

The time step is set to 1 fs, and the total simulated duration is 2.5 ns. To focus on the mechanics of the decohesion and avoid thermal activation effects, the simulation is performed on the NVT ensemble at the prescribed temperature of 0.1 K through a Nose–Hoover thermostat. Before the loading, the atoms in the bicrystal beam are properly relaxed to cancel out any initial internal stress which may have been introduced by the design of the bicrystal and its GB, which fixes exactly the position of all the atoms. When the prescribed velocity is applied, the mode I crack propagates in a controlled way along the GB.

The crack tip velocity is on the order of 25 m/s, which may be compared to the velocity under uncontrolled steady-state dynamic propagation observed in MD simulations, at 1500 m/s, which is about one seventh of the Rayleigh wave speed.⁴⁵ The difference between these values is about two orders of magnitude, which justifies the assumption that the controlled crack propagation is quasi-static.

We define N cohesive zone volume elements (CZVE) of dimensions L_x by L_y along the GB. These CZVE are 2D-volumes that are called volumes to highlight the fact that they are used to consider atoms that are within it. More precisely, atoms inside these areas in the reference configuration are definitively assigned to the CZVE (Fig. 7(b)). The stresses in these regions are computed from the 2D-adapted definition of the atomic stress, based on the Virial theorem³⁶ and defined by

where i and j are the Cartesian coordinates, $mα$ denotes the mass of atom α, and $vα$ is its velocity relative to the macroscopic motion. $Ωα$ is the surface area of atom α in the present configuration, which may be approximated, when possible, to the one of an atom in the honeycomb lattice in the reference configuration: $Ω0hon=334d2$ with d as the interatomic distance. $rαβ=rα−rβ$ is the displacement vector between atoms α and β, and $fαβ$ is the interatomic force exerted on atom α by atom β. The sum is performed over the set of atoms $N(α)$⁠, neighbors of atom α that are in its interaction range. The 2D character of the above defined atomic stress lies in the division, in Eq. (1), by the atomic surface area instead of the atomic volume as it is done in the usual 3D definition of the atomic stress.

The 2D-stress in the kth CZVE is then deduced through a time and volume average of the atomic stress over the CZVE

where $Nτ=2500$ is the number of time steps t in the interval τ = 2.5 ps over which the average is performed and A^k is the set of atoms belonging to the kth CZVE. $Ωk$ denotes the surface area of the kth CZVE and is computed as the surface area $Ω0hon$ times the number of atoms in the CZVE. This measure of $Ωk$ neglects the inhomogeneities in atomic surfaces introduced by the pentagonal and heptagonal defects (mesoscopically compensated) as well as the small local strain around the crack tip.

Note that the atomic stress as defined in Eq. (1) is of practical interest and provides qualitative information regarding the stress field at the atomic level but shall not be considered to be a macroscopic stress measure.³⁶ In contrast, space and time average of the atomic stress is related to the continuum notion of Cauchy stress to which it converges in the thermodynamic limit.^34–36 The tension perpendicular to the GB, which is of interest for the mode I TSL, is defined as the yy component of the 2D-stress measure: $σyyk$⁠.

Since the crack is atomically sharp, the opening displacement in the kth CZVE is defined as the average distance between the atom lines forming the crack edges in this particular CZVE (the opening is defined as zero in the relaxed configuration). Like for the stress, the opening displacement in each CZVE is averaged over 2500 time steps.

The selection of the size of the CZVE over which the stress is averaged, and our choice of measure of the opening displacement at the crack edges is developed and justified in Sec. IV C.

We focus our study on two high angle grain boundaries. The first one denoted $(3,1)|(3,1)$ is a symmetric GB of tilt angle 27.8° and the second one denoted $(7,0)|(4,4)$ is an asymmetric GB of tilt angle 30.0° (Figs. 2(b) and 2(c)). Fracture simulations of low angle GBs show a competition between crack propagation along the GB and deviation of the crack inside one of the grains. This may be due to the structure of low angle GBs, where the atomic defects (pentagons and heptagons) do not form an almost continuous line but are separated by domains of undisturbed lattice (hexagons).

As a reference, we also determine the TSL for a crack that propagates in the bulk of pristine graphene along the $[11¯00]$ direction (or armchair), i.e., the direction y shown in Fig. 2(a).

For each one of the 55 CZVEs dispersed along the GB, we superpose in Fig. 8 all the traction-separation points measured from the beginning of the simulation to the time 1.25 ns. In fact, each CZVE goes through the steps where the local stress progressively increases until the point where the crack goes through the CZVE and the stress drops. These different steps appear in Fig. 8 for all the 55 CZVE, and it can be seen that the traction-separation relation propagates in a self-similar way as the crack propagates through the sample.

For the purpose of implementation in FEM, we need to derive from the set of tractions and separations obtained in MD an analytical expression of the TSL that includes physically significant parameters, namely, the peak traction t_m and the work of separation w_sep. Fig. 8 suggests the use of a bilinear TSL such as the one proposed by Geubelle and Baylor¹³ that is entirely defined by three parameters: the peak traction t_m, the corresponding separation δ_m for the ascending part, and the critical separation δ_c at which the traction vanishes for the decreasing part. The work of separation is then expressed by

One notices that the value of δ_m does not appear in the expression of w_sep, given the approximate triangular shape chosen for representing the TSL. Fig. 8 shows the bilinear TSL fitted to the simulation points for the GB $(3,1)|(3,1)$⁠, and Table I reports the parameters that characterize the bilinear TSL for armchair bulk graphene, the GB $(3,1)|(3,1)$⁠, and the GB $(7,0)|(4,4)$⁠. The comparison of the parameters of the three derived TSL shows that the main characteristics of the TSL differ by less than 10%. In particular, the work of separation that corresponds to the fracture toughness of the GB shows little difference with pristine graphene. Quantitatively, the fracture toughness of the two studied GBs is less than 7% lower than the one of pristine graphene with a $[11¯00]$ crack.

The CZM approach has to be consistent with the energetic approach of crack propagation developed by Griffith⁴⁷ for virgin homogeneous solids and generalized by Irwin and Orowan⁴⁸ to fracture occurring between unlike bodies or at singular interfaces.

This consistency lies in the equality between the work of separation w_sep of the CZM and the fracture toughness G_c of the material or the GB.^10,49 For a demonstration of the equivalence between w_sep and G_c based on the J-integral, see Lawn.⁵⁰ On the one hand, the work of separation corresponds to the mechanical energy absorbed by the cohesive forces per unit length of fully opened crack (δ ≥ δ_c). It is a characteristic of the TSL that describes the cohesive traction. On the other hand, the fracture toughness, also called crack-resistance energy, expresses the ability of a material or an interface to resist to crack propagation. It gives a criterion for an equilibrium crack propagation called as the Griffith energy balance concept⁵⁰ 

where G is the mechanical energy release rate defined by

with Π as the mechanical potential energy and a as the crack length.

We are thus able to investigate the internal consistency of the TSL in terms of energy by comparing the values of w_sep and G_c.

First, the calculation of w_sep from TSL parameters is straightforward by applying Eq. (3). We can estimate that there is an uncertainty of about 10% on the value of w_sep that comes from an uncertainty of ≈5% on both t_m and δ_c.

Second, to derive G_c, we use the equilibrium crack propagation condition Eq. (4). Indeed, in our simulations, the crack propagation is controlled in the sense that to a continuously varying prescribed displacement corresponds a quasi-continuous variation of the crack length. As may be seen in Fig. 9, the crack moves by successive jumps that correspond to the periodic pattern of the GB and this characterizes bond trapping, a phenomenon that we discuss later. Notwithstanding this atomic scale specificity, from the continuum point of view, we observe on average controlled crack propagation unlike the typical catastrophic failure observed in uniform tensile loading experiments.

With the displacement boundary conditions applied in our simulation, the mechanical potential energy Π involved in the definition of G Eq. (5) is equal to the elastic energy U_el.

U_el appears as a component of the internal energy U_int that is a primary output of MD defined as the sum of the potential U_pot and kinetic U_kin energies introduced in Sec. III B. Note that in our simulations, because of the thermostat at 0.1 K and the low velocity of the prescribed displacement, the kinetic energy which is three orders of magnitude less than the potential energy is negligible, hence, U_int = U_pot.

To compute U_el, we note that

where $γeg$ is a generalized measure of the edge energy per unit of length, that is

• For a crack in the bulk, $2γeg=2γe$ with γ_e as the edge energy per unit of length of the newly created edge.

• For a crack along a GB, $2γeg=2γe−γGB$ with γ_e as the edge energy per unit of length of the newly created edge and γ_GB as the energy recovered from the destruction of the GB. Note, however, that this decomposition is in the case of graphene rather artificial. Indeed, some pentagon and heptagon defects may remain, after fracture, on either side of the crack, making thus the decomposition of $2γeg$ in surface energy created and grain boundary energy recovered meaningless.

In both cases, $2γeg$ is defined and obtained from the difference in potential energy U_pot between the broken and non-broken relaxed configurations with periodic boundary conditions (PBCs) shown in Fig. 10 divided by the crack length. Indeed, from Eq. (6), the difference in internal energy (that is potential energy in the absence of kinetic energy) between the two configurations of Fig. 10 corresponds to $2γega$⁠.

We then compute the evolution of U_el—and thus Π—by inverting Eq. (6).

We describe the state of the system at time t with two variables: q(t) is the prescribed displacement (see Sec. III D) and a(t) is the crack length; hence, $Π(q,a)=Uel(q,a)$ is a function of these two variables that are linked, at equilibrium crack propagation, by Eq. (4).

We use a simple Bernoulli–Euler beam model of the DCB simulation with an isotropic linear elastic constitutive law to relate G to the total derivative of $Π(q,a)$ with respect to a, taking into account the relation q(a) satisfied during crack propagation. The linear isotropic model is reasonable given that strains in the DCB are, except near the crack tip, less than 0.5% (see Fig. 4). Denoting E as the 2D-Young's modulus with dimension force per length and h as the width of DCB's arms (Fig. 6), the mechanical potential energy can be written⁵¹ 

Writing the equilibrium crack propagation condition Eq. (4) with G computed with Eqs. (5) and (7), we obtain the relation between a and q

That allows us to write Π as a function of a only

Hence, G_c may be deduced from the total derivative of $Π̃(a)$ that is computed from the molecular dynamics outputs

In Fig. 11, we plot the function $Π̃(a)$ as computed from MD. We can see that notwithstanding the saw-tooth pattern in $Π̃(a)$ related to atomic scale phenomena, the overall evolution of the function that is to be considered in the continuum approach is linear in accordance with Eq. (9), except in the crack initiation region where the beam model is not appropriate. The slope of the linear fit to the function $Π̃(a)$ and Eq. (10) yields G_c. For the GB $(3,1)|(3,1)$⁠, we obtain G_c = 2.7 eV/Å.

Table II presents the values of the work of separation for the modeled TSLs, the fracture toughness, and the edge energy for the three studied configurations. Their comparison shows that the work of separation of the derived TSLs falls within 10% of the fracture toughness for the configurations studied. Therefore, the TSLs proposed for FEM implementation are energetically consistent with the energy based approach of Griffith for crack propagation. The values of fracture toughness obtained for the grain boundaries and the bulk graphene are consistent with the empirical measures of the fracture toughness of polycrystalline graphene performed by Zhang et al.³⁰ They have measured the critical fracture stress of nanocracked polycrystalline graphene sheets and obtained an empirical value of the fracture toughness of 15.9 J/m² or 3.3 eV/Å. This measure which represents a microscopic homogenized measure of the fracture toughness is expected to lie in the range of the nanoscopic fracture toughnesses of bulk graphene and its grain boundaries. Indeed, our computation of nanoscopic fracture toughnesses reported in Table II is of the same order as the microscopic fracture toughness measured by Zhang et al.³⁰ 

Table II also reports the values of $2γeg$ that correspond to $2γe$ for the bulk and $2γe−γGB$ for the GBs. In all the cases, the fracture toughness is about 35% higher than $2γeg$ which is consistent with the fact that $2γeg$ is a lower bound of the energy per unit length required to propagate a crack. The difference between $2γeg$ and fracture toughness is likely to be attributed to the phenomenon of bond trapping.^52,53 Whereas it is often assumed that the fracture toughness corresponds to the edge energy created decreased by the grain boundary energy, bond trapping at the atomic scale may require greater amount of energy to propagate the crack. The concept of bond trapping comes from the fact that fracture propagation is ultimately determined by the strength of interatomic bonds. The crack may be caught in a metastable state where an energy barrier needs to be crossed to allow the propagation.^54,55 The extra energy is supposed to be transported in the lattice vibrations or, for the thermostated case, in the heat bath itself.

This thermal dissipation appears clearly when we check the energy balance of the global system. Thermodynamic quantities that characterize energy transfers are the internal energy U_int, the work W received by the system from external forces, and the heat Q received by the system from the thermostat. Fig. 12 shows that during the delamination process, the work provided by the external reaction forces is larger than the change in internal energy, and concurrently, heat is released by the system to the thermostat (Q < 0). Since we check that $Uint−W−Q$ is constant with time (see Fig. 12, numerically, the deviation in $Uint−W−Q$ from zero remains less than 3% of W), the first law of thermodynamics $ΔUint=W+Q$ is satisfied, which is to say, in MD terms, that the energy is conserved. The difference between U_int and W indicates that not all the work provided by the external forces to the system is transmitted to the internal energy. This thermal dissipation is consistent with the bond trapping that we have previously identified.

Because of the high gradients of stress at the tip of the crack, we expect the traction-separation response to depend on the size of the CZVE. An interesting discussion on the dependence of the traction-separation relation on the exact definition of the opening displacement and the traction at the atomic scale can be found in Gall et al.⁴² We explain here that the choice of the CZVE size in which stress and displacement are computed is justified by quantitative analysis of the TSL. Generally speaking, the CZVE should be small enough to resolve the stress gradient near the crack tip but large enough to ensure that the atomic scale variations of the atomic stress (i.e., the value of the Virial stress for individual atoms) are smoothened. The last condition ensures that the Virial stress averaged over the CZVE converges sufficiently close to the Cauchy stress. Nevertheless, since the TSL is itself an interpolation of the stress-separation measures, a perfect convergence of the CZVE average to the continuum definition of the Cauchy stress is not necessary.

First, the width L_y in Fig. 7 of the CZVE should be large enough to encompass the pre-stress field, due to the pentagon-heptagon defects, whose size characterizes the “width of the GB.” The lower bound for L_y thus lies about 8 Å. The choice of L_y should also resolve the stress gradient in y-direction that spans over about 30 Å (30 Å is the length over which the atomic stress is reduced by four). Within this range, we choose L_y = 10 Å, noting that any value between 8 Å and 20 Å provides TSL's parameters that differ by less than 5%.

The second parameter is the length L_x of the CZVE. Yamakov et al.¹⁴ noticed that the length of the CZVE gives a length scale inherent to the computed TSL. This length scale will be used in the FEM as the interface element length. Thus, considering the mesh size of the FEM in which the TSL will be implemented provides an enlightening insight for the choice of L_x.

There exists an intrinsic length scale along the crack path when implementing the CZM which is the cohesive zone length. The cohesive zone length is a characteristic length of the process zone over which traction is exerted. It is defined by Hillerborg et al.¹⁰ as

where E is the Young's modulus of the material.

For the FEM implementation of the CZM to converge, the FE mesh along the cohesive surface elements must resolve the stress gradient in the cohesive zone; practically, the cohesive zone length should contain at least two to three elements.⁵⁶ In the case of graphene, taking $E=340 N/m$ (Ref. 5) and $Gc=3.0 eV/Å, tm=25 N/m$ as average values for the studied GBs, the cohesive zone length is: $lcz=2.6 nm$⁠. Therefore, in order to have three elements or more within the cohesive zone, the length of the cohesive zone elements should be less than or equal to l_e = 0.9 nm. This first analysis provides an upper bound for the length L_x of the CZVE. On the other hand, L_x should be large enough to give a good estimate of the Cauchy stress and, if possible, should be a multiple of the repeating length of the pattern of defects in order to smooth out the atomic scale stress variations.

Giving these constraints, the length of the CZVE is taken close to 9 Å, the exact values being adjusted with the repeating length of the GB that are 8.9 Å for the GB $(3,1)|(3,1)$ and 17.0 Å for GB $(7,0)|(4,4)$⁠. We therefore choose for the CZVE length: L_x = 8.9; 8.5; 9 Å, for the GB $(3,1)|(3,1)$⁠, the GB $(7,0)|(4,4)$⁠, and bulk graphene respectively.

The modification of L_x on the TSL is such that an increase in L_x reduces t_m while keeping the work of separation w_sep constant. Therefore, the comparison of the maximum traction t_m of the TSL with the expected value provides another way to set the value of L_x. Indeed, t_m corresponds physically to the maximum stress sustainable by the grain boundary and gives the criterion for crack initiation. By performing equibiaxial tensile simulations on flawless bicrystals (that correspond to the crack tip stress state for which the TSL is derived), we have determined the expected value of t_m. We have performed simulations on bicrystal samples of dimensions 60 nm × 30 nm containing two GBs to enforce periodic boundary conditions in the two in-plane directions. See Fig. 13 for a schematic view of the performed tensile simulations. We have obtained, under equibiaxial stress state, the strength of the grain boundary and the bulk by computing the maximum normal component of the stress in the x-direction, reached before failure. For the three configurations, GB $(3,1)|(3,1)$⁠, GB $(7,0)|(4,4)$⁠, and bulk graphene, we compute a strength under equibiaxial stress state of (24 ± 0.5) N/m. The maximum stresses t_m of the TSLs given in Table I correspond to this reference value of strength, further confirming the consistency of the choice of 9 Å for L_x.

Lastly, we discuss the definition of the opening displacement. Yamakov et al.¹⁴ defines the separation as the distance between the CZVE averaged over all the atoms in the CZVE while Krull and Yuan¹⁶ only take the pair of closer atoms on both sides of the crack to define the separation. It can be understood that the location where the opening displacement is measured only affects the value of the separation δ_m at which the maximum stress is reached. This is because the opening displacement is, regardless of where it is measured, defined relative to the uncracked state for which δ = 0; whether it is defined between the crack edges or between areas encompassing the edges does not change the value of δ when the crystals are relaxed. Thus, the critical displacement δ_c at which the traction falls to zero is independent of the exact definition of the separation and it is only δ_m, and consequently, the initial slope of the TSL that is affected by that choice. The dependence of the TSL on the definition of the separation, as described above, has been numerically observed by Zhou et al.¹⁵ Given that the modification of the definition of the separation amounts to changing the initial slope of the cohesive zone while keeping the other parameters t_m and δ_c unchanged, we examine the meaning of the initial slope with respect to the FEM implementation of the TSL. The slope of the increasing part of the TSL is referred as the interface stiffness; however, it does not really correspond to a physical concept since it depends on the size of the volume that is considered to be the interface. Instead, the value of the interface stiffness is rather based on numerical considerations in the FEM implementation. The interface stiffness should be large enough to make sure that the contribution of the cohesive elements to the global compliance remains small but not too high to avoid numerical problems such as stress oscillations.^56,57 Given the fact that it is mainly numerical aspects related to FEM that govern the choice of the interface stiffness used in the implementation of the CZM, we provide herein a TSL based on the physical definition of the separation for an atomically sharp crack. This definition corresponds to a measure of the separation exactly at the atoms forming the crack edges. A measure further away from the crack would include a contribution of the bulk to the interface description. For the FEM implementation, the interface stiffness would have to be properly adapted to the numerical context.

As discussed in Sec. IV C, for a successful implementation of the CZM in the FEM, the cohesive zone length should contain at least three elements, bringing the element size down to l_e = 0.9 nm. To reduce the computational cost of such a FEM simulation, it is of practical interest to be able to implement the CZM with coarser mesh sizes that exceed the cohesive zone length. Such a method has been proposed by Turon et al.⁵⁶ who showed that delaminations were well predicted by FEM simulations with a mesh size up to three times the cohesive zone length provided that the TSL was properly scaled. In short, the scaling of the TSL consists, when the length l_e of a cohesive interface element does not resolve the cohesive zone length any longer, in reducing the peak traction as $1/le$ while keeping the fracture toughness constant.

More precisely, the scaling procedure proposed by Turon et al.⁵⁶ to accurately model a delamination process with CZM implemented in a FEM with elements larger than the cohesive zone length l_cz is based on the following idea: Although it is the peak traction of the TSL t_m that governs the crack initiation, the crack propagation is essentially controlled by the work of separation w_sep. Further, it was observed numerically that changing the peak traction does not significantly alter the FEM predictions of a delamination process but that lowering that peak traction improves the convergence of the results with respect to the element size l_e. To be understood, these observations may be put in perspective with the need to resolve the cohesive zone length with the element size. The cohesive zone length is related to the size of the elements that discretize it through the formula

Then, by combining Equations (11) and (12), one gets the relation between the peak traction and the size of one element

Eq. (13) is a way to express the finite element size l_e given the peak traction t_m of the TSL and by fixing the number of elements N_e in the cohesive zone length to a value large enough in order to discretize the stress profile accurately enough.

The strategy proposed consists in scaling the TSL in such a way that, first, N_e is kept to the chosen value—N_e = 3 is a minimum for a resolution of the cohesive zone—and second, Eq. (13) is satisfied. Practically, one may fix l_e to an arbitrary value $le¯$ larger than l_cz/N_e and scale the peak traction of the TSL to the artificial value $tm¯$ that results from Eq. (13) in which $le¯$ has been substituted to l_e. This amounts to defining an artificial larger cohesive zone length $lcz¯=Nele¯$ that is properly discretized. Moreover, the scaling of the TSL should be done at a constant work of separation in order to conserve a good mechanical energy balance. The critical opening displacement δ_c will be increased accordingly.

Turon et al.⁵⁶ have shown numerically that this procedure allows to accurately predict the delamination process for an element size $le¯$ set up to three times the cohesive zone length l_cz. The accuracy of the results can be seen in the estimate of the macroscopic load versus displacement. However, with such a scaled TSL, the stress concentration at the crack tip will not be well predicted.

We have observed a quite similar rescaling behavior when changing the CZVE length L_x in the MD extraction of the TSL. Indeed, Fig. 14 and the TSL parameters reported in Table III show that when the length of the CZVE exceeds 0.9 nm (i.e., one third of the cohesive zone length), that is to say, it does not resolve the cohesive zone length, the peak traction of the TSL is reduced as $1/Lx$ while conserving w_sep. Saether⁵⁸ has suggested that the length of the CZVE L_x should be the same as the cohesive surface element size l_e in the FEM mesh, and our numerical results put in perspective with the work of Turon et al.⁵⁶ tend to confirm this idea.

On the basis of the work of Turon et al.,⁵⁶ we expect that a FEM simulation with a cohesive surface element size of 5 nm embedding the corresponding scaled-TSL will properly model the crack propagation for a lower computational cost, although the stress concentration at the crack tip will not be resolved. Table IV gives the scaled-TSL parameters for the two GBs and the bulk corresponding to an element size l_e = 5 nm. However, this scaling procedure has been validated for the delamination process only and it may be, since the crack initiation is governed by the value of t_m, that the prediction of the crack initiation will not be accurate when using the elements that are larger than 0.9 nm. The proposed scaled-TSL may be a way to model a phenomenon at a reasonable computational cost, but will miss other features for which a finer mesh will remain necessary.

For the purpose of developing a multiscale model of nanoindentation experiments performed on polycrystalline graphene, we have developed cohesive zone models of intergranular fracture that may be implemented in FEM. We have derived, from molecular dynamics, traction-separation laws that characterize high tilt angle grain boundaries. Grain boundaries of tilt angle 27.8°, 30.0°, and pristine graphene have been investigated. We have performed MD simulations that reproduce the biaxial stress state to which the graphene is subjected during the indentation experiments in order to derive TSLs that correspond to the framework in which they will be implemented.

The MD simulations are designed to ensure a controlled crack propagation that allows an energetic validation of the TSL based on a comparison of the work of separation of the TSL and the fracture toughness of the sample. The energy absorbed by the cohesive tractions corresponds to the fracture toughness measured through the energy release rate during the equilibrium crack propagation. Our computation agrees with the experimental measures of the fracture toughness of polycrystalline graphene performed by Zhang et al.³⁰ at the microscale. However, the fracture toughness turns out to be greater by about 35% than the energy thermodynamically necessary to create the two new edges: $2γe−γGB$⁠. We attribute this difference to the phenomenon of bond trapping that specifically appears in atomic scale fracture.

The consistency of the TSL is discussed in relation to the mesh size effects introduced by the averaging procedure over the cohesive zone volume elements. When the scale over which atomic quantities are averaged is too large to resolve the traction distribution in the cohesive zone, in our case larger than 0.9 nm, the parameters of the TSL depend on the length of the volume elements. This dependence is such that the maximum traction t_m of the TSL is lower for coarser meshes while conserving the work of separation. When considering the FEM implementation of the CZM, the optimization of the computational cost encourages us to provide a TSL that could be implemented with coarser meshes. Therefore, we propose particular scaled-TSLs adapted to finite element surface elements ranging from 0.9 to 5 nm.

We have focused our study on two high angle grain boundaries and one type of intragranular crack under conditions of biaxial traction at the crack tip. The derived TSLs show that the TSL parameters, and in particular, the fracture toughness of the high angle grain boundaries, remain comparable to the fracture toughness of the bulk. This provides a new point of view on the impact of the grain boundaries on the mechanical properties of graphene. Further work will be needed to extend the cohesive zone model to low angle grain boundaries. Moreover, we have noticed a dependence of the TSL on the stress state applied at the tip of the crack, the TSLs are currently derived for quasi-equibiaxial stress states that correspond to what a graphene sheet experiences under the indentation tip. Supplementary investigations would be necessary to extend the TSL to other stress states.

We thank J. J. Marigo for fruitful discussions and L. Pastewka for his support regarding the REBO2+S interatomic potential. J.W.K. gratefully acknowledges support from the National Science Foundation Grant Nos. CMMI-1363093 and CMMI-1437450.