Study on Lubrication Characteristics of C4-Alkane and Nanoparticle during Boundary Friction by Molecular Dynamics Simulation

: Lubricant has been widely applied to reduce wear and friction between the contact surfaces when they are in relative motion. In the current study, a nonequilibrium molecular dynamics (NEMD) simulation was speciﬁcally established to conduct a comprehensive investigation on the dynamic contact between two iron surfaces in a boundary friction system considering the mixed C4-alkane and nanoparticles as lubricant. The main research objective was to explore the effects of ﬂuid and nanoparticles addition on the surface contact and friction force. It was found that nanoparticles acted like ball bearings between the contact surfaces, leading to a change of sliding friction mode to rolling friction mode. Under normal loads, plastic deformation occurred at the top surface because nanoparticles were mainly supporting the normal load. By increasing the number of C4-alkane molecules between two contact surfaces, the contact condition has been changed from partial to full lubrication. In addition, an attractive force from the solid–liquid LJ interaction between C4-alkane and surfaces was observed at the early stage of sliding, due to the large space formed by wall surfaces and nanoparticles. The ﬁndings in this paper would be beneﬁcial for understanding the frictional behavior of a simple lubricant with or without nanoparticles addition in a small conﬁnement. on surfaces contact and friction force evolution. This study provides a fundamental understanding on the roles of mixed ﬂuid and nanoparticles in the boundary lubrication system that will be beneﬁcial indicates the system height evolution against the simulation time, and a comparison between the ﬂat and nanoparticle models is performed. It is evident from the ﬁgure that there is no obvious system height change for the ﬂat model under all normal loads. This is due to the high stiffness of pure Fe. During the subsequent sliding process, a slight ﬂuctuation of system height is observed, which is due to the same surface structures of both upper and lower walls. In contrast, much more stable system heights are observed in the nanoparticle models when the normal load varies between 0.25 and 0.75 GPa. However, when the normal load is increased up to 1.0 GPa, there is an obvious drop in the system height. Moreover, the drop becomes more severe with the sliding distance, which indicates that plastic deformation may occur at the top surfaces of both upper and lower walls. More details will be discussed in the later section of this paper.


Introduction
Recent experiments have found that adding nanoparticles to a lubricant has a significant influence on reducing wear and friction. The types of nanoparticles include metallic [1,2], metallic oxides [3][4][5] and non-metallic [6][7][8][9][10]. A comprehension discussion and summary on the roles of nanoparticles in oil lubrication has been carried out recently by Dai et al. [11]. They have compared the effects of chemical and physical properties of nanoparticles as well as their morphology and size on the lubricating performances. Due to the chemical composition of nanoparticles, a tribofilm was formed to protect the friction surface [12]. The size of nanoparticles showed visible effects on both friction and wear. The ideal particle size mainly depends on the working condition. Small nanoparticles can easily enter the contact interface, yielding a relatively low friction coefficient [13]. The main morphology of nanoparticles is spherical, followed by sheet, granular, onion and tube, which determines the lubrication conditions which are rolling, sliding, exfoliation, or a combination of all [14]. Moreover, nanoparticles can also affect lubricant flow behavior. A so-called plug flow occurring with nanoparticles and localized shear at surfaces were observed by Sperka and co-authors [15]. In addition, Ghaednia et al. [16] and Jackson et al. [17] have also studied the interactions between nanoparticles with the lubricant or sur-faces, and they both found a friction reduction mechanism separately when nanoparticles were present in the lubricants.
However, these mechanisms cannot be directly observed by experiments. Molecular dynamic (MD) simulations provide one of the most cost-effective and precise approaches to investigate molecular or atomistic level phenomenon [18][19][20]. To date, several simulations based on the molecular dynamics method have been conducted to investigate surface contact and the frictional characteristics of various lubricants and nanoparticles. There are models including single asperity [21], amorphous [22] or 2D rough surface [23], 3D rough surface contact [24], thin film model with flat surface [25] and 3D rough surface [26,27]. It should be pointed out that all these above-mentioned simulations with sliding walls are nonequilibrium molecular dynamics (NEMD) simulations in tribology [28]. According to previous studies [29][30][31][32][33][34][35], there is no doubt that NEMD simulations provide unique insights into understanding the nanoscale friction and lubrication mechanisms. For example, Lee et al. [29] studied the rolling resistance of a rigid Ni sphere of Cu substrate. The distribution of contact pressure during indentation was found to be mainly affected by the atomic terraces of the sphere. Their simulations also indicated that smaller nanoparticles and surface with adhesion showed a larger friction coefficient. Joly-Pottuz et al. [30] simulated a rolling-sliding behavior of carbon onions between DLC surfaces, which agreed well with TEM and STM experimental observations. Bucholz et al. [31] found that the interfacial bonds during friction determined the rolling and sliding condition. The NEMD simulations by Eder et al. [32] were about the abrasion process considering a rough iron surface having multiple hard abrasive particles. Ewen et al. [33] developed a model containing two different carbon nanoparticles between iron surfaces and studied the influence of nanoparticle type and coverage. They also discussed the influences of local pressure and velocity on friction. Shi et al. [34] studied the effect of nanoparticle shape on its frictional condition. Their simulation results revealed that the movement pattern turned into rolling from sliding when the nanoparticles were closer to a sphere. Su et al. [35] mixed water and nanodiamonds that had charges as lubricant confined by two gold surfaces. They claimed that the nanodiamonds with negative charges dramatically increased friction force due to the combination of rolling and cutting movement patterns.
However, the reports about the anti-wear mechanisms of fluid containing nanoparticles are still limited. Lv et al. [36] carried out MD simulations of Cu-Ar nanofluids confined within two solid surfaces. It has been found that the nanoparticles exerted a supporting force at the plates and further reduced the contact of two solid surfaces, weakening the friction effect. In addition, with an increase in pressure, nanoparticles were cut and absorbed into the solid surface that possessed a potential influence to fill for rough plates. Ji et al. [37] performed MD simulations to investigate the roles of nanoparticles on the sliding friction process. According to their simulation results, nanoparticles had a filling effect to smooth the contact surface and further to reduce the friction force. Hu et al. [38,39] studied the effect of nanoparticles by considering interactions with the base fluid. They observed a higher transition pressure due to the presence of nanoparticles. Moreover, they also compared the effect of the diamond and SiO 2 nanoparticles during sliding. Specifically, hard nanoparticles were found to polish the friction surfaces.
In the boundary and mixed lubrication conditions, direct metal/metal contacts often occur due to less liquid lubricants [40]. Therefore, it is very critical to explore the effect of fluid combining with nanoparticles on preventing direct surfaces contact. In our previous reports [25][26][27], we have systematically studied the influence of surface roughness during dry or lubricated contact conditions. The current work is an extension of these studies. NEMD simulations considering the lubricating conditions between two iron wall surfaces using the mixed C4-alkane and nanoparticles is specifically performed to investigate the influences of fluid and nanoparticle and associated mechanisms on surfaces contact and friction force evolution. This study provides a fundamental understanding on the roles of mixed fluid and nanoparticles in the boundary lubrication system that will be beneficial for understanding the frictional behavior of a simple lubricant with or without adding nanoparticles in a small confinement.

