Concurrent Multiscale Simulations of Rough Lubricated Contact of Aluminum Single Crystal

: In this paper, a concurrent multiscale simulation strategy coupling atomistic and continuum models was proposed to investigate the three-dimensional contact responses of aluminum single crystal under both dry and lubricated conditions. The Hertz contact is performed by using both the multiscale and full molecular dynamics (MD) simulations for validation. From the contact area, kinetic energy and stress continuity aspects, the multiscale model shows good accuracy. It can also save at least ﬁve times the computational time compared with the full MD simulations for the same domain size. Furthermore, the results of lubricated contact show that the lubricant molecules could e ﬀ ectively cover the contact surfaces; thereby separating the aluminum surfaces and bearing the support loads. Moreover, the surface topography could be protected by the thin ﬁlm formed by the lubricant molecules. It has been found that the contact area decreases obviously with increasing the magnitude of load under both dry and lubricated contacts. Besides, a decrease in contact area is also seen when the number of lubricant molecules increases. The present study has conﬁrmed that the dimension of lubricated contacts could be greatly expanded during the simulation using the proposed multiscale method without sacriﬁcing too much computational time and accuracy.


Introduction
Tribology is an interdisciplinary field of research that studies the adhesion, friction, lubrication and wear of contacts from the nano-scale to macro-scale, which is quite complicated, as the processes are usually governed by physical mechanisms at different length scales and can be affected by many factors, such as load, temperature, velocity, environment and so on [1][2][3][4][5][6][7]. For contact mechanics, in particular, not only the long-range elastic and plastic deformation play roles, but also the atomic scale roughness at the interface, reported by Luan and Robbins, has critical effects on the contact area and stress [4]. They revealed the breakdown of continuum models, and intrigued new investigations for a number of contacting problems at the nanoscale, which are still open for in-depth analyses. Over the last 40 years, the molecular dynamics (MD) simulations were widely used to explore the mechanisms of contact, friction, wear and lubrication at the nanoscale for various materials [8][9][10][11][12]. As a complement to conventional experiments, MD simulation can produce some detailed observations which cannot be obtained by experimental instruments. Despite the advancing development of computer science, MD still suffers from the temporal and spatial limitations. It is still very challenging to slow the timescale of MD down to realistic vales, although many efforts have been devoted to approaching such a target, such as accelerated molecular dynamics [13] and multiscale simulations [14][15][16]. The finite element method (FEM), based on the framework of continuum mechanics, has also been widely using the random midpoint displacement (RMD) algorithm, following the previous MD simulations [34,35]. The corresponding RMS roughness of the generated atomic surface was 0.204, 0.492 and 0.771 nm, respectively. The atomic region of the substrate consists of three layers: (2) the deformable layer, (4) the thermostat layer and (3) the pad layer. The contact occurs among the rigid plane, lubricant and deformable layer. The thermostatting technique was applied to the atoms in the thermostat layer to maintain the constant temperature of 300 K for the system. Along the x and y directions, a periodic boundary condition was applied to both the atomic and FEM regions. The bottom boundary of the FEM region was fixed while the top surface of the deformable layer was free. The FEM region consists of 33,600 elements and 6615 nodes, corresponding to 1,433,600 atoms. A constitutive anisotropic law describes the relation between stress and strain for each integration point. The elastic constants C11 = 114.0, C12 = 61.6, C44 = 31.6 (unit: GPa) were taken from Mishin's work [36]. The key ideas of the coupling strategy for the HSM is described here. The pad layer as a part of the FEM region provides a physical boundary for the atomic region. The displacements of each atom (e.g., uβ) in the pad layer are interpolated from the nodal displacements [26] as shown in Equation (1): where N is the shape function of atom β corresponding to the nodes in each element. mnodes is the number of the nodes in each element. As can be seen in Figure 1, the top boundary of the FEM region lies in the thermostat layer. The nodal displacements on this boundary (e.g., di) were obtained by The key ideas of the coupling strategy for the HSM is described here. The pad layer as a part of the FEM region provides a physical boundary for the atomic region. The displacements of each atom (e.g., u β ) in the pad layer are interpolated from the nodal displacements [26] as shown in Equation (1): where N is the shape function of atom β corresponding to the nodes in each element. m nodes is the number of the nodes in each element. As can be seen in Figure 1, the top boundary of the FEM region lies in the thermostat layer. The nodal displacements on this boundary (e.g., d i ) were obtained by averaging the displacements of atoms within the average spheres, as shown in Equation (2). The radius of each sphere is 4a, corresponding to the smallest size of the element edge. The four-node tetrahedron elements were used to divide the continuum region: where n atoms is the number of atoms in the average circle for node i. The weight ω of each atom within the averaging circle with a radius of r av can be calculated by a weighting function [37]: where d = r/r av and r is the distance between the atom within the average circle and the node. The timestep is 2 fs. Before loading, the substrate is allowed to relax for equilibrium. The system has 0.4 ns to equilibrate after applying a load. To study the effects of different pressure on the lubricated contact, 0.05, 0.15 and 0.25 GPa were applied on the rigid plane, respectively. During the loading, at an interval of 200 timesteps, the forces and stresses were recorded. The calculation time of the multiscale simulation ranges from 74 to 98 h due to the different number of atoms. For the different cases here, the number of FEM nodes are not changed.

