In this paper, we present an approach for characterizing the interfacial region using the molecular dynamics (MD) simulations and the shear deformation model (SDM). The bulk-level mechanical properties of graphene-reinforced nanocomposites strongly depend on the interfacial region between the graphene and epoxy matrix, whose thickness is about 6.8–10.0 Å. Because it is a challenge to experimentally investigate mechanical properties of this thin region, computational MD simulations have been widely employed. By pulling out graphene from the graphene/epoxy system, pull-out force and atomic displacement of the interfacial region are calculated to characterize the interfacial shear modulus. The same processes are applied to 3% grafted hydroxyl and carboxyl functionalized graphene (OH-FG and COOH-FG)/epoxy (diglycidyl ether of bisphenol F (DGEBF)/triethylenetetramine (TETA)) systems, and influences of the functionalization on the mechanical properties of the interfacial region are studied. Our key finding is that, by functionalizing graphene, the pull-out force moderately increases and the interfacial shear modulus considerably decreases. We demonstrate our results by comparing them with literature values and findings from experimental papers.

## Introduction

Thermoset epoxy resins are lightweight and exhibit remarkable mechanical performance compared to thermoplastic polymers due to the strong and entangled covalent bonding between the epoxy resins and hardener. To date, the epoxy resins have been adopted as polymer matrix of nanoscale fillers such as nanoparticles, carbon nanotubes (CNTs), and graphene in various studies [1–6]. Among them, graphene nanofillers—a thin sheet of sp^{2} bonded carbon—are known to bring enhancements on various physical and chemical properties such as mechanical stiffness, strength, electrical and thermal conductivity, and glass transition temperature at relatively low graphene concentration in the composites [5–8]. Compared to its fierce rival, CNT, the graphene-based nanofillers are cost-effective due to its possibility of mass production and might outperform CNT fillers regarding enhancements in mechanical, thermal, and electrical properties [8–12].

One of the major factors that influence bulk-level mechanical properties of nanocomposites is the interfacial region, a thin area in the polymer adjacent to the nanofillers, which is structurally different from that of the bulk polymer as revealed from atomistic simulations and transmission electron microscopy (TEM) image [1,5,13]. Because it is a challenge to characterize mechanical properties of this thin region in experimental ways, whose thickness is about 6.8–10.0 Å, computational molecular dynamics (MD) simulations have been widely used [1,5,13]. However, the mechanical properties of the interfacial region differ from characterization methods; the inverse calculation method expects higher elastic properties [1] whereas the interfacial energy differentiation method expects lower elastic properties than that of the bulk matrix [14]. Furthermore, even though there have been numerous studies on the functionalized graphene, which is known to enhance filler dispersion and interfacial bonding, there has been little research on characterizing mechanical properties of the interfacial region in the functionalized graphene (FG)/epoxy systems.

In this paper, we obtained the pull-out force and shear modulus of the interfacial region through the pull-out simulations, and conducted the same simulated tests on 3% graphited hydroxyl and carboxyl functionalized graphene (OH-FG and COOH-FG)/epoxy systems to analyze influences of the functionalization on mechanical properties of the interfacial region. Our key finding is that, by functionalizing graphene, the pull-out force increases moderately and the interfacial shear modulus considerably decreases. Although it is widely known that graphene functionalization enhances the effective stiffness of nanocomposites, at nanoscale, the stiffness of the interfacial region decreased in the current study; this is elucidated through findings from experimental literatures. Section 2 describes MD simulation models and methods. In Sec. 3, we propose a shear deformation model (SDM) that is combined with MD simulation results to characterize the interfacial region. Section 4 presents results and discussions. Finally, conclusions are made in Sec. 5.

## Molecular Dynamics Simulations

In order to construct nanoscale molecular structures, calculate the elemental strain and atomistic potential energy, and conduct the pull-out simulations, commercial software for molecular dynamics simulations, materialsstudio 2017 package [15], was adopted. COMPASSII (condensed-phase optimized molecular potentials for atomistic simulation studies) force field was applied to describe interatomic interactions as it has been widely employed to simulate systems containing amorphous polymers and carbon nanofillers [16,17].

### Preparation of Model Structures.