Materials and Model Setup
The MD simulation is performed through a large scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [41], which has been proven to be a very powerful MD simulator. As displayed in Figure 1, the geometry model consists of two <001> oriented BCC iron upper and lower walls, randomly distributed C4 molecules and 16 iron rigid nanoparticles with radius 8 Å, which leads to about 19% surface coverage by the nanoparticles. It should be noted that the use of iron for walls and nanoparticles is a simplification in this study. The iron wall is divided further into three regions: rigid region, thermostat region and free deformable region. Periodic boundary conditions have been defined along the X/Y lateral directions. In contrast, the Z normal direction is variable when the system height fluctuates; the system is compressed and sheared under 4 normal loads of 0.25, 0.5, 0.75 and 1.0 GPa. The lower and upper iron walls slide against each other with a velocity of 20 m/s in the X direction, which results in the lubricant shearing at 40 m/s. It should be stated that the relative motion velocity of 40 m/s is quite high compared to the practical applications and real experiments; however, it is mainly used in this study to model several sliding cycles within a reasonable simulation time. In the current investigation, different amounts of C4-alkane molecular are used to construct lubricant flow, which could provide a different lubricated condition, such as less or full lubrication. friction force evolution. This study provides a fundamental understanding on the roles of mixed fluid and nanoparticles in the boundary lubrication system that will be beneficial for understanding the frictional behavior of a simple lubricant with or without adding nanoparticles in a small confinement.

Materials and Model Setup
The MD simulation is performed through a large scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) [41], which has been proven to be a very powerful MD simulator. As displayed in Figure 1, the geometry model consists of two <001> oriented BCC iron upper and lower walls, randomly distributed C4 molecules and 16 iron rigid nanoparticles with radius 8 Å, which leads to about 19% surface coverage by the nanoparticles. It should be noted that the use of iron for walls and nanoparticles is a simplification in this study. The iron wall is divided further into three regions: rigid region, thermostat region and free deformable region. Periodic boundary conditions have been defined along the X/Y lateral directions. In contrast, the Z normal direction is variable when the system height fluctuates; the system is compressed and sheared under 4 normal loads of 0.25, 0.5, 0.75 and 1.0 GPa. The lower and upper iron walls slide against each other with a velocity of 20 m/s in the X direction, which results in the lubricant shearing at 40 m/s. It should be stated that the relative motion velocity of 40 m/s is quite high compared to the practical applications and real experiments; however, it is mainly used in this study to model several sliding cycles within a reasonable simulation time. In the current investigation, different amounts of C4-alkane molecular are used to construct lubricant flow, which could provide a different lubricated condition, such as less or full lubrication.