Force Field
As the plastic deformation on the rough surface with the atoms reshuffling or dislocations usually involves during compression and shearing between the rough surfaces, the EAM potential [38] was used to model the interactions between the aluminum atoms in the substrate. The atomic interactions between the rigid plane and the substrate were calculated by the Lennard-Jones (LJ) potential, which is characterized by the distance and energy parameters σ and ε, respectively. The LJ potential was used to simulate the aluminum oxide surface by reducing ε [34]. A united-atom (UA) model [39] was used to model the hexadecanes molecules. The UA model simplified all the CH X groups with pseudo carbon atoms. Linear chain molecules of hexadecane with the chemical formula C 16 H 34 were chosen as the lubricant liquid. A chain of hexadecane is made of connected CH 3 and CH 2 groups as shown in Figure 2a. The optimization of the conformation of the hexadecane molecule in Figure 2b is attributed to the intra-molecular interactions, which include bond stretching, angle bending and the dihedral angle torsion given in Equations (1)- (3). To study the effects of the amount of lubricant on the lubricated contact, the different number of C 16 H 34 molecules was added (0-4000) between the aluminum tribo-surfaces: averaging the displacements of atoms within the average spheres, as shown in Equation (2). The radius of each sphere is 4a, corresponding to the smallest size of the element edge. The four-node tetrahedron elements were used to divide the continuum region: where natoms is the number of atoms in the average circle for node i. The weight ω of each atom within the averaging circle with a radius of rav can be calculated by a weighting function [37]: where d = r/rav and r is the distance between the atom within the average circle and the node. The timestep is 2 fs. Before loading, the substrate is allowed to relax for equilibrium. The system has 0.4 ns to equilibrate after applying a load. To study the effects of different pressure on the lubricated contact, 0.05, 0.15 and 0.25 GPa were applied on the rigid plane, respectively. During the loading, at an interval of 200 timesteps, the forces and stresses were recorded. The calculation time of the multiscale simulation ranges from 74 to 98 h due to the different number of atoms. For the different cases here, the number of FEM nodes are not changed.

Force Field
As the plastic deformation on the rough surface with the atoms reshuffling or dislocations usually involves during compression and shearing between the rough surfaces, the EAM potential [38] was used to model the interactions between the aluminum atoms in the substrate. The atomic interactions between the rigid plane and the substrate were calculated by the Lennard-Jones (LJ) potential, which is characterized by the distance and energy parameters σ and ε, respectively. The LJ potential was used to simulate the aluminum oxide surface by reducing ε [34]. A united-atom (UA) model [39] was used to model the hexadecanes molecules. The UA model simplified all the CHX groups with pseudo carbon atoms. Linear chain molecules of hexadecane with the chemical formula C16H34 were chosen as the lubricant liquid. A chain of hexadecane is made of connected CH3 and CH2 groups as shown in Figure  2a. The optimization of the conformation of the hexadecane molecule in Figure 2b is attributed to the intra-molecular interactions, which include bond stretching, angle bending and the dihedral angle torsion given in Equations (1)- (3). To study the effects of the amount of lubricant on the lubricated contact, the different number of C16H34 molecules was added (0-4000) between the aluminum tribosurfaces: Metals 2020, 10, 965

