A molecular dynamics study on the interaction between epoxy and functionalized graphene sheets

The interaction between graphene and epoxy resin was studied using molecular dynamics simulations. The interfacial shear strength and pull out force were calculated for functionalised graphene layers (carboxyl, carbonyl, and hydroxyl) and epoxy composites interfaces. The influence of functional groups, as well as their distribution and coverage density on the graphene sheets were also analysed through the determination of the Young's modulus. Functionalisation proved to be detrimental to the mechanical properties, nonetheless according to interfacial studies the interaction between graphene and epoxy resin increases.


Introduction
Engineering structures are increasingly comprised of carbon fibre reinforced polymer composites due to their low weight to stiffness ratio, resistance to environmental impact, and processability. Understanding the interaction between the fibres and the polymer matrix at different length scales is crucial for developing improved carbon fibre reinforced polymers. At the molecular level carbon fibres consist of stacked layers of graphene. This is a two dimensional layer with carbon-carbon sp 2 hybridized covalent bonds that result in outstanding mechanical properties. However, the crystalline and smooth surface of the graphene hinders the wetting by the polymer matrix. Recent reports suggest that the covalent bonding of functional groups to graphene can overcome this problem [1]. Nevertheless, it has an impact on the mechanical properties of the graphene due to a change in hybridization and introduction of defects. The chemisorption changes the hybridization of the grafted carbon atom from sp 2 to sp 3 resulting in an elongation of the bonds between the sp 3 carbon and its neighbouring sp 2 hybridized atoms. This elongation shortens the angles between them, and the grafted carbon buckles off plane in course of pyramidalization [2]. The mechanical properties of graphene and in particular the influence of functional groups have been addressed by several authors applying Molecular Dynamics (MD), in particular grafting density, orientation and/or location of hydroxyl, methyl, and carboxyl groups [3][4][5].
On the other hand the adhesion between polymer and graphene is related to the chemistry of both surfaces. As mentioned previously, by covalently bonding of chemical groups on the graphene the interaction between graphene and the polymer matrix can be increased improving the efficiency of load transfer [6]. Hence, a compromise needs to be made between the increase in load transfer capability and the decrease of the mechanical properties of the graphene due to functionalisation. MD simulations are advantageous in addressing nanoscale interfacial issues as it calculates microscopic interaction between molecules that can be scaled up to the macro level through statistical mechanics, obtaining predictions of bulk properties. Interfacial properties of CNT and polymer composites have been examined through MD by applying the pull out method [4,[7][8][9]. More recently graphite and graphene have been added to this circle as reported by Jang et al. [10] using vinyl ester resin as matrix, and Liu et al. [11] who compared two different matrices, polyethylene and poly(methyl methacrylate).
In this work we employed MD to calculate the interfacial shear strength (ISS) between an epoxy resin and functionalised graphene layers. Each layer has a single type of chemical group -carboxyl, carbonyl, and hydroxyl -and are dispersed on the graphene with a 10 % coverage density. Their mechanical and structural properties, the influence of coverage density (3.0 and 5.0 %) and placement of hydroxyl groups (10 % grafting as a large single cluster and as smaller clusters) were also assessed.

Molecular models
2.1.1. Epoxy. Epoxy resin was built using the 3D atomistic tool from Materials Studio 6.0 [12]. The polymer was based on a 50 % crosslinked bisphenol A diglycidyl ether (DGEBA) with triethylenetetramine (TETA) consisting of 744 atoms. The crosslinked molecule was replicated resulting in a periodic cell with 70 packed molecules and 1 configuration which originated 12 250 atoms. Cell dimensions are 5.2 x 4.0 x 4.7 nm and cell angles α = β = γ = 90˚. The density was set to 1.15 g.cm -3 , an average value for epoxy resin density as reported in literature [13]. The final structure has a density of 2.24 g.cm -3 . Functional groups (COOH, C=O, and OH) were individually and randomly grafted onto the 96 carbon atoms from the graphene layer resulting in 10 % functionalization.

Simulation Setup
The COMPASS force field [14] was used to model the interactions between atoms. In calculating the non-bonded interactions, van der Waals and electrostatic interactions are described by Lennard Jones 6-9 and Coulombic functions, respectively. The atom-based summation method with cut-off radii of 9 Å and long range corrections was used to calculate non-bond interactions. Dynamics were run with time step of 1 fs and at 298K in every case, unless stated otherwise. To equilibrate the epoxy system successive MD simulations were run. First, dynamic simulations with the canonical ensemble (NVT) and Andersen thermostat were conducted for 70 ps before annealing. Next, simulated annealing was applied in order to gradually minimise the energy of the structure to prevent it from being trapped in a conformation that represents a local minimum. Thus higher temperatures allow the simulation to overcome energy barriers to move into other low energy areas hence, the temperature is periodically increased from an initial temperature to a mid-cycle temperature and back. At the end of each complete temperature cycle, the lowest energy structure produced by that cycle is output to the trajectory file. The annealing schedule was performed with an initial temperature (300 K), a mid-cycle temperature (600 K), and four heating ramps per cycle and with 500 simulation steps in each dynamics step. The annealing was performed in an isothermal-isobaric ensemble (NPT) with Andersen thermostat and Parrinello-Rahman barostat, respectively. This allowed the shape and volume of the cell to change during annealing. Following annealing a final equilibration dynamic simulation was performed for 110 ps using the NPT ensemble with Andersen thermostat and Parrinello-Rahman barostat. After equilibrium the angles of the simulation cell were corrected again to 90˚ as they have been changed during NPT. This step is necessary in order to combine the epoxy system with the graphene molecular system. For graphene systems dynamics were performed with a NVT ensemble using Berendsen thermostat during 100 ps to ensure equilibration.

Interaction at the interface.
The interface between the epoxy and graphene was prepared by layering the crosslinked epoxy cell on top of the energy minimized graphene cell. For calculations only non-covalent bond interactions were considered. Epoxy was set to match the graphene layer and follow lattice orientation from the latter as the former is more flexible and adjusts better to the more rigid graphene. The layer was gradually compressed in the Z coordinate until it reached a density close to the theoretical value of epoxy, 1.15 g.cm -3 . Each compression cycle was followed by a dynamics NVT for 50 ps. To avoid large energy increases to the system, this compression was performed stepwise [15]. Finally, NPT dynamics was performed with Berendsen as both thermostat and barostat for 100 ps in order to equilibrate the final structure.

Mechanical properties of graphene.
Each molecular model was energy minimized using a combination of the steepest decent and the conjugate gradient methods. The graphene was strained in the x-axis direction with steps of 0.5 %. A dynamics run using NVT ensemble was performed after each step with Berendsen thermostat for 100 ps. Atomistic stress calculations were based on the virial expression [16]. The stress-strain plots were used to calculate the Young's modulus. The influence of grafting density on the mechanical properties of graphene was studied by randomly grafting hydroxyl groups to the graphene structure in different percentages, 3.0 % and 5.0 %. The effect of the distribution of hydroxyl groups as a single cluster or several small clusters was assessed by comparing these models with dispersed hydroxyl groups on the graphene. This latter study was performed for 10 % functionalised layers.

Pull out simulations
For the calculation of the interaction energies and the interfacial shear strength the periodic boundary conditions were removed and hydrogen atoms were added to the graphene edges and fixed. The graphene was pulled in the x-axis direction with increments of 5 Å until complete separation of graphene and epoxy, as shown in figure 2. Geometry optimization was performed after each pull out.

Graphene mechanical properties
The Young's modulus for pristine and functionalised graphene are presented on figure 3. Close to what is reported in the literature for MD simulations [16][17][18], the pristine graphene obtained the highest Young's modulus, 1.35 TPa. This is attributed to its sp 2 hybridised structure. A decrease on the Young's modulus is verified with the grafting of functional groups onto the graphene layer. Both -COOH and -OH groups decrease the Young's modulus of graphene 42 % and 47 %, respectively. The largest drop occurs for carbonylated graphene, 53 %. With the grafting of -COOH and -OH groups the carboncarbon bond hybridisation changes to sp 3 in order to promote the chemical bond (figures 5 and 6). This change of in-plane sp 2 -hybridisation to out-of-plane sp 3 allows a larger mobility of the carbon due to unsupported sp 3 bonds [19,20]. Several simulations of functionalised graphene layers have reported a decrease on the mechanical properties of graphene [3][4][5]. With a different type of grafting -C=O promotes the complete breaking of the sp 2 bond between one neighbour atom of the grafted carbon creating a vacancy defect, which has been previously studied with AIREBO potential in MD simulations [21], and proved a decrease in strength of graphene with 2.0 % vacancy concentration.  The structure of graphene was studied by calculating the bond length distribution for systems before and after equilibration. Pristine graphene has a homogeneous length distribution of 1.42 Å prior to equilibration [16]. Some minor out-of-plane displacements occurred after equilibration as a consequence of ring strain relief and thermal vibrations. This leads to a small change in bond length distribution of about 0.20 Å (figure 5). The influence of functional groups is however more noticeable right before equilibration. The grafting of carboxylic groups on the graphene and consequent conversion of local carbon displacements from in plane sp 2 to out-of-plane sp 3 hybridisation, while at the same time the grafted carbon is pulled off plane in a pyramidal structure (figure 4), increases the bond length between grafted carbon and structural carbons to approximately 1.53 Å (figure 5). Thus decreasing the surrounding bonds to 1.40 Å. After equilibration the bond length distribution broadens. Out-of-plane displacements are even clearer. Hydroxyl groups have a similar behaviour as carboxylated surfaces except for out-of-plane displacements after equilibration which are less severe ( figure 6).
Unlike the previous functional groups, carbonyl promote the breaking of a sp 2 hybridised bond between the grafted carbon and one neighbouring carbon in order to covalently bond an oxygen atom by sp 2 hybridisation. This creates a defect on the graphene surface, resulting in an out-of-plane displacement similar to carboxyl groups with broader bond length distribution. After equilibrium the bond length distribution increases ( figure 6). Such behaviour has been reported in the literature [2,19]. Calculations showed a decrease of the Young's modulus with increasing functionalisation, from 0.77 to 0.74 to 0.72 TPa for 3.0 %, 5.0%, and 10 % grafting, respectively. The major effect is however verified for the smallest degree of functionalisation as with increasing grafting percentages the influence on the Young's modulus is not so noticeable. This seems to suggest that a smaller degree of funtionalisation is enough to decrease the mechanical properties of graphene, and an increase of that percentage will not have a much larger effect. These results come in agreement with a study by Qin and Buehler [3] who verified the influence of functional groups density on the fluctuations of the graphene sheets and consequently their mechanical properties. Comparisons between 10 % grafted groups distributed randomly, as a single cluster, and as several clusters (figure 8) were made. The single cluster obtained the lowest Young's modulus, 0.68 TPa, which can be understood by a larger interaction between concentrated hydroxyl groups together with a higher amount of sp 3 bonds confined in a reduced area. The smaller clusters model on the other hand have less impact on the elastic modulus, comparable to randomly dispersed groups, between 0.72 -0.75 TPa. This study proves that density and distribution of functional groups influence the mechanical properties of graphene, as shown by Pei et al. [4].
Regarding their structural properties, graphene with single cluster and several clusters have an average bond length distribution of 1.42 Å. Out-of-plane displacements are more concentrated for the single cluster, whilst for small clusters and randomly distributed groups they are dispersed and less pronounced.

Interfacial characterisation
The interaction energy was estimated from the energy difference (∆E) between the total potential energy of the composite and the sum of the potential energies of individual molecules as shown in equation (1): Where Etotal is the total potential energy of the system, Egraphene the potential energy of the graphene, and Eepoxy the potential energy of the epoxy resin. The interaction energy is negative during each stage of the pull-out showing the existence of attractive interaction between polymer and reinforcement. The ISS, through van der Waals interactions, increases linearly with applied strains, and tend to stabilize to a constant value close to zero after complete separation of the surfaces as a result of the lack of interaction.
The interfacial shear strength can be evaluated by analysing the changes in potential energy during the pull out of the graphene from the epoxy and following the relations on equations (2) and (3): Where W is the width of the graphene (W = 49 Å) and L the length (L = 51 Å). As there are no periodic boundary conditions applied and we have a two-layer system, the graphene only interacts with one side of the epoxy, thus we have simply w. Furthermore, even though L is assumed, figure 9 shows that a distance of 60 to 65 Å (depending on the system) is a good value taking into account that after that length there are no considerable interactions.
The pull out force was obtained through differentiation of the potential energy curve and based on equation (4): Where ∆X is the total displacement. The results displayed on figure 10 show that the force required for pulling the pristine graphene from the polymer is extremely lower than for functionalised layers thus suggesting very little interaction between the two. On opposite sides, carboxylated and carbonylated graphene seem to interact very well with the epoxy has the forces required for pulling the graphene layers are the highest of the four models studied (14 nN), despite lowering considerably after the initial stage. Even the ISS are similar for both, 113 and 111 MPa, respectively. Graphene-OH is the one among the functionalised layers who presents less interaction with the epoxy. It reaches maximum values of 38.8 MPa for ISS and 2.5 nN of force required for the pull out, during the initial stages.

Conclusion
The results obtained show that functionalisation of graphene sheets improve considerably the ISS and pull out force for epoxy resin composites. Pristine graphene seems to have almost no interaction with the polymer whereas carboxyl and carbonyl groups increase considerably this interaction. Hydroxylated graphene shows an improvement of the pull out force and ISS compared to pristine graphene however these values are much lower than the other functionalised graphene models studied. Despite an increase in interface interactions, the determination of the Young's modulus showed a decrease of the mechanical properties of functionalised layers, with the smallest drop occurring for carboxylated graphene, 42 %, and the highest for carbonylated, 57 %. A higher coverage density also proved a decrease of the Young's modulus. This is expected due to the change of hybridisation of grafted carbons, which also implies a broadening of the bond length distribution. The organisation of these groups in a single cluster also means a greater negative impact on the mechanical properties. Overall the results suggest a good commitment between loss of mechanical properties and improvement of interfacial strength should be achieved with carboxyl functional groups.