*Article* **Current and Future Trends of Resource Misallocation in the Construction Industry: A Bibliometric Review with Grounded Theory**

**Jingxiao Zhang 1,\*, Fangyu Dong 1, Pablo Ballesteros-Pérez 2, Hui Li 3,\* and Martin Skitmore 4**

**Abstract:** Resource misallocation (RM) refers to the existence of marginal output inequalities between different industries or companies in an economy. Prior studies of RM have mostly focused on effect analysis, construction industry structure upgrades, and organization management. However, these studies have been fragmented and unrelated. This paper analyzes the status quo, consequences, and emerging trends of RM research at the macroscopic level based on current problems and with the aim of exploring potential solutions. Drawing on grounded theory, a qualitative analysis using text-mining is used to analyze the characteristics of 124 RM-related papers. The results more comprehensively and systematically reveal that current RM research encompasses four major dimensions of sources and concepts, misallocation degree measurement and characterization, focused issues (field), and RM research deficiencies. Methods for measuring RM have also been developed from the simple proportional method to current mainstream methods (e.g., growth rate decomposition and variant substitution). We conclude that, in order for this discipline to thrive and effectively reduce RM, future research into RM should focus on core categories, especially the reform of market-oriented factors, transformation of government functions, construction industrial structure adjustment, and methods of income distribution. This systematic review provides a discipline oversight and uncovers necessary and potential research directions.

**Keywords:** resource misallocation; grounded theory; market distortion; total factor productivity


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 tribo-surfaces: 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:

**Figure 2.** Schematic of one hexadecane molecule's (**a**) initial structure and (**b**) after the optimization of conformation. **Figure 2.** Schematic of one hexadecane molecule's (**a**) initial structure and (**b**) after the optimization of conformation. The pink and yellow atoms indicate the united atoms of CH2 and CH3 .

The pink and yellow atoms indicate the united atoms of CH2 and CH3.

Etorsion = c0 + c

1

$$\mathbf{E}\_{\rm bond} = \mathbf{k}\_{\rm b} \left( \mathbf{r} - \mathbf{r}\_{0} \right)^{2},\tag{4}$$

$$\mathbf{E}\_{\text{bend}} = \mathbf{k}\_{\theta} \left(\theta - \theta\_0\right)^2 \,\text{\,\,\,\,\,}\tag{5}$$

[1 + cos(3ϕ)], (6)

$$\mathbf{E}\_{\text{torsion}} = \mathbf{c}\_0 + \mathbf{c}\_1[1 + \cos\phi] + \mathbf{c}\_2[1 - \cos(2\phi)] + \mathbf{c}\_3[1 + \cos(3\phi)],\tag{6}$$

[1 - cos(2ϕ)] + c3

[1 + cosϕ] + c2

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).

**Table 1.** Parameters for the Lennard–Jones (LJ) potentials for the non-bonded interactions.

a taken from Ref. [34]. b taken from Ref. [39]. c calculated from the Lorentz–Berthelot rules.

#### *2.3. 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.

#### *2.4. 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.

#### **3. Results and Discussion 3. Results and Discussion**

#### *3.1. 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. *3.1. 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

*Metals* **2020**, *10*, x FOR PEER REVIEW 6 of 14

good computational cost and accuracy, the MU model was here extended to study lubricated contacts.

From the contact area, kinetic energy and stress continuity aspects, our MU simulations show good

of C–C bonds—thereby affecting the system heights in Figure 3.