of 14
The corresponding parameters in Equations (4)-(6) are taken from Reference [39]. The inter-molecular and liquid-solid interactions are described by the LJ potential with a cutoff distance of 1.4 nm. The LJ parameters for the liquid-solid interactions were estimated via the Lorentz-Berthelot mixing rules (Table 1).

Contact Area
Calculations for the real contact area are of importance at the atomic scale, especially when the dimension of the devices is reduced down to the nanometer. The real contact area has an important effect on heat conduction, friction and adhesion. In continuum mechanics, the contact area can be obtained based on the edge of the contact zone. It is difficult to accurately define the atomic area and even the atomic contact. The atomic contact is determined by either repulsive force [28,40] or distance [35,41]. Here, the contact is established when the distance of two atoms is less than or equals to 0.5 nm [30], and the real contact area is estimated from a projected rectangular cuboid method [34]. In this method, the n × n grids as a bottom face of the rectangular cuboids are plotted on the x-y plane, on which all the atoms of the deformable layer are projected. Those atoms are divided into the corresponding grids according to their x and y coordinates. The length of each cuboid is the maximum z coordinate values of the atoms within it. Similar cuboids can be plotted for the rigid plane. The top faces of the deformable layer cuboids and the bottom faces of the rigid plane cuboids are used to determine the atomic contact.

Multiscale Strategy Valication
Here, we validate our MU model by performing a 3D Hertz contact, benchmarked with full MD simulations, following the the Anciaux and Molinali's work [32]. Details of the model setup for the Hertz contact is shown in Figure S1 and demonstrated in the supplementary materials. The contact area, kinetic energy and stress continuality were considered. Figure S2a shows the relationship between the contact area and the normalized indent displacement. The contact region consists of the atoms on the stepped tip and is close to the area of a circle. The contact area was determined by averaging the distances between the center of the circle and the outermost atoms. In terms of the contact area, the MU predicts the MD well. It is found from Figure S2a that the results between the two models are very close. From the energy point of view, the MU works well based on the results in Figure S2b. Although there are a few deviating points, the kinetic energies of MU and MD are close to each other. In the MU model, the atomic, or the virial stress was calculated for the MD region while the Cauchy stress for the FEM region. The stress component σ zz of both the MU model and the MD model when normalized indent displacement is 1.31a was shown in Figure S2c. First and foremost, the HSM can ensure stress continuity after comparing σ zz of both the FEM region and the MD region (dotted square). Another important result shown in Figure S2c is that the MU model can predict the full MD model well.
From the contact area, kinetic energy and stress continuity aspects, our MU simulations show good accuracy compared with the full MD simulations. More importantly, it can save at least 5 times the computational time compared with the full MD simulations for the same domain size. With relatively good computational cost and accuracy, the MU model was here extended to study lubricated contacts.

