Temperature-Dependent Mechanical Properties of Graphene/Cu Nanocomposites with In-Plane Negative Poisson's Ratios

Negative Poisson's ratio (NPR), also known as “auxetic”, is a highly desired property in a wide range of future industry applications. By employing molecular dynamics (MD) simulation, metal matrix nanocomposites reinforced by graphene sheets are studied in this paper. In the simulation, single crystal copper with crystal orientation [1 1 0] is selected as the matrix and an embedded-atom method (EAM) potential is used to describe the interaction of copper atoms. An aligned graphene sheet is selected as reinforcement, and a hybrid potential, namely, the Erhart-Albe potential, is used for the interaction between a pair of carbon atoms. The interaction between the carbon atom and copper atom is approximated by the Lennard-Jones (L-J) potential. The simulation results showed that both graphene and copper matrix possess in-plane NPRs. The temperature-dependent mechanical properties of graphene/copper nanocomposites with in-plane NPRs are obtained for the first time.

Negative Poisson's ratio (NPR), also known as "auxetic", is a phenomenon associated with a material that experiences transverse expansion when a tensile load is applied in the longitudinal direction. On the other hand, if the material is under a compressive load, it contracts transversely. The NPR materials are a novel class of materials, which have advantages that conventional materials may not possess [22]. Materials with NPR have triggered significant interest in the scientific community because of their unique physical properties [23][24][25][26][27] including high shear modulus, high indentation resistance [28], large plane strain fracture toughness [29], and excellent shock and sound absorption abilities [30,31], and all of which are highly desirable properties for engineering applications. Specifically, the NPR materials can be used in microelectromechanical systems (MEMS) and fastening devices due to their auxetic characteristics, and they can be also added into engineering structures (e.g., sandwich or laminated panels in the aerospace industry) as energy absorption components. The natural materials found with NPR are polymers with microporous characteristics. Such polymer consists of a randomly or periodically distributed 3D array of particles interconnected by fibrils in the microscopic scale. For instance, microporous polytetrafluoroethylene (PTFE) film has a network of individual fibrils which can be extended by several microns in both directions when under uniaxial tension [32,33]. On the other hand, auxetic materials could be achieved via designing and tailoring their microscopic compositions and structures. Such man-made materials consist of a large number of unit cells, which are arranged into certain periodical truss patterns such as reentrant type, chiral type, or rotating units [32]. Recently, it has been reported that laminated composites reinforced with microfibers can gain out-of-plane NPR [34] through design and optimization of stacking sequences. However, there are still some limitations in engineering applications for the above NPR materials. For example, the majority of materials fabricated only have NPR in certain directions where the stiffness is low in those directions. Hence, such materials may only be used as the core for sandwiches which need to be enhanced by composite face sheets. In addition, NPR is usually achieved in the out-of-plane direction of a structural element rather than in the in-plane direction [35][36][37]. The research on composites with in-plane NPR is still in the conceptual stage [38,39] due to the unavailability of such type of NPR materials and their associated material properties for structural applications. Nevertheless, as hinted in the literature [40], there is a possible approach to obtain in-plane auxetic composites by using the matrix and the reinforcing phase that both possess the in-plane auxetic property.
In most polycrystalline metal materials, the value of Poisson's ratio is in the vicinity of 1/3. However, it is different for metal materials with single crystals where NPR may exist in certain crystallographic directions [41]. Early research [42] suggested that some single crystal metals including copper with a face-centered cubic (FCC) lattice might have an auxetic phenomenon in the ½1 1 0 direction when subjected to uniaxial ½1 1 0 loading. To date, the NPR for single-layer graphene has been observed by researchers using different approaches, including the Monte Carlo method [43], molecular dynamics and statics [44,45], molecular structural mechanics [46], first principle theory [47], and membrane theory [48].
In this study, the classical MD method using the LAMMPS package [49] is employed to construct the single crystal copper and the software package OVITO [50] is utilized to generate graphic presentations from the MD simulation results. For graphene, discrepancies exist in the MD studies on Poisson's ratio due to the selection of different potentials for the C-C bond. From existing literature [51], the absolute NPR (-0.158) can be obtained by MD simulations with the Tersoff potential [52] but cannot be achieved with the potentials based on a reactive empirical bond order (REBO) potential [53]. Hence, a hybrid potential [54] based on both REBO and Tersoff potentials is used to address this issue. Finally, the MD models of graphenereinforced copper composites are established and used for the prediction of the temperature-dependent NPRs and other material properties for the composites.

