Large-Scale Molecular Dynamics Elucidates the Mechanics of Reinforcement in Graphene-Based Composites

Using very large-scale classical molecular dynamics, the mechanics of nano-reinforcement of graphene-based nanocomposites are examined. Simulations show that signiﬁcant quantities of large, defect-free, and predominantly ﬂat graphene ﬂakes are required for successful enhancement of materials properties in excellent agreement with experimental and proposed continuum shear-lag theories. The critical lengths for enhancement are approximately 500 nm for graphene and 300 nm and for graphene oxide (GO). The reduction of Young’s modulus in GO results in a much smaller enhancement of the composite’s Young’s modulus. The simulations reveal that the ﬂakes should be aligned and planar for optimal reinforcement. Undulations substantially degrade the enhancement of materials properties.


Introduction
To improve the mechanical performance of polymer composites, the inclusion of nanoscale 2D materials such as graphene has attracted considerable interest. [1,2]Graphene possesses exceptional mechanical properties: atomic force microscopy (AFM) indentation has shown that graphene has a Young's modulus on the order of one TPa and intrinsic strength of around 130 GPa. [3,4]OI: 10.1002/adma.202302237   Transferring these impressive properties to the composite context is key to producing polymer composites with vastly improved mechanical properties.[7] In addition, graphene under an external stress may crumple or bend, thereby reducing the transfer of stress onto the embedded graphene and furnishing little reinforcement.
Raman spectroscopy is an important tool to examine the strain of graphene flakes when embedded in a polymer matrix.The sensitivity of chemical bonds to local strain conditions results in a shifting of the Raman vibrational bands. [6,8,9]The use of Raman spectroscopy to measure the stress/strain characteristics of fillers in composite materials was pioneered by Galiotis and co-workers for fibers, [10] such as carbon and aramid. [11]They showed that Raman spectroscopy could measure the fiber strain distribution and subsequently be converted into an interfacial shear stress distribution. [12,13]Such strain distributions were also successfully determined by Raman spectroscopy for 1D fillers with a nanoscale radius, such as single-and double-walled carbon nanotubes. [14]n the study by Gong et al., [15] the authors showed that highresolution Raman spectroscopy can be used to visualize the strain field within an embedded 2D graphene sheet.Gong et al. [15] found that embedded graphene in a polymer matrix, when subjected to external tensile stress, manifests considerable transfer of stress, indicating the viability of using graphene for improving mechanical properties.The stress builds up in the graphene sheet through a process that is well described by continuummechanics-based shear-lag models where the flake deforms the most in the central region, while deforming the least at the edges, leading to high interfacial shear stress at the edge. [4,16]This shearlag process had previously been found to be the stress-transfer mechanism in 1D fibers, [17] confirming the theoretical models of mechanical enhancement going back to the 1950s. [18]The study of Gong et al. [15] showed that shear-lag models can be successfully modified for 2D fillers as well.
An important characteristic of the flake is the stress-transfer length, which can be determined from Raman spectroscopy or estimated using shear-lag models.The stress transfer length is defined as the distance from the edge of the flake to the point in the flake where the stress has reached its maximum value.The critical flake length is the minimum length of flake required to reach total stress build up and is twice the stress transfer length.
The study of Gong et al. shows that macroscopic theories of stress transfer for embedded plates are broadly applicable at the nanoscale. [15]From the strain profile, the authors derived a stress transfer length of approximately three micrometers; graphene flakes with a dimension less than this were predicted to have reduced reinforcement properties.The accurate identification of the stress transfer length is highly important as typical scalable techniques to create graphene polymer composites, such as liquid shear exfoliation of graphite, produce small few-layer graphene flakes, which may be only a few micrometers in lateral dimensions, and therefore comparable to the predicted transfer length of graphene. [19,20]However, the resolution of the technique meant that the strain distribution near the edges of the flake, i.e., less than 1 μm from the edge, was not fully resolved. [21]n recent studies, it was found that the strain field at the edge deviated from classical shear-lag theory, which was attributed to unintentional chemical doping and edge effects. [22]The authors proposed two domains for the graphene flake: near the edges of the flake (<2 μm) the stress transfer is reduced, while in the central region stress transfer occurs due to an elastic shear-lag processes.Recently, Manikas et al. [23] examined the stress transfer of a well-defined graphene microribbon aligned with the loading direction using Raman spectroscopy; they determined that the transfer lengths in flakes of monolayer graphene depend on the applied strain and are lower than previously considered, within a range from 500 to 1000 nm.
There is, therefore, a need to unambiguously determine the transfer length for graphene by understanding the stress transfer mechanism at the atomistic level.To do so, we have used very large-scale molecular dynamics simulations to establish how the efficiency of stress-transfer depends upon the size of the graphene sheet, the interfacial attraction with the polymer matrix and the extent of planarity of the sheet.In this work, we have systematically varied the size of the reinforcing sheet up to a maximum graphene sheet length of 2 μm and examined the subsequent reinforcement of the composite.To examine the effect of altering the interfacial attraction between the polymer and the flake, we have also considered graphene oxide (GO) flakes.MD simulations allow us to visualize the strain field within the graphene/GO sheet.Both the reinforcement as a function of sheet size and the strain profile along the sheet can be directly compared to shear-lag models, allowing us to determine the validity of the model for the crucial micrometer and sub-micrometer sized flakes and unambiguously determine the stress transfer length for monolayer graphene.
Determining the validity of shear-lag models is important as these models do not consider edge effects or the interphase region where the polymer density on the surface of the sheet is greater than in the bulk and causes considerably different behavior to that of the bulk polymer.Similarly, the effect of wrinkles or undulations of the graphene flake is ignored in shear-lag models, but will emerge in molecular dynamics simulation due to thermal motion of the atoms.

