Determination of elastic constants of functionalized graphene-based epoxy nanocomposites: a molecular modeling and MD simulation study

Functionalization of graphene is the best way to create a high degree of dispersion and bonding to polymer matrix in order to obtain high performance composites. The effects of carboxyl (−COOH) functionalized graphene (FG) on the mechanical properties of its epoxy-based nanocomposites have been examined by molecular dynamics (MD) simulations. Simulations cells of nanocomposites with varying wt% of FG (1, 2, and 3 wt%) were constructed using Material Studio 6.0. The MD simulation findings of nanocomposites reveal that they have better mechanical properties such as elastic modulus, bulk modulus, shear modulus, and the Poisson’s ratio than pure epoxy. Furthermore, the computational results of nanocomposites have been effectively confirmed with available experimental data. Therefore, the current MD simulation shows a decent computational sign for the existing experimental and simulation outcomes on mechanical properties of FG/epoxy nanocomposites.


Introduction
Epoxies, also known as thermosetting polymers, are made up of three-dimensional (3D) cross-linked structures. Epoxies having 3D cross-linked structures exhibit improved characteristics, making the materials more appealing for use in a variety of applications such as electronic packaging, adhesives, coatings, and composite materials. Among the properties connected with polymers studied by various researchers, little effort has been made to examine the cross-linking process and structure-property of cross-linked polymers and their composites. The cross-linked polymers have permanent covalent bonds and cannot be re-shaped and re-melted. Several researchers examined epoxy polymers made from EPON 862 resin that had been cross-linked with triethylenetetramine (TETA) [1,2] curing agents and diethyl toluene diamine (DETDA) [3][4][5][6]. Researchers analyzed the cross-linking of diglycidyl Ether of Bisphenol A (DGEBA) resin with methylene glycol di-p-aminobenzoate (TMAB) [7], isophorone diamine (IPD) [8], diamine, diamino diphenyl sulfone (DDS) [9], tri DETDA [10], poly-oxy propylene (POP) diamines [11], and methylene di-aniline (MDA) [12].
Tsige and Stevens [13] investigated the crack behavior of strongly cross-linked polymer networks in terms of crosslinking functionality and interfacial bond density. Hölck et al. [14] estimated the mechanical properties of nanocomposites formed by cross-linking phenol novolac epoxy resin with bisphenol-A hardener. Using MD simulations, they investigated the interfacial characteristics of cross-linking of epoxy phenol novolac (EPN) and SiO 2 layered nanocomposites. Also, various researchers examined the performance of cross-linked epoxy nanocomposites reinforced with SiC [1], Al 2 O 3 [1] nanoparticles, and carbon nanotubes (CNTs) [15,16] using MD simulations. Several other works including MD simulations and the modeling of cross-linked nanocomposites were discovered in the literature [17][18][19].
Graphene is an allotrope of carbon one-atom-thick consisting of a single layer of carbon atom arranged in two dimensional (2-D) honeycomb lattice which was discovered by Novoselov et al. [20]. Graphene can be manufactured by the exfoliation of graphite as a single layer, double layers, and multilayers. Many researchers are attracted towards graphene due to its outstanding characteristics such as high strength (1.0 TPa) [21], stiffness (0.7 TPa), thermal conductivity 3500 W m −1 K −1 and surface area (2630 m 2 /g) [22,23], optical transmittance (97.7%) [24], and electrical conductivity with the recent density of 108 A/cm 2 and able to carrying densities of electrical higher than copper [25]. In recent times, functionalized graphene (FG) and its derivatives are used as a reinforced filler in place of graphene to overcome the problems of dispersion as well as an agglomeration of graphene in polymer composites and improve the mechanical properties of polymer nanocomposites followed as potentially durable, tough, and lightweight materials [26][27][28]. Researchers belong to the field of polymer composites that preferred FG as a reinforcing element due to its outstanding mechanical and thermal properties [29].
Liang et al. [28] investigated the mechanical properties of graphene oxide (GO)/poly (vinyl alcohol) (PVA) nanocomposites fabricated using a water solution processing method. They found that Young's modulus and tensile strength of nanocomposites increased by 62% and 72%, respectively, by the addition of 0.7 wt% of GO. Likewise, Qi et al. [30] fabricated liquid crystalline epoxy grafted GO infused epoxy composites. Results showed that the reinforcement of 1.5 wt% of fillers enhances the tensile strength of composites by 48% (from 55.43 to 80.85 MPa) as compared to pure epoxy.
Nowadays, molecular modeling and simulation of materials are effective methods to solve real-world problems at an optimized cost and offer to explore the research up to a certain extent. Nayebi et al. [31] studied the mechanical properties of graphene/polythiophene (PT) composites by MD simulations. Computed Young's modulus of graphene/ PT composites was 99.6 GPa, which is higher than the calculated value of pure PT, i.e., 1.29 GPa, as well as the experimental value of 1.5 GPa. In another research, Lee et al. [32] examined the mechanical properties of FG reinforced epoxy nanocomposites using MD simulation. They have found that graphene containing both −COOH and −NH 2 functional groups exhibits improved tensile strength, modulus, and toughness of nanocomposites compared to pristine graphene. Based on MD simulation study, Lin et al. [33] reported that the increasing weight fraction of graphene in epoxy increases the Young's and shear moduli of graphene/ polymethyl methacrylate (PMMA) nanocomposites; however, the temperature rises deteriorate its overall properties. Shen et al. [34] examined the cryogenic mechanical properties of the nanocomposite at room temperature. They obtained enhanced Young's modulus and tensile strength of nanocomposites compared to pure epoxy at graphene loading of 0.1 wt%. Yu [35] reported an improvement of 4.563 GPa in Young's modulus of aminated graphene (AG)/epoxy nanocomposites, i.e., 51% higher value compared to pure epoxy.
To the best of the author's knowledge, there is a less analogous study on molecular simulation of LY556 epoxy resin cross-linked with DETA curing agent to estimate the mechanical properties of −COOH-FG/epoxy nanocomposites with different wt% of graphene using MD simulation methods. In this paper, the authors examined the effect of −COOH-FG on the mechanical properties of cross-linked epoxy-based nanocomposites. Key investigation parameters such as bulk moduli (K), Young's moduli (E), Poisson's ratio (ν), and shear moduli (G) of epoxy nanocomposites have been measured for three different wt% of −COOH-FG. MD simulation results were compared with the investigational results as well as in good agreement with the experimental work.