For the composite matrix, we adopted the crosslinked epoxy system, where diglycidyl ether of bisphenol F (DGEBF) was used as the epoxy resin, and triethylenetetramine (TETA) was used as the hardener, as shown in Fig. 1. The crosslink can be generated in two different ways in MD, which are the dynamic crosslinking simulation and the representative molecule method [18]. Although both methods yield reasonable coincidence with experimental data [18], the dynamic crosslinking simulation was adopted in this paper to describe a complex configuration between the epoxy matrix and the graphene. It is known that the thermoset epoxies develop high crosslinking density [18–20]. However, according to Khare and Khare [21], a large cutoff distance such as 40 Å and extensive computational power are needed to achieve almost 100% crosslinking density. Therefore, a conversion ratio of 80% with 10 Å cutoff distance was set to be our target crosslinking density in order to cut down the computational power while maintaining high crosslinking density. Changes in mechanical properties by different crosslinking densities are beyond the scope of the present study. In order to examine mechanical properties of the pure epoxy matrix of 80% target conversion ratio, an amorphous cell containing 360 molecules of DGEBF and 120 molecules of TETA was developed. Mechanical properties were then calculated from the amorphous cell through the constant strain energy method [22].

### Interfacial Thickness Calculation.

A vertical density distribution profile was obtained from the dynamically crosslinked graphene-epoxy nanocomposite to define the thickness of the interfacial region. Detailed simulation procedures are explained in this section.

A three-layered amorphous cell, which consists of epoxy/hardener mixture (matrix1), graphene monolayer, and epoxy/hardener mixture (matrix2), was constructed as shown in Fig. 2. Before the dynamic crosslinking simulation, the cell went through constant number of particles, temperature, and pressure (NPT) ensemble dynamics with Parrinello barostat at 1 atm and 500 K for 50 ps to minimize the total potential energy and compress the layered amorphous cell. Reactive crosslinking atoms in matrix1 and matrix2 were designated differently during the dynamic crosslinking simulation even though both layers were made up of the same epoxy resin and hardener (DGEBF and TETA). This was because if they are considered as one family on both matrices, the reactive atoms over and under the graphene can be crosslinked, which is physically impossible. Therefore, the dynamic crosslinking simulation was conducted separately in matrix1 and matrix2, with the target conversion of 80%. Table 1 shows compositional details and the achieved crosslinking density of matrix1 and matrix2. Although the target conversion was set as 80%, the resulting crosslinking density could not be achieved exactly 80% owing to the process of the dynamic crosslinking simulation, where new crosslinks were developed at every cutoff distance. Therefore, the resulting crosslinking density was slightly different from each simulation. The crosslinking densities of graphene/epoxy and FG/epoxy systems are shown in Table 1. As the deviation is not significant, the different crosslinking densities on each system were not taken into account in this study. Periodic boundary conditions were assumed in the repeated unit cell for the crosslinking simulations.

The amorphous cell then went through a series of energy minimization processes. At first, the geometric structure of the amorphous cell was optimized through the conjugate gradient method. After stabilizing the geometric configurations, the repeated unit cell was subjected to a series of ensemble dynamics to minimize the total potential energy, while it avoids a local energy minimum at the same time. First, an constant number of particles, volume, and temperature (NVT) ensemble with Andersen thermostat was performed at 900 K for 9 ps with 0.3 fs time-step. At this stage, enough thermal energy was given to each atom to avoid being trapped to a local energy minimum. Next, another NVT ensemble was performed at 600 K for 15 ps with 0.5 fs time-step, and an NPT (constant pressure and temperature) ensemble with Berendsen barostat was performed at 298 K and atmospheric pressure for 50 ps, 1 fs time-step. The lattice size of this energy minimized amorphous cell was about $60\xd760\xd760\xc5$. The density distribution profile was then obtained along the normal direction of the graphene, which is *z*-direction.

### Pull-Out Simulations.

In order to comprehend interfacial dynamics and obtain mechanical properties of the interfacial region, pull-out simulations were performed through analogous methods suggested by several studies [23–27]. The *x*-direction lattice parameter of the amorphous cell was elongated to 75.5 Å, which is a summation of the graphene length (60 Å) and van der Waals cut off distance (15.5 Å), in order to remove the nonbonded interaction between successive cells along the *x*-direction [23,24]. The top and bottom regions of the graphene/epoxy system in the *z*-direction, marked by rounded rectangles in Fig. 3, were fixed [23,26,27]. At every step, the graphene was pulled out to 1.5–3 Å, and the cell went through energy minimization by the conjugate gradient method [23–25]. The ensemble equilibration was not performed in this step to avoid huge computational burden; in fact, there is no clear evidence that the ensemble equilibration guarantees results, which are compatible with experiments. Successive steps were repeated until the graphene was segregated from the epoxy matrix. Figure 4 shows the scheme of the pull-out simulation at different pull-out displacements. At each step, *x*-coordinate of the representative atom sets and the interfacial energy, which will be discussed in Sec. 3, were recorded in a study table.

