Molecular Dynamics Investigation of Graphene Nanoplate Diffusion Behavior in Poly-α-Olefin Lubricating Oil

Graphene as a type of novel additive significantly enhanced the tribological performance of blended lubricating oil. However, the dispersibility of graphene with long-term stability in lubricating oil is still a challenge. Chemical modification for graphene, rather than using surfactants, provided a better method to improve the dispersibility of graphene in lubricants. In this study, the equilibrium molecular dynamics (EMD) simulations were carried out to investigate the diffusion behavior of graphene nanoplates in poly-α-olefin (PAO) lubricating oil. The effects of graphene-size, edge-functionalization, temperature, and pressure on the diffusion coefficient were studied. In order to understand the influence of edge-functionalization, three different functional groups were grafted to the edge of graphene nanoplates: COOH, COON(CH3)2, CONH(CH2)8CH3 (termed GO, MG, and AG, respectively). The EMD simulations results demonstrated that the relationships between diffusion coefficient and graphene-size and number of functional groups were linear while the temperature and pressure had a nonlinear influence on the diffusion coefficient. It was found that the larger dimension and more functional groups provided the lower diffusion coefficient. AG with eight CONH(CH2)8CH3 groups exhibited the lowest diffusion coefficient. Furthermore, the experimental results and radial distribution function for graphene-PAO illustrated that the diffusion coefficient reflected the dispersibility of nanoparticles in nanofluids to some degree. To our best knowledge, this study is the first time the diffusion behavior of graphene in PAO lubricating oil was investigated using EMD simulations.

Over the last decades, graphene has attracted significant attention owing to its remarkable physical and chemical properties, including excellent thermal conductivity, high strength, and superior chemical inertness [16,17].Moreover, as a novel lubricating oil additive with good environmental compatibility, graphene has proved to enhance dramatically the tribological properties of the lubricating oils [18,19].The research studied by Eswaraish et al. [20] demonstrated that extreme pressure capacity, frictional coefficient and wear resistance were respectively improved by 40%, 80%, and 33% after adding graphene into the engine oil.Wu et al. [21] investigated the tribological performance of 4010 aviation lubricants containing graphene nanoplates.Compared to base oil, the friction coefficient and wear scar diameter were reduced by 27% and 43% with the addition of graphene.In addition, Mao et al. [22] carried out experiments to study the influence of the micromorphology of graphene on lubrication performance.The results illustrated that a graphene lubricant with regular edges has the best lubrication properties as friction reduction and anti-wear ability were improved by 27.9% and 14.1% of those of base oil.
However, as observed in our previous study [21], graphene was easy to form aggregates in lubricating oils, which was a challenge for the further applications of graphene as additives.In this case, Zheng et al. [23] added surfactant (span-80, C 24 H 44 O 6 ) into PAO4 lubricating oil to improve the dispersion stability of graphene.In other ways, the methods of chemical modification of graphene have been widely used to obtain stable blended nanofluids [24,25].Wu et al. [26] attached red phosphorus to graphene plates using a mechanical synthesis method and the mixed lubricating oil was uniformly dispersed for two months.Zhang et al. [27] developed a novel boehmite/graphene oxide nanohybrid and was observed to keep stable dispersion in lubricant base oil (VHVI8) at least 48 h.
Previous research on the dispersion stability of graphene focused on macroscopic experiments.However, limited efforts have been devoted to the investigations on the dispersion behavior of graphene from the microscopic perspective.In order to understand the microscopic behavior and obtain details that are inaccessible from experiments, molecular dynamics (MD) simulations were selected by many researchers [28][29][30][31][32][33][34].Saathoff et al. [35] explored how three different edge terminations affected the dispersion stability of graphene in n-methyl-pyrrolidone (NMP) via MD simulations.Konatham et al. [36] reported that functional groups grafted on the graphene plates were effective for the stable dispersion of graphene in liquid organic alkanes.Chen et al. [37] investigated the aggregation rate of two parallel graphene in dodecane.The results illustrated that graphene aggregated when parts of two graphene are in contact with each other.However, only a few efforts have been devoted to investigating the dispersibility of graphene in lubricating oils.Its detailed process from a microscopic perspective is still unclear, which limited the wide application of graphene as additives in the tribology area.In this study, the dispersibility of graphene was evaluated using a diffusion coefficient, similar to the method provided by Cha et al. [38].The EMD simulations were employed to calculate the diffusion coefficient of graphene-in-PAO4 systems to reveal the effects of graphene-size, functional groups, temperature, and pressure on dispersibility of graphene in PAO lubricant.