Single Crystal Copper.
Copper is an ideal substrate that supports growing graphene layers. In this study, we create a MD model of single crystal copper along the ½1 1 0, ½1 1 0, and ½0 0 1 crystal orientations as shown in Figure 1(a). The simulation box containing 5632 copper atoms is of 40:89 Å × 40:89 Å × 39:76 Å in dimensions. The boundary conditions are set as periodical in the three directions. An embedded-atom method (EAM) potential [55] developed by Dawn and Baskes [56] is employed to create the interaction between each pair of two copper atoms. The mathematical form of the general total energy for a pair of copper atoms is where F i ðρÞ is the embedding energy of atom i which is a function of the atomic electron density ρ, ϕ is the pair potential interaction, and α and β are the element types of atoms I and J. In the case of this study, α and β are the same. The lattice constant and the atomic weight of the FCC copper are selected as 3.615 Å and 63.55, respectively.
In the simulation, the tensile and shear behaviors controlled by strain are applied to this MD model to study its in-plane and out-of-plane mechanical properties. The temperature varying from 300 K to 1000 K is also considered in the simulation to observe the effect of temperature on the mechanical properties of the single crystal copper.

Single-Layer Graphene.
Graphene is a two-dimensional material. However, we cannot ignore its effective thickness in the mechanical analysis of Young's modulus and shear modulus. The MD model of the monolayer graphene with 201:68 Å × 200:21 Å is built, which is, respectively, referred to as the zig-zag direction and the armchair direction. A hybrid empirical bond order potential [54] for the C-C interaction based on REBO [53] and Tersoff [52] potentials is employed herein. The cohesive energy of this potential is written as a sum over individual bond energies: where V R and V Λ are the pairwise attractive and repulsive contributions, respectively, f C is the cutoff function, b is the bond order, and r is the distance between the nearest two carbon atoms. The detailed definitions of these physical quantities can be found in the literature [54]. The C-C bond length is well known as 1.42 Å, and the atom weight is 12.011. For the sake of convenience, this hybrid potential is named as the Erhart and Albe potential (EAP).
To determine the in-plane properties, such as E 11 , E 22 , G 12 , ν 12 , and ν 21 , the simulation tests in three directions are carried out. The force f is only applied on the carbon atoms closest to the boundary for each simulation, as shown in Figure 1(b). If we assume that there are M atoms in the zigzag direction and N atoms in the armchair direction, then we can first obtain stretching rigidities (C 11 and C 22 ) and shear rigidity (C 66 ) as: or where L x and L y are the origin lengths of the graphene in the zig-zag direction and the armchair direction, respectively, and Δ denotes the deformation caused by the force. It is worth noting that Equations (5) and (6) are valid when ΔL x or ΔL y is very small. There is divergence on the effective thickness of graphene which affects the values of its moduli. Usually, the graphite thickness is used to estimate the effective thickness of graphene. However, this method is not accurate for a mechanical concept. In this study, we treat the graphene sheet as a thin plate in the framework of mechanics and the effective thickness of the graphene sheet will be obtained after the central deflection of graphene is measured under a uniform load. Shen [57] has derived the relationship between the load q and the central deflection W based on the classical plate theory: where D 11 is the flexural rigidity of graphene in the zig- W are parameters and defined as where superscript G represents graphene. Finally, the effective thickness of the graphene sheet (h e ) can be derived from Equation (7): The detailed derivation for the load-deflection relationship in Equation (7) is presented in the appendix. In the simulation, the boundary atoms are fixed and a small transverse force is applied to the rest of the atoms.
Once the effective thickness is determined by MD simulations, the effective Young's modulus and shear modulus of the graphene sheet can be obtained uniquely from Equations (3) to (6). The TECs at different temperatures can also be obtained from the relaxation process. For the computation of the TEC, the standard relationship is used [58,59]: where α is the TEC with unit K -1 . T and T 0 are the simulation temperature and reference temperature, respectively. L x0 and L y0 are the lengths of graphene in the zig-zag direction and the armchair direction, respectively, corresponding to the reference temperature T 0 .