RMS Roughness
Before compression, the rigid plane was above the substrate with a distance larger than the cut-off radius of the interaction potential. As mentioned above, the original RMS roughness of the generated atomic surface was 0.204, 0.492 and 0.771 nm, respectively. The atoms at the surface lack the full neighbors, and are attracted by the atoms at the subsurface. Therefore, the relaxing process at a temperature of 300 K modified the surface roughness with RMS = 0.198, 0.484 and 0.760 nm, respectively. Although the RMS roughness of the surfaces changes, it is still close to the required values. After relaxation, the top plane was forced to move downwards under constant pressure. Without any resisting force, the top plane will impact on the substrate with a high speed. This is not acceptable, since it will cause significant plastic deformation on the substrate. Therefore, the top plane was forced to move at a very low speed until it touched the deformable layer. The system heights were calculated during the whole compression. In Figure 3, the relationship between the system height and simulation time for dry and lubricated contact is shown. It indicates that all the simulations reached equilibrium after 0.4 ns. The heights in the lubricated cases were higher than the dry counterparts, as shown in Figure 3a,b. The thickness of the thin film formed by the hexadecane molecules contributed to the system height. In Figure 3c, two thousand hexadecane molecules had no effect on the system height when the surface roughness was 0.8 nm (indicated by the dashed circle in the Figure), since there was not enough quantity to fill the cavity between the two contacting surfaces. Meanwhile, four thousand hexadecane molecules lifted the top plane by a distance of 0.4 nm as indicated by the dashed arrow in Figure 3c. The system height is also related to the chain length of the lubricant molecules [42]. The effective separation for the contacting surfaces is dependent on the added amount of lubricant and surface roughness. A few lubricant molecules can work well for relatively flat surfaces, while much more molecules are required to fill a cavity in rougher surfaces and effectively separate them. The volume of the hexadecane molecule is also related to the rotation of C-C bonds-thereby affecting the system heights in Figure 3. From the contact area, kinetic energy and stress continuity aspects, our MU simulations show good accuracy compared with the full MD simulations. More importantly, it can save at least 5 times the computational time compared with the full MD simulations for the same domain size. With relatively good computational cost and accuracy, the MU model was here extended to study lubricated contacts.