### Functionalized Graphene.

The same process was applied to the FG/epoxy systems: hydroxyl functionalized graphene (OH-FG)/epoxy system and carboxyl functionalized graphene (COOH-FG)/epoxy system. Both OH and COOH functional groups were grafted to pristine graphene with the same grafting percentage of 3%. For example, the pristine graphene used in this study consists of 1312 carbon atoms. Three percentage grafting means that about 41 out of 1312 carbon atoms are bonded to the functional groups. Figure 5 shows structural configurations of the pristine graphene, OH-FG, COOH-FG, and COOH-FG/epoxy system. The FG/epoxy systems were built in the same process discussed in Sec. 2.2, where the FGs replaced the pristine graphene.

## Shear Deformation Model for Determination of Interfacial Properties

The basic concept of the SDM was applied to the pull-out simulation to obtain mechanical properties of the interfacial region. Figure 6 represents a configuration of the SDM, where $ug$ and $ui$ are the displacements of the matrix adjacent to the graphene and outer matrix, respectively. $T$ and $t$ are the thickness of the total interfacial region and graphene layer each, $L$ is the *x*-direction length of the matrix, $s$ is the coordinate axis starting in the middle of the embedded graphene, parallel to the *x*-direction.

where $Fi$ is the pull-out force. According to the shear lag model, the interfacial shear stress follows a distribution shown in Fig. 7 [28]. The pull-out force can be obtained from the interfacial shear stress integrated with the contact area. However, it also can be easily obtained from the interfacial energy under the assumption of linear elasticity.

where $x$ is the pull-out displacement. The left side of Eq. (6) then can be numerically calculated from the interfacial energy profile.

*m*1….

*m*24 at the upper boundary and

*m*25…

*m*48 at the lower boundary of the interfacial region) where they were arranged at equal intervals, as shown in Fig. 6. By charting

*x*-coordinates of the representative atoms at every pull-out step, $ug\u2212ui$ can be calculated numerically so that we can obtain the integral at the right side of Eq. (6). Now, the only unknown in Eq. (6) is $Gi$. Hence, the interfacial shear modulus $Gi$ is expressed as

## Results and Discussions

### Crosslinking Simulation.

The elastic moduli of the pure crosslinked epoxy model were compared with those listed in the literature, as shown in Table 2. Both Young's modulus and shear modulus of the present study are about 10–20% higher than the literature values. For example, according to Kim et al. [18], Young's modulus and shear modulus are obtained as 3.76 and 1.37$GPa$, respectively on 62.5% crosslinking density model, which are 10% and 11% lower than the present results of 80% crosslinking density. However, results in the present study are reasonable considering that the values obtained from MD simulations in the literature show enhanced mechanical properties at a high crosslinking density [1,18,29].

### Interfacial Thickness Definition.

Figure 8 is the density distribution profile along the *z*-axis. The peak in the middle represents the relatively high density of the graphene compared to that of the bulk matrix. There are slightly higher density regions on both sides of the peak in the middle, which are noted by dashed boxes in Fig. 8. The thickness of the interfacial region was assumed to be 8.5 Å each, which is bigger than the red boxes for two reasons: (1) to match the thickness with the literature values (6.8–10.0 Å) [1,5] and (2) to clarify the displacement difference between upper and lower end of the interfacial region.

### Pull-Out Simulation.

The purpose of the pull-out simulation is to calculate the right side of Eq. (9) numerically so that we can obtain the interfacial shear modulus $Gi$. The pull-out force $Fi$ on the numerator was calculated from the energy increment divided by the pull-out displacement while the integration on the denominator was obtained from the displacement of the representative atom sets. Before that, we have to demonstrate the deformation of the matrix is within the linear elastic zone. Therefore, results of the pull-out simulation will be presented in following five sections: pull-out force calculation, demonstration of the linear elasticity, integration of the interfacial region displacement, characterization of the interfacial shear modulus and comparison of graphene/epoxy and FG/epoxy systems.

#### Pull-Out Force Calculation.