Graphene-Reinforced Copper
Nanocomposite. The MD model of the graphene-reinforced copper nanocomposite is illustrated in Figure 1(c). The EAM and EAP potentials are also used in the simulation of composite materials for the interactions of Cu-Cu and C-C. The Lennard-Jones (L-J) potential which is theoretically a good molecular dynamics model for long and short distances for neutral atoms and molecule is adapted for the C-Cu interaction. Moreover, the C-Cu interaction has already been demonstrated to meet the prediction from the Lennard-Jones potential in previous investigations [60,61]. This potential has also been widely used in existing literature [62][63][64] on the MD study of carbon/copper molecular structures. To estimate the interaction between copper and carbon atoms, the L-J potential is used: In Equation (18), r is the distance between two atoms. Obviously, this potential depends on the depth ε and the equilibrium interatomic distance σ. These parameters are determined to be ε = 0:01996 eV and σ = 3:225 Å based on existing studies [65]. On the other hand, the common divisor of the two atoms' lattice parameters is considered in the geometric design of the MD models for the composites to reduce the effect of mismatch, which will lead to internal stress in the MD models during the phases of relaxation and loading.
In the simulation, five different weight fractions (2.8%, 3.7%, 4.8%, 6.1%, and 7.1%) of the graphene sheet are taken into account. The temperature effect is also considered in the simulation. The time step is set as 1 fs, and the thermal relaxation consumes 40000 steps within the context of an isothermal-isobaric (NPT) ensemble. The tensile deformation of the simulation box is controlled by a strain rate of 10 -6 /fs, and the whole load process takes 100000 steps while it takes only 10000 steps for the shear deformation.