Force Field
In this work, there might be a severe plastic deformation that intensively depends on the proper force field, which could not be handled correctly through simple potentials including the Morse and Lennard-Jones potentials. Therefore, the embedded-atom method (EAM) force field designed for solid Fe was used in this work [42].

Force Field
In this work, there might be a severe plastic deformation that intensively depends on the proper force field, which could not be handled correctly through simple potentials including the Morse and Lennard-Jones potentials. Therefore, the embedded-atom method (EAM) force field designed for solid Fe was used in this work [42].
The C4-alkane is simulated by the TraPPE United Atom force field which utilizes pseudo-atoms to represent the CHx groups [43]. It should be noted that, TraPPE United Atom for field is not the best one to accurately describe the material properties compared to other advanced force fields, but it is still widely used from the consideration of computational efficiency [25][26][27][43][44][45]. Thus, we have applied the TraPPE United Atom force filed for atomistic simulations in order to reduce the computational cost, which is very important for NEMD simulations when there are a huge number of cases to run. The bond stretching of E bond , angle bending of E bend and torsion angle of E torsion have been calculated by Equations (3)-(5), respectively. The Lennard-Jones potential is used to model the interactions between the segment-pairs that are separated by more than three bonds as displayed in Equation (1). The interactions between surface-particle, liquid-particle, and particle-particle are also calculated by LJ potential. The parameters of unlike interactions were calculated by the Lorentz-Berthelot combining laws; see Equation (2). All the TraPPE-UA and LJ potential parameters are listed in Table 1.
39.0279464 1.54 * -LJ parameters of Fe are for surface-particle, liquid-particle and particle-particle interactions.

MD Simulation Procedures
The detailed model geometry has been introduced in the previous section. Each MD simulation is divided into three steps. The first step is dynamic relaxation with 200 ps. The atoms of both lower and upper rigid layers are constrained to maintain the system height, which allows these lubricant molecules to fully relax and the model to reach equilibrium with no variation in the energy of the system. In the second step, uniformly distributed loads of 0.25, 0.5, 0.75 and 1.0 GPa are added on the atoms belonging to the upper rigid layer. On the other hand, the atoms of the lower rigid layer are fixed along the normal direction, and therefore the lubricant is compressed. However, the rigid layers are limited within a maximum distance of 0.00025 Å every time step during the limited compression to achieve gentile touch between deformable layers and rigid particles. Within this stage, in order to maintain the temperature at 300 K, the Nose-Hoover thermostat [48] is applied on the atoms of thermostat layers with a damping coefficient of 100 fs. After the limited compression for 600 ps, the rigid layers are updated freely at each time step for the subsequent 200 ps. The last step is a sliding or shearing process. The atoms of upper rigid layer are sliding at a velocity of 20 m/s along the X direction for 3000 ps. The same sliding velocity is applied to the atoms of lower rigid layer but in an opposite direction, which results in the lower and upper walls sliding and lubricant shearing. During this step, the energy is dissipated via thermostat regions. In order to conduct the numerical integration of atomic classical equations of motions, the Velocity-Verlet algorithm is used in the present study, and the time step is set to 2 fs. The total simulation time was 4000 ps. For the post-processing of simulation data, both friction force and system height are saved every 4 ps, and a smooth window with size of 20 is used for figure plotting.

Results and Discussion
During the dynamic relaxation step, different numbers of C4-alkane from 500 to 4000 molecules are used as lubricants. Figure 2 exhibits a typical example that 1500 C4alkane molecules and 16 rigid particles are fully relaxed during the first 200 ps. All the nanoparticles are close to the lower wall. This is because their initial positions are closer to the lower wall than those of the upper wall. For simplicity, flat model, particle model and lubricated model are used in the subsequent sections to present the model without nanoparticles and C4-alkane molecules, the model only lubricated with nanoparticles, and the model lubricated with a mixture of nanoparticles and C4-alkane molecules, respectively.
Metals 2021, 11, x FOR PEER REVIEW 5 of 1 direction, and therefore the lubricant is compressed. However, the rigid layers are limite within a maximum distance of 0.00025 Å every time step during the limited compressio to achieve gentile touch between deformable layers and rigid particles. Within this stage in order to maintain the temperature at 300 K, the Nose-Hoover thermostat [48] is applie on the atoms of thermostat layers with a damping coefficient of 100 fs. After the limite compression for 600 ps, the rigid layers are updated freely at each time step for the sub sequent 200 ps. The last step is a sliding or shearing process. The atoms of upper rigi layer are sliding at a velocity of 20 m/s along the X direction for 3000 ps. The same slidin velocity is applied to the atoms of lower rigid layer but in an opposite direction, whic results in the lower and upper walls sliding and lubricant shearing. During this step, th energy is dissipated via thermostat regions. In order to conduct the numerical integratio of atomic classical equations of motions, the Velocity-Verlet algorithm is used in the pre sent study, and the time step is set to 2 fs. The total simulation time was 4000 ps. For th post-processing of simulation data, both friction force and system height are saved ever 4 ps, and a smooth window with size of 20 is used for figure plotting.