Results and Discussion
Our graphene-nanocomposite model consists of a single narrow graphene nanoribbon in a poly(vinyl alcohol) (PVA) polymer matrix (Figure 1).We consider both graphene andGO composites with nanoribbons of lengths ranging from 100 to 2000 nm, with the orthogonal direction kept constant at 10 nm.The longest nanoribbon (2 μm) in our study is comparable to the experimental study of May et al. (2.3 μm).The sheets are discontinuous, that is they do not cross periodic boundaries in the simulation cell.The volume fraction is kept constant at below 2%, consistent with graphene-polymer composites, with a polymer thickness of approximately 10 nm.An illustration of the model setup is shown in Figure 1a and a list of models is given in Table S1 (Supporting Information).The ribbons were chosen to have zigzag edges.The distance between the edges in the x direction of periodic replicas is 20% of the length of the flake, except for the 2000 nm flake systems, where it is initially set at 200 nm.The size of the systems and their atom count are given in Table S1 (Supporting Information).It should be emphasized that our models are fully atomistic and therefore chemically specific for the graphene/GO-PVA composite.The largest models contain approximately 30 million atoms.It is clear that sampling equilibria with MD requires ensembles of simulations, [24] which leads to increased reliability, reproducibility, and a tighter control of standard uncertainty.Therefore, we created an ensemble for each system containing seven replicas.This number was chosen as it is a reasonable point at which there are diminishing returns concerning the confidence in the average Young's modulus for a test system (300 nm graphene flake), see Figure S1 (Supporting Information).The replicas were given a different random seed for the generation of atomic velocities drawn from a Maxwell-Boltzmann distribution.After equilibration by simulated annealing and cooling to a glassy state, the composites were subjected to uniaxial tensile strain in the x direction (the direction of the sheets) at a constant strain rate (1 × 10 8 s −1 ) with a zero pressure condition imposed on the lateral faces (see Computational Methods section and Supporting Information for more details).Although this strain rate is much higher than laboratory experiments due to the timescales possible using molecular dynamics, for very small strains (0.2%) in the elastic regime, the results should be comparable with experiment.
PVA was chosen as the matrix polymer to directly compare with the experimental results of May et al., where the authors separated ultrasonication-exfoliated graphene flakes into two different sizes by selective centrifugation [25] before dispersing in PVA.Significant reinforcement of approximately 70% of the theoretical maximum was observed for flakes of mean size 2.3 μm, while much lower reinforcement was observed for flakes of mean size 1.1 μm.This very high level of reinforcement observed by May et al. for the larger flake size demonstrates there is good interfacial stress transfer between PVA and graphene, and that PVA is an excellent model polymer with which to study graphene nanocomposites.
In Figure 2, we show the Young's modulus for our simulated composites of graphene monolayer flakes embedded in a PVA matrix (E c ) as a function of flake size, calculated using a linear regression of the computed stress-strain curve between strains of 0 to 0.2% averaged over the replicas in each ensemble.We find The initial lattice dimension in the z direction is 10 nm.c) Top-down view of the graphene/GO ribbon-PVA composite.The simulation supercell lattice dimension is initially 20 nm in the y direction.The ribbon is 10 nm wide, leading to approximately 10 nm between the edges of periodic replicas in the y direction.The ribbon is aligned in the x direction with lengths varying from 100 to 2000 nm.d) The edges of the ribbon are terminated by hydrogen atoms.Our simulations contain ribbons with zigzag edges along the length of the ribbon.e) The shear-lag model for deformation under uniaxial strain/stress: the deformation in the stiffer graphene/GO flake is greatest in the central region, while deforming the least at the edges.
that the Young's modulus is over double that of bulk PVA for the smallest 100 nm flake (7.5 GPa, with bulk PVA 3 GPa) before sharply rising between 100-400 nm.After 400 nm, the rate of increase in E c decreases markedly, with E c of the 500 nm flake (15.5 GPa) being comparable to that of the largest flake size in our study (2000 nm) of 17 GPa.This latter E c is over five times higher than that of pure PVA (3 GPa).Our simulations therefore identify a threshold for the size of a pristine monolayer graphene flake where stress-transfer to the stiffer flake becomes efficient and reinforcement occurs.
In the experimental study of May et al., [25] the authors observed highly efficient reinforcement for their largest graphene flake (2.3 μm), confirming the results of our simulations.The authors also observed a significant increase in Young's modulus when L/t was increased from 1000 to 2000, where L is the length of the flake and t is the thickness (assumed to be 1 nm from their TEM studies) for graphene flakes embedded in PVA, [25] indicating that efficient stress transfer occurs when the flake length is greater than 1 μm, which is longer than that found in our simulations.It should be noted, however, that the critical length has been shown to be dependent on the number of layers of graphene in the flakes. [26]1 nm thickness corresponds to few-layer graphene of approximately three flakes, while our simulations contain monolayer graphene.
We can also compare with critical lengths derived from the Raman spectra of large graphene flakes, defined as the distance from the edge to the plateau region in the strain profile.Our observed critical length for high reinforcement efficiency (400-500 nm) is similar to the recent studies of Manikas et al. [23] who used shifts in Raman spectra to estimate the transfer length of highly aligned graphene ribbons of 70 μm length embedded in a poly(methyl methacrylate) (PMMA) matrix.The transfer length Figure 2. a) Illustration of the simulated unaxial deformation of the composite: the lattice dimension of the simulation supercell in the x direction is increased at every timestep, with the atomic coordinates correspondingly scaled.The strain is ((L x ′ − L x )∕L x , where L x ′ and L x are the current lattice and initial lattice dimension in the x direction, respectively.The Young's modulus is calculated as the gradient of the stress response with strain between strain of 0-0.2%.b) The computed Young's modulus for the graphene-composites as a function of ribbon length (blue) and the comparison with shearlag theory (red).The corresponding length efficiency factor is shown in the inset.c) The computed Young's modulus for GO-composites as a function of ribbon length (blue).The shear-lag predicted values are also shown for different thickness of the GO flake (and the corresponding volume fraction).In both cases, the computed Young's modulus is below that predicted by the shear-lag model.
as found to depend on the applied strain and was in the range from 500 to 1000 nm, thereby indicating critical length sizes of 1-2 μm, depending on the strain.
Although the study of Manikas et al. used PMMA rather than PVA, as the interaction between both polymers and graphene is purely via van der Waals interactions, we would expect the interfacial energy and the stress transfer lengths to be comparable.
We can compare our results to the predictions of shear-lag theory, which has been shown to successfully capture the mechanics of stress transfer from the polymer matrix to the reinforcing graphene flakes. [4,8,25]Here we follow the approach of May et al., [25] where we first consider the effect of uniform strain of a perfectly aligned graphene flake by using the simple "rule of mixtures" (ROM), modified by shear-lag considerations.The modulus of a composite E c is given by the ROM equation as: where E f is the Young's modulus of the graphene filler, V f is the volume fraction of the graphene and E m is the Young's modulus of the polymer matrix.In our simulations, the volume fraction of graphene flakes is approximately 1.5%, and assuming a Young's modulus of graphene and PVA of 1 TPa [3] and 3 GPa respectively (the latter from our simulations), gives a ROM result of 17.95 GPa.We can therefore see in Figure 2 that the Young's modulus of the composite containing our largest 2000 nm-long flake is close to this value (17 GPa), illustrating that graphene flakes, when perfectly aligned and of size greater than the critical length, can effectively reinforce the polymer matrix to close to its maximum ROM value.Following May et al., [25] detailed analysis of shear-lag theory produces a modified ROM that includes a length efficiency factor reflecting the dependence of reinforcement on flake length and increases from 0 to 1 with increasing length.
In Figure 2b, we have converted our E c results into the length efficiency factor  le .As shown by May et al., a detailed analysis based on shear-lag theory finds that  le is given by ref. [25]: where L is the length of the graphene flake, t is its thickness, and G m is the shear modulus of the polymer, which we assume is 1GPa. [25]Using a thickness of 0.335 nm for monolayer graphene, we can compute the  le for our composite systems using Equation (3) and compare them to our results from simulation, shown in Figure 2b.As we can see, the agreement is very good, indicating that the shear-lag model is an excellent approximation for The resulting strain within the 500 and 2000 nm graphene flakes, with an imposed strain of 0.2% on the composite.All profiles show shear-lag behavior, with a strain-plateau in the center.The differences in the x component of the bonds matches or exceeds the global strain, while the strain in bond lengths is less, indicating a corresponding reduction in the y and z components of the bond.The stress transfer length, i.e., the distance to reach maximum strain, is longer for the 2000 nm flake (≈500 nm) than for the 500 nm flake (≈250 nm).
the behavior of pristine monolayer graphene.This agreement shows that the interphase region does not need to be considered for aligned monolayer graphene.The difference in critical length between our results and those determined experimentally by May et al. is also captured by shear-lag models due to the difference in the thickness of the flakes; for our simulations, the sheets are unambigously monolayers with a thickness of 0.335 nm.Experimentally, the flakes are multilayer (≈1 nm, corresponding to three layers), [25] which significantly extends the critical flake length in Equation ( 3).Molecular simulation gives information on the behavior of individual atoms.Using atom positions during the stress-strain simulations, we can analyze the strain field that evolves when the composite is under uniaxial tension.This can be directly compared with the strain fields determined experimentally by the shift in Raman vibration bands, such as those of Manikas et al. [23] In Figure 3 we show the local strain across the graphene sheet, computed via the change in bond distances for carbon bonds that predominantly point in the axial direction of strain, divided by their original length (Figure 3a). Figure 3b shows the strain profile for our largest flake.We see the local strain increasing from zero at the edges to a plateau region at about 500 nm, where the local strain is approximately constant, before decreasing again at 1500 nm (500 nm from the edge), again to zero at the far edge.This behavior is expected of a flake undergoing shear-lag behavior where the stress is zero at the edge and increases to a maximum value.It is also consistent with the experimental results of Manikas et al. [23] The strain field was calculated by taking the snapshot at an applied strain of 0.2% and continuing the simulation at this fixed strain state, calculating the difference in bond lengths and the corresponding x component of the bond compared to the zero strain state.The local strain computed from the bond lengths is approximately 0.11% in the plateau region, which is less than the global strain, while the local strain computed from the change in the x component matches or even exceeds the global strain of 0.2%.We conclude that, under unaxial strain, the bonds extend in the strain direction and also become more aligned with the strain direction.The changes in carboncarbon bond distances and angles within the sheet, combined with their spring constants, determines the elastic constants of the graphene sheet. [27]In Figure 3b, we show the strain field for the 500 nm flake.This flake also reaches a plateau region with an x bond component strain matching that of the global strain.Interestingly, the stress transfer length is shorter than for the larger flake (≈250 nm), indicating that the stress transfer length is length dependent and smaller flakes may be more efficient than estimated using stress transfer lengths derived from larger flakes.
We can also estimate the strain profile using shear-lag theory. [15]For a given level of strain applied to the polymer matrix, ɛ m , the variation of strain in the graphene flake, ɛ f , with position x, where x = 0 corresponds to the centre of the flake, will be of the form: where and G m is the matrix shear modulus, E f is the Young's modulus of the graphene flake, l is the length of the graphene flake in the x (axial) direction, t is the thickness of the graphene, T is the total thickness (which is equivalent to the distance between periodic replicas in the z direction in our simulations, L z ) and s is the aspect ratio of the graphene (l/t) in the x direction. [15]Equation (5)  shows that the strain within the flake increases very quickly under the applied strain (to ≈91% by 100 nm), unlike the slow increase we observe in Figure 3.Our simulations therefore suggest that the efficiency of strain transfer is less than that predicted by shear-lag theory.This needs to be considered when designing composite materials.Tension in the longitudinal direction has been shown to cause orthogonal buckling to occur in the transverse direction due to the Poisson contraction. [28,29]Orthogonal buckling can weaken the polymer-graphene interface, which will affect the stresstransfer length.[32] The high value for graphene embedded in polymer is to be expected, as the polymer will resist out of plane deformations, assuming there is no significant roughness in the polymer surface.
Our model contains a perfectly aligned graphene flake of regular rectangular geometry.Manikas et al. showed experimentally when such a graphene flake undergoes uniaxial strain, no orthogonal buckling occurs. [23]In Figure S7 (Supporting Information), we show the average amplitude of undulations in the transverse direction to uniaxial strain direction at zero strain and 0.2% strain.There is almost no detectable change in the out-of-plane amplitudes, indicating that no lateral buckling has occurred.[32] At 0.2% strain, we find a lateral compressive strain in the polymer matrix of −0.06%, significantly below the critical buckling strain.The strain within the flake calculated in the y direction shows that there is no shearlag behavior, as expected due to the 10 nm width of the flake, and compression within the graphene sheet is almost entirely in-plane (Figure S8, Supporting Information).We can therefore conclude that at low strains (<0.2%), there is no lateral buckling and consequently no weakening of the polymer-graphene interface, despite the very high aspect ratio of our graphene flake.
Molecular simulation allows us to examine the behavior of the polymer-graphene interface.We have calculated the density of polymer atoms perpendicular to the graphene surface.In Figure S4 (Supporting Information), we observe a layer close to the surface at a distance of 3.8 Å.The peak in this density profile is unchanged when the flake is extended from zero strain to 0.2% (Figure S6, Supporting Information), again indicating that there is no weakening of the polymer-graphene interfacial bonding at a strain of 0.2%.
In Figure 4, we show the behavior of the graphene flakes when the strain is fixed.We observe that the stress response does not relax to a converged value as expected, but instead oscillates around a constant value.The maxima and minima of the stress response correspond to a vibrational mode that increases and decreases the length of the flake in the axial direction.This is the longitudinal vibrational mode shown in Figure 4a.The perfect correlation between the stress response and the length of the sheet (i.e., the distance in the x direction) as a function of time shown Figure 4b, indicating that this behavior is due to longitudinal vibrations.This vibrational mode has recently been observed as a low-energy mode in Raman spectra for armchair graphene nanoribbons. [33]ere, we also observe the vibrational mode for zigzag graphene nanoribbons.The period of vibration is dependent on the length of the flake and can significantly modulate the stress response of the flake.We can see in Figure 4c that the period of vibration is proportional to the length of the flake, which was also observed by Overbeck et al. [33] We have also addressed the role of surface functionalization through studying the reinforcement behavior of GO.The oxygen containing functional groups on GO will increase the interaction with the polar polymer matrix, which should allow for better stress transfer and reduce the critical length.A visualization of the GO surface is shown in Figure S3 (Supporting Information).Nevertheless, the functional groups tend to damage the graphene lattice, thereby reducing the stiffness; the Young's modulus of GO is approximately 300 GPa, compared to one TPa for pristine graphene.The thickness of a GO flake is considerably more than that of pristine graphene, with an interlayer separation of 0.8 nm, [34] compared to 0.35nm for graphene.Our simulations are able to determine whether these factors increase or decrease the reinforcement.To illustrate the greater interfacial interaction between the polymer and GO compared to graphene, in Figure S5 (Supporting Information) we show the radial distribution functions between hydroxyl and epoxy groups on the GO surface and hydroxyl groups on the polymer.Figure S5 (Supporting Information) shows that all functional groups on the GO surface are fully participating in hydrogen bonding with the polymer.
In Figure 2c we show the Young's modulus of our GO reinforced composites as a function of GO flake size.We see that the critical length for enhancement of E c is shorter than for graphene (≈300 nm), as previously observed experimentally for GO flakes embedded in PMMA. [35]However, the composite Young's modulus is much reduced compared to that of graphene (≈10 GPa vs 17 GPa).Using a thickness for a GO flake of 0.8 nm, a correspondingly larger volume fraction of 3.5%, and a Young's modulus of 300 GPa, we can compare to the predictions of shear-lag modified ROM equation, shown in Figure 2c.However, the thickness of a GO sheet is not uniform; for example, graphitic regions will have a thickness of 0.335 nm.If we consider the thickness of a GO sheet to be 0.6 nm, the shear-lag predictions are decreased compared to a thickness of 0.8 nm, but it is still greater than our simulated Young's moduli.It is clear that, irrespective of our choice of GO thickness, the efficiency of stress transfer is reduced in our simulations compared to that predicted by shearlag models, unlike graphene flakes.This discrepancy may be due to greater undulations of the GO flake which reduce the Young's modulus of the GO sheet.
To further investigate the effect of undulations on the reinforcement efficiency of graphene and GO, we have conducted stress-strain simulations on undulating flakes.To create an undulatory flake, we equilibrated our systems while allowing the graphene flake to move when the polymer was a low density (≈0.7 g mL −1 ), rather than keeping the flake fixed.This resulted in flakes with many more undulations, as shown in a snapshot from simulation in Figure 4d and the corresponding Fourier transform of the height function of the surface, shown in Figure 4e.The stress response of the 500 nm flake (green) with a fixed global strain of 0.2% with the corresponding length of the flake (black) lines.The stress response is clearly shown to be modulated by the longitudinal vibrational mode, with greater stress when the flake is longer.The period of the vibrational mode is approximately 0.16 ns.c) Linear relationship between the length of the flake and the period of the longitudinal vibrational mode.d-f) Illustration of the reduction of the Young's modulus of the composite with undulations of the flake.d) A snapshot from simulation of a portion of a graphene flake allowed to undulate by equilibriating in an initially low-density polymer matrix.e) Fourier transform of the height of the undulatory graphene (black) and planar graphene (red).The undulatory graphene shows much greater intensity for both long and short wavelength undulations.f) Corresponding stress-strain behavior of the undulatory and planar graphene-PVA composites.The Young's modulus is reduced from approximately 17 GPa to 6.5 GPa.
The Fourier transform shows increased amplitude for short wavelengths (10-60 nm) and at long wavelengths (500-2000 nm).The amplitude of the largest undulations is approximately 2 Å, less than the 1 nm intrinsic wrinkles observed for suspended monolayer graphene caused by thermal fluctuations, [36] but comparable to that of previous atomistic simulations of suspended graphene (1 Å). [37] It should be noted that the flake still lies predominantly in the xy plane.In Figure S4 (Supporting Information) we compare the density profile of polymer carbon atoms perpendicular to the graphene surface for planar and undulatory graphene.We see there is very little difference in the two profiles and the total amount of polymer atoms in the first peak closest to the graphene sheet is unchanged.In Figure 4f we find that the Young's modulus is significantly reduced for the undulatory flake compared to the flatter flake (6.5 ± 1.5 GPa as opposed to 16.8 ± 1.4 GPa).Similar reductions are seen for GO composites (4.9 ± 0.2 GPa as opposed to 10.6 ± 0.2 GPa).This is due to the much reduced Young's modulus of undulatory graphene/GO for small strain values, where the strain merely dampens some of the undulations in the graphene flake rather than activating the in-plane stiffness.This is a significant observation as it shows that undulations of the graphene flake results in a substantial reduction in the Young's modulus of the composite, not due to any debonding or delamination of the surface and subsequent reduction in interfacial stress transfer, but from an intrinsic reduction of the Young's modulus of the graphene flake from undulations.