Molecular structure of epoxy resin (DGEBA) and curing agent (DETA)
Epoxy resin is used as a thermoset polymer matrix which has many epoxide groups for each molecule and displays outstanding high elastic modulus, thermal conductivity, and corrosion resistance. Due to these exceptional properties, epoxy resin as a polymer matrix was used to construct polymer nanocomposites in different engineering areas.   Fig. 1(b) shows that the DETA consists of five reactive sites with primary as well as secondary amine groups and hence is a multifunctional (five-fold-functional) reactant. Due to this, each DETA molecule can be reacted with five DGEBA molecules, also capable to react with another DETA molecule by its second epoxide head. Hence, DGEBA and DETA can produce 3D cross-linked epoxy polymers.

Cross-linking process
In epoxy resin, the C-O bond of each epoxide group is reformed to form a reactive -CH 2 site (Fig. 2); thus, it can be cross-linked with a curing agent. Figure 3(a) represents the molecular model of LY556. Hydrogen atoms present in the amine (−NH 2 ) groups of the curing agent (di-ethylene toluene di-amine −DETDA) react with the epoxide groups of epoxy resin molecules ( Fig. 3(b)). Furthermore, result cured   (Fig. 3(c)) react with the epoxy molecule at the site of HN and NH 2 and one curing agent molecule reacts at the site of another epoxide group. Hence, the crosslinking was established between the curing agent and epoxy resin represented in Fig. 3(d). Finally, the effect remains and produces more cross-link motions in all directions between the curing agent and epoxy resin [36,37].