Results and Discussion
During the dynamic relaxation step, different numbers of C4-alkane from 500 to 400 molecules are used as lubricants. Figure 2 exhibits a typical example that 1500 C4-alkan molecules and 16 rigid particles are fully relaxed during the first 200 ps. All the nanopar ticles are close to the lower wall. This is because their initial positions are closer to th lower wall than those of the upper wall. For simplicity, flat model, particle model an lubricated model are used in the subsequent sections to present the model without nano particles and C4-alkane molecules, the model only lubricated with nanoparticles, and th model lubricated with a mixture of nanoparticles and C4-alkane molecules, respectively

Dry Contact
Aiming to study the influence of nanoparticles during boundary friction, studies o the flat and particle models are conducted first. After relaxation, different normal load varying from 0.25 to 1.0 GPa are applied on iron atoms of the upper rigid layer. Figure   Figure 2. Lubricated model with 1500 C4-alkane molecules in the relaxation step.

Dry Contact
Aiming to study the influence of nanoparticles during boundary friction, studies on the flat and particle models are conducted first. After relaxation, different normal loads varying from 0.25 to 1.0 GPa are applied on iron atoms of the upper rigid layer. Figure 3 Metals 2021, 11, 1464 6 of 16 indicates the system height evolution against the simulation time, and a comparison between the flat and nanoparticle models is performed. It is evident from the figure that there is no obvious system height change for the flat model under all normal loads. This is due to the high stiffness of pure Fe. During the subsequent sliding process, a slight fluctuation of system height is observed, which is due to the same surface structures of both upper and lower walls. In contrast, much more stable system heights are observed in the nanoparticle models when the normal load varies between 0.25 and 0.75 GPa. However, when the normal load is increased up to 1.0 GPa, there is an obvious drop in the system height. Moreover, the drop becomes more severe with the sliding distance, which indicates that plastic deformation may occur at the top surfaces of both upper and lower walls. More details will be discussed in the later section of this paper.
Metals 2021, 11, x FOR PEER REVIEW 6 of 16 indicates the system height evolution against the simulation time, and a comparison between the flat and nanoparticle models is performed. It is evident from the figure that there is no obvious system height change for the flat model under all normal loads. This is due to the high stiffness of pure Fe. During the subsequent sliding process, a slight fluctuation of system height is observed, which is due to the same surface structures of both upper and lower walls. In contrast, much more stable system heights are observed in the nanoparticle models when the normal load varies between 0.25 and 0.75 GPa. However, when the normal load is increased up to 1.0 GPa, there is an obvious drop in the system height. Moreover, the drop becomes more severe with the sliding distance, which indicates that plastic deformation may occur at the top surfaces of both upper and lower walls. More details will be discussed in the later section of this paper.  . Comparison of the system height evolution history against the simulation time under various normal loads between the flat and particle models.
Although the system height predicted by the flat models in Figure 3 remains stable, the corresponding simulated friction force is relatively large and fluctuates dramatically during the sliding process. Figure 4a,b reveals the influence of normal load on the evolution history of friction force against the simulation time obtained in different simulation systems. For flat models (Figure 4a), the upper and lower surfaces contact directly. There is an unusual friction force change when the simulation time is increased from 250 to 500 ps, as marked by the red dashed circle in the figure. This phenomenon is mainly due to the alignment pattern between the upper and lower walls. Similar jump in the friction force was also observed by Stephan et al. [49], in which the simulated tangential force increased as the chip started to form at the beginning of the lateral movement. On the other hand, An and co-authors [50] have successfully proposed a quantitative relationship of the nanoscale friction force and coefficient between liquid and solid based on their atomic force microscopy experimental characterizations. To better illustrate this effect, a typical system slice along the vertical direction is shown in Figure 5. At 300 ps, the upper wall normally approaches the lower wall. Due to the repulsive force between two deformable layers, the force Fx is increased. However, the force decreases immediately when two Figure 3. Comparison of the system height evolution history against the simulation time under various normal loads between the flat and particle models.
Although the system height predicted by the flat models in Figure 3 remains stable, the corresponding simulated friction force is relatively large and fluctuates dramatically during the sliding process. Figure 4a,b reveals the influence of normal load on the evolution history of friction force against the simulation time obtained in different simulation systems. For flat models (Figure 4a), the upper and lower surfaces contact directly. There is an unusual friction force change when the simulation time is increased from 250 to 500 ps, as marked by the red dashed circle in the figure. This phenomenon is mainly due to the alignment pattern between the upper and lower walls. Similar jump in the friction force was also observed by Stephan et al. [49], in which the simulated tangential force increased as the chip started to form at the beginning of the lateral movement. On the other hand, An and co-authors [50] have successfully proposed a quantitative relationship of the nanoscale friction force and coefficient between liquid and solid based on their atomic force microscopy experimental characterizations. To better illustrate this effect, a typical system slice along the vertical direction is shown in Figure 5. At 300 ps, the upper wall normally approaches the lower wall. Due to the repulsive force between two deformable layers, the force Fx is increased. However, the force decreases immediately when two surfaces have the shortest distance because the walls can move freely along both X and Y directions.
surfaces have the shortest distance because the walls can move freely along both X and Y directions.    In comparison with the results shown in Figure 4a by flat models, the particle models predict much smaller friction forces as displayed in Figure 4b for all the studied normal loads. The wall blocks separating from each other are observed when the hard nanoparticles are added in the boundary friction system. This agrees well with the observations by Tao et al. [51] that nanoparticles are able to divide the rubbing planes and thereby avoid direct surface contact. Specifically, while the normal load is ≤0.75 GPa, the friction force drops nearly to zero. This is likely owing to the rolling effect of nanoparticles when the normal load is very low [39]. Figure 6 reveals that the rolling effect of nanoparticles helps to decrease the friction force. For the normal load of 0.25 GPa, the nanoparticles are supporting the upper surface against the normal load but also smoothly rolling between the surfaces during the sliding process. The shape of atoms in gray color remains unchanged through the whole sliding process, suggesting the absence of plastic deformation in the friction pairs. When the normal load rises to 1.0 GPa, there is an obvious rise in the friction force as can be seen from Figure 4b. It means that the nanoparticles induce the plastic deformation in the top surface and thus lead to the increase of friction force. This deformation is shown in Figure 7. It should be noted that no nanoparticles are shown in the figure for clearness. The atoms of the top surface are deformed to accommodate the rigid nanoparticles, which is to ensure the rolling effect between surfaces. In comparison with the results shown in Figure 4a by flat models, the particle models predict much smaller friction forces as displayed in Figure 4b for all the studied normal loads. The wall blocks separating from each other are observed when the hard nanoparticles are added in the boundary friction system. This agrees well with the observations by Tao et al. [51] that nanoparticles are able to divide the rubbing planes and thereby avoid direct surface contact. Specifically, while the normal load is ≤0.75 GPa, the friction force drops nearly to zero. This is likely owing to the rolling effect of nanoparticles when the normal load is very low [39]. Figure 6 reveals that the rolling effect of nanoparticles helps to decrease the friction force. For the normal load of 0.25 GPa, the nanoparticles are supporting the upper surface against the normal load but also smoothly rolling between the surfaces during the sliding process. The shape of atoms in gray color remains unchanged through the whole sliding process, suggesting the absence of plastic deformation in the friction pairs. When the normal load rises to 1.0 GPa, there is an obvious rise in the friction force as can be seen from Figure 4b. It means that the nanoparticles induce the plastic deformation in the top surface and thus lead to the increase of friction force. This deformation is shown in Figure 7. It should be noted that no nanoparticles are shown in the figure for clearness. The atoms of the top surface are deformed to accommodate the rigid nanoparticles, which is to ensure the rolling effect between surfaces. In order to gain a better understanding on the plastic deformation under boundary friction condition, particle models with different numbers of nanoparticles are also carried out in this study. Specifically, four cases with 4, 8, 12 and 16 nanoparticles placed between the friction surfaces are considered. Moreover, only the normal load of 1.0 GPa is applied in these cases to ensure that the nanoparticles are able to penetrate the wall blocks. Compared with the flat model, the addition of 4 nanoparticles results in a slight increase of system height, but it remains stable during sliding as can be seen from Figure 8. Although additions of 12 and 16 nanoparticles lead to almost the same system height during compression step, the system height of the model containing 12 nanoparticles drops more quickly once the sliding starts. It has been found that fewer nanoparticles share larger compression stress and thus penetrate deeper into the wall blocks. This phenomenon becomes more obvious when the model contains 8 nanoparticles. In order to understand the influence of the number of nanoparticles, the friction force obtained in the model containing 4, 8, 12 or 16 nanoparticles is plotted in Figure 9, respectively. A comparison clearly indicates that more nanoparticles ensure the rolling effect and further reduce the friction force. In order to gain a better understanding on the plastic deformation under boundary friction condition, particle models with different numbers of nanoparticles are also carried out in this study. Specifically, four cases with 4, 8, 12 and 16 nanoparticles placed between the friction surfaces are considered. Moreover, only the normal load of 1.0 GPa is applied in these cases to ensure that the nanoparticles are able to penetrate the wall blocks. Compared with the flat model, the addition of 4 nanoparticles results in a slight increase of system height, but it remains stable during sliding as can be seen from Figure 8. Although additions of 12 and 16 nanoparticles lead to almost the same system height during compression step, the system height of the model containing 12 nanoparticles drops more quickly once the sliding starts. It has been found that fewer nanoparticles share larger compression stress and thus penetrate deeper into the wall blocks. This phenomenon becomes more obvious when the model contains 8 nanoparticles. In order to understand the influence of the number of nanoparticles, the friction force obtained in the model containing 4, 8, 12 or 16 nanoparticles is plotted in Figure 9, respectively. A comparison clearly indicates that more nanoparticles ensure the rolling effect and further reduce the friction force. In order to gain a better understanding on the plastic deformation under boundary friction condition, particle models with different numbers of nanoparticles are also carried out in this study. Specifically, four cases with 4, 8, 12 and 16 nanoparticles placed between the friction surfaces are considered. Moreover, only the normal load of 1.0 GPa is applied in these cases to ensure that the nanoparticles are able to penetrate the wall blocks. Compared with the flat model, the addition of 4 nanoparticles results in a slight increase of system height, but it remains stable during sliding as can be seen from Figure 8. Although additions of 12 and 16 nanoparticles lead to almost the same system height during compression step, the system height of the model containing 12 nanoparticles drops more quickly once the sliding starts. It has been found that fewer nanoparticles share larger compression stress and thus penetrate deeper into the wall blocks. This phenomenon becomes more obvious when the model contains 8 nanoparticles. In order to understand the influence of the number of nanoparticles, the friction force obtained in the model containing 4, 8, 12 or 16 nanoparticles is plotted in Figure 9, respectively. A comparison clearly indicates that more nanoparticles ensure the rolling effect and further reduce the friction force.   Figure 10 shows the surface morphology of the lower wall after sliding process. It is evident that grooves form during friction, which indicates that there is a cutting action from nanoparticles on surfaces. Such kind of grooves or scratches formation on the worn surface has often been experimentally observed in either nanoscale or microscale [52,53]. It has to be mentioned that only atoms of the lower wall with Z > 25 Å are shown in Figure  10 to clearly demonstrate the cutting effect. The corresponding color bar indicates the Z   Figure 10 shows the surface morphology of the lower wall after sliding process. It is evident that grooves form during friction, which indicates that there is a cutting action from nanoparticles on surfaces. Such kind of grooves or scratches formation on the worn surface has often been experimentally observed in either nanoscale or microscale [52,53] It has to be mentioned that only atoms of the lower wall with Z > 25 Å are shown in Figure  10 to clearly demonstrate the cutting effect. The corresponding color bar indicates the Z Figure 9. Comparison of the evolution of friction force against the simulation time between the flat model and particle models containing 4, 8, 12 and 16 nanoparticles. Figure 10 shows the surface morphology of the lower wall after sliding process. It is evident that grooves form during friction, which indicates that there is a cutting action from nanoparticles on surfaces. Such kind of grooves or scratches formation on the worn surface has often been experimentally observed in either nanoscale or microscale [52,53]. It has to be mentioned that only atoms of the lower wall with Z > 25 Å are shown in Figure 10 to clearly demonstrate the cutting effect. The corresponding color bar indicates the Z position of the atoms of the lower wall. As displayed in Figure 10a, all 4 nanoparticles deeply penetrate into the wall blacks under the normal load of 1.0 GPa. In addition, the atoms next to nanoparticles tend to accumulate and act like a cutting tool to remove the materials during the sliding process. Very similar results are also found in the model with 8 nanoparticles, as displayed in Figure 10b. By increasing the number of nanoparticles to 12 and 16, more grooves are formed with smaller height difference at block surfaces, as can be seen from Figure 10c,d, respectively. position of the atoms of the lower wall. As displayed in Figure 10a, all 4 nanoparticles deeply penetrate into the wall blacks under the normal load of 1.0 GPa. In addition, the atoms next to nanoparticles tend to accumulate and act like a cutting tool to remove the materials during the sliding process. Very similar results are also found in the model with 8 nanoparticles, as displayed in Figure 10b. By increasing the number of nanoparticles to 12 and 16, more grooves are formed with smaller height difference at block surfaces, as can be seen from Figure 10c,d, respectively.