RMS Roughness
Before compression, the rigid plane was above the substrate with a distance larger than the cutoff radius of the interaction potential. As mentioned above, the original RMS roughness of the generated atomic surface was 0.204, 0.492 and 0.771 nm, respectively. The atoms at the surface lack the full neighbors, and are attracted by the atoms at the subsurface. Therefore, the relaxing process at a temperature of 300 K modified the surface roughness with RMS = 0.198, 0.484 and 0.760 nm, respectively. Although the RMS roughness of the surfaces changes, it is still close to the required values. After relaxation, the top plane was forced to move downwards under constant pressure. Without any resisting force, the top plane will impact on the substrate with a high speed. This is not acceptable, since it will cause significant plastic deformation on the substrate. Therefore, the top plane was forced to move at a very low speed until it touched the deformable layer. The system heights were calculated during the whole compression. In Figure 3, the relationship between the system height and simulation time for dry and lubricated contact is shown. It indicates that all the simulations reached equilibrium after 0.4 ns. The heights in the lubricated cases were higher than the dry counterparts, as shown in Figure 3a, b. The thickness of the thin film formed by the hexadecane molecules contributed to the system height. In Figure 3c, two thousand hexadecane molecules had no effect on the system height when the surface roughness was 0.8 nm (indicated by the dashed circle in the Figure), since there was not enough quantity to fill the cavity between the two contacting surfaces. Meanwhile, four thousand hexadecane molecules lifted the top plane by a distance of 0.4 nm as indicated by the dashed arrow in Figure 3c. The system height is also related to the chain length of the lubricant molecules [42]. The effective separation for the contacting surfaces is dependent on the added amount of lubricant and surface roughness. A few lubricant molecules can work well for relatively flat surfaces, while much more molecules are required to fill a cavity in rougher surfaces and effectively separate them. The volume of the hexadecane molecule is also related to the rotation of C-C bonds-thereby affecting the system heights in Figure 3. When the contact occurs, the top plane will flatten the substrate asperities as shown in Figure  4b, c. The RMS roughness of flattened surface changes as the load increases. The RMS of the flattened surface is denoted as RMSf, while that of the initial surface RMSu. The change of RMS, therefore, is defined as ∆RMS = RMSf − RMSu. Firstly, it indicates that RMSf decreases with the increasing loads for three surfaces, and this relationship is independent of RMSu. However, the percentage of change in RMS roughness is related to RMSu. There was a much smaller change for the rougher surfaces When the contact occurs, the top plane will flatten the substrate asperities as shown in Figure 4b,c. The RMS roughness of flattened surface changes as the load increases. The RMS of the flattened surface is denoted as RMS f , while that of the initial surface RMS u . The change of RMS, therefore, is defined as ∆RMS = RMS f − RMS u . Firstly, it indicates that RMS f decreases with the increasing loads for three surfaces, and this relationship is independent of RMS u . However, the percentage of change in RMS roughness is related to RMS u . There was a much smaller change for the rougher surfaces (solid triangle and square in Figure 4a), because the contact for those rougher surfaces occurred on a few high asperities which support the loads. It is the rest part on the surface that dominated the roughness. However, the surface RMS changed more for RMS = 0.2 nm, because much more asperities on the substrate contacted with the top plane were simultaneously flattened. The peak-to-valley (PV) distances [34] indicate roughly the number of the contacting asperities. The PV distance here is determined by the difference between the maximum peak height to the maximum valley depth of the atomistic surface. The PV of the atomic surfaces for RMS = 0.2, 0.5 and 0.8 nm are 0.172, 0.3285 and 0.514 nm, respectively. The height of the highest asperity for RMS = 0.2 nm is approximately 0.86 nm. Therefore, many asperities support the top plane simultaneously. Figure 5 shows the effects of the molecule number on the RMS change. Zheng et al. [42] demonstrated that the appropriate amount of lubricant is dependent on the RMS roughness. The increase in the molecule number reduces the change in RMS roughness or asperity flattening as shown in Figure 5. After compression, the lubricant molecules filled the valleys and covered the surface as shown in the insets. With 4000 molecules, the surface could be entirely covered. The asperities were covered by the lubricant, and were not contacted by the top plane. Consequently, the RMS changed only by about 1.8%. In addition, Wu et al. [43] indicated that the lubricant viscosity has an important effect on the changes in the surface roughness. In their work, the smaller reduction in the surface roughness was obtained for the lubricant with a higher viscosity. The viscosity in Wu's work plays a similar role as the lubricant amount here. A well-formed lubricant film reduced the direct metal-metal contact, and protected the surfaces. For the lubricated contact, the new surface was formed by the lubricant molecules and substrate atoms. The RMS roughness of the new surface was 0.33 and 0.06 nm for 3000 and 4000 molecules, respectively. They are both much lower than the original roughness of 0.8 nm.
Metals 2020, 10, x FOR PEER REVIEW 7 of 14 (solid triangle and square in Figure 4a), because the contact for those rougher surfaces occurred on a few high asperities which support the loads. It is the rest part on the surface that dominated the roughness. However, the surface RMS changed more for RMS = 0.2 nm, because much more asperities on the substrate contacted with the top plane were simultaneously flattened. The peak-to-valley (PV) distances [34] indicate roughly the number of the contacting asperities. The PV distance here is determined by the difference between the maximum peak height to the maximum valley depth of the atomistic surface. The PV of the atomic surfaces for RMS = 0.2, 0.5 and 0.8 nm are 0.172, 0.3285 and 0.514 nm, respectively. The height of the highest asperity for RMS = 0.2 nm is approximately 0.86 nm. Therefore, many asperities support the top plane simultaneously. Figure 5 shows the effects of the molecule number on the RMS change. Zheng et al. [42] demonstrated that the appropriate amount of lubricant is dependent on the RMS roughness. The increase in the molecule number reduces the change in RMS roughness or asperity flattening as shown in Figure 5. After compression, the lubricant molecules filled the valleys and covered the surface as shown in the insets. With 4000 molecules, the surface could be entirely covered. The asperities were covered by the lubricant, and were not contacted by the top plane. Consequently, the RMS changed only by about 1.8%. In addition, Wu et al. [43] indicated that the lubricant viscosity has an important effect on the changes in the surface roughness. In their work, the smaller reduction in the surface roughness was obtained for the lubricant with a higher viscosity. The viscosity in Wu's work plays a similar role as the lubricant amount here. A well-formed lubricant film reduced the direct metal-metal contact, and protected the surfaces. For the lubricated contact, the new surface was formed by the lubricant molecules and substrate atoms. The RMS roughness of the new surface was 0.33 and 0.06 nm for 3000 and 4000 molecules, respectively. They are both much lower than the original roughness of 0.8 nm.

