Molecular dynamics simulation of the mechanical properties of multilayer graphene oxide nanosheets

Multilayer graphene oxide (GO) is an attractive candidate for new applications in nanoelectromechanical materials and structural reinforcement nanocomposites due to its strong mechanical properties. In this study, the mechanical properties and failure mechanism of multilayer GO nanosheets were studied by non-equilibrium molecular dynamics simulation. The simulated Young's modulus, fracture stresses, and fracture strains were found to be consistent with the experimentally measured values. The effects of the surface oxidation content of GO and the stacking layer number on these mechanical properties were investigated. The oxidation content has a larger influence on the mechanical properties compared with the layer number. The failure of multilayer GO nanosheets undergoes a relatively slow cracking process due to the existence of functional groups and the stacking layers. There appears to be different twodimensional stress distributions on multilayer GO sheets from the outer layer to inner layer. The Young's modulus and the fracture strength of the middle layer are generally larger than those in the outer layer. The fracture of the outside GO sheet begins first, and then the failure of the inner GO sheet occurs with a delayering process. The simulation result is expected to improve understanding of the mechanical behavior of multilayer GO nanosheets.


Introduction
Graphene oxide (GO) inherits the two-dimensional structure of graphene, but contains surface functional groups, mainly including hydroxyl and epoxy groups in the plane, as well as carbonyl and carboxyl groups on the edge. 1 Different from monolayer GO nanosheets 2-4 and GO paper, GO nanosheets 5-7 exhibit unique properties due to their intermediate transition structure and have promising applications in a wide range of elds, such as in structural reinforced composites, [8][9][10] super capacitors, 11 exible electronic materials, 12,13 and biology. 14 For GO nanosheets, excellent mechanical properties such as high strength and strong toughness are generally essential. To further promote the development of GO in engineering applications by regulating the structure of GO, it is necessary to have a good understanding of its mechanical behavior.
At present, a considerable number of experiments have been conducted to measure the mechanical properties of monolayer GOs and GO paper. 2,4,15,16 There is a signicant difference in the mechanical properties of GO materials on various thickness scales. For example, when the thickness of GO was reduced, the elastic modulus increases rapidly. Different fracture mechanisms have been proposed for monolayer GO and GO papers. It is generally accepted that the fracture of monolayer GO shows an intraplanar structure cleavage, 16 whereas GO paper typically fails by an interplanar delamination and shear deformation. 4 Molecular simulation 17,18 has been applied to study the mechanical behavior of single-layer GO sheets. It has been found that the relative concentration of hydroxyl and epoxide groups has larger effect on the deformation and failure behavior of single-layer GO. Mechanical properties of single-layer GO sheet with different types of functional groups have been simulated, respectively, to characterize the effect of functionalization. However, these functional groups can not represent the structure and composition on actual GO sheets. Currently, very limited works have been reported for the mechanical properties of multilayered GO nanosheets. It is highly interest to recognize the different mechanical properties and fracture mechanisms between single-layer and multilayer GO sheets. Wang 19 et al. provided the preliminary study on fracture behavior of multilayer GO, in which, only single stress-strain curve was characterized with insufficient fracture mechanism. In a recent study, 20 the overall strength and fracture behavior of GO nanosheets were experimentally examined. This study demonstrated the inherent defects and precracked layer determined the fracture behavior of actual GO materials. However, these actual defects might disguise the essential mechanical behavior. Currently, the mechanical properties and fracture mechanism of defect-free GO nanosheets still remains largely unexplored.
In this paper, Reaxff-based molecular dynamic (MD) simulations were performed to comprehensively investigate the mechanical properties of multilayer GO nanosheets. The fracture behavior of each single layer in the multilayer GO has been separately revealed. The different fracturing performance for individual single layer in the GO nanosheets has been displayed. Various thicknesses of GO nanosheets were considered and the effects of both functional groups and layer number of GO sheet were investigated.