Conclusion
Using very large MD simulations, we have simulated the elastic response of graphene/GO-PVA nanocomposites.We have shown that shear-lag models provide a very good description of the behavior of graphene flakes, accurately predicting the Young's modulus and the transfer length for monolayer graphene perfectly aligned with the direction of strain.The only inputs to the shear-lag model are the length, thickness, and Young's modulus of graphene and the Young's and shear moduli of the polymer, indicating that the presence of a interphase region of high polymer density near the graphene surface does not affect the efficiency of reinforcement.This important result shows that macroscopic theories can be successfully used for reinforcing graphene flakes, even down to very small length scales (<500 nm).Our simulations show that perfectly aligned graphene sheets with a dimension of 2000 nm provide reinforcement close to that predicted by the ROM.It is clear that this is an upper bound and there are several reasons why even pristine, highly aligned monolayer graphene would be less efficient at reinforcement.Even relatively small undulations, with amplitudes of less than 1 nm, can significantly reduce the reinforcement.Small strains only serve to dampen the undulations in the sheet.This is especially the case for GO, where the oxidized regions create very flexible sheets.As a result, our MD simulations show that even highly aligned GO has a lower Young's modulus than that predicted by shear-lag models.We have also shown that longitudinal vibration modes can also have a significant effect on the stress response of graphene flakes.
Summarizing, we have shown that while flakes over 500 nm in length can significantly increase the elastic properties through a shear-lag process, the most important factor determining the efficiency of stress transfer is the undulations of the flake.Processing conditions such as melt compounding or solvent processing often produce flakes with bent, wrinkled, or crumpled shapes, caused by stresses or defects. [2,38]Reducing the undulations may be achieved by controlling the residual strain of the embedded graphene flake, as highly strained graphene flakes will have dampened undulations.However, large uniaxial strains may cause buckling in the transverse direction due to the Poisson contraction, thereby reducing any potential reinforcement.These phenomena will be the subject of future investigations.It is clear that, if graphene is to reach its full potential for creating exceptionally stiff and strong materials, research effort should be directed toward producing planar graphene/GO fillers in composites.