Lubricated Contact
When lubricant is present in the boundary friction system, the wall surfaces are subjected to a combined influence of nanoparticles, adhesion as well as the lubricant, particularly while the average thickness of the lubricant is comparable to the surface roughness and the size of nanoparticles [23]. In this section, MD simulations considering the C4-alkane as the lubricant are carried out to investigate its influence on the boundary friction. Figure 11 displays the system height evolution for the lubricated models with C4alkane and 16 nanoparticles. Each system has an initial system height, which is apparently dependent on the number of C4-alkane molecules. However, all systems have the same height reduction rate because of the limited distance movement at each time step. This step could ensure the gentile touch between two wall surfaces and nanoparticles. Before the sliding process, it is important to guarantee the initial contact surface. Otherwise, the sharp contact could result in severe indentation on the wall surface, remaining in the following sliding process. From Figure 11, it has been found that C4-alkane molecules less than 1000 are unable to support the upper wall as the height does not change compared

Lubricated Contact
When lubricant is present in the boundary friction system, the wall surfaces are subjected to a combined influence of nanoparticles, adhesion as well as the lubricant, particularly while the average thickness of the lubricant is comparable to the surface roughness and the size of nanoparticles [23]. In this section, MD simulations considering the C4-alkane as the lubricant are carried out to investigate its influence on the boundary friction. Figure 11 displays the system height evolution for the lubricated models with C4alkane and 16 nanoparticles. Each system has an initial system height, which is apparently dependent on the number of C4-alkane molecules. However, all systems have the same height reduction rate because of the limited distance movement at each time step. This step could ensure the gentile touch between two wall surfaces and nanoparticles. Before the sliding process, it is important to guarantee the initial contact surface. Otherwise, the sharp contact could result in severe indentation on the wall surface, remaining in the following sliding process. From Figure 11, it has been found that C4-alkane molecules less than 1000 are unable to support the upper wall as the height does not change compared with the particle model. However, 2000 molecules of C4-alkane start supporting against the normal load as the height increases. More obvious increases could be observed in the systems which have 3000 or 4000 molecules. Normal load: 1.0 GPa Figure 11. Comparison of the system height evolution history against the simulation time between the particle and lubricated models under the normal load of 1.0 GPa. The number of C4-alkane molecules varies from 500 to 4000. Similar to the system height evolution, the friction force also varies when the number of lubricant molecules changes. From the results shown in Figure 12, it can be found that C4-alkane with 500 and 1000 molecules has nearly no contributions as a lubricant to decrease the friction force. However, when the number of molecules is greater than 1500 or above, an obvious change in the friction force is visible. It is due to the flow of lubricant starting to support the upper wall against the normal load. Therefore, the particles are subjected to less compression. It seems that the surfaces are only lubricated by fluid, and similar results were also observed in our previous study [25]. It is worth noting that, in addition to the system height, the nanoparticle-surface indentation depth is also a very important parameter which makes it possible to link and compare the simulation and experimental results [54,55]. Such kind of analysis will be conducted in our future work, where a series of new NEMD simulations considering different types of nanoparticles with 3D rough surfaces will be developed. Figure 11. Comparison of the system height evolution history against the simulation time between the particle and lubricated models under the normal load of 1.0 GPa. The number of C4-alkane molecules varies from 500 to 4000. Similar to the system height evolution, the friction force also varies when the number of lubricant molecules changes. From the results shown in Figure 12, it can be found that C4-alkane with 500 and 1000 molecules has nearly no contributions as a lubricant to decrease the friction force. However, when the number of molecules is greater than 1500 or above, an obvious change in the friction force is visible. It is due to the flow of lubricant starting to support the upper wall against the normal load. Therefore, the particles are subjected to less compression. It seems that the surfaces are only lubricated by fluid, and similar results were also observed in our previous study [25]. It is worth noting that, in addition to the system height, the nanoparticle-surface indentation depth is also a very important parameter which makes it possible to link and compare the simulation and experimental results [54,55]. Such kind of analysis will be conducted in our future work, where a series of new NEMD simulations considering different types of nanoparticles with 3D rough surfaces will be developed. To better explore the effect and mechanism of the lubricant fluid in the boundary lubrication system, the component of friction force corresponding to the systems having 500, 1000, 1500 and 2000 molecules is plotted in Figure 13, respectively. When the system only has 500 C4-alkane molecules mixed with nanoparticles, the supporting force is mainly contributed by the nanoparticles, as can be seen from Figure 13a. There is no direct supporting force from the lower wall as the gap distance is larger than the LJ cutoff distance. It is interesting that the force Z from C4-alkane decreases with the sliding process and even becomes negative after 1500 ps. It means that C4-alkane helps to compress the upper and lower walls. This is likely due to the large space formed by wall surfaces and nanoparticles. In addition, C4-alkane is free to flow and naturally the molecules tend to accumulate. The adsorption between C4-alkane and surfaces finally results in an attractive force between them. However, the negative force from C4-alkane disappears quickly when the flow is under compression or there is less space, which could be observed in Figure 13b. During the compression step and early stage of sliding, the supporting force from C4-alkane is near zero. It starts to increase at 2500 ps as the system height decreases slightly from then on. When 1500 molecules are mixed with nanoparticles, the force components remain stable during the sliding process. However, on the other hand, 2000 molecules exhibit a completely different influence on the supporting force. Figure 13d shows that nanoparticles exert an attractive force on wall surfaces, which results in a larger compression force on the lubricant flow. In this condition, nanoparticles could move freely within the lubricant flow. However, the nanoparticle atoms are still in the cut-off distance with surface atoms, which finally results in an attractive force. Therefore, these results reveal a great influence of the number of C4-alkane molecules on its lubricity within a boundary friction system and a transition from partial to full lubrication. To better explore the effect and mechanism of the lubricant fluid in the boundary lubrication system, the component of friction force corresponding to the systems having 500, 1000, 1500 and 2000 molecules is plotted in Figure 13, respectively. When the system only has 500 C4-alkane molecules mixed with nanoparticles, the supporting force is mainly contributed by the nanoparticles, as can be seen from Figure 13a. There is no direct supporting force from the lower wall as the gap distance is larger than the LJ cutoff distance. It is interesting that the force Z from C4-alkane decreases with the sliding process and even becomes negative after 1500 ps. It means that C4-alkane helps to compress the upper and lower walls. This is likely due to the large space formed by wall surfaces and nanoparticles. In addition, C4-alkane is free to flow and naturally the molecules tend to accumulate. The adsorption between C4-alkane and surfaces finally results in an attractive force between them. However, the negative force from C4-alkane disappears quickly when the flow is under compression or there is less space, which could be observed in Figure 13b. During the compression step and early stage of sliding, the supporting force from C4-alkane is near zero. It starts to increase at 2500 ps as the system height decreases slightly from then on. When 1500 molecules are mixed with nanoparticles, the force components remain stable during the sliding process. However, on the other hand, 2000 molecules exhibit a completely different influence on the supporting force. Figure 13d shows that nanoparticles exert an attractive force on wall surfaces, which results in a larger compression force on the lubricant flow. In this condition, nanoparticles could move freely within the lubricant flow. However, the nanoparticle atoms are still in the cut-off distance with surface atoms, which finally results in an attractive force. Therefore, these results reveal a great influence of the number of C4-alkane molecules on its lubricity within a boundary friction system and a transition from partial to full lubrication.

