Characterizing the fracture parameters of a graphene sheet using atomistic
simulation and continuum mechanics
Jia-Lin Tsai
*, Shi-Hua Tzeng, Yu-Jen Tzou
Department of Mechanical Engineering, National Chiao Tung University, Hsinchu 300, Taiwan
a r t i c l e
i n f o
Article history:Received 22 August 2009
Received in revised form 11 October 2009 Available online 22 October 2009 Keywords:
Graphene sheet Fracture
Strain energy release rate Molecular dynamics simulation
a b s t r a c t
The fracture behavior of a graphene sheet, containing a center crack (length of 2a) was characterized based on the atomistic simulation and the concept of continuum mechanics. Two failure modes, i.e., opening mode (Mode I) and sliding mode (Mode II), were considered by applying remote tensile and shear loading, respectively, on the graphene sheet. In the atomistic simulation, the equilibrium configu-rations of the cracked graphene subjected to applied loadings, before and after the crack extension of 2Da, were determined through molecular dynamics (MD) simulation, from which the variation of the potential energy and the strain energy release rate of the discrete graphene sheet because of crack extension was calculated accordingly. It is noted that because of the discrete attribute, there is no stress singularity near the crack tip, and thus, the concept of stress intensity factor that is generally employed in the continuum mechanics may not be suitable for modeling the crack behavior in the atomistic structures. For the sake of comparison, the continuum finite element model with the same geometric parameters and material properties as the atomistic graphene sheet was constructed, and the corresponding strain energy release rate was calculated from the crack closure method. Results indicated that the strain energy release rates obtained from the continuum model exhibit good agreement with those derived from discrete atomistic model. Therefore, it is suggested that the strain energy release rate is an appropriate parameter, which can be employed in the atomistic model and the continuum model for describing the fracture of cova-lently bonded graphene sheet.
Ó 2009 Elsevier Ltd. All rights reserved.
1. Introduction
In the theory of linear elastic fracture mechanics (LEFM), stress intensity factor (SIF) and strain energy release rate (SERR) are the two parameters that are generally introduced in modeling the frac-ture behaviors of a continuum media. The stress intensity factor basically defines the amplitude of stress singularity at the crack tip, providing essential information, such as stresses, strain, and displacement fields, near the crack tip. It was postulated that the onset of crack extension is governed by the local singular stress field near the crack tip (so-called near tip stress field), and thus, the fracture behavior can be described appropriately using SIF. On the other hand, the strain energy release rate is regarded as the measurement of change in strain energy with an infinitesimal crack extension. It was assumed that when the SERR reaches a crit-ical value, the crack begins to propagate. Although the two quanti-ties were developed based on different physical concepts (one is from the local stress/displacement filed, and the other is from en-ergy variation), their applications to fracture mechanics are
essen-tially equivalent in a linear elastic continuum (Anderson, 1995). In other words, the two fracture parameter can be adopted alterna-tively in predicting the initiation of the crack embedded in a con-tinuum solid. Nevertheless, if the scale is deduced to a nano range, the material system examined is no longer a continuum so-lid, but a discrete molecular structure. Under this nano scale, it is interesting to explore if the two fracture parameters developed originally based on continuum concept can be extended and appli-cable to the discrete atomistic structure.
Jin and Yuan (2005b) investigated the macroscopic fracture parameters from both atomistic simulation and the continuum model, indicating that the near-tip stress field calculated from atomistic simulation agrees well with the continuum stress. A sim-ilar observation that the local stress near the crack tip can be well described by 1=pffiffiffir dependence was also addressed by Ome-ltchenko et al. (1997).Miller et al. (1998)examined the Mode I deformation near the crack tip in a single crystal nickel and then computed the critical stress intensity factor, indicating that the atomistic stress beyond the first three lattice spacings from the crack tip is consistent with the elastic 1=pffiffiffircurve. Because of the atomically blunt crack tip, the atomistic stress in the region near the crack tip deviated from the elastic solution (singular stress 0020-7683/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved.
doi:10.1016/j.ijsolstr.2009.10.017
* Corresponding author. Tel.: +886 3 5731608; fax: +886 3 5720634. E-mail address:jialin@mail.nctu.edu.tw(J.-L. Tsai).
Contents lists available atScienceDirect
International Journal of Solids and Structures
j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / i j s o l s t rfield). The mixed mode crack in a single crystal nano-plate was characterized byZhou et al. (2009)using quasi-continuum simula-tion technique. They addressed that the atomic stress fields neigh-boring the crack tip are also singular and controlled by the atomic stress intensity factor. If the near-tip stress singularity can be val-idated in the atomistic simulation, the concept of stress intensity factor may be employed to describe the fracture behavior of dis-crete atomistic structures. In addition to the stress intensity factor, the strain energy release rate for atomistic graphene sheet with center crack was calculated using J-integral approaching with the assumption of nonlinear constitutive relation (Jin and Yuan, 2005a). The critical value of J-integral Jcwas introduced to measure
the onset of crack extension of graphene sheet.
In this study, the atomistic structures as well as the continuum model of graphene sheet with central crack were established, respectively. In the discrete model, the local stress field near the crack tip was evaluated by using the hardy stress formula (Hardy, 1982; Zimmerman et al., 2004), and the suitability of the stress intensity factor employed to the atomistic structure was thor-oughly examined. In addition, the strain energy release rate in the atomistic model was determined using the global energy meth-od and the crack closure methmeth-od. For the sake of comparison, the continuum graphene model was established, and the correspond-ing near tip stress field and strain energy release rate were evalu-ated from the finite element analysis. The fracture parameters suitable for characterizing the atomistic structure and continuum model were of concern.
2. Atomistic simulation
2.1. Construction of atomistic structures of graphene sheet
Graphite structure is constructed by the carbon layers where the carbon atoms are arranged in a hexagonal pattern. The intera-tomistic distance between the adjacent carbon atoms is 1.42 Å, and the associated atomistic interaction is covalently bonded by SP2
hybridized electrons, the bond angle of which is 120° to each other (Cho et al., 2007). In order to investigate the mechanical properties of the graphite, the atomistic structures have to be constructed in conjunction with the appropriately specified atomistic interaction. In the description of graphite structure, two kinds of atomistic interactions are normally taken in account: one is bonded interac-tion, such as the covalent bond, and the other is the non-bonded interaction, i.e., van der Waals and electrostatic forces. Among the atomistic interactions, the covalent bond between two neigh-boring carbon atoms that provides the building block of the pri-mary structure of the graphite may play an essential role in the mechanical responses. Such bonded interaction can be described using the potential energy that consists of bond stretching, bond angle bending, torsion, and inversion (Rappe and Casewit, 1997). As a result, the total potential energy of the graphite contributed from the covalent bond can be written explicitly as
Ugraphite¼ X Urþ X Uhþ X U/þ X Ux ð1Þ
where Uris a bond stretching potential; Uhis a bond angle bending
potential; U/ is a dihedral angle torsional potential; and Uxis an inversion potential. For graphite structures under in-plane defor-mation, the atomistic interaction is mainly governed by the bond stretching and bond angle bending; therefore, the dihedral torsion and inversion potentials that are related to the out-of-plane defor-mation were disregarded in the modeling. The explicit form for the bond stretching and bond angle bending can be approximated in terms of elastic springs as (Li and Chou, 2003).
Ur¼ 1 2krðr roÞ 2 ð2Þ Uh¼ 1 2khðh hoÞ 2 ð3Þ where kr and kh are the bond stretching force constant and angle
bending force constant, respectively. The constants kr¼ 93; 800mole nmkcal 2and kh¼ 126mole radkcal 2selected from AMBER force field for carbon–carbon atomic-interaction (Cornell et al., 1995) was employed in our molecular simulation. The parameters r0and
h0represent bond length and bond angle in equilibrium position,
which are assumed to be 1.42 Å and 120°, respectively, for the graph-ite atomistic structures. In general, in addition to the bonded interac-tion, the non-bonded interacinterac-tion, i.e., van der Waals and electrostatic forces, has to be considered in the atomistic simulation. However, as compared to the covalent bond, the non-bonded interaction exhibits less contribution on the mechanical properties of graphene sheet; therefore, it was neglected in our current atomistic model.
The center crack of length 2a in the graphene sheet was created by eradicating the associated covalent bond so that for the atom pairs across the crack surface, there is no atomistic interaction tak-ing place. The dimension of the graphene sheet with center crack is illustrated inFig. 1. In order to construct the atomistic structures of cracked graphene sheet with a stress-free state, the modified NPT ensemble with the characteristics of varying simulation box in shape and size (Melchionna et al., 1993) was conducted in the MD simulation at the time increment of 1 fs for 100 ps.Figs. 2 and 3illustrate the potential energy and stress histories during the simulation, respectively. It can be seen that after the minimum potential energy is attained, the stress components are also approaching zero. Subsequently, the loading, i.e., uniaxial tensile stress and pure shear stress, was applied respectively on the boundary of graphene sheet to simulate the Mode I and Mode II fracture behaviors. After the modified NPT ensemble (Melchionna et al., 1993) with time increment at 1 fs for 100 ps was performed, the atomistic configuration of graphene sheet with the applied loading was determined. It should be mentioned that in the mod-ified NPT, the stress components could be applied independently on the simulation box; however, in the original NPT, only hydro-static stress can be implemented on the surfaces of the simulation box. The stress history in the modified NPT ensemble is presented inFig. 4where tensile stress (2 GPa) is applied on the graphene sheet to simulate the Mode I fracture. It is noted that in the study, the DL-POLY package originally developed by Daresbury Labora-tory (Smith and Forester, 2001) was modified to simulate the uni-axial loading/pure shear loading on the graphene sheet.
2.2. Local stress field near the crack tip
It is well-known that in the continuum media, the square root stress singularity occurs near the crack tip, and thus, the concept
Fig. 1. The dimension of the graphene sheet with center crack where L = 1022.4 Å, 2w = 1475.8 Å and 2a = 100.8 Å.
of stress intensity factor can be defined accordingly and utilized to characterize the fracture behavior. Based on the linear elastic frac-ture mechanics concept, we attempted to measure the local stress field in the discrete atomistic model and to understand if there is singular stress field existing in the vicinity of the crack tip. There have been various definitions of local stress in atomistic simulation such as virial stress, BDT (atomic) stress (Basinski et al., 1971), Lut-sko stress (Cormier et al., 2001), and Hardy stress (Hardy, 1982). Numerical results demonstrated that the Hardy is appropriate for both homogeneous and inhomogeneous deformation; moreover, it is very accurate and robust, and it is superior to BDT and Lutsko stress formulation in the local stress calculation (Zimmerman et al., 2004; Webb et al., 2008). Thus, the Hardy stress formulation ex-pressed as rðr; tÞ ¼ X N a¼1 mað
v
av
Þ ðv
av
Þwðra rÞ þ1 2 XN a¼1 XN b–a rab FabBab ðrÞ ( ) ð4Þwas employed to evaluate the local stress field near the crack tip of the graphene sheet, where ma;
v
aand rarepresent the mass, veloc-ity and position of atoma
, respectively;v
and w represent the con-tinuum velocity and localization function. In addition, rab¼ ra rband raand rbrepresent the position of atom
a
and b, respectively.Fabis the force between atoms
a
and b. By considering thelocaliza-tion funclocaliza-tion w as a Gaussian funclocaliza-tion
wðra rÞ ¼ 1 ðpffiffiffiffi
p
hÞ3exp ðra rÞ2 h2 " # ð5Þ bond function Bab ðrÞ is then expressed as Bab ðrÞ ¼ Z 1 0 1 ðpffiffiffiffip
hÞ3exp ðkrabþ ra rÞ2 h2 " # dk ð6Þwhere h is the smoothing length and is equal to 1.9 at the present simulation. It is noted that the kinetic term in Eq.(4)would equal zero when the analysis is performed at zero temperature. In addi-tion, the size of the localization volume used in the stress calcula-tion is 10 Å. During the employment of Hardy formulacalcula-tion, the angle bond interaction given in Eq.(3)needs to be converted into the atomic forces acting on the atoms, the detail procedure can be found elsewhere (Machida, 1999).
2.3. Calculation of strain energy release rate
In addition to the stress intensity factor, the strain energy re-lease rate is the other fracture parameter that has been frequently used in the linear elastic fracture mechanics. In order to apply the concept of the strain energy release rate in the discrete atomistic model, the energy difference before and after the crack extension has to be calculated appropriately. Two methods, global energy method and crack closure method, are introduced to determine the energy difference.
2.3.1. Global energy method
In the atomistic simulation, the total energy is contributed from the potential energy because of the bonded interaction and the ki-netic energy caused by the atomistic thermal oscillation. Because the simulation is conducted at 0 K, the effect of kinetic energy can be neglected, and the potential energy is regarded as the total energy of the discrete graphene system. When the atomistic con-figuration of the cracked graphene sheet under remote loading
Time (ps) T o ta l P o ten ti al E n er g y (k c a l/ m o l) 25 50 75 100 -2E-07 -1E-07 0 1E-07 2E-07
Fig. 2. Total potential energy histories during the NPT ensemble in MD simulation.
Time(ps) St re s s (G Pa ) 0 25 50 75 100 -1E-08 -8E-09 -6E-09 -4E-09 -2E-09 0 2E-09 4E-09 6E-09 8E-09 1E-08 σ xx σyy τxy
Fig. 3. Stress histories during the NPT ensemble in MD simulation.
Time (ps) Stress (GPa) 25 50 75 100 -0.5 -0.25 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 2.25 2.5 2.75 3 σxx σyy τxy
was obtained through MD simulation, the potential energy of the system with crack length equal to 2a was calculated based on the bonded potential function given in Eqs. (2) and (3). Subse-quently, the center crack length was extended from 2a to 2a þ 2Da by removing the covalent bonds in the crack tips. Because of the characteristics of the discrete structure, the crack extension lengthDa is equal to the lattice size of graphene, 2.45 Å. Berendsen NVT ensemble (Berendsen et al., 1984) was carried out on the new system (crack length equal to2a þ 2Da) for 100 ps to achieve the total energy minimization; meanwhile, the minimized value of the total energy was directly calculated from the MD simulation. It is noted in the NVT ensemble that the whole volume of the sys-tem is remaining fixed, and no external work was supplied into the system. As a result, the total energy variation being considered as the energy loss during the crack extension can be utilized to calcu-late the strain energy release rate of the atomistic graphene struc-tures as follows:
G ¼ U2aþ2Da U2a
2Dat ð7Þ
where U2aþ2Da is the potential energy of the graphene sheet with
crack length 2a þ 2Da, and U2ais the potential energy of the
graph-ene sheet with crack length 2a. In addition, ‘‘t” is the thickness of the graphene that is assumed to be 3.4 Å (Iijima et al., 1992).
2.3.2. Crack closure method
The energy difference during the crack extension in the atomis-tic structure can also be calculated from the crack closure concept where the amount of energy variation during the crack extension is equal to the work required to close the crack back to its original length ð2aÞ. This concept, originally proposed by Irwin (1957), was normally adopted to measure the strain energy release rate in continuum media. The work needed to close the crack for a Mode I fracture can be estimated as
DW ¼
1 2X3
i¼1
Fiyðuþiyþ uiyÞ ð8Þ
where Fiyis the concentrated force acting on atoms Aiand Biwhen
crack length is 2a, and uþ
iyand uiyare the corresponding opening
dis-placement for the atoms Aiand Bias shown inFig. 5when the crack
extends 2Da. In the MD simulation, the force Fiyis obtained by
arti-ficially removing the corresponding covalent bonds (one stretching bond and four angle bonds) and then calculating the unbalanced forces acting on atoms Ai and Bi. It is noted that the unbalanced
forces are only present on the atoms Ai and Bi. ði ¼ 1—3Þ when
the corresponding covalent bonds (one stretching bond and four an-gle bonds) are eradicated from the graphene structure, and thus the work required to close the crack can be evaluated in terms of the unbalanced forces and the corresponding displacement of the atoms. Because of symmetric geometry and loading conditions, the calculated unbalanced forces are almost the same, but they are acting in the opposite directions as illustrated inFig. 5. More-over, the displacement of atoms Aiand Biafter the crack extension
can be simply obtained by measuring the corresponding position vectors in the MD simulation. It is noted that in the continuum media, the work estimated from Eq. (8) can be directly used in the calculation of stain energy release rate. However, for the dis-crete atomistic system, in addition to the work done, the bonded energy that was eradicated during the crack extension has to be recovered in order for the discrete system to be completely back to its original state. In other words, the covalently bonded energy contributed from the stretching bond and the angle bending bond has to be evaluated and accounted for in the crack closure method. In each crack tip, the potential energy needed to be retrieved con-tains one stretching bond (bond a) and four angle bending bonds
(bonds b, c, d and e) as shown inFig. 6.Table 1illustrates the poten-tial energy associated with these interatomistic bonds. Thus, the recaptured potential energy in each crack tip can be written as Ebonds¼ ðEaþ Ebþ Ecþ Edþ EeÞ ð9Þ
A combination of Eqs.(8) and (9) yields the total energy release associated with the crack extension ofDa, from which the Mode I strain energy release rate for the graphene sheet is deduced as
GI¼ 1 2Datð2DW þ 2EbondsÞ ð10Þ
a
F
2yF
2yA
1 − 1yu
+ 1yu
+ 2yu
− 2yu
+ 3yu
− 3yu
F
1yF
1yF
3yF
3yA
2A
3B
1B
2B
3a
a
F
2yF
2yA
1 − 1yu
+ 1yu
+ 2yu
− 2yu
+ 3yu
− 3yu
F
1yF
1yF
3yF
3yA
2A
3B
1B
2B
3a
Fig. 5. Atomic force and opening displacement in Mode I fracture for graphene sheet.
a
A2
B3
A3
B1
A1
B2
a
b
c
e
d
a
A2
B3
A3
B1
A1
B2
a
b
c
e
d
A2
B3
A3
B1
A1
B2
a
b
c
e
d
Fig. 6. The bonded energy should be recovered in crack closure method for the atomistic graphene structure.
In the same manner, the work needed to close the crack for Mode II fracture can be estimated as
DW ¼
1 2 X3 i¼1 Fixðuþixþ u ixÞ ð11Þwhere Fixis a concentrated force acting on atoms Aiand Biin the x
direction. In addition, uþ
ixand uixare the sliding displacement for the
atoms Ai and Bi, respectively, when the crack extend in sliding
mode. Again, the strain energy release rate for the Mode II fracture in the atomistic graphite sheet can be expressed in terms of the crack closure method as
GII¼
1
2Datð2DW þ 2EbondsÞ ð12Þ
where EBondsdenotes the potential energy should be recovered from
the covalent bonds associated with the shear deformation of the graphite sheet.
3. Continuum model
In the previous section, the graphene was modeled as a discrete atomistic structure, and the energy difference during crack exten-sion was calculated from the MD simulation. On the other hand, when the graphene sheet was characterized as an isotropic homo-geneous solid, the fracture parameters, i.e., stress intensity factor and strain energy release rate, can be determined simply from the LEFM. It is noted that in the continuum model, the atomically blunt crack was assumed to be a line sharp crack in the LEFM anal-ysis. It is an interesting task to determine if the fracture parameters calculated from the LEFM would be compatible to those obtained from the atomistic graphene sheet model. To understand if there is a fracture parameter being able to bridge the gaps between the discrete and continuum models would be the main concern of the study. In order to resolve the aforementioned questions, the fracture parameters of the continuum model, i.e., strain energy
re-lease rate as well as the stress intensity factor, were calculated based on LEFM theory, and the results were compared to the cor-responding atomistic simulation.
For a finite plate with a center crack, the stress intensity factor can be approximated in a polynomial form as (Anderson, 1995)
K ¼
r
pffiffiffiffiffiffip
asecp
a 2W h i1=2 1 0:025 a W 2 þ 0:06 a W 4 ð13Þ where a is the half crack length, W is the half graphene width, andr
is the remote loading acting on the graphite sheet. Because the graphene sheet is pretty thin, the plane stress condition was as-sumed in the modeling, and thus the strain energy release rate can be expressed in terms of stress intensity factor as
G ¼K
2
E ð14Þ
where E is the Young’s modulus of the graphene sheet. It is noted that the material properties utilized in the continuum model was calculated based on the molecular dynamics simulation with the interatomistic energy described earlier, and the results are pre-sented inTable 2.
In addition to the analytical LEFM solution, the finite element method is an alternative technique to evaluate the fracture param-eters in continuum solids.Fig. 7illustrates the mesh of the graph-ene sheet. It is noted that the fine meshes in the vicinity of crack tips were generated in order to precisely simulate the singular stress field. For the comparison purpose, the same geometric con-figuration and loading condition as used in the atomistic simula-tion was introduced in the FEM model. The stress intensity factor of the continuum solid can be determined from the near-tip stress field (Sanford, 2003) as KI¼ lim r!0
r
yy ffiffiffiffiffiffiffiffiffi 2p
r p ð15Þ KII¼ lim r!0r
xy ffiffiffiffiffiffiffiffiffi 2p
r p ð16Þ where KI and KII are Mode I and Mode II stress intensity factor,respectively, and
r
yyandr
xyare the tensile and shear stresscompo-nents near the crack tip. In addition to the stress intensity factor, the strain energy release rate of the continuum graphene sheet can be determined from the finite element method using the crack closure integral (Sun and Jih, 1987) as
GI¼ lim da!0 1 2daf b yðu a y ua0yÞ ð17Þ GII¼ lim da!0 1 2daf b xðu a x ua0xÞ ð18Þ where fb
x and fyb denotes the nodal force at node b in the x and y
direction. ua
x and uay are the displacement components at node a,
and ua0
x and ua0y are the displacement components at node a’ as
shown inFig. 8. Table 1
Bonded energy for discrete graphite structure under mode I loading.
Stretch bond (a) Angle bond (b) (d) Angle bond (c) (e)
0.111 0.029 0.000027
Unit: kcal/mol.
Table 2
Material properties of graphene sheet obtained from MD simulation.
Young’s modulus (GPa) Shear modulus (GPa) Poisson’s ratio
790.7 311 0.27
4. Results and discussion
Fig. 9illustrates the stress distribution near the crack tip ob-tained from the continuum model (FEM solution and LEFM solu-tion) as well as the atomistic model associated with the same remote tensile loading. Apparently, the stress field calculated based on the LEFM solution and the FEM solution demonstrates 1=pffiffiffixstress singularity near the crack tip. On a contrary, in the dis-crete atomistic model, the stress field begins to deviate from the 1=pffiffiffix singularity solution in one lattice distance from the crack tip; therefore, there is no apparent stress singularity observed near the crack tip. This non-singular stress field in the atomistic model could be attributed to two possible reasons. One reason is that near to the crack tip where the deformation gradient is high, the expres-sion of Hardy stress is expected to be non-local (Eringen et al., 1977). The other one could be the atomically blunt crack tip in the atomistic model, which may be responsible for the non-singu-lar stress. No matter which one caused the non-singunon-singu-lar stress field, it is evident that adopting the stress intensity factor concept in the atomistic model may make the fracture parameter
ambigu-ous. In addition, it was found that when the position is beyond one lattice spacing from the crack tip, the stress field calculated from the discrete model coincides with the FEM solution well. The same attributes of the stress field was also observed in the Mode II frac-ture when the graphene sheet is subjected to remote shear loading. Because of the characteristic of the non-singular stress at the crack tip, it is suggested that the stress intensity factor defined based on singular stress assumption may not be suitable for describing the fracture in the discrete model.
The calculation of total energy variations obtained from the glo-bal energy method and crack closure method for Mode I crack extension in molecular dynamics simulation are illustrated in Ta-ble 3. Moreover, the strain energy release rates for the Mode I frac-ture calculated based on the discrete atomistic model and its continuum counterpart are presented inTable 4. It can be seen that the strain energy release rates obtained from the continuum model is in good agreement with those in discrete model. Moreover, both the global energy method and the crack closure method provide al-most the same strain energy release rate in the discrete graphene model. Since the global energy method in the MD simulation is rel-atively simple, it was suggested that it be utilized in the calculation of the strain energy release rate in the molecular simulation. Sim-ilar tendency for the Mode II fracture of graphene sheet was also observed, and the corresponding results are presented inTables 5
b
a
a’
a
δ
δ
a
b
y
f
b
x
f
b
a
a’
a
δ
δ
a
b
y
f
b
x
f
Fig. 8. Schematic of crack closure technique used in FEM analysis.
Distance(Å) σyy (GPa) 0 2 4 6 8 10 2 4 6 8 10 12 LEFMSolution FEM HardyStress(MD)
Fig. 9. Stress distribution near the crack tip obtained from the continuum model (FEM solution and LEFM solution) and atomistic model.
Table 3
Comparison of total energy variation calculated based on global energy method and crack closure method during crack extension in MD simulation (Mode I case).
Global energy method Crack closure method
Total energy variation DW Ebonds
1.95 1.61 0.34
Unit: kcal/mol.
Table 6
Comparison of strain energy release rate obtained from continuum model and discrete atomistic model (Mode II case).
Continuum model Atomistic model LEFM
(finite plate)
FEM Global energy method Crack closure method Strain energy release rate ðJ=m2Þ 0.0201 0.0203 0.021 0.0208 Table 5
Comparison of total energy variation calculated based on global energy method and crack closure method during crack extension in MD simulation (Mode II case).
Global energy method Crack closure method Change in total potential energy DW Ebonds
0.504 0.351 0.151
Unit: kcal/mol. Table 4
Comparison of strain energy release rate obtained from continuum model and discrete atomistic model (Mode I case).
Continuum model Atomistic model LEFM
(finite plate)
FEM Global energy method Crack closure method Strain energy release rate ðJ=m2Þ 0.081 0.0794 0.0814 0.0809
and 6. In light of the forgoing, it is revealed that the strain energy release rate can be regarded as a physical quantity that is capable of establishing connections between the atomistic simulation and the continuum media for characterizing the fracture of covalently bonded graphene sheet.
5. Conclusion
The fracture behavior of a graphene sheet with a center crack was characterized using atomistic simulation and linear elastic fracture mechanics (LEFM). In the atomistic simulation, the graph-ene was regarded as an atomistic structure, containing discrete carbon atoms; nevertheless, in the LEFM, it was modeled as an iso-tropic homogeneous media. Results from atomistic simulation indicated that because of the discrete attribute, there is no stress singularity near the crack tip. Therefore, the concept of stress intensity factor, which is generally employed in the continuum mechanics, may not be suitable for modeling the crack behavior in the atomistic graphene sheet. In order to validate the strain en-ergy release rate concept, the enen-ergy variation before and after the crack extension were evaluated in both continuum and atomistic model. For the discrete atomistic model, two methods, i.e., the glo-bal energy method and the crack closure method, were employed to compute the energy variation as well as the strain energy re-lease rate. On the other hand, the finite element analysis was per-formed in the continuum graphene sheet for the estimation of the strain energy release rate. It was denoted that the strain energy re-lease rate calculated based on the global energy method and crack closure method are almost the same. Furthermore, it was demon-strated that the global energy method is a simple manner in the atomistic simulation for the calculation of the strain energy release rate. A comparison of atomistic simulation with finite element re-sults illustrated that the strain energy release rates obtained from continuum model coincided with that in the discrete model asso-ciated with the same loading condition. As a result, the concept of strain energy release rate is regarded as a physical quantity that can establish connections between the atomistic simulation and continuum modeling for modeling the fracture of covalently bonded graphene sheet.
References
Anderson, T.L., 1995. Fracture Mechanics: Fundamentals and Applications, second ed. CRC Press, Inc., Boca Raton.
Basinski, Z.S., Duesbery, M.S., Taylor, R., 1971. Influence of shear stress on screw dislocations in a model sodium lattice. Canadian Journal of Physics 49, 2160– 2180.
Berendsen, H.J.C., Postma, J.P.M., van Gunsteren, W.F., DiNola, A., Haak, J.R., 1984. Molecular dynamics with coupling to an external bath. Journal of Chemical Physics 81, 3684–3690.
Cho, J., Luo, J.J., Daniel, I.M., 2007. Mechanical characterization of graphite/epoxy nanocomposites by multi-scale analysis. Composites Science and Technology 67, 2399–2407.
Cormier, J., Rickman, J.M., Delph, T.J., 2001. Stress calculation in atomistic simulations of perfect and imperfect solids. Journal of Applied Physics 89, 99– 104.
Cornell, W.D., Cieplak, P., Bayly, C.I., Gould, I.R., Merz Jr., K.M., Ferguson, D.M., Spellmeyer, D.C., Fox, T., Caldwell, J.W., Kollman, P.A., 1995. A second generation force field for the simulation of proteins, nucleic acids, and organic molecules. Journal of the American Chemical Society 117, 5179–5197.
Eringen, A.C., Speziale, C.G., Kim, B.S., 1977. Crack-tip problem in non-local elasticity. Journal of the Mechanics and Physics of Solids 25, 339–355. Hardy, R.J., 1982. Formulas for determining local properties in molecular-dynamics
simulations: shock waves. Journal of Chemical Physics 76, 622–628. Iijima, S., Ajayan, P.M., Ichihashi, T., 1992. Growth model for carbon nanotubes.
Physical Review Letters 69, 3100–3103.
Irwin, G.R., 1957. Analysis of stresses and strains near the end of a crack traversing a plate. Journal of Applied Mechanics 24, 361–364.
Jin, Y., Yuan, F.G., 2005a. Atomistic simulations of J-integral in 2D graphene nanosystems. Journal of Nanoscience and Nanotechnology 5, 1–9.
Jin, Y., Yuan, F.G., 2005b. Nanoscopic modeling of fracture of 2D graphene systems. Journal of Nanoscience and Nanotechnology 5, 601–608.
Li, C., Chou, T.W., 2003. A structural mechanics approach for the analysis of carbon nanotubes. International Journal of Solids and Structures 40, 2487–2499. Machida, M., 1999. Principles of Molecular Mechanics. John Wiley & Sons, Inc., New
York, USA.
Melchionna, S., Ciccotti, G., Holian, B.L., 1993. Hoover NPT dynamics for systems varying in shape and size. Molecular Physics 78, 533–544.
Miller, R., Tadmor, E.B., Phillips, R., Ortiz, M., 1998. Quasicontinuum simulation of fracture at the atomic scale. Modelling and Simulation in Materials Science and Engineering 6, 607–638.
Omeltchenko, A., Yu, J., Kalia, R.K., Vashishta, P., 1997. Crack front propagation and fracture in a graphite sheet: a molecular-dynamics study on parallel computers. Physical Review Letters 78, 2148–2151.
Rappe, A.K., Casewit, C.J., 1997. Molecular Mechanics Across Chemistry. University Science Books, Sausalito, Calif.
Sanford, R.J., 2003. Principles of Fracture Mechanics. Prentice Hall, Inc., Upper Saddle River, NJ.
Smith, W., Forester, T.R., 2001. The DL_POLY User Manual, version 2.13. Daresbury Laboratory.
Sun, C.T., Jih, C.J., 1987. On strain energy release rates for interfacial cracks in bi-material media. Engineering Fracture Mechanics 28, 13–20.
Webb III, E.B., Zimmerman, J.A., Seel, S.C., 2008. Reconsideration of continuum thermomechanical quantities in atomic scale simulations. Mathematics and Mechanics of Solids 13, 221–266.
Zhou, T., Yang, X., Chen, C., 2009. Quasicontinuum simulation of single crystal nano-plate with a mixed-mode crack. International Journal of Solids and Structures 46, 1975–1980.
Zimmerman, J.A., Webb III, E.B., Hoyt, J.J., Jones, R.E., Klein, P.A., Bammann, D.J., 2004. Calculation of stress in atomistic simulation. Modelling and Simulation in Materials Science and Engineering 12, S319–S332.