Figure 9 shows the interfacial energy along the pull-out distance. The interfacial energy increases linearly from pull-out displacement 10 Å to 50 Å. This tendency coincides with CNT pull-out simulations reported by Li et al. [23], where they divided the graph into three stages: (1) initial ascent stage, (2) subsequent platform stage, and (3) final descent stage [23]. The length of both first and third stages is similar to van der Waals cutoff distance, which is also coincident with the literature [23]. According to Eq. (8), the slope of stage 2 is equal to a pull-out force $2.24nN$.

#### Demonstration of the Linear Elasticity.

We assumed that the deformation of the matrix is within the linear elastic zone until the pull-out displacement reaches to 35 Å. In Sec. 4.3.5, the pull-out force was calculated on stage 2, within a domain of 10–50 Å pull-out displacement. However, to demonstrate the linear elasticity of the interfacial region, we conservatively set the boundary of linear elasticity as 35 Å. The graphene in the graphene/epoxy system was removed not to apply any force to the matrix, as shown in Fig. 10(b). Then NPT ensemble was conducted to the system for 500 ps, and the average of displacements of the 48 representative atom sets was investigated, as shown in Fig. 10(c). The average of displacements decreases as the NPT ensemble proceeds, converging to minus six around 300 ps. Note that the negative convergence value is caused by the cell shrinkage during the NPT simulation. Because the graphene was removed, the cell shrinkage is unavoidable during the NPT simulation. As shown in Fig. 10(d), a remarkable drop in *x*-length from 70 Å to 50 Å is observed while little increase is observed in *y*- and *z*-length. Since the *x*-length of the cell decreased to 50 Å, the average of displacements converged to a negative value.

As the graphene is the only medium that applies the shear force to the matrix, removing the graphene means no more shear force is applied to the matrix. The NPT ensemble dynamics to the system at 1 atm is analogous to a situation where the shear force is just removed from the system. Because the linear elasticity means that the deformed shape caused by an external force fully recovers once the force is removed, the convergence of the average of displacements implies that the system is within the linear elastic zone.

#### Integration of Interfacial Region Displacement.

#### Characterization of Interfacial Shear Modulus.

Now, everything needed for characterizing $Gi$ in Eq. (9) has been acquired. By substituting 20.4 Å to *T*, 3.4 Å to *t*, and 60 Å to *b*, $Gi$ was obtained at each pull-out step. Figure 12 represents the interfacial shear modulus from the pull-out displacement 15 Å to 35 Å included in stage 2. The dotted line in Fig. 12 indicates the average value, 2.34$GPa$. There are fluctuating $Gi$ values over the average up to the pull-out displacement of 20 Å. These points may be attributed to a nonuniformity of the displacements (i.e., $ugandui$) of the epoxy matrix and small values of the displacement integration in Eq. (9). The numerator in Eq. (9), pull-out force, remained constant during the stage 2 pull-out steps. However, the displacement integration in the denominator of Eq. (9) fluctuated around 2–15 Å^{2} with an average of 9 Å^{2}, which led to relatively large fluctuation on the interfacial shear modulus. These fluctuations diminished as the displacement of the interfacial region (or the displacement integration) increased; this will be discussed in Sec. 4.3.5.

#### Comparison of Pristine Graphene/Epoxy and Functionalized Graphene/Epoxy Systems.

The same process was applied to 3% grafted OH-FG and COOH-FG/epoxy systems. As shown in Fig. 13, FG/epoxy systems showed higher pull-out forces than the pristine graphene/epoxy system as the functional groups enhance van der Waals interaction between FGs and matrix. Among the FGs, the pull-out force of the COOH-FG/epoxy system was higher than that of the OH-FG/epoxy system. Figure 14 shows how the functional groups develop an interlaced structure with the epoxy matrix. The COOH-group develops a more entangled structure with the matrix than the OH-group, which enables higher van der Waals interaction in the COOH-FG system. In the same manner, because there are no functional groups in the graphene/epoxy system, the smallest van der Waals interaction is expected for the graphene/epoxy system among the others.

A sharp decrease of the interfacial shear modulus was observed on FG/epoxy systems in Fig. 15(c). This is because of an increment on the displacement of the interfacial region in FG/epoxy systems, as shown in Fig. 16, where more deformation of the interfacial region was observed in FG/epoxy systems than the graphene/epoxy system. The displacement of the interfacial region can be expressed numerically as the displacement integration, the denominator of Eq. (9). Average of the displacement integration is about 9 Å^{2} for graphene/epoxy, 27 Å^{2} for OH-FG/epoxy and 194 Å^{2} for COOH-FG/epoxy system.