Molecular structure of functionalized graphene
Chemical functionalization was formed on the covalent linkage of functional groups onto the carbon domain of graphene. Graphene is used as an efficient reinforcement in polymer composites. Dispersion of graphene is considered the most common difficulty that needs to be focused upon to use graphene as an effective reinforcement in polymer composites. Various organic molecules, like −NH n , −CH n , −OH, and −COOH, are formed by chemical covalent bonds among the polymer chains.
Here, different functional groups are discussed, out of which the carboxyl group has high reactivity and could be instantly absorbed into the epoxy resin, as described by Lee et al. [37]. It can be authorized to overcome the problem of dispersion and quantity of graphene in the epoxy matrix. Thus, the chemical covalent bond mechanism of functional groups on the surface of graphene helps maintain the stable bonds between graphene and epoxy. In this work, graphene was functionalized with carboxyl functional groups, which are randomly grafted to the carbon atoms of the graphene sheet, as shown in Fig. 4(b). Moreover, hydrogen atoms were attached at the end of the graphene sheet to avoid their unsaturated boundary effects [38] represented in Fig. 4(a). Chemical functionalization did obtain a robust adhesive interface and quality dispersion between the adjacent polymer chains and graphene. All the constructed models were relaxed through energy minimization processes using the MD simulation package in Material Studio 6.0.

Software and force field
A force field used to model the interatomic interactions is the condensed-phase optimized molecular potential for atomistic simulation studies (COMPASS). It is the first ab initio force field that enables accurate and simultaneous prediction of gas-phase properties (structural, conformational, vibrational, etc.) and condensed-phase properties (equation of state, cohesive energies, etc.) for a broad range of molecules and polymers. Polymer consistent force-field (Pcff30) force field belongs to ab initio consistent force field family with the same functional form. Pcff30 force field has been effectively applied to simulate the cross-linked DGEBA with DETA curing agent, nanocomposites with different wt% of FG. It can also be used for cohesive energies, mechanical properties, heat capacities, and elastic constants. It handles electron delocalization in aromatic rings utilizing a charging library rather than bond increments. Pcff30 has been used in the entire dynamic direction and many molecules as an ab initio force field. Ab initio force field is used to determine the properties of many molecules like polymers, most repeated organics, and minor inorganic molecules [39,40]. Pcff30 force field contains three different terms of energy such as cross-terms, bonded energy terms, and non-bonded energy terms. Covalent bond stretching, torsion angle rotation energy, and bond angle bending of the polymer chain represent the bonded energy terms and van der walls forces, coulombic electrostatic forces represents cross-terms and non-bonded energy terms correspondingly, here, the Lennard-Jones potential equation is used to obtain van der Wall interaction. Various functional forms in the Pcff30 force field have been discussed by H. Sun [41] and are as follows in Eq. 1: The valence energy term consisted of the covalent bondstretching energy, E bond ; the bond-angle bending energy, E angle ; the torsion-angle rotation energy of the polymer (1) E total = E valence + E cross−term + E non−bond chain, E torsion ; and the out-of-plane energy, E oop , as discussed in Eq. 2.
The terms for cross-interactions included the dynamic variations among the bond stretching interactions between two bonds, E bond−bond ; bond-bond interactions between two valence angles associated with a common vertex atom, E angle−angle ; stretch bond interactions between a two-bond angle and one of its bonds, E bond−angle ; stretch-torsion interactions between a dihedral angle and one its middle bond, E middle−bond−torsion ; stretch-torsion interactions between a dihedral angle and one of its end bonds, E end−bond−torsion ; bend torsion interactions between a dihedral angle and one of its valence angles, E angle−torsion ; and bend torsion interactions between a dihedral angle and its two valence angles, E angle−angle−torsion , as discussed in Eq. 3.
The non-bond energy terms represent the Coulomb electrostatic force, E Columb , and the vdW force, E vdW , as discussed in Eqs. 4 and 5.
where q is the atomic charge, ɛ is the dielectric constant, and r ij is the distance between atom pairs. "b" is the internal coordinates of the bond, "is the length of the other O-H bond, "θ" is the bond angle, (2) E valence = E bond + E angle + E torsion + E oop ij "φ" is the torsion angle obtained by twisting the atoms about the bond axis, and x is the out of plane angle.
, B ij are fitted from quantum mechanics calculations and are implemented into the Discover module of Materials Studio.