Contact Area
The effects of roughness and lubricant on the contact area calculated by the project method will be studied in this section. Figure 4a indicates that the contact area goes up with the increasing loads for the surfaces with three different RMS. Under the load of 0.05 GPa, the contact area of RMS = 0.2 nm is close to 50%, which is much higher than 4% and 3.5% of the RMS = 0.5 and 0.8 nm. As mentioned in the previous section, the heights of the asperities on the surface of RMS = 0.2 nm are relatively low, so a number of them were in contact simultaneously. For the rougher surface, however, only few asperities supported the load applied from the top plane. With the increasing loads, the contacting asperity deformed and their heights decreased while the lower non-contacting asperities were joined to support the top plane. Therefore, the contact area of the three surfaces was increased in Figure 4a. Considering the different topography due to the different random seed when generating the selfaffine surfaces [35] and adhesive strength [34], the results could deviate a little. The relationship between the contact area and the loads obtained in this work qualitatively agrees with the published literature [44].
For the lubricated contacts, the contact area still increased with the loads at the range of 0.05 GPa (red diamond) to 0.25 GPa (purple square) as shown in Figure 6a. However, the contact area reduced with the increasing number of molecules. The lubricants could separate the substrate from the top plane. The separating gap depends on the lubricant amount. It indicates that under 0.25 GPa the contact area of 0 and 1000 molecules is almost identical while the contact area of 4000 molecules was reduced by 75%. However, it was not enough for the 2000 molecules to fill the valleys on the rough surface. From the point of view of normal force, it could also provide some clues. The normal force on the top plane of lubricated contact was divided into two parts: normal force FS from the substrate and normal force FL from the lubricant. The sum of FS and FL is equal to the normal force F0 on the top plane of dry contact (no lubricant molecules). Figure 6b shows the normal force on the top plane supported by dry and lubricated contacts as a function of molecule number. In the figure, the support force FS/F0 from the asperities roughly decreased with the larger number of molecules, which is consistent with the trend in the contact area. At the same time, the normal force FL/F0 from the

Contact Area
The effects of roughness and lubricant on the contact area calculated by the project method will be studied in this section. Figure 4a indicates that the contact area goes up with the increasing loads for the surfaces with three different RMS. Under the load of 0.05 GPa, the contact area of RMS = 0.2 nm is close to 50%, which is much higher than 4% and 3.5% of the RMS = 0.5 and 0.8 nm. As mentioned in the previous section, the heights of the asperities on the surface of RMS = 0.2 nm are relatively low, so a number of them were in contact simultaneously. For the rougher surface, however, only few asperities supported the load applied from the top plane. With the increasing loads, the contacting asperity deformed and their heights decreased while the lower non-contacting asperities were joined to support the top plane. Therefore, the contact area of the three surfaces was increased in Figure 4a. Considering the different topography due to the different random seed when generating the self-affine surfaces [35] and adhesive strength [34], the results could deviate a little. The relationship between the contact area and the loads obtained in this work qualitatively agrees with the published literature [44].
For the lubricated contacts, the contact area still increased with the loads at the range of 0.05 GPa (red diamond) to 0.25 GPa (purple square) as shown in Figure 6a. However, the contact area reduced with the increasing number of molecules. The lubricants could separate the substrate from the top plane. The separating gap depends on the lubricant amount. It indicates that under 0.25 GPa the contact area of 0 and 1000 molecules is almost identical while the contact area of 4000 molecules was reduced by 75%. However, it was not enough for the 2000 molecules to fill the valleys on the rough surface. From the point of view of normal force, it could also provide some clues. The normal force on the top plane of lubricated contact was divided into two parts: normal force F S from the substrate and normal force F L from the lubricant. The sum of F S and F L is equal to the normal force F 0 on the top plane of dry contact (no lubricant molecules). Figure 6b shows the normal force on the top plane supported by dry and lubricated contacts as a function of molecule number. In the figure, the support force F S /F 0 from the asperities roughly decreased with the larger number of molecules, which is consistent with the trend in the contact area. At the same time, the normal force F L /F 0 from the lubricant increased with the higher number of molecules. The decrease in the normal force F S /F 0 under each load indicates that the number of the contacting asperities dropped. At 0.25 GPa, more lubricant endured the larger force. When 4000 hexadecane molecules were used, F L reached 87.17 nN under 0.25 GPa, which corresponds to 85.45% of F S in dry contact. The project method can be used to calculate the volume occupied by the hexadecane molecules-thereby estimating its density. For RMS = 0.8 nm with 4000 molecules under 0.15 and 0.25 GPa, the hexadecane has the density of 6.9 and 7.2 g/cm 3 , respectively. The estimations are very close to 7.7 g/cm 3 obtained from the previous experiments [45].
lubricant increased with the higher number of molecules. The decrease in the normal force FS/F0 under each load indicates that the number of the contacting asperities dropped. At 0.25 GPa, more lubricant endured the larger force. When 4000 hexadecane molecules were used, FL reached 87.17 nN under 0.25 GPa, which corresponds to 85.45% of FS in dry contact. The project method can be used to calculate the volume occupied by the hexadecane molecules-thereby estimating its density. For RMS = 0.8 nm with 4000 molecules under 0.15 and 0.25 GPa, the hexadecane has the density of 6.9 and 7.2 g/cm 3 , respectively. The estimations are very close to 7.7 g/cm 3 obtained from the previous experiments [45].