**Figure 3**. Variation of the system heights with the time for the root mean square (RMS) = (**a**) 0.2 nm, (**b**) 0.5 nm, (**c**) 0.8 nm under 0.25 MPa. The red lines indicate dry contact while the green and blue lines indicate lubricated contact with 2000 and 4000 hexadecane molecules, respectively. **Figure 3.** Variation of the system heights with the time for the root mean square (RMS) = (**a**) 0.2 nm, (**b**) 0.5 nm, (**c**) 0.8 nm under 0.25 MPa. The red lines indicate dry contact while the green and blue lines indicate lubricated contact with 2000 and 4000 hexadecane molecules, respectively.

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 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 (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.

**Figure 4.** (**a**) RMS change ∆RMS and contact area A/A0 of dry contact as a function of pressure. The circle, square and triangle indicate RMS = 0.2, 0.5 and 0.8 nm, respectively. Changes in the topography before (**b**) and after (**c**) loading for the surfaces with RMS = 0.2, 0.5 and 0.8 nm, respectively.

**Figure 4.** (**a**) RMS change ∆RMS and contact area A/A0 of dry contact as a function of pressure. The circle, square

loading for the surfaces with RMS=0.2, 0.5 and 0.8 nm, respectively.

**Figure 5.** Variations of the RMS change ∆RMS with the molecule number under 0.25 GPa. The RMS roughness of the surface is 0.8 nm. Insets show the surface topography and the hexadecane molecule distribution after loading. **Figure 5.** Variations of the RMS change ∆RMS with the molecule number under 0.25 GPa. The RMS roughness of the surface is 0.8 nm. Insets show the surface topography and the hexadecane molecule distribution after loading.

#### *3.2. Contact Area 3.2. 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 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].

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 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 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

experiments [45].

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/cm3 , respectively. The estimations are very close to 7.7 g/cm3 obtained from the previous experiments [45]. 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/cm3 , respectively. The estimations are very close to 7.7 g/cm3 obtained from the previous

lubricant endured the larger force. When 4000 hexadecane molecules were used, FL reached 87.17

*Metals* **2020**, *10*, x FOR PEER REVIEW 9 of 14

**Figure 6.** Contact area (**a**) and the normalized normal force (**b**) as a function of the lubricant molecule number for surface with a roughness RMS = 0.8 nm. The molecule number of zero indicates dry contact. The load is 0.25 GPa in (b). **Figure 6.** Contact area (**a**) and the normalized normal force (**b**) as a function of the lubricant molecule number for surface with a roughness RMS = 0.8 nm. The molecule number of zero indicates dry contact. The load is 0.25 GPa in (**b**).

#### *3.3. Surface Pressure 3.3. 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 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 surfacepressure 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.

loads.

from the contact zone.

deformation occurred in the FEM region, since the asperities were flattened in the corresponding area

**Figure 7.** Deformation of the FEM region (**a**) before compression and (**b**) after compression for dry contact. **Figure 7.** Deformation of the FEM region (**a**) before compression and (**b**) after compression for dry contact.

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 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.

*Metals* **2020**, *10*, x FOR PEER REVIEW 11 of 14

**Figure 8.** Pressure distribution in the dry contact (**a**) and 4000 molecules lubricated contact (**b**) under 0.05 GPa (red line), 0.15 GPa (green line) and 0.25 GPa (blue line). **Figure 8.** Pressure distribution in the dry contact (**a**) and 4000 molecules lubricated contact (**b**) under 0.05 GPa (red line), 0.15 GPa (green line) and 0.25 GPa (blue line).

#### **4. Conclusions 4. Conclusions**

lubricant molecules to the load-bearing.

published version of the manuscript.

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 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 **Supplementary Materials:** The following are available online at http://www.mdpi.com/2075-4701/10/7/965/s1, Figure S1: Geometry of the multiscale (a) and full MD (b) simulations, Figure S2: Comparisons between MU and full MD simulations: (a) Normalized contact area as a function of normalized indent displacement. (b) Kinetic energy near the substrate surface as a function normalized indent displacement. (c) Changes of stress component σzz in a line along the z direction when normalized indent displacement is 1.31a. Indent displacement and Z are normalized by lattice constant of a.

**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.

**Funding:** This research was funded by the Fundamental Research Funds for the Central Universities (grant number 2020JBM059), the National Natural Science Foundation of China (grant number 51905296) and Australian Research Council (ARC) Discovery Early Career Researcher Award (DECRA) project (DE180100124).

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

#### **References**

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Metals* Editorial Office E-mail<EMAIL_ADDRESS>www.mdpi.com/journal/metals

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34

www.mdpi.com ISBN 978-3-0365-7342-7