Materials and Simulation Method
In this work, all models were established using Materials Studio 8.0 (Accelrys Software Inc., San Diego, CA, USA).The consistent valence forcefield (CVFF) implemented in the software package Materials Studio was used for both PAO molecules and graphene.For the blended nanofluids, we used 1-decene trimer, the major component of PAO 4 [39], as the base lubricating oil.As shown in Figure 1, the graphene-in-PAO4 system consisted of one cubic cell of 6.12 nm in length (containing 20,300, 20,306, 20,322, and 20,362 atoms for using graphene (G), GO, MG and AG with side length of 1 nm, respectively) and contained a single graphene and 220 1-decene trimer molecules.Periodic boundary conditions were employed in the three directions.Four different graphene sizes with 1, 2, 3 and 4 nm were selected in the extraction, which resulted in the 0.57, 1.84, 3.15, and 7.54 wt.% concentrations of graphene in nanofluids.Graphene (G), along with its modified alternatives GO (G-COOH) [40], MG (G-COON(CH3) 2 ) [41], and AG (G-CONH(CH 2 ) 8 CH 3 ) [42] were employed to explore their diffusion behavior in PAO 4 lubricating oil, as displayed in Figure 2. Graphenes with 4 nm side length were used for investigating the effect of the type and number of functional groups, temperature and pressure on diffusion coefficient.The graphene was gradually modified to study the effect of number of functional groups (2, 4, 6, and 8 corresponding to 2.9, 5.8, 8.7, and 11.6%, respectively).
Crystals 2018, 8, x FOR PEER REVIEW 3 of 11 temperature and pressure on diffusion coefficient.The graphene was gradually modified to study the effect of number of functional groups (2, 4, 6, and 8 corresponding to 2.9, 5.8, 8.7, and 11.6%, respectively).All the EMD simulations were performed via Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [43] in this study.N (number of molecules), P (pressure), and T (temperature) ensemble for all temperatures were selected with a Nose-Hoover [44,45] style to retain the system temperature.The non-bonded interactions contained both electrostatic and van der Waals (vdW) interactions in the EMD simulations, where the time step length was 1 fs.The vdW interaction with temperature and pressure on diffusion coefficient.The graphene was gradually modified to study the effect of number of functional groups (2, 4, 6, and 8 corresponding to 2.9, 5.8, 8.7, and 11.6%, respectively).All the EMD simulations were performed via Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [43] in this study.N (number of molecules), P (pressure), and T (temperature) ensemble for all temperatures were selected with a Nose-Hoover [44,45] style to retain the system All the EMD simulations were performed via Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [43] in this study.N (number of molecules), P (pressure), and T (temperature) ensemble for all temperatures were selected with a Nose-Hoover [44,45] style to retain the system temperature.The non-bonded interactions contained both electrostatic and van der Waals (vdW) interactions in the EMD simulations, where the time step length was 1 fs.The vdW interaction with a 10 Å cut off was explained by the standard 12/6 Lennard-Jones potential.
where i is the depth of potential well, r represents the distance between two atoms, and σ is the distance where the interaction between atoms is equal to zero.The short range electrostatic interactions up to a distance of 10 Å were accounted for by applying the Coulombic law Here C is the energy conversion constant, and q i , q j represent the charges of two atoms, respectively.The long range electrostatic interactions beyond 10 Å were processed using the particle-mesh Ewald (PME) method with a precision of 10 −5 .In order to investigate the influence of temperature on diffusion behavior of graphene in PAO 4 lubricant, the simulations were performed over a temperature range of 300-400 K with a degree interval of 20 K. Energy minimization was conducted before the EMD simulations began.The system was equilibrated for 10 ns at an N (number of molecules), P (pressure), and T (temperature) ensemble until its energy fluctuations become stable, as shown in Figure 3.This provided the initial point for all simulations.After equilibration, the production steps of 10 ns were conducted at an N (number of molecules), V (volume), and T (temperature) ensemble and the time step was 1 fs.
stals 2018, 8, x FOR PEER REVIEW 4 of ere i is the depth of potential well, r represents the distance between two atoms, and σ distance where the interaction between atoms is equal to zero.The short range electrosta eractions up to a distance of 10 Å were accounted for by applying the Coulombic law Here C is the energy conversion constant, and i q , j q represent the charges of two atom pectively.The long range electrostatic interactions beyond 10 Å were processed using the particl sh Ewald (PME) method with a precision of 10 −5 .In order to investigate the influence perature on diffusion behavior of graphene in PAO 4 lubricant, the simulations were performe r a temperature range of 300-400 K with a degree interval of 20 K. Energy minimization w ducted before the EMD simulations began.The system was equilibrated for 10 ns at an N (numb olecules), P (pressure), and T (temperature) ensemble until its energy fluctuations become stab shown in Figure 3.This provided the initial point for all simulations.After equilibration, t duction steps of 10 ns were conducted at an N (number of molecules), V (volume), and perature) ensemble and the time step was 1 fs.