Mechanical Properties and NPR of Single Crystal
Copper. As a typical FCC metal, single crystal copper shows an in-plane auxetic property when the load is applied in the ½1 1 0 direction. Table 1 shows the inplane properties of the single crystal copper. From the table, in-plane NPRs (ν 12 and ν 21 ) are obtained through linear fitting MD simulation results and it is found that auxeticity is more significant with the increase in temperature. However, the elastic moduli (E 11 and E 22 ) and shear modulus (G 12 ) are decreased as temperature rises. We also find that the shear modulus is much smaller than the elastic moduli. The thermal expansion coefficient is also obtained and listed in Table 1. The out-of-plane properties of the single crystal copper, including E 33 , G 13 , and G 23 , are listed in Table 2. It is found that the out-of-plane elastic modulus is nearly half of the in-plane elastic modulus. The out-ofplane shear moduli G 13 and G 23 are almost four times of the in-plane shear modulus.  3.2. In-Plane Properties and Out-Of-Plane Bending Behavior of the Graphene Sheet. The variation of the stretching rigidities, the shear rigidity, the in-plane Poisson's ratios, and the thermal expansion coefficients of the graphene from 300 K to 1000 K is shown in Table 3. It can be seen that the stretching rigidities are decreased slightly as temperature increases. However, the trend of the shear rigidity is opposite to that of the stretching rigidities, which is in agreement with the MD results of Lin et al. [7]. Like the single crystal copper, the NPR of the graphene sheet is increased when temperature rises. However, the thermal expansion coefficients of the graphene sheet are negative, which means that graphene contracts when temperature increases. The effect of different loading directions on NPR has been investigated in our previous work [45] where the changes of NPRs with an increased load in the zigzag and armchair directions were also studied. We found that the gap between the NPRs in the zigzag and armchair directions is remarkable when the load is small but the gap is gradually reduced with the increase in the load. It was believed that the geometric configuration in different chiral directions plays a key role, because the distance between neighboring carbon atoms is increased once deformation of the graphene is occurred, which causes the reduction of the C-C bonds. With the increase in the interatomic distance, the effect of geometric configuration becomes weak, which results in the reduction of the gap between the NPRs in those two directions. The bending behavior is then performed to determine the effective thickness of graphene. As mentioned above, the small force is set to be 10 -4 nN and 1 million steps are taken for the balance of the whole system. The center deflection of graphene is calculated as the average value among transverse displacements of the nearest carbon atoms to the center point. In our study, the strain load of 10 -6 /fs is applied on the MD model and the averaged lattice parameters for the time step are converged to within 0.0001 nm. The strain load 10 -6 /fs with unit fs -1 means that we apply about 2 × 10 −5 nm as an incremental displacement per femtosecond on the MD model of graphene. The equilibrium state depends on the difference of lattice parameters within a time step. If the difference is lower than 0.0001 nm, the whole model can be regarded to be in an equilibrium state. The graphene is an atomicscale two-dimensional hexagonal lattice made by carbon atoms. In a dimensionally fixed graphene, the geometric center point must be in a hexagon cell. In our computa-tion, the transverse displacements of the six carbon atoms of the hexagon cell in the center area of graphene are averaged to obtain the deflection of the graphene for each time point. Applying Equation (15), we obtain the effective thickness of the graphene being 0.159 nm with the deflection of 5.784 Å at room temperature when the equivalent applied pressure q is 3:75 × 10 −3 GPa. Then, the density of the graphene can be obtained, which is 4717 kg/m 3 . Finally, the effective moduli of graphene can also be obtained based on Equations (3)-(6). Table 4 presents the elastic moduli E 11 and E 22 as well as the shear modulus G 12 for the graphene sheet.