Simulation models
In this paper, the epoxy resin cross-linked with a curing agent has been used as a polymer matrix. A single-layer graphene sheet with 164 carbon atoms was selected as nanofiller in the nanocomposites. Figure 5 represents the simulation model of pure cross-linked epoxy resin. Figure 6(a-c) shows the simulation cells of epoxy nanocomposites reinforced by three different wt% (1, 2, and 3 wt%) of −COOH-FG. In simulation models, the atom has been constructed via MD simulation method to examine the elastic properties of nanocomposites by assuming the cut-off distance (5 Å). The geometry optimization was done at the maximal displacement of 0.005 Å, atomic force 0.05 eV, the total energy 1 × 10 −6 eV/atom, stress 0.005 GPa, and cut-off energy 270 eV. Furthermore, the amorphous cell was constructed at room temperature with the help of a module tool at an initial density of 0.5 g/cm 3 .

Simulation procedure for graphene reinforced epoxy nanocomposites
In the present work, mechanical properties of FG/epoxy nanocomposites are predicted by MD simulations. Firstly, the models are exposed to several steps of equilibration procedure followed by a series of energy minimization. This was accomplished via Discover Minimization module using the Steepest Descent technique with 10,000 iterations with a convergence rate of 0.1 Kcal/mol/A. The systems were equilibrated in the NVT (isothermal) ensemble for a short period before NPT (isothermal-isobaric) ensemble to minimize the total energy of the initial system [42]. During the equilibration process, the volume change has been neglected because the volume fraction of the particles varies a little from the initial unit cell. Figure 7 depicts the FG/epoxy nanocomposite model after equilibration, which contains an amorphous cell-sized 44.8 × 44.8 × 44.8 Å 3 , with a low density of 0.5 g/cm 3 formed at room temperature under periodic boundary conditions. To achieve the stable structure of a model, the potential, lattice parameters, and non-bond energy have been minimized using MD simulations. Later, a step-by-step process of MD simulation has been completed with a time of 1 fs and 1 atm at 298 K to obtain equilibrium density and a balanced system. The inter-layer separation distance of graphene is assumed to be 3.4 in the assigned simulation [43].