esults and Discussion
Graphene was found to easily aggregate and form sediment in lubricating oil.This proper ited the application of graphene as lubricants additives, although it was proved to significant

Results and Discussion
Graphene was found to easily aggregate and form sediment in lubricating oil.This property limited the application of graphene as lubricants additives, although it was proved to significantly enhance the tribological performance.In order to improve the dispersibility of graphene in PAO lubricating oil, it is necessary to quantitatively describe the influences of graphene-size, edge-functionalization, temperature and pressure on the dynamic properties of graphene in lubricating oil.The diffusion coefficient of graphene was calculated to describe its dynamic behavior in the nano range.The mean square displacements (MSDs) was another critical factor to describe the microscopic dynamics of nanoparticles, which reflects the rheological properties of nanofluids.Therefore, the motion of graphenes was analyzed by calculating the MSDs and diffusion coefficient (D) using Equations ( 3) and ( 4) [46,47].
Here N is the total number of atoms, r i (t) is the position of atom i at time t and r i (t 0 ) represents the original location of atom.The appropriate size of additives corresponded to its excellent dispersion stability.In this study, different dimensions of graphene with 1, 2, 3, and 4 nm side length were used in the simulations to study the effect of graphene-size on the diffusion coefficient.
Figure 4 showed the MSDs and diffusion coefficients of single graphene in the PAO 4 lubricant.The fluctuations of MSDs were very small, and indicated that the equilibrium time of 10 ns was a suitable time period, as shown in Figure 4a.The results showed that the smaller graphene moves more than the larger one.Figure 4b exhibited that the diffusion coefficient decreased with increasing graphene-size.Compared to 1 nm graphene-size, the diffusion coefficient of 4 nm was reduced by 88.18%.Similar results were also found in another study, which illustrated that the larger nanoparticles were more stable than that the smaller ones [38].
Crystals 2018, 8, x FOR PEER REVIEW 5 of 11 ( ) Here N is the total number of atoms, ( ) i r t is the position of atom i at time t and ( ) represents the original location of atom.The appropriate size of additives corresponded to its excellent dispersion stability.In this study, different dimensions of graphene with 1, 2, 3, and 4 nm side length were used in the simulations to study the effect of graphene-size on the diffusion coefficient.
Figure 4 showed the MSDs and diffusion coefficients of single graphene in the PAO 4 lubricant.The fluctuations of MSDs were very small, and indicated that the equilibrium time of 10 ns was a suitable time period, as shown in Figure 4a.The results showed that the smaller graphene moves more than the larger one.Figure 4b exhibited that the diffusion coefficient decreased with increasing graphene-size.Compared to 1 nm graphene-size, the diffusion coefficient of 4 nm was reduced by 88.18%.Similar results were also found in another study, which illustrated that the larger nanoparticles were more stable than that the smaller ones [38].Figure 5 exhibited the diffusion coefficients of four types of graphene containing two functional groups, which were shown in Figure 2, the original graphene (G), GO (G-COOH), MG (G-COON(CH3)2), and AG (G-CONH(CH2)8CH3) at 300 K and 0.1 MPa.It was clear that the diffusion coefficients decreased as the graphene-size became larger.When the graphene-size increased to 4 nm from 1 nm, the diffusion coefficients for GO, MG and AG were reduced by 83.75%, 84.91% and 84.04%, respectively.In addition, it was found that functional groups had a significant effect on the diffusion coefficient.Compared to G, the diffusion coefficients of GO, MG and AG were reduced by 48.18%, 70.91% and 91.45%, respectively when the graphene-size was 1 nm.This result demonstrated that AG, showing the smallest diffusion coefficient, was most stable in base oil.Similarly, Cha et al. [36] found that low speeds of graphene in liquids had less opportunity to aggregate.The diffusion Figure 5 exhibited the diffusion coefficients of four types of graphene containing two functional groups, which were shown in Figure 2, the original graphene (G), GO (G-COOH), MG (G-COON(CH 3 ) 2 ), and AG (G-CONH(CH 2 ) 8 CH 3 ) at 300 K and 0.1 MPa.It was clear that the diffusion coefficients decreased as the graphene-size became larger.When the graphene-size increased to 4 nm from 1 nm, the diffusion coefficients for GO, MG and AG were reduced by 83.75%, 84.91% and 84.04%, respectively.In addition, it was found that functional groups had a significant effect on the diffusion coefficient.Compared to G, the diffusion coefficients of GO, MG and AG were reduced by 48.18%, 70.91% and 91.45%, respectively when the graphene-size was 1 nm.This result demonstrated that AG, showing the smallest diffusion coefficient, was most stable in base oil.Similarly, Cha et al. [36] found that low speeds of graphene in liquids had less opportunity to aggregate.The diffusion coefficient of GO was smaller than that of G, which estimated that GO had the better dispersibility in lubricating oil and is consistent with the experimental result.As depicted in Figure 6, GO kept uniformly dispersed in lubricating oil for one week while the G exhibited complete sedimentation after setting for 100 h.In order to further understand the relationship between the diffusion coefficient and dispersibility of graphene in lubricating oil, the radial distribution functions (RDFs) for graphene-PAO molecules were computed at 300 K and 0.  In order to investigate the interaction between graphenes and PAO molecules, the RDFs between the centers of mass of graphenes and PAO molecules were calculated using graphenes with 4 nm side length.The RDFs correspond to the data obtained during 10 ns after equilibrium under NPT ensemble.As depicted in Figure 7, there are five sets of peaks when the distances between graphenes  In order to investigate the interaction between graphenes and PAO molecules, the RDFs betwe centers of mass of graphenes and PAO molecules were calculated using graphenes with 4 n In order to investigate the interaction between graphenes and PAO molecules, the RDFs between the centers of mass of graphenes and PAO molecules were calculated using graphenes with 4 nm side length.The RDFs correspond to the data obtained during 10 ns after equilibrium under NPT ensemble.As depicted in Figure 7, there are five sets of peaks when the distances between graphenes and PAO molecules are 1.325, 2.325, 2.575, 3.525 and 4.024 Å.The intensity becomes lower with the longer distance between graphene and PAO molecules, and the first set of peaks are the most important to evaluate the association between graphenes and PAO molecules.The higher peak at shorter separation was representative of the stronger association between graphene and PAO molecules; namely, stronger association represented better dispersibility [36].It was clear that AG exhibited the best dispersibility while G showed the weakest.Four types of graphene, G, GO, MG and AG were ranked according to their dispersibility in lubricating oil as: G < GO < MG < AG.
Crystals 2018, 8, x FOR PEER REVIEW 7 of 11 functional groups provided more intensive association between graphenes and PAO molecules, which allowed graphene to stably dispersed in oil.The influence of numbers of functional groups on diffusion coefficient was investigated using graphene with 1 nm of side length.The numbers of 2, 4, 6 and 8 of functional groups were grafted to the edge of graphene and the original graphene provided a baseline.Due to limitation of our computation ability, graphenes with side length of 4 nm were not selected in this study, although larger ones could reflect the dynamic properties better.Figure 8 shows that the diffusion coefficients of GO, MG, and AG linearly decreased with respect to increasing number of functional groups.In comparison with the original graphene, the reduction of diffusion coefficients for various graphenes containing different number of functional groups are listed in Table 1.The diffusion coefficients of GO, MG and AG containing eight functional groups were reduced by 87.27%, 96.36%, and 96.73%, respectively.The results illustrated that more functional groups resulted in better dispersibility for graphenes in lubricating oil.The reason might be attributed to the fact that graphene with more functional groups provided more intensive association between graphenes and PAO molecules, which allowed graphene to stably dispersed in oil.Different graphenes with side length of 4 nm were employed to investigate the effect of temperature and pressure on the diffusion coefficient of the four types of graphenes in oil.In our previous tribology tests [21], we found that the flash temperatures were higher than 300 but lower than 345 K when tribo-pairs were lubricated by oil containing graphene.In this case, temperatures ranging from 300 to 400 K were selected to study the diffusion coefficients of graphenes.Figure 9 showed that diffusion coefficients of graphenes increased with respect to increasing temperature, and  Different graphenes with side length of 4 nm were employed to investigate the effect of temperature and pressure on the diffusion coefficient of the four types of graphenes in oil.In our previous tribology tests [21], we found that the flash temperatures were higher than 300 but lower than 345 K when tribo-pairs were lubricated by oil containing graphene.In this case, temperatures ranging from 300 to 400 K were selected to study the diffusion coefficients of graphenes.Figure 9 showed that diffusion coefficients of graphenes increased with respect to increasing temperature, and AG exhibited the smallest diffusion coefficients at all simulated temperatures.In addition, the effect of temperature on diffusion coefficient was nonlinear, which was different from that of graphene-size and number of functional groups.The diffusion coefficients for G, GO, MG, and AG were improved by 310%, 326%, 383% and 422% when the temperature increased to 400 K from 300 K, which illustrated that AG was the most easily affected by temperature, and then MG, GO, and G.The results could be explained by functional groups being more easily affected by temperature than original graphene with stable structure.
The relationship between diffusion coefficient of graphenes and pressure was investigated under pressures ranging from 0.1 to 500 MPa, as shown in Figure 10.The effects of pressure on diffusion coefficients of graphenes were nonlinear, which was similar to temperature.However, the diffusion coefficients of graphenes decreased as the high pressure was applied in simulations.A similar result was obtained by Michalis et al. [48].Moreover, the original graphene exhibited the largest diffusion coefficients under all simulated pressures while AG the smallest.When the pressure was raised up to 500 MPa the diffusion coefficients for G, GO, MG, and AG were reduced by 98.7%, 97.4%, 95.1%, and 93.6%, respectively.
proved by 310%, 326%, 383% and 422% when the temperature increased to 400 K from 300 K hich illustrated that AG was the most easily affected by temperature, and then MG, GO, and G.The esults could be explained by functional groups being more easily affected by temperature than riginal graphene with stable structure.The relationship between diffusion coefficient of graphenes and pressure was investigated nder pressures ranging from 0.1 to 500 MPa, as shown in Figure 10.The effects of pressure on iffusion coefficients of graphenes were nonlinear, which was similar to temperature.However, the iffusion coefficients of graphenes decreased as the high pressure was applied in simulations.A imilar result was obtained by Michalis et al. [48].Moreover, the original graphene exhibited the rgest diffusion coefficients under all simulated pressures while AG the smallest.When the pressure as raised up to 500 MPa the diffusion coefficients for G, GO, MG, and AG were reduced by 98.7% 7.4%, 95.1%, and 93.6%, respectively.The relationship between diffusion coefficient of graphenes and pressure was investigated nder pressures ranging from 0.1 to 500 MPa, as shown in Figure 10.The effects of pressure on iffusion coefficients of graphenes were nonlinear, which was similar to temperature.However, the iffusion coefficients of graphenes decreased as the high pressure was applied in simulations.A imilar result was obtained by Michalis et al. [48].Moreover, the original graphene exhibited the argest diffusion coefficients under all simulated pressures while AG the smallest.When the pressure as raised up to 500 MPa the diffusion coefficients for G, GO, MG, and AG were reduced by 98.7% 7.4%, 95.1%, and 93.6%, respectively.