Surface Pressure
Surface pressure was calculated based on the projection method in section 2.3. There are other approaches for calculating surface pressure, but for the comparisons the projection method was used.
The virial stress of each atom was calculated at every 20 ps. After the cuboids were obtained, the surface pressure was estimated by averaging the virial stress along the z direction of all atoms within each cuboid. Before that, we show in Figure 7 the deformation of the FEM region due to the applied loads from the top plane on the substrate. Before compression, there is almost no deformation at the coarse region as shown in Figure 7a. With the compression of the top plane, the deformation is small in this region as shown in Figure 7b, as the FEM region is far from the surface. However, more

Surface Pressure
Surface pressure was calculated based on the projection method in Section 2.3. There are other approaches for calculating surface pressure, but for the comparisons the projection method was used. The virial stress of each atom was calculated at every 20 ps. After the cuboids were obtained, the surface pressure was estimated by averaging the virial stress along the z direction of all atoms within each cuboid. Before that, we show in Figure 7 the deformation of the FEM region due to the applied loads from the top plane on the substrate. Before compression, there is almost no deformation at the coarse region as shown in Figure 7a. With the compression of the top plane, the deformation is small in this region as shown in Figure 7b, as the FEM region is far from the surface. However, more deformation occurred in the FEM region, since the asperities were flattened in the corresponding area on the surface, indicating that the coarse region can properly accommodate the elastic deformation from the contact zone. deformation occurred in the FEM region, since the asperities were flattened in the corresponding area on the surface, indicating that the coarse region can properly accommodate the elastic deformation from the contact zone.  Figure 8 shows the pressure distribution for the unlubricated and lubricated contacts. To obtain the distribution histogram, the surface pressure in each cuboid was binned with an interval of 0.05 GPa. The shape of the probability distribution in this work was in agreement with the previous works [30,34]. It seems that the system size and the roughness of the top surface have little effect on the shape of the probability distribution. In Figure 8a, the sum of the negative pressure distribution is 88.1%, 80.0% and 66.79% for 0.05, 0.15 and 0.25 GPa, respectively. Meanwhile, in Figure 8b, the obtained values are 89.2%, 72.8% and 63.4%, respectively. More than half of the surface area mainly suffered the attractive forces. In addition, for both dry and lubricated surfaces, the probability sum of the positive pressure goes up as the loads increase as shown in Figure 8. Further compression caused the increase in the repulsive forces between the top plane, lubricants and substrate. As demonstrated by Cheng et al. [46], the lubricant molecules spread the pressure over a larger contact area. The shifting of probability with the increasing load for the lubricated contacts is much more than that for the dry contacts, indicating the significant role of lubricant molecules in bearing the loads.  Figure 8 shows the pressure distribution for the unlubricated and lubricated contacts. To obtain the distribution histogram, the surface pressure in each cuboid was binned with an interval of 0.05 GPa. The shape of the probability distribution in this work was in agreement with the previous works [30,34]. It seems that the system size and the roughness of the top surface have little effect on the shape of the probability distribution. In Figure 8a, the sum of the negative pressure distribution is 88.1%, 80.0% and 66.79% for 0.05, 0.15 and 0.25 GPa, respectively. Meanwhile, in Figure 8b, the obtained values are 89.2%, 72.8% and 63.4%, respectively. More than half of the surface area mainly suffered the attractive forces. In addition, for both dry and lubricated surfaces, the probability sum of the positive pressure goes up as the loads increase as shown in Figure 8. Further compression caused the increase in the repulsive forces between the top plane, lubricants and substrate. As demonstrated by Cheng et al. [46], the lubricant molecules spread the pressure over a larger contact area. The shifting of probability with the increasing load for the lubricated contacts is much more than that for the dry contacts, indicating the significant role of lubricant molecules in bearing the loads.