Computational Methods
As reported in previous studies investigating the structures formed by graphene and GO flakes in polymeric media [39] and aqueous environments, [40] molecular simulation is a powerful method by which to study these systems because the composition and structure of graphene and GO can be controlled.To reach length scales on the order of micrometers of a embedded graphene sheet with atomistic resolution requires a very thin graphene ribbon to be tractable.Therefore, systems were created with one dimension that varies from 100 to 2000 nm, with the orthogonal direction kept constant at 10 nm.The sheets are discontinuous, i.e. they do not cross periodic boundaries in the simulation cell and are surrounded by polymer molecules.The volume fraction was kept constant at below 2%, consistent with graphene-polymer composites, with a polymer thickness of approximately 10 nm (Figure 1).
Model Description: The graphene/GO systems was created by following the procedure.A rectangular graphene or GO sheet was created with dimensions y = 10 nm, x = 100 nm, 200 nm, 300 nm, 400 nm, 500 nm, and 2000 nm.A corresponding PVA polymer system was created by first constructing a small polymer system containing a dummy graphene sheet.The simulation box was 20 nm × 20 nm × 20 nm.The dummy sheet was in the center of the simulation box lying the xy plane, with dimensions 20 nm × 10 nm (i.e., periodic in the x direction, but containing edges in the y direction).Around the sheet was filled with PVA polymer molecules of 100 monomer units with a density of 0.7 g mL −1 .To add the polymer to the simulation box we used a Monte Carlo procedure to produce a low energy state, drawing on the methods of Theodorou and Suter [41] and the lookahead procedure developed by Meirovitch. [42]The details of this method have been previously published. [43]This method generated an amorphous polymer system without high energy overlaps.The Lennard-Jones parameters of the dummy graphene sheet were altered to repel the polymer molecules significantly further away than graphene ( = 10 Å).This system was subsequently relaxed by molecular dynamics at 500 K for 1 ns, before reducing the temperature to 300 K over 2 ns.The dummy graphene sheet was not allowed to move.
The dummy graphene sheet was then removed, and the system replicated in the x direction to the size of the corresponding graphene/GO sheet plus 20% of the length of the flake (10% for the 2000 nm flakes).The void created by the removal of the dummy graphene sheet was then filled by the graphene/GO sheet, that is the sheet possesses edges in the x direction (see Figure 1 for a schematic diagram).The system was then relaxed, with the polymer molecules filling the voids at the ends of the graphene/GO sheets.The largest simulations in this study, containing graphene/GO sheets of 2 μm in length, contain approximately 30 million atoms.A list of the systems studied, including the number of atoms, the initial and final lattice dimensions is given in Table S1 (Supporting Information).
To reach thermodynamic equilibrium, a simulated annealing approach was used to overcome energy barriers.The graphene/GO sheets were required to remain as much in-plane as possible, while at the same time allowing the sheet to relax due to interactions with polymer.The simulated annealing approach therefore concentrated on the polymer molecules, while leaving the sheets fixed, before relaxing the whole system at 300 K. Initially, the simulations were run at 300 K for 0.5 ns at 300 K, with only polymer molecules allowed to move for the first 0.25 ns and both polymer and sheet allowed to move for final 0.25 ns.The final 0.25 ns were computed in a NpT ensemble, allowing the simulation box to relax.The simulation box was then fixed (i.e., a NVT ensemble), and the simulation was run with only polymer molecules allowed to move at 900 K, for 0.3 ns.The system was cooled down to 300 K over 0.2 ns.The local density across the simulation cell was checked in the xz plane to ensure there were no voids, and that the density of polymer across the simulation box was constant.The simulation box was then allowed to relax for 0.2 ns at 300 K, allowing both polymer and sheet to move.The potential energy and volume evolution was checked to ensure any drift was small.Any drift in potential energy observed was less than 0.3% in the final 0.1 ns of simulation at 300 K. To simulate the stress-strain behavior, the simulation was quenched down to 200 K, to ensure it was below the glass transition temperature of the polymer.This reduction in temperature was performed over 0.2 ns, with all molecules allowed to move in a NpT ensemble.To create the undulatory graphene flakes, the flakes were not fixed in the initial equilibriation at low polymer density, but were allowed to move.The flakes were subsequently fixed for the simulated annealing part of the equilibriation process described above.No bond breaking or formation is possible in classical atomistic molecular dynamics simulations, therefore no reactions were simulated at higher temperatures; instead configurational space of the polymers was explored before cooling down to 300 K. To simulate the stress-strain behavior, we took the final timestep from the simulations at 200 K. Uniaxial tensile strain simulations in the x direction were then performed at a constant strain rate (1 × 10 8 s −1 ) with a zero pressure condition imposed on the lateral faces.The stress components were determined from the pressure tensor, calculated via the virial stress, to give the stress-strain behavior.We therefore plot the instantaneous stress while the strain increases.The Young's modulus was calculated by a linear regression of the stress response between strains of 0 to 0.2%.The stress-strain profiles for graphene and GO are shown in Figures S9 and S10 (Supporting Information), respectively.
Graphene Oxide Models: The GO systems have functional groups on the basal plane and the free edges.Graphene edges were terminated with hydrogen atoms.All of the GO models contain carboxylic acid groups on the edges of the flakes.These groups were protonated.GO contains equal numbers of epoxy and hydroxyl groups on the basal surface.A schematic diagram of the atomic structure of the GO surface is shown in Figure S2 (Supporting Information).
To generate the distribution of oxidation sites on the basal surfaces of the GO flake, a recently developed algorithm was used, which employs quantum mechanical simulations of reactive intermediates [44] to determine the progress of oxidation on the flake surface. [45]The algorithm creates GO structures with distinct oxidized and graphene domains, as visualized in experiments.Figure S3 (Supporting Information) shows the snapshots from our simulations, which illustrate the distinct domains of sp 2 and sp 3 phases as created by the GO model builder.The GO models created here had a carbon oxygen ratio of 2.5, corresponding to GO commonly produced using the modified Hummers method. [46]nsembles and Averages: Ensembles of simulations were required with MD simulations, leading to increased reliability, reproducibility and estimates of uncertainty.[24] The optimal ensemble size was evaluated by measuring confidence in our quantity of interest (the Young's modulus).In the Supporting Information, the 95% confidence interval for predictions of the Young's modulus for the graphene 300 nm flake system are shown.There is little decrease in the confidence interval with ensembles greater than seven (Figure S1, Supporting Information).8] All potentials used the OPLS forcefield parameter set, [49] which had previously been successfully used to capture the interactions between graphene and graphene oxide with solvent and polymer.[50][51][52][53] PVA is also well represented by OPLS, exhibiting good agreement with experimental values such the glass-transition temperature and radial distribution functions, illustrating the ability of OPLS to capture the hydrogen-bonding of PVA.[54] A list of the OPLS parameters is given in the Supporting Information.