Note that both the pull-out force and the displacement integration increased when graphene was functionalized. The pull-out force is proportional to the interfacial shear modulus because it is in the numerator of Eq. (9), whereas the displacement integration is inversely proportional to the interfacial shear modulus because it is the denominator of Eq. (9). Therefore, the functionalization has both increasing effect and decreasing effect on the interfacial shear modulus. However, the pull-out force increased less than two times whereas the displacement integration increased even more than 20 times at COOH-FG/epoxy system. This led to the overall decrease of the interfacial shear modulus of the FG/epoxy systems. More decrease of the interfacial shear modulus in the COOH-FG/epoxy system compared to the OH-FG/epoxy system can be explained in the same way.

Comparing the pristine graphene/epoxy (Fig. 12), OH-FG/epoxy, and COOH-FG/epoxy systems (Figs. 15(a) and 15(b)), the fluctuation of the $Gi$ values decreased in the order of graphene/epoxy, OH-FG/epoxy, and COOH-FG/epoxy systems. It is due to the relatively larger displacement integration. Although the displacement integration still fluctuated around 170–220 Å^{2} in the COOH-FG/epoxy system, the fluctuation of the interfacial shear modulus was minimal because the displacement integration was the highest among the three systems.

Current simulation results are compared with literature values, as shown in Table 3. Yu et al. defined the thickness of the interfacial region between nanoparticle and epoxy matrix as 6.8 Å and inversely derived a stiffness matrix of the interfacial region using Mori-Tanaka micromechanics [1]. Because they used a nonfunctionalized nanoparticle in their atomistic simulations, their result can be compared to graphene/epoxy system in this study; they show moderate consistency. Arash et al. [14] obtained Young's modulus of the interfacial region between CNT and poly (methyl methacrylate) matrix from the double differentiation of interfacial energy. In their study, information of the interfacial shear modulus is not given, but we can estimate it by applying two assumptions: (1) the interfacial region is isotropic and (2) Poisson's ratio is 0.3. With these assumptions, we can obtain the interfacial shear modulus as 0.28–0.9 $GPa$, which is lower than the results in this study. A possible reason for the inconsistency is that they used a thermoplastic polymer that has lower mechanical properties than the epoxy. Also, the definition of the interfacial region is ambiguous because they obtained the interfacial Young's modulus only from the interfacial energy.

The decrease of the interfacial shear modulus in FG/epoxy systems does not imply that the total shear modulus would decrease or the nanocomposites would become less stiff. In fact, the interfacial region appears to be small compared to the whole nanocomposite considering that usually less than a few weight percent of the graphene is added to a polymer matrix, which means that a decreasing effect on the interfacial shear modulus is marginal. Instead, the functionalization still has positive effects on the stiffness such as better dispersion, which may lead to an overall increase of stiffness parameters of the nanocomposites. Furthermore, the less stiffened interfacial region of the FG/epoxy system plays a role as a bumper between the graphene and the bulk matrix, which may lead to a remarkable increase in fracture toughness. This tendency is found from an experiment done by Cha et al. [4]. They conducted mechanical testing on melamine functionalized CNT (M-CNT)/epoxy composite. They reported a remarkable enhancement on fracture toughness (95–100%) of the M-CNT/epoxy composite compared to that of neat epoxy, while Young's modulus and tensile strength increased only 64% and 22% [4]. This result is consistent with the interpretation of this study. As a result, graphene functionalization enhances toughness and ductility of the nanocomposites while it degrades stiffness of the interfacial region.

## Conclusions

A combined MD-SDM approach for characterizing the shear modulus of the interfacial region was developed. The pull-out force was calculated for the pristine graphene/epoxy, OH-FG/epoxy, and COOH-FG/epoxy systems. The pull-out force of OH-FG/epoxy system and COOH-FG/epoxy system were 5% and 54% higher than that of the graphene/epoxy system. The interfacial shear modulus of OH-FG/epoxy system and COOH-FG/epoxy system decreased 56% and 85%, respectively, compared to that of graphene/epoxy system. This implies that the graphene functionalization makes the interfacial region less stiff. The proposed combined MD-SDM method can be applied for various interfacial regions in nanocomposites between nanofillers and bulk matrix.

## Acknowledgment

This work was supported by National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP) and SNU Undergraduate Research Program. Authors are grateful for their support.

## Funding Data

National Research Foundation of Korea (2017R1A2B4004996).

Seoul National University (Undergraduate Research).