Conclusions
In the present study, a comprehensive NEMD simulation model considering two iron wall surfaces, mixed C4-alkane and nanoparticles as lubricant has been established to explore the influences of nanoparticles and fluid on the surface contact and friction force during boundary friction system. The simulation results revealed that nanoparticles were acting like ball bearings between contact surfaces, leading to a change from the sliding friction to rolling friction. When a large normal load greater than 0.75 GPa was applied, there was severe plastic deformation at the top surface. Moreover, C4-alkane with a number of molecules lower than 1500 was unable to support the upper wall, which had a similar friction behavior to that of the cases with only nanoparticles. In addition, there was an attractive force from C4-alkane molecules to surfaces, which resulted in an extra compression. This was due to the large space formed by wall surfaces and nanoparticles during the sliding process. However, with the increase of the number of molecules, the surfaces were gradually separated by the lubricant fluid from the mixed lubrication to totally thin film lubrication, and the friction force dropped dramatically.

Conclusions
In the present study, a comprehensive NEMD simulation model considering two iron wall surfaces, mixed C4-alkane and nanoparticles as lubricant has been established to explore the influences of nanoparticles and fluid on the surface contact and friction force during boundary friction system. The simulation results revealed that nanoparticles were acting like ball bearings between contact surfaces, leading to a change from the sliding friction to rolling friction. When a large normal load greater than 0.75 GPa was applied, there was severe plastic deformation at the top surface. Moreover, C4-alkane with a number of molecules lower than 1500 was unable to support the upper wall, which had a similar friction behavior to that of the cases with only nanoparticles. In addition, there was an attractive force from C4-alkane molecules to surfaces, which resulted in an extra compression. This was due to the large space formed by wall surfaces and nanoparticles during the sliding process. However, with the increase of the number of molecules, the surfaces were gradually separated by the lubricant fluid from the mixed lubrication to totally thin film lubrication, and the friction force dropped dramatically.