Figure 1 .
Figure 1.a) 3D visualization of the simulation setup illustrating the graphene/GO ribbon orientation (dark green) in the PVA polymer matrix.Periodic boundary conditions are shown as blue lines.The ribbon lies in the xy plane, does not cross any periodic boundaries and is completely surrounded by polymer molecules.b) Side view of the graphene/GO ribbon-PVA composite (xz plane).The initial lattice dimension in the z direction is 10 nm.c) Top-down view of the graphene/GO ribbon-PVA composite.The simulation supercell lattice dimension is initially 20 nm in the y direction.The ribbon is 10 nm wide, leading to approximately 10 nm between the edges of periodic replicas in the y direction.The ribbon is aligned in the x direction with lengths varying from 100 to 2000 nm.d) The edges of the ribbon are terminated by hydrogen atoms.Our simulations contain ribbons with zigzag edges along the length of the ribbon.e) The shear-lag model for deformation under uniaxial strain/stress: the deformation in the stiffer graphene/GO flake is greatest in the central region, while deforming the least at the edges.

Figure 3 .
Figure3.a) The strain field within the flake is calculated within the graphene flake by considering changes in the bond distances between the nonstrained flake (dark grey lines) and the strained atom positions (green in the ball-and-stick representation).Only bonds which are predominantly in the x direction are considered (for example, the bond highlighted in (a)).Differences in bond distance (r ′ − r) and the difference in the x component of the bond are calculated, and converted to a local bond strain by dividing by r and the x component of the unstrained bond, respectively.The local bond strains are binned along the x direction of the flake into 300 strips.b) The resulting strain within the 500 and 2000 nm graphene flakes, with an imposed strain of 0.2% on the composite.All profiles show shear-lag behavior, with a strain-plateau in the center.The differences in the x component of the bonds matches or exceeds the global strain, while the strain in bond lengths is less, indicating a corresponding reduction in the y and z components of the bond.The stress transfer length, i.e., the distance to reach maximum strain, is longer for the 2000 nm flake (≈500 nm) than for the 500 nm flake (≈250 nm).

Figure 4 .
Figure 4. a) An illustration of the atomic displacement involved in the longitudinal vibrational mode of a zigzag edge graphene ribbon.This mode has been observed as a length-dependent low-energy peak in Raman spectra.b)The stress response of the 500 nm flake (green) with a fixed global strain of 0.2% with the corresponding length of the flake (black) lines.The stress response is clearly shown to be modulated by the longitudinal vibrational mode, with greater stress when the flake is longer.The period of the vibrational mode is approximately 0.16 ns.c) Linear relationship between the length of the flake and the period of the longitudinal vibrational mode.d-f) Illustration of the reduction of the Young's modulus of the composite with undulations of the flake.d) A snapshot from simulation of a portion of a graphene flake allowed to undulate by equilibriating in an initially low-density polymer matrix.e) Fourier transform of the height of the undulatory graphene (black) and planar graphene (red).The undulatory graphene shows much greater intensity for both long and short wavelength undulations.f) Corresponding stress-strain behavior of the undulatory and planar graphene-PVA composites.The Young's modulus is reduced from approximately 17 GPa to 6.5 GPa.