Results and discussion
There are three main methods to evaluate the mechanical properties of nanocomposites by MD simulations, fluctuation formula [44], dynamics (constant-stress molecular dynamics), and static (constant-strain minimization) [45].
In the present work, the static method was used to calculate the mechanical properties of −COOH-FG-based crosslinked epoxy systems. Stress-strain behavior for linear elastic materials can be described by Hooke's law (Eq. 6): Where, i; j = 1, 2, 3, 4, 5, and 6; σ i ; and ε j are the 6-dimensional stress-strain vectors, and C ij is the 6 × 6 stiffness matrix.
The strain amplitude was set to 0.005 to obtain the elastic constants. The stress components were calculated using the so-called viral expression (Eq. 7) [46].
Where V is the volume; m k and u k stand for the mass and velocity of the kth particle, respectively; rkl is the distance between kth and lth particles; and flk represents the force exerted on lth particle by kth particle.
The Lamé coefficients, λ and μ, can be calculated using any two of the following equations (Eqs. 8-13). Where, K, E, ν, and G are bulk moduli, Young's moduli, Poisson's ratio, and shear moduli, respectively. Table 1 represents the elastic stiffness matrix of pure cross-linked epoxy obtained by the Forcite calculation available in the module tool of Material Studio software to calculate its mechanical properties. Similarly, Tables 2, 3, and 4 elucidate the elastic stiffness matrices of epoxy nanocomposites reinforced with 1, 2, and 3 wt% of −COOH-FG, respectively. In the nanocomposite unit cell, graphene grafted with carboxyl functional groups acts as nano reinforcement, while cross-linked epoxy resin acts as a matrix.
(13) G = Table 5 depicts and compares the estimated mechanical properties of −COOH-FG cross-linked epoxy nanocomposites to available research data. As the loading of the FG increases from 1 to 2 or from 2 to 3 wt% at a constant volume, the number of molecules increases in epoxy system. Furthermore, the molecules in the epoxy system have enough time to find the optimal position that minimizes total energy for a well-balanced unit cell structure. After carrying out numerous simulations, it was discovered that functionalization is beneficial to the elastic properties of the resulting nanocomposites. Simulation results prove that the functionalized graphene-containing epoxy system shows enhanced mechanical properties   compared to the pure epoxy system. It has been noticed that the FG containing epoxy resin system with 3 wt% of fillers loading showed higher elastic moduli compared to the pure epoxy system. It is concluded that covalently grafted −COOH functional groups on graphene surfaces increased elastic moduli values, which is illuminated by hybridization theory. Pristine graphene has an out-of-plane sp 2 -hybridized geometrical configuration, and −COOH-functionalization transforms the hybridization to sp 3 , which adds to an increase in covalent bonding between the polymer chains and the FG in the nanocomposites. Thus, the chemical covalent bond mechanism of functional groups on graphene's surface contributes to the stability of the bonds between FG and epoxy resin. Surprisingly, the calculated results agree well with the results obtained by Lee et al. [32], who estimated the improved mechanical properties of epoxy nanocomposites infused with graphene containing both −NH 2 and −COOH functional groups. Moreover, the MD simulation results revealed that the elastic modulus of −COOH + −NH 2 FG reinforced epoxy nanocomposites increases from 3.1 to 3.45 GPa compared to −COOH-FG reinforced nanocomposites. Sun et al. [47] noticed that Young's modulus of graphene/epoxy nanocomposites increases by 2.74 GPa with filler wt% of 0.88. In addition, Rafiee et al. [48] observed a ~30% improvement in Young's modulus of nanocomposites with ~0.3% weight fraction graphene nanoribbon in an experimental study. The examined value of the elastic modulus of −COOH-FG/epoxy nanocomposites is higher than the other experimental values, as shown in Table 5. It could be because simulation systems have no imperfections and thus represent an ideal situation, whereas physical models contain a variety of imperfections that result in compact calculated properties.

Conclusions
Molecular models of epoxy nanocomposites with varying wt% of −COOH-FG have been constructed to determine the elastic constants such as Young's modulus, bulk modulus, shear modulus, and Poisson's ratio. In this study, the MD simulation approach with Pcff30 force field has been effectively applied to simulate the cross-linked epoxy resin with DETA curing agent. The constant strain approach has been adopted to calculate the mechanical properties of nanocomposites. The MD results showed that the mechanical properties of nanocomposites such as moduli, Young's, bulk, shear moduli, Poisson's ratio, and shear moduli significantly improved by adding −COOH-FG to cross-linked epoxy. The best results were obtained at 3.0 wt% of FG, which was in good agreement with experimental as well as related MD simulation work. The MD simulation showed more interaction energy among the graphene and surrounding matrix for higher wt% of graphene due to −COOH functional groups on the surface of graphene. Therefore, −COOH-FG is used as a filler material to improve the elastic constants of the FG/epoxy nanocomposites. Results of the present work have been compared with available experimental and MD simulation data as illustrated in Table 5, recognizing that the current MD simulation shows a decent computational sign for the existing experimental and simulation outcomes on mechanical properties of FG/epoxy nanocomposites.
Author contribution Aman Yadav: conceptualization, modeling and simulation, methodology, formal analysis, investigation, visualization, writing -original draft, and writing -review and editing; Amit Kumar: supervision and writing -review and editing. Kamal Sharma: writing -review and editing. A.K. Pandey: writing -review and editing.

Data availability NA.
Code availability NA.