Conclusions
In the present study, a concurrent multiscale model that combines the MD and FEM methods was proposed and successfully applied to study the three-dimensional rough lubricated contact of aluminum single crystal. We performed the Hertz contact to compare the MU model with the full MD model. Our MU simulations show good accuracy from the contact area, kinetic energy and stress continuity aspects, and it can save at least five times the computational time compared with the full MD simulations for the same domain size. With relatively good computational cost and accuracy, the MU model here has been extended to study lubricated contacts. A systematic investigation considering the effects of the lubricant amount, surface roughness and loads during contact in mixed lubrication at the atomic scale was conducted. It was found that the sufficient lubricant molecules could fill the cavity and cover the large portion of contacting surfaces and effectively separate the tribo-pair. Meanwhile, the insufficient lubricant molecules can work well for relatively flat surfaces. In addition, the results show that the increasing number of lubricant molecules leads to smaller changes in surface RMS under the studied loads, varying from 0.05 to 0.25 GPa, indicating that the thin film formed by the lubricant molecules is beneficial to maintain the initial surface topography. For both the lubricated and unlubricated contacts, the contact area decreases with the increasing loads. As the tribo-pair gap increases with the increasing lubricant amount, the contact area decreases, and the larger portion of loads is supported by the lubricant molecules. The changes in pressure distribution for the different loads in the lubricated contacts confirm the contributions of the lubricant molecules to the load-bearing.
Author Contributions: Conceptualization, software, visualization, writing-original draft, J.Z.; methodology, J.Z. and L.S.; validation, writing-review and editing, L.S. and Z.W. All authors have read and agreed to the published version of the manuscript.

Conclusions
In the present study, a concurrent multiscale model that combines the MD and FEM methods was proposed and successfully applied to study the three-dimensional rough lubricated contact of aluminum single crystal. We performed the Hertz contact to compare the MU model with the full MD model. Our MU simulations show good accuracy from the contact area, kinetic energy and stress continuity aspects, and it can save at least five times the computational time compared with the full MD simulations for the same domain size. With relatively good computational cost and accuracy, the MU model here has been extended to study lubricated contacts. A systematic investigation considering the effects of the lubricant amount, surface roughness and loads during contact in mixed lubrication at the atomic scale was conducted. It was found that the sufficient lubricant molecules could fill the cavity and cover the large portion of contacting surfaces and effectively separate the tribo-pair. Meanwhile, the insufficient lubricant molecules can work well for relatively flat surfaces. In addition, the results show that the increasing number of lubricant molecules leads to smaller changes in surface RMS under the studied loads, varying from 0.05 to 0.25 GPa, indicating that the thin film formed by the lubricant molecules is beneficial to maintain the initial surface topography. For both the lubricated and unlubricated contacts, the contact area decreases with the increasing loads. As the tribo-pair gap increases with the increasing lubricant amount, the contact area decreases, and the larger portion of loads is supported by the lubricant molecules. The changes in pressure distribution for the different loads in the lubricated contacts confirm the contributions of the lubricant molecules to the load-bearing.

Conflicts of Interest:
The authors declare no conflict of interest.