Material Properties of Graphene/Cu Nanocomposites.
In this section, five different graphene weight fractions for the graphene-reinforced copper composites will be taken into consideration. The calculation of the graphene weight fractions can be achieved based on the number of carbon and copper atoms and their respective atom weight. Then, the volume fractions can be directly obtained by [66] where V and w denote the volume fraction and the weight fraction, respectively. ρ is the mass density. The superscripts or subscripts G and C represent graphene and copper, respectively. For example, one of our MD models for composites contains 112064 atoms in all, of which 14720 atoms are carbon and the rest are copper atoms. Accordingly, the obtained graphene weight fraction is 2.8%. Taking 8960 kg/m 3 as the mass density of copper, the volume fraction of graphene is then calculated to be 5%. The considered five weight fractions and their corresponding volume fractions are listed in Table 5.   Figure 2 illustrates the effect of graphene volume fraction on the in-plane Young's moduli and shear modulus of the composites at room temperature (300 K). It is worth noting that the 11 direction is along the zig-zag direction of the graphene as well as the ½1 1 0 crystal direction of the copper while the 22 direction is along the armchair direction of the graphene as well as the ½1 1 0 crystal direction of the copper. We found that E 11 is slightly larger than E 22 and they are both improved as the volume fraction of graphene is increased. We also notice that the rising trend of the moduli becomes slower when the graphene volume is higher than 7%. The out-of-plane moduli, E 33 , G 13 , and G 23 , with different V G are depicted in Figure 3. Although the Young's modulus in the thickness direction is enhanced by increasing the volume fraction of graphene, the incremental rate is relatively smaller when compared with that of the in-plane Young's moduli. It can be found that E 33 is only improved by 7% with 13% volume fraction of graphene when compared with that of the pure copper. Unlike the in-plane shear modulus, the out-of-plane shear moduli, G 13 and G 23 , are actually weakened by adding graphene reinforcement. We observe that the out-of-plane shear properties are degraded by as much as half when the composite contains only 5% volume fraction of graphene. The detailed moduli of the composites with different volume fractions of graphene under various temperatures are presented in Table 6.
In order to study the auxetic behavior of the graphenereinforced copper composites, the NPRs of the composite with 9% graphene volume fraction and its component materials are all shown in Figure 4. As mentioned above, the graphene sheet and single crystal copper are both inplane auxetic. Unlike in [21] where out-of-plane NPRs were obtained for graphene/Cu composites, in the current study, the in-plane NPRs of the resulting composite are obtained. From Figure 4, it can be found that the inplane NPRs (ν 12 and ν 21 ) of the composite are between those of graphene and copper but closer to the graphene NPRs. With increase in temperature, the NPRs of the composite are decreased and the NPRs can reach about -0.09 at 1000 K. The effect of graphene volume fraction on ν 12 of the composite is shown in Figure 5. Unlike the effect of temperature, the increase in graphene content may weaken auxetic characteristics of the composite. The TECs of the copper and graphene as well as the composite varying from 300 K to 1000 K are shown in Figure 6. Due to α 11 and α 22 for all three materials being very close, only their α 11 is compared in Figure 6. Obviously, the TEC of graphene is negative but the TEC of copper is positive. The TEC of the composite is still positive but lower than that of the copper. In other words, the addition of graphene with negative TECs decreases the TECs of the  copper based composite. Hence, internal stress may exist in the graphene-reinforced copper composites once temperature is varied. It is observed that the TEC trend of the composite is slightly different from that of graphene or copper due to the coupling effect. Table 7 shows detailed NPRs and TECs of composites from 300 K to Table 6: Temperature-dependent properties of graphene-reinforced copper composites (11 and 22 directions are, respectively, referred to as zig-zag and armchair directions, and the 33 direction is referred to as the thickness direction).  Figure 4: The effect of temperature on NPRs of the composite (V G = 9%) and its component materials. 1000 K. It is worth noting that the interaction between the carbon atom and copper atom is approximated by the L-J potential in the present MD models of graphenereinforced composites. Hence, the influence of graphene reinforcement on the mechanical properties and Poison's ratio of the composites mainly depends on the mechanical loading transfer at the interface of the graphene sheet and copper matrix.

Conclusion
This study not only presents novel MD results for MMCs reinforced by graphene sheets but also provides a new approach for the design of auxetic materials and structures. The graphene sheet and single crystal copper both showing in-plane auxeticity are selected as reinforcements and matrix, and their detailed material properties are obtained through MD simulations. To determine the effective material properties of graphene-reinforced metal matrix composites (GRMMCs), the MD simulations are carried out under different graphene volume fractions and temperature ranges from 300 K to 1000 K. The results showed that the in-plane NPRs of GRMMCs are decreased as the temperature rises or the graphene volume fraction increases and the maximum NPRs may reach -0.11 at 1000 K when the volume fraction of graphene V G is 5%. A graphene sheet with length a, width b, and effective thickness h e under a distributed transverse pressure q is considered. When the effective thickness is thin enough, the classical thin plate theory is adapted and the displacement field of an arbitrary point (X, Y) on the graphene can be expressed as where U, V, and W are midplane displacements in a Cartesian coordinate (X, Y, and Z) and t is the time. Then, the corresponding strains in the von Kármán type can be written as : The stress-strain relationship in a plane stress state can be expressed as   where Ω is the area of the plate. Using the variational principle, the von Kármán-type nonlinear equations for the thin plate can be derived as [57]  and Lð Þ contains the geometric nonlinearity terms in the von Kármán sense and can be given as

Conflicts of Interest
The authors declare that there is no conflict of interest.