Characterizing mechanical properties of graphite using molecular dynamics simulation

(1)

Characterizing mechanical properties of graphite using molecular

dynamics simulation

Jia-Lin Tsai

*

_{, Jie-Feng Tu}

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 19 March 2009 Accepted 19 June 2009 Available online 23 June 2009 Keywords:

Nano materials (A) Mechanical (E) Atomic structure (F)

a b s t r a c t

The mechanical properties of graphite in the forms of single graphene layer and graphite flakes (contain-ing several graphene layers) were investigated us(contain-ing molecular dynamics (MD) simulation. The in-plane properties, Young’s modulus, Poisson’s ratio, and shear modulus, were measured, respectively, by apply-ing axial tensile stress and in-plane shear stress on the simulation box through the modified NPT ensem-ble. In order to validate the results, the conventional NVT ensemble with the applied uniform strain filed in the simulation box was adopted in the MD simulation. Results indicated that the modified NPT ensem-ble is capaensem-ble of characterizing the material properties of atomistic structures with accuracy. In addition, it was found the graphene layers exhibit higher moduli than the graphite flakes; thus, it was suggested that the graphite flakes have to be expanded and exfoliated into numbers of single graphene layers in order to provide better reinforcement effect in nanocomposites.

1. Introduction

With the characteristics of high strength and stiffness, the graphite has been used as reinforcements in composite materials

[1]. The natural graphite is constructed by numbers of graphene layers with interlayer spacing of around 3.4 Å. Through chemical oxidation in the environment of sulfuric and nitric acid, the acid intercalant can be intercalated into the graphite galleries to form an intercalated graphite compound. Subsequently, by applying ra-pid heating because of the vaporization of the acid intercalant in the graphite galleries, the interacted graphite was significantly expanded along the thickness direction and converted into the expanded graphite (EG). After a mechanical mixer together with sonication process, the expanded graphite was dispersed and exfoliated into the polymer matrix to form graphite-reinforced nanocomposites. The synthesizing process for manufacturing the nanocomposites was discussed in detail in the literatures [2,3]. Recently, Si and Samulski [4] employed platinum nanoparticles adhered to the graphene to prevent the aggregation of isolated graphene sheets during the drying process. However, the stacked graphene structures (so called graphite flakes) are commonly observed in TEM micrographs and XRD examination[5], and it is a challenging task to fully exfoliate the aggregated graphene sheets. In fact, graphite flakes together with graphene layers are commonly observed in graphite nanocomposites and it is impor-tant to clarify if the two atomistic configurations of the graphite,

i.e., graphite ﬂakes and single graphene layer, would have the same mechanical properties. Moreover, in order to accurately character-ize the mechanical properties of the graphite-reinforced nanocom-posites, an exploration of the fundamental properties of the graphite associated with different microstructures is required.

Cho et al.[6]performed a molecular structural analysis to calcu-late the graphite’s elastic constants. The in-plane properties of graphite were derived by considering the geometric deformation of a single graphene sheet subjected to in-plane loading. However, for the out-of-plane properties, they modeled the graphic as gra-phic flake with multi graphene layers, the non-bonded atomistic interactions of which were described using Lennard-Jones poten-tial function. Scarpa et al.[7]proposed truss-type analytical mod-els in conjunction with cellular material mechanics theory to describe the in-plane elastic properties of single layer graphene sheets. It was found that the analytical results and the numerical results obtained from finite element shows good agreement with existing numerical values. Hemmasizadeh et al.[8]who correlated the force-depth results obtained from MD simulation with the large deflection formulation of circular plates loaded at center to evaluate the effective Young’s modulus of graphene sheet and the corresponding wall thickness of the single sheet. By using MD simulation, Bao et al.[9]investigated the variations of Young’s modulus of graphite, which contains different numbers of graph-ene layers (one to five layers). Results indicated that there is no considerable difference in Young’s modulus between the single layer of graphene and graphite flakes with five layers of graphene. Reddy et al.[10]modeled the elastic properties of a finite-sized graphene sheet using continuum mechanics approach based on

*Corresponding author. Tel.: +886 3 5731608; fax: +886 3 5720634. E-mail address:[email protected](J.-L. Tsai).

Contents lists available atScienceDirect

Materials and Design

(2)

Brenner’s potential [11]. The computed elastic constants of the graphene sheet are found to follow the orthotropic material behav-iors. In addition to the fundamental material properties, the vibra-tional responses of single layer graphene sheets were investigated through a molecular structural mechanics approach [12]. Both mode shapes and natural frequency of single graphene sheet were considered in their analysis. In light of the forgoing, most studies characterize the elastic properties of the graphite based on the behavior of a single graphene sheet; the mutual influences of the adjacent graphene layers on the mechanical responses in the graphite flakes are rarely taken into account. As previously mentioned, both the graphene layer and the graphite flakes were commonly found in the nanocomposites, so it is not adequate to utilize the properties of the graphene sheet instead of the graphite flakes in the modeling of graphite-reinforced nanocomposites.

In this study, the mechanical properties of the graphite flakes and the graphene were systematically characterized using MD simulation. Both bonded and non-bonded interactions were ac-counted for in the description of the atomistic graphite structures. By applying uniaxial tensile loading on the atomistic graphite structures, the Young’s modulus and Poisson’s ratio were deter-mined from the strain field in the deformed configuration. In the same manner, the shear modulus was predicted from the shear deformation associated with the applied shear stress. The proper-ties of the single graphene layer were then compared to those of the graphite flakes with multi-layers of graphene.

2. Molecular dynamics simulation

2.1. Construction of atomistic structures of graphite

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. In naturally occurring or high quality synthetic graphite, the car-bon layers are attacked along the thickness direction in AB type se-quence with interlayer spacing of approximately 3.4 Å as shown in

Fig. 1. Hereafter, the graphite with several carbon layers lumped together is referred to as graphite ﬂakes. Because the adjacent car-bon layers are held together by the weak van der Waals force, after proper processing[2,3], the stacked carbon layers can be dispersed and separated into a single layer that is usually called graphene sheet or graphene layer.

In order to investigate the mechanical properties of the graphite ﬂakes and the graphene layer, the atomistic structures have to be constructed in conjunction with the appropriately speciﬁed

atom-istic interaction. In the description of graphite structure, two kinds of atomistic interactions are normally taken in account; one is bonded interaction, 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 be-tween two neighboring carbon atoms that provides the building block of the primary 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 stretch-ing, bond angle bendstretch-ing, torsion, and inversion as illustrated in

Fig. 2 [13]. Therefore, the total potential energy of the graphite con-tributed from the covalent bond is given 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[14]:

Ur¼ 1 2krðr r0Þ 2 ð2Þ Uh¼ 1 2khðh h0Þ2 ð3Þ

where krand khare the bond stretching force constant and angle

bending force constant, respectively. The constants kr¼

93; 800 kcal

mole nm2 and kh¼ 126_{mole rad}kcal 2 selected from AMBER force ﬁeld for carbon–carbon atomic-interaction [15]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 graphite atomistic structures.

In addition to the bonded interaction, the non-bonded interac-tion between the carbon atoms was regarded as the van der Waals force, which can be characterized using the Lennard-Jones (L-J) po-tential as UvdW¼ 4u r0 rij 12 r0 rij 6 " # ð4Þ

where rijis the distance between the non-bonded pair of atoms. For

the hexagonal graphite, the parameters u = 0.0556 kcal/mole and r0= 3.40 Å suggested in the literature [16] were adopted in the

modeling. Moreover, the cutoff distance for the van der Waals force

1.42 Å 3.4 Å

A

B

1.42 Å 3.4 Å

A

B

Fig. 1. Schematic of graphite structures.

i

j

i n k j

(b)

(a)

(c)

(d)

i n k j

ω

j

i

k

i

j

i

j

i n k j i n k j i n k j

ω

i n k j

ω

j

i

k

j

i

k

φ θ

Fig. 2. A schematic of the inter-atomic potential: (a) bond stretch; (b) bond angle bending; (c) dihedral angle torsion; (d) inversion

(3)

is assigned to be 10 Å, which means that beyond this distance, there are no more van der Waals interactions taking place.

In order to model the material properties of graphite ﬂakes and the graphene sheet, the simulation box suitable for representing the corresponding atomistic structures has to be established.

Fig. 3shows the schematic of the simulation box for graphite flakes and the graphene sheet as well. A periodic boundary condition was implemented on all surfaces to demonstrate the infinite graphite structures. It is noted in the graphene sheet that the dimension of the simulation box in the thickness direction is set to be large enough that the van der Waals interaction between the neighbor-ing layers can not be attained. This especial design of the simula-tion box is intended to simulate the exfoliated graphene sheets. The atomistic structures with stress-free configuration were ob-tained by performing an NPT ensemble with time increments at 1 fs for 100 ps (the total iteration steps are 100,000) until the potential energy accomplished a stable value. The NPT ensemble represents that the pressure and temperature of the system may approach to the specified values during the MD simulation. This process is accomplished through the Berendsen thermostat

[17,18]by scaling the velocities and positions of atoms at each step to push the temperature and pressure toward the desired value (P = 0 and T = 0 K).Figs. 4 and 5demonstrate the potential energy

history and the temperature variation, respectively, for the graphene sheet during the NPT ensemble. It seems that the potential en-ergy attains a stable value rapidly, and the temperature also approaches to 0 K. Based on the observations, it was suggested that the current atomistic structure is in the equilibrium condition and suitable for the characterization of the material properties. The MD simulation was carried out under the DL-POLY package originally developed by Daresbury Laboratory[19]in conjunction with the homemade subroutine for post-processing.

2.2. Characterizing the Young’s modulus and Poisson’s ratio of the graphite

The methodology developed to evaluate the mechanical proper-ties of the atomistic structures was motivated from the technique commonly used in the continuum solid. For continuum solids, the Young’s modulus and Poisson’s ratio are measured from the simple tension test. The same concept was extended and applied to the atomistic structures by means of a modiﬁed NPT ensemble in MD simulation with the characteristics of varying a simulation box in shape and size[20]. In other words, axial stresses can be implemented on both sides of the simulation box with other faces’ being traction free as shown inFig. 6. Again, after the energy min-imization process, the equilibrated graphite atomistic structure under axial loading was obtained, and the Young’s modulus and Poisson’s ratio was deﬁned in the continuum manner as

27.2 Å A B A B A B A B 3.4 Å 34.0 Å 3.4 Å

(a)

(b)

Fig. 3. Schematic of atomistic model in the MD simulation for: (a) graphene sheet and (b) graphite ﬂakes.

time (ps) T o ta le ne rgy (k ca l/m ol ) 25 50 75 100 1194.1928 1194.1930 1194.1933 1194.1935 1194.1938 1194.1940 1194.1943 1194.1945 1194.1948 1194.1950

Fig. 4. Variation of potential energy in the modiﬁed NPT ensemble for a graphene sheet. time (ps) Temperature (K) 25 50 75 100 -1E-06 0 1E-06 2E-06 3E-06 4E-06 5E-06 6E-06 7E-06 8E-06 9E-06 1E-05 1.1E-05 1.2E-05 1.3E-05

Fig. 5. Variation of temperature in the modiﬁed NPT ensemble for a graphene sheet.

x

1

x

2

x

3

(a)

(b)

(4)

E ¼

r

e

1 ð5Þ

m

12¼

e

2

e

1 ð6Þ

where e1is the strain component measured in the loading direction,

and e2is the strain component measured in the lateral direction. As

is noted in Eq.(5),

r

should be the stress directly acting on the graphite structure. However, for the case of graphene sheet as shown inFig. 6a, because the dimension of the graphene sheet in thickness direction is not compatible to the size of the simulation box, the stress in the graphene sheet has to be converted from the stress acting on the simulation box,

r

box, in terms of the

geo-metric parameters as

r

¼

r

boxh

t ð7Þ

where h is the height of the simulation box, and t is the thickness of the graphene sheet, which is equal to 3.4 Å. The Young’s modulus and Poisson’s ratio obtained from MD simulation for the graphite ﬂakes and graphene sheet are presented, respectively, inTable 1.

Instead of the application of axial stress, the graphite properties can be calculated alternatively by applying a small amount of strain on the simulation box with the other strain components remaining at zero. Because this approach is commonly employed in MD simulation to determine the Young’s modulus of atomistic structures[21–23], it was utilized to verify the results presented in the previous section. When the graphite in the stress-free state was obtained through MD simulation, a small amount of axial strain component was applied on the simulation box while the other strain components remain at zero as shown inFig. 7. After an energy minimization process in the NVT ensemble where the

volume and temperature are ﬁxed during the simulation, the asso-ciated stress in the deformed conﬁguration of the atomistic struc-ture was calculated through the virial theorem[18]as

r

ij¼ 1 V0 X i<j rijfijT ! ð8Þ

In Eq.(8), rijand fijdenotes the atomic distance and the

corre-sponding interaction force between any two atoms. Vorepresents

the total volume of the simulation box. In the same manner, the stress component on the simulation box has to be converted into the one on the graphene sheet by using Eq.(7). It is noted that in Eq. (8), because the atomistic model was simulated at 0 K, the velocity term caused by temperature was neglected in the stress computation. For the graphite layer structures, the constitutive relation under plane stress assumption can be written in terms of the stiffness matrix Cijas

r

11

r

22

s

12 8 > < > : 9 > = > ;¼ C11 C12 0 C21 C22 0 0 0 C66 2 6 4 3 7 5

e

11

e

22

c

12 8 > < > : 9 > = > ; ð9Þ

and the values of the entries in the Cijmatrix can be calculated from

the measured stress components corresponding to the applied strain ﬁled in the NVT ensemble. Once the values were determined, the Young’s modulus and Poisson’s ratio of the graphite can be de-duced as[24] E ¼C11C22 C12C21 C22 ð10Þ

m

12¼ C12 C22 ð11Þ

Table 2 shows the values of the elements in the stiffness matrix with regard to the applied strain of 0.001. In addition, the calculated modulus and Poisson’s ratio are also included in the same table. It was found that the modulus and Poisson’s ratio derived based on Eqs.(10) and (11)are in agreement with the predictions obtained Table 1

Comparison of in-plane elastic constants of graphite ﬂakes and graphene sheet obtained from the modiﬁed NPT ensemble.

E (TPa) m G12(TPa)

Graphene 0.912 0.261 0.358

Graphite ﬂakes 0.795 0.272 0.318

(b)

(a)

Fig. 7. Axial deformation applied in: (a) graphene sheet and (b) graphite ﬂakes (the solid line denotes the undeformed conﬁguration and the dashed line indicates the deformed shape).

Table 2

The elements in the stiffness matrix for the graphene sheet and graphite ﬂakes calculated based on axial strain deformation in conventional NVT ensemble.

C11(GPa) C12(GPa) C21(GPa) C22(GPa) E (TPa) m

Graphene 977.91 254.46 254.74 978.19 0.912 0.26

Graphite ﬂakes 864.29 235.16 235.21 864.34 0.8 0.272

(a)

(b)

(5)

from the uniaxial stress method. This conformity clearly reveals that the modiﬁed NPT ensemble with the technique of applied load-ing is suitable for characterizload-ing material properties of the atomistic structure. In the next section, the modiﬁed NPT ensemble will be adopted to characterize the shear modulus of the graphite with re-spect to the applied shear loading.

2.3. Characterization of shear modulus G12

By following the same technique used in the early section, the shear modulus of the atomistic structure can be evaluated via the application of in-plane shear stress on the simulation box as shown inFig. 8. This process is accomplished by conducting the modified NPT ensemble in MD simulation. After the energy mini-mization process, the deformed configuration of the simulation box was calculated from which the shear strain associated with the applied shear stress was determined. If the deformation is small, the shear modulus of the graphite can be defined based on the theory of linear elasticity as

G ¼

s

c

ð12Þ

where

s

is the applied shear stress, and

c

is the corresponding shear strain determined from MD simulation. The shear moduli calculated with Eq.(12)for the graphite ﬂakes and graphene sheet are also listed inTable 1.

3. Results and discussion

Results presented inTable 1indicate that the single graphene sheet demonstrates higher Young’s modulus and shear modulus than the graphite flakes. Thus, it was suggested that to achieve bet-ter mechanical properties of nanocomposites, the aggregated graphite flakes need to be exfoliated in the form of graphene sheets and uniformly dispersed into the matrix systems. Moreover, according to the relationship between Young’s modulus, shear modulus, and Poisson’s ratio, it was found that both graphite flakes and graphene demonstrate isotropic in-plane properties. This iso-tropic property could be attributed to the hexagonal array of the carbon atoms.

For the purposes of comparing, the calculated material proper-ties of the graphite are listed together with other published predic-tions inTable 3. It was revealed that the moduli obtained from the current model are a little less than those listed in the literature although the discrepancy is not much. This difference could be resulting from the different potential functions employed in the modeling of the atomistic interaction of the carbon atoms. On the other hand, it should be indicated that most of the published values are calculated based on the graphene sheet except the one addressed by Bao et al.[9]who investigated the Young’s modulus of graphite with numbers of graphene layers (from one layer up to five layers). In their investigation, there is no significant difference in Young’s modulus between the single graphene layer and the graphite flake with five-layer graphene. It is possible that the dis-similarity may not be considerable just by comparing the single layer graphene with the five-layer graphene. On the contrary, our prediction considers the periodic boundary condition in the

thick-ness direction and would be close to the behavior of the graphite ﬂakes with numbers of graphene layers. This is the reason why in our simulation, the graphic ﬂakes would exhibit different material properties from the single graphene sheet. In addition, for the sake of comparison, the experimental values for the graphite structures provided by Blakslee et al.[25]were added in

Table 3. It shows that the values predicted based on the graphene sheet model have a better agreement with the experimental data.

4. Conclusions

The in-plane properties of graphene sheet and the graphite flakes were investigated using MD simulation. Two approaches were introduced to calculate the Young’s modulus and Poisson’s ra-tio: one is the applied axial stress on the simulation box under the modified NPT ensemble, and the other is the applied axial strain on the simulation box through the NVT ensemble. The values calcu-lated based on the two approaches have good agreement with each other, and the modified NPT ensemble, similar to the continuum mechanics approach, was regarded as an effective way to character-ize the material properties of atomistic structures. Furthermore, the shear modulus was evaluated by taking the ratio of the applied in-plane shear stress to the corresponding shear strain. Because of the hexagonal array of the carbon atoms, the in-plane shear modulus, Young’s modulus, and Poisson’s ratio of the graphite flakes and graphene sheet satisfy the isotropic properties. A comparison of in-plane properties of the graphene sheet and graphite flakes re-veals that the single graphene sheet exhibits higher modulus than the graphite flakes; therefore, the exfoliation of the graphite flakes into graphene layers is essential in order to have better mechanical properties of graphite-reinforced nanocomposites.

References

[1] Fukushima H, Drzal LT. Graphite nanoplatelets as reinforcements for polymers: structural and electrical properties. In: Proceedings of the American society for composites, The 17th technical conference; October 2002. p. 21–23. [2] Li J, Kim JK, Sham ML. Conductive graphite nanoplatelet/epoxy

nanocomposites: effect of exfoliation and UV/ozone treatment of graphite. Scripta Mater 2005;53:235–40.

[3] Fukushima H, Drzal LT, Rook BP, Rich MJ. Thermal Conductivity of Exfoliated Graphite Nanocomposites. J Therm Anal Calorim 2006;85(1):235–8. [4] Si Y, Samulski ET. Exfoliated graphene separated by platinum nanoparticles.

Chem Mater 2008;20:6792–7.

[5] Yasmin A, Luo JJ, Daniel IM. Processing of expanded graphite reinforced polymer nanocomposites. Compos Sci Technol 2006;66:1182–9.

[6] Cho J, Luo JJ, Daniel IM. Mechanical characterization of graphite/epoxy nanocomposites by multi-scale analysis. Compos Sci Technol 2007;67: 2399–407.

[7] Scarpa F, Adhikari S, Phani AS. Effective elastic mechanical properties of single layer graphene sheets. Nanotechnology 2009;20:065709.

[8] Hemmasizadeh A, Mahzoon M, Hadi E, Khandan R. A method for developing the equivalent continuum model of a single layer graphene sheet. Thin Solid Films 2008;516:7636–40.

[9] Bao WX, Zhu CC, Cui WZ. Simulation of Young’s Modulus of single-walled carbon nanotubes by molecular dynamics. Physica B 2004;352:156–63. [10] Reddy CD, Rajendran S, Liew KM. Equilibrium conﬁguration and continuum

elastic properties of ﬁnite sized graphene. Nanotechnology 2006;17:864–70. [11] Brenner DW. Empirical potential for hydrocarbons for use simulating the

chemical vapor deposition of diamond ﬁlms. Phys Rev B 1990;42(15): 9458–71.

[12] Sakhaee-Pour A, Ahmadian MT, Naghdabadi R. Vibrational analysis of single-layered graphene sheets. Nanotechnology 2008;19:085702.

[13] Rappe AK, Casewit CJ. Molecular mechanics across chemistry. Sausalito (CA): University Science Books; 1997.

Table 3

Comparison of the predicted values with others listed in the literatures.

Graphene Graphite ﬂakes Cho et al.[6] Bao et al.[9] Reddy et al.[10] Blakslee[25]

E (TPa) 0.912 0.795 1.153 1.026 0.671 1.020

m 0.261 0.272 0.195 – 0.428 0.160

(6)

[14] Li C, Chou TW. A structural mechanics approach for the analysis of carbon nanotubes. Int J Solids Struct 2003;40(10):2487–99.

[15] Cornell WD, Cieplak P, Bayly CI, Gould IR, Merz Jr KM, Ferguson DM, et al. A second generation force ﬁeld for the simulation of proteins, nucleic acids, and organic molecules. J Am Chem Soc 1995;117(19):5179–97.

[16] Battezzatti L, Pisani C, Ricca F. Equilibrium conformation and surface motion of hydrocarbon moleculars physisorbed on graphite. J Chem Soc 1975;71: 1629–39.

[17] Berendsen HJC, Postma JPM, van Gunsteren WF, DiNola A, Haak JR. Molecular dynamics with coupling to an external bath. J Chem Phys 1984;81(8): 3684–90.

[18] Allen MP, Tildesley DJ. Computer simulation of liquids. Oxford: Clarendon Press; 1990.

[19] Smith W, Forester TR. The DL_POLY user manual, version 2.13. Daresbury (Warrington, England): Daresbury Laboratory, CCLRC; 2001.

[20] Melchionna S, Ciccotti G, Holian BL. Hoover NPT dynamics for systems varying in shape and size. Mol Phys 1993;78(3):533–44.

[21] Agrawal PM, Sudalayandi BS, Raff LM, Komanduri R. A comparison of different method of Young’s Modulus determination for single-wall carbon nanotubes (SWCNT) using molecular dynamics (MD) simulations. Comput Mater Sci 2006;38:271–81.

[22] Suzuki K, Nomura S. On elastic properties of single-walled carbon nanotubes as composite reinforcing ﬁllers. J Compos Mater 2007;41(9): 1123–35.

[23] Adnan A, Sun CT, Mahfuz H. A molecular dynamics simulation study to investigate the effect of ﬁller size on elastic properties of polymer nanocomposites. Compos Sci Technol 2007;67(3–4):348–56.

[24] Gere JM. Mechanics of materials. Cheltenham (UK): Nelson Thornes; 2001. [25] Blakslee OL, Proctor DG, Seldin EJ, Spence GB, Weng T. Elastic constants of