Methodology and simulation details
The Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) 21 was utilized to perform the MD simulation. We describe the bonded and nonbonded interaction between different atoms in the GO systems with ReaxFF potential. [22][23][24] The potential energy in Reaxff eld includes several different components as following: where E bond are the bond energies; E over , the energy to penalize over coordination of atoms; E under , the energy to stabilize undercoordination of atoms; E val , the valence-angle energies; E pen , the penalty energy of over/under coordination in central atom; E tors , the torsion-angle energies; E conj , the energy of the conjugation effects to the molecular energy; E Coulomb , the nonbonded Coulomb energies; and E vdWalls , the nonbonded van der Waals interactions. ReaxFF includes a bond order/bond energy relationship, which allows for bond formation and dissociation during molecular dynamics simulations. Representative GO nanosheets were chosen with a planar dimension of 4.7 nm by 4.7 nm and various thicknesses, depending on the number of stacking layers. The hydroxyl and epoxy functional groups randomly and uniformly distributed on the both sides of the basal plane of GO, based on the Lerf-Klinowski model, 25 in which the ratio of hydroxyl groups and epoxy groups was set to be 1 : 1. The carbonyl and carboxyl groups on the edge of GO plane were not considered in this study. These groups have weak effects on the mechanical properties because they only represent small portion of atoms with dangling bonds only attached to plane edge. 26 In the present work, we did not consider the interlayer water molecules, because the main goal of this work is to study the intrinsic mechanical behavior of multilayer GO by comparing with single-layer GO. It has been found that the elastic modulus determined by the interlamellar hydrogen bond network, which was mediated by functional groups and water molecule, is relatively weak compared with the intrinsic stiffness of the GO sheets. 7 We constructed 5 structure models stacking with one to ve layers of GO nanosheets, named with GO-I, GO-II, GO-III, GO-IV, GO-V, respectively. We used ve different oxidation contents from 20% to 40%, which is dened as the oxygen-to-carbon ratio 26 for the different stacking layers of GO nanosheets.
Prior to the uniaxial tensile loading of the GO nanosheets, we rstly made the initial GO model to a local minimum conguration with a conjugate-gradient algorithm. The system was then heated to 300 K with the Nose-Hoover thermostat and relaxed at this temperature for 25 ps to bring the system into an equilibrated conguration. As shown in the Fig. 1, during the heating and relaxation process, the periodic boundary was performed along the in-plane of GO nanosheets, and the outplane of GO nanosheets was free. 20 The boundary condition was selected to eliminate large undulations associated with free planar edges. 26 A time step of 0.25 fs was used for all MD simulation. As shown in the Fig. 1, uniaxial tensile loading was applied by dilating the GO nanosheets under an equal affine transformation to atom position with dilating the simulation cell along the loading direction. 20 This deformation methodology is equivalent to uniaxial tensile loading along the loading direction. We carried out the uniaxial tensile loading of GO nanosheet with a constant strain rate of 0.0001 ps À1 . In the uniaxial tensile simulation, the stress was calculated as 27 where s ab is the average atomic-level stressor tensor. F b ij is the interatomic forces including bonding and nonbonding forces between atom i and atom j. r a ij is the distance between atom i and atom j, V 0 is the initial volume of the structure, and the interlayer spacing of the multilayer GO models was set to be Loading is performed along the armchair direction which is normalized against the assumed interlayer spacing. 0.7 nm, which is the effective layer thickness that has been generally used in the previous study of GO lms. 15,20,28 Result and discussion At rst, the uniaxial tensile simulation of a series of monolayer GOs with different oxidation contents was carried out along the directions of armchair and zigzag. We obtained the corresponding Young's modulus, fracture strength, and fracture strain respectively. As shown in Fig. 2, the Young's modulus of monolayer GO generally decreases with the increase of oxidation content, decreasing from $450 GPa at 20% of oxidation content to about $200 GPa at 40%. The results agree well with the experimentally measured values of 250 AE 150 GPa by Gomez-Navarro et al. 29 and 207.6 AE 23.4 GPa by Suk et al., 15 as well as the previously simulated values. 17,30 This monotonous decrease of Young modulus with the oxidation content can be explained as the fact that the addition of oxygen-containing groups makes the strong sp 2 carbon bonds transform to weak hybrid sp 3 carbon bonds, thereby reducing the energy stability of GO basal plane. 17,18,30 As shown in Fig. 2(a), the Young's modulus in armchair direction is slightly less than that in zigzag direction, which is probably due to different random arrangement of epoxide groups on GO plane. For example, if the epoxy groups, attached on the C-C bonds parallel to the armchair direction, have larger proportion than those attached on other C-C bonds, this arrangement of epoxy functional groups could lead to the increase of fracture strain along the armchair direction. However, the discrepancy of Young's modulus between the two directions can be negligible, comparing to the effect of oxidation content. This insensibility of chirality for Young's modulus of monolayer GO sheet has been similarly observed in previous researches. 17,19 As expected, we nd that the fracture strength ( Fig. 2(b)) of monolayer GO decreases with the oxidation content. This is in good agreement with the previous results. 18 Fig. 2(c) shows the fracture strain of GOs under various contents of oxidation is in the range of 10-16%. This result conforms to the simulated results of 10-15%, 19,30 as well as the measured values of 10-13%. 19 In short, the simulated mechanical properties of monolayer GO in this work are in good agreement with the previous reference results.
In this work, we mainly simulated the mechanical properties of multilayer GOs. Herein, we only conducted uniaxial tensile simulation along the armchair direction. Fig. 3(a) shows the typical stress-strain curve under uniaxial loading for the 5-layer stacked GO nanosheets (GO-V) with the oxidation content of 25%. The stress-strain curve can visually capture the changes in the stress with the strain during the tensile process. In the stress-strain curve, point A, B, and C represent the characteristic points on the tensile process. The stress-strain curve before point A shows the linear elastic behavior of GO with small strain, and the carbon structure on the basal plane was not destroyed yet. As the strain increases, the plastic deformation appeared with the crack and defects of the GO nanosheets aer the point A. The maximum stress of the multilayer GO   The Young's modulus was determined by tting the linear region of the loading stress-stress curve. 20 Meanwhile, the fracture strength and fracture strain were generally dened as the point when the tensile stress reaches the maximum value. 31 For the 5-layer stacked GO nanosheets under the oxidation content of 25%, the simulated Young's modulus is $314 GPa, and the fracture strength is $31 GPa. Cao et al. 20 performed a tensile test on multilayer GO nanosheets by using SEM in microeletromechanical systems (MEMS). The Young's modulus and the fracture strength for ve GO samples with thickness ranging from 75 AE 13 nm to 24 AE 4 nm were found to be 103-291 GPa and 4-11 GPa, respectively. Although our simulated mechanical properties of GO nanosheets are generally larger than the experimental values, both of them exhibit the same order of magnitude. The main reason for this discrepancy is the existing structural defect in real GO nanosheets, 13,14 which can cause attenuation of GO strength. Meanwhile, the jagged stacking feature of actual GO nanosheets also causes a decrease of mechanical properties.
Upon reaching the maximum tensile stress, the stress-strain curve of GO-V demonstrates slow and smooth decline aer the failure point, which is obviously different from that observed in monolayer GO. Generally, monolayer GO undergoes drastic and brittle failure near maximum loading 31 for hydroxyl-rich GOs and ductile deformation for epoxide-rich GOs, which is due to the epoxide-to-ether transformation, thus leading to the wavelike fracture behavior. This distinguished tensile performance of multilayer GO nanosheets is probably caused by the mutual interference between the neighboring layers. In addition, it is noted that the wave-like fracture behavior of the simulated multilayer GO stress-strain curve in this work is also distinct from the previous reported result, 19 in which the large brittle behavior has been observed. We can conclude that the increased percentage of hydroxyl groups used in the previous work 19 might cause GO to behave as a brittle material. Fig. 3(b) shows a snapshot of the overall structural evolution of 5-layer GO nanosheets under different strains. It can be observed that failure occurs along each layer plane without interlayer cleavage. This is in agreement with the experimental observation that the mechanical failure of GO nanosheets occurred by intraplanar facture of individual layers. 20 In addition, the fracture feature of each layer in the multilayered GO nanosheet is not synchronous. From 10% to 15% in the fracture strain, the outer layer of GO nanosheets rst breaks. As shown in Fig. 3(b), when stretching to the strain of 10%, the outer layer carbon (red circles) in the basal plane began to break and then the broken area increased. The obvious breakage of inner layer appeared when stretching to 12.5%. When the strain is up to 17.5%, the overall GO nanosheets fractured. This asynchronous behavior is probably due to interlayer interaction in the multilayer stacked structure, in which the inner layers possess relative stronger cooperative interactions with neighboring layers. In addition, during the tensile process, functional groups were found to reorganize and aggregate, leading to distortion in the gallery space of GO nanosheets. This has been similarly observed in the previous work. 32 In order to further explore the effect of stacking level number on the tensile behavior of GO nanosheets, a series of multilayer GO sheets with various thicknesses of GO sheets, from 1 to 5 stacking layers, denoted as GO-I, GO-II, GO-III, GO-IV, and GO-V respectively, have been simulated. Fig. 4 shows the comparison of mechanical properties for these GO nanosheets. Overall, for all the multilayer GO nanosheets, the stiffness and strength decrease with the oxidation content. The fracture strain shows certain oscillation with the oxidation content. These mechanical properties are similar with those observed in monolayer GO sheet. 17,19,30 As seen from Fig. 4(a), the Young's modulus of GO nanosheets decreases to some extent with the number of layers increasing. This indicates that enhanced number of layers could produce inhibitive effect on the Young's modulus of GO nanosheets. Comparatively, as shown in Fig. 4(b), the fracture strain and fracture strength do not show obvious changing trend with the layer number. Under the uniaxial tension, the fracture stress of the multilayer GO nanosheets has similar strength to monolayer GOs. Meanwhile, the fracture strain of GO nanosheets also maintains certain uctuating level from 11% under the lower oxidation content (20%) to 17% under the higher oxidation content (40%). This fracture strain of multilayer GO nanosheets is in good agreement with the experimental values (8-15%) by Cao et al. 20 Fig . 5 shows the corresponding comparison of mechanical properties of each sheet in the 5-layer stacked GO nanosheets (GO-V). On the whole, the Young's modulus and fracture strength for each sheet of the ve-layer GO decrease as the oxidation content increases, which are similar to those for monolayer GO. 17,30 This result means that the effect of stacking on the mechanical properties of GO is relatively weak, as compared with the oxidation content.
Further inspection of Fig. 5(a) shows that the Young's modulus of middle layer is slightly larger than that in the outer layer. Similar observation has been found for the fracture strain of individual GO sheet. This is consistent with the preceding observation that the fracture evolution process of each layer in multilayer GO nanosheets is not synchronized and the outer sheet of multilayer GO breaks rstly. This also conrmed that the fracture evolution process of multilayer GO is different from monolayer GO.
We further computed the stress-strain curves of each single layer in the 5-layer stacked GO nanosheets and analyzed the change of surface atomic structure. Fig. 6-8 show the stressstrain curves and the corresponding stress change during the tensile processes for the rst layer (layer-1), second layer (layer-2), and third layer (layer-3), respectively, in the 5-layer stacked GO (GO-V). For comparison, we also labeled the three characteristic pints (A, B, C) in the stress-strain curves, which have been used in Fig. 3(a) in the tensile process of GO-V. Fig. 6(a) shows the strain-strain curve of the layer-1, in which the fracture occurs under the overall strain of $12.3%, which is before the breaking point B of the GO-V. However, the layer-2 started to crack aer the point B ( Fig. 7(a)) with the fracture strain $13%. Furthermore, in Fig. 8(a), the stress of layer-3 reaches the rupture point when loaded to $14.8%, which is further behind the point B. As a result, we can observe that, when the layer-1 reaches the maximum stress point ($30 GPa), the stress in the layer-2 and the layer-3 continues to increase. When the layer-2 attains the fracture point ($34 GPa), the stress of layer-3 continues to increase. The fracture behavior of individual layer in the 5-layer GO is developed from the outer sheet (layer-1) to the inner sheet (layer-3), displaying delayed breakage for each layer.
The upper panel in Fig. 6(b) further shows the surface twodimensional (2-D) stress distribution of the layer-1 during the tensile process. At the whole strain of 10% for the GO-V, which is before the breakage strain, larger yellow and red regions appear in the 2-D stress diagram, suggesting the existence of higher stress on the layer-1 sheet. This is because, at this stage, there appears to be a large portion of sp 2 bonds on this plane, which could withstand applied stress loading. From the   corresponding conguration (low panel in Fig. 6(b)), it was observed that for the rst GO layer in the early tensile process, epoxide group can transform into ether group (C-O-C) by cleavage of the C-C bond associated with epoxide rings. At this stage, the skeleton structure on the basal plane of GO has not been broken yet, so the stress continues to increase. The exible ether bond allows GO to tolerate greater deformation and slows down the fracture of GO basal plane. At the strain of 12.8% of GO-V, large stress also appears in the 2-D stress distribution of layer-1. At this strain stage, the carbon structure on the layer-1 plane was destroyed with the appearance of C-C chain. This is consistent with the fact that at the strain stage, layer-1 has achieved the fracture point. At point C, corresponding to the strain of 15% of GO-V, the green color distribution becomes obvious, and the surface stress is reduced. The basal plane structure of the layer-1 is completely destroyed, corresponding to a decline in the stress-strain curve. From Fig. 6, we can nd that the fracture is rstly occurred in the region where the functional groups are located because the sp 3 bond is the weakest connection of GO structure and is the rst place to fracture. 33 Similar stress behavior has been observed for the second GO layer (Fig. 7(b)). However, the red color regions in the layer-2 are more broadly distributed than that in layer-1 ( Fig. 6(b)), showing an enhanced stress. This implies a stress-transfer from outer to inner layers between the adjacent GO layers. From the local conguration at the C point of Fig. 7(b), we found that the small crack appeared in the edge of GO, which did not show in Fig. 6(b), and it spreads to the center of the carbon basal plane slowly. In Fig. 8(b), we obverse that there appear even large stress distributions in the middle layer under the same strain condition, as compared with the rst and second layer. From the local conguration snapshots, only breaking of scattered carbon bonds and small vacancy can be formed in the third layer. Thus, under the same overall strain, the breakage of layer-3 is slower than the layer-1 and layer-2 with bearing more tough stress. The different stress distributions demonstrate the surface stress could transfer from outer GO layers to inner GO layers. This behavior has been reported in previous experimental study for the deformation and failure mechanism in GO paper. 34 In short, for the non-defective multilayer GO nanosheets, GO sheet begins to break from the weakest sp 3 C-C bond in the basal plane, and meanwhile the fracture of multilayer GO nanosheets originates from the outer layers to inner layer. The fracture process in the middle layer of multilayer GO sheets is generally delayed.

Conclusion
In this work, comprehensive MD simulations were performed to investigate the mechanical properties and fracture behavior of multilayer GO nanosheets. The Young's modulus, fracture stress, and the fracture strain were found to be 175-425 GPa, 22-41 GPa, and 10-16%, respectively. These simulated mechanical properties are in good agreement with the experimental values. We observed that the mechanical properties of multilayer GO nanosheets degrade with the increase of oxidation content except the fracture strain. With the increase of thickness, the Young's modulus of multilayer GO decreases slightly, but the fracture strain and strength are usually insensitive to the GO thickness. Meanwhile, the failure mechanism of multilayer GO nanosheets was revealed through the analysis of stress-strain behavior of individual layer in the multilayer GO nanosheets. The fracture of the outside GO layer begins rstly, and then the failure of inner GO sheet occurs with delayered process. It was shown that the middle GO layer possesses stronger toughness, as compared with the outer layer. The simulation results provide additional fracture mechanism of multilayer GO nanosheets and highlight the potential to exploit GO nanosheets as structure reinforcing materials by tuning the structures of GO nanosheets.

Conflicts of interest
There are no conicts to declare.