Figure 4 .
Figure 4. (a) MSD curves and (b) diffusion coefficients of graphene with various sizes at 300 K and 0.1 MPa.

Figure 4 .
Figure 4. (a) MSD curves and (b) diffusion coefficients of graphene with various sizes at 300 K and 0.1 MPa.

1 11 Figure 5 .
Figure 5. Diffusion coefficients of graphene with various functional groups calculated under NVT ensemble at 300 K and 0.1 MPa.

Figure 5 .Figure 5 .
Figure 5. Diffusion coefficients of graphene with various functional groups calculated under NVT ensemble at 300 K and 0.1 MPa.

Figure 7 .
Figure 7. Radial distribution function for graphene-PAO molecules at 300 K and 0.1 MPa from EMD simulations.

Figure 7 .
Figure 7. Radial distribution function for graphene-PAO molecules at 300 K and 0.1 MPa from EMD simulations.

Figure 8 .
Figure 8. Diffusion coefficients of graphene containing various number of functional groups calculated at NVT ensemble at 300 K and 1 atm.

Figure 8 .
Figure 8. Diffusion coefficients of graphene containing various number of functional groups calculated at NVT ensemble at 300 K and 1 atm.

Figure 9 .
Figure 9. Diffusion coefficients of graphene containing two functional groups with respect to increasing temperature at 0.1 MPa.

Figure 9 .
Figure 9. Diffusion coefficients of graphene containing two functional groups with respect to increasing temperature at 0.1 MPa.

Figure 9 .
Figure 9. Diffusion coefficients of graphene containing two functional groups with respect to increasing temperature at 0.1 MPa.

Figure 10 .
Figure 10.Diffusion coefficients of graphene containing two functional groups with respect to increasing pressures at 300 K.

Figure 10 .
Figure 10.Diffusion coefficients of graphene containing two functional groups with respect to increasing pressures at 300 K.

Table 1 .
The reductions in diffusion coefficients for various graphene containing different numbers of functional groups at 300 K and 0.1 MPa (%).

Table 1 .
The reductions in diffusion coefficients for various graphene containing different numbers of functional groups at 300 K and 0.1 MPa (%).