Plastic–elastic Model for Water-based Lubrication Considering Surface Force

Water-based lubrication is an effective method to achieve superlubricity, which implies a friction coefficient in the order of 10−3 or lower. Recent numerical, analytical, and experimental studies confirm that the surface force effect is crucial for realizing water-based superlubricity. To enhance the contribution of the surface force, soft and plastic materials can be utilized as friction pair materials because of their effect in increasing the contact area. A new numerical model of water-based lubrication that considers the surface force between plastic and elastic materials is developed in this study to investigate the effect of plastic flow in water-based lubrication. Considering the complexity of residual stress accumulation in lubrication problems, a simplified plastic model is proposed, which merely calculates the result of the dry contact solution and avoids repeated calculations of the plastic flow. The results of the two models show good agreement. Plastic deformation reduces the local contact pressure and enhances the function of the surface force, thus resulting in a lower friction coefficient.


Introduction
Friction is ubiquitous in daily life and industrial production. Recent statistics show that approximately 23% of global energy consumption is associated with friction, and this ratio is gradually declining owing to progress in tribology [1]. In the last three decades, researchers have focused on reducing friction and obtaining superlubricity, which implies a friction coefficient in the order of 10 −3 or lower [2,3]. Previous studies confirm that utilizing an aqueous salt solution as a lubricant is an effective method to realize liquid superlubricity under high loads [4][5][6]. The aqueous superlubricity mechanism can be ascribed to the combined action of hydration and hydrodynamic effects [4,7]. The hydration effect was initially proposed as strong repulsion between charged surfaces across high-concentration electrolytes during the measurement of surface force at the microscale [8,9]. Hydration repulsion originates from the steric force among hydrated ions, which are composed of central ion and surrounding water molecules. This effect was further developed to achieve excellent lubricating performance at both the micro and macro scales [4][5][6][10][11][12]. Hydrodynamic effects have been investigated since the nineteenth century, when the well-known Navier-Stokes and Reynolds equation were proposed. These equations consistently provide the foundation for hydrodynamic lubrication studies.
When the relative sliding velocity of a friction pair decreases, the lubrication regime changes from full-film hydrodynamic lubrication to mixed lubrication and then to boundary lubrication, as the hydrodynamic effect weakens gradually. Conventional hydrodynamic lubrication barely contribute to friction reduction at lowvelocity conditions but can achieve low friction under this condition via hydration lubrication. To improve the hydration effect, a friction system where the nominal contact areas are in close proximity with each other is preferred because the hydration force is particularly effective at extremely small surface spacings [8,10,12]. Therefore, materials such as ultrahigh-molecular-weight polyethylene (UHMWPE) and poly-ether-ether-ketone  35:117 (PEEK) are suitable options as friction pairs in waterbased lubrication [6,[13][14][15]. Compared with ceramics or sapphire, these soft plastic materials are more susceptible to deformation and contact under compression, which is conducive to improving the contribution of the surface force. They are widely used in artificial cartilages and water-lubricated marine stern tube bearings. These frictional components typically operate in an aqueous environments and propagate at low speeds. It is difficult to form lubricating film under such circumstances. The materials mentioned above, which exhibit excellent tribological performance, can be utilized to avoid severe abrasion or loud noise caused by friction. Although many researchers have experimentally demonstrated the water-based lubricating performance of these materials, they are unable to directly observe the mechanism at the microscale. Whether the hydration effect is applicable to these friction systems and the extent to which the material properties should be modified to realize superlubricity are difficult to determine using conventional methods. Hence, numerical analysis may be a feasible alternative method.
The numerical model of plastic-elastic materials is significantly more complicated than that of purely elastic materials. The numerical plastic model at the early stage merely limits the local contact pressure below the yield stress of the materials [16,17]. This method does not involve strict mechanical derivation. Moreover, the surface residual displacement or subsurface stress fields are not available. Jacq et al. proposed a semi-analytical plastic-elastic contact model for theoretically smooth and rough surfaces with a single asperity [18]. Subsequently, it was further developed using a more complicated rough surface and more efficient calculations [19,20]. However, the numerical procedure involved was laborious. Wang et al. proposed a new method for calculating the residual stress and residual displacement [21]. First, they obtained the strain field in half-space induced by the plastic strain of a small cubic element using the superposition method. The half-space solution was derived from adding an infinite-space solution from two mirror elements and a half-space solution with normal stress only [22]. This superposition method successfully eliminated the shear stress on the free surface of the half-space. Subsequently, the residual stress was derived from a strain induced using Hooke's law. Finally, the residual displacement was calculated from the normal stress mentioned above using the deformation-pressure influence coefficient [23]. The plastic-elastic contact model based on this calculation procedure was validated by comparing the results with those obtained using the finite element method. To further improve efficiency, Wang et al. [24] developed a plastic-elastic contact model via the three-dimensional fast Fourier transform method and a parallel computational strategy by adopting the influence coefficient between plastic strain and residual stress proposed by Liu et al. [25,26]. Compared with the previous model, this model does not require the calculation of induced strain. It is currently the most efficient plastic-elastic model.
The method to calculate the residual stress and residual displacement mentioned above is a general method that is applicable to various circumstances. Combining it with a contact model based on the conjugate gradient method [27] resulted in a plastic-elastic contact model was realized [21]. Meanwhile, combining it with a full mixed lubrication model [28] resulted in a plasto-elastohydrodynamic lubrication model [29]. Based on a waterbased lubrication model that considers the surface force proposed previously by the authors [7], a plastic-elastic model of water-based lubrication was established in the current study. Owing to the sensitivity of the plastic flow and surface force calculation, a mixed lubrication model with high accuracy was considered [30]. By determining the lubrication behavior between epoxy surfaces across a 0.9% sodium chloride solution, which is the same condition as in a previous experimental study [6], the effect of plastic flow on water-based lubrication was investigated using the new model.

Basic Equations
The elastic deformation v is calculated via the Boussinesq integration in the entire computational domain, as Eq. (1) [16]: where p tot is the total pressure, which comprises the surface force pressure p sf and hydrodynamic pressure or direct-contact pressure p; E' is the effective elastic modulus of the friction pair [7,30].
After discretizing the integration, the summation form of Eq. (1) is as Eq. (2): where C is the deformation-pressure influence coefficient [23]; α, β, ζ, and ξ are the discrete coordinates. The upper-case letters indicate the dimensionless form of the corresponding parameters [30].
The surface force, which constitutes the total pressure, is particularly effective at extremely low surface spacings (< 2 nm), and includes the van der Waals, electric double layer, and hydration forces [7,8,10].  35:117 where A H is the Hamaker constant; h is the surface spacing or film thickness; ε 0 and ε are the vacuum and relative dielectric constants, respectively; κ is the reciprocal of the Debye length; ψ 1 and ψ 2 are the surface potentials of each surface; p 0 and λ 0 are the amplitude and decay length of the hydration force, respectively; h 0 is the surface spacing at which the hydration force reaches its maximum.
For the other constituent of the total pressure, the entire computational domain is partitioned into a hydrodynamic lubrication region and a direct-contact region. The hydrodynamic pressure in the lubrication region is governed by the steady-state Reynolds equation.
where p is the local hydrodynamic pressure; U is the entrainment speed; η and ρ are the viscosity and density of the lubricant, respectively.
In an actual ball-on-disk experiment, the upper surface is typically fixed, and the lower disk rotates at a specified speed, indicating a pure sliding regime. The x-coordinate coincides with the sliding direction of the lower surface, which is also the movement direction of the lubricant.
Based on the Barus relationship, the viscosity of a water-based solution varies with pressure [31].
where η 0 is the ambient viscosity of water.
In this study, the viscosity-pressure coefficient α was set as 0.8 GPa -1 to fit the experimental viscosity of water at room temperature [32]. The change in density of waterbased lubrication with pressure obeys the following relationship proposed by Dowson and Higginson [33]: where ρ 0 is the ambient density of water.
In the direct-contact region, the direct-contact pressure is calculated using Eq. (7) [30]: where the superscript "old" implies the result of the previous iteration.
In Eqs. (3), (4), and (7), the film thickness comprises five components, as Eq. (8): (3) where the first term represents a general spherical profile, whereas the second to fifth terms represent the surface roughness profile δ, rigid-body displacement h r , elastic deformation v, and plastic residual displacement v p , respectively.
A series of calculations must be conducted to obtain the residual displacement. First, the subsurface stress field is calculated based on the total pressure distribution [34].
where γ corresponds to the coordinate along the depth direction; ij denotes the stress direction; D N and D S are the influence coefficients for calculating the subsurface stress from the normal pressure and surface shear stress, respectively; s is the shear stress on the surface, namely the friction stress, which is calculated by Newton's law of viscosity in the lubrication region. The surface force does not affect the calculation of the friction stress [7]. In the direct-contact region, the boundary lubrication friction coefficient was set as 0.2 for the epoxy friction pair [6].
Next, the plastic flow is investigated based on the J2-flow theory with a known subsurface stress field. The plastic strain ε * can be calculated using the radial return method [21,35]. Subsequently, the residual stress σ * is computed using the superposition method mentioned above [24][25][26].
where ij and kl denotes the direction of stress and strain, respectively; T is the influence coefficient between the plastic strain and residual stress; ε * and ε' * are the plastic strain matrices of the local and mirror domains, respectively; τ zz is the normal stress on the surface caused by one domain.
Finally, the plastic residual displacement is obtained by multiplying the normal stress by 2.

Numerical Procedure
The calculation steps are illustrated in Figure 1. After parameter initialization, three influence coefficient matrices were calculated, including the elastic deformation-pressure coefficient C in Eq. (2), subsurface stresspressure coefficients D N and D S in Eq. (9), and residual stress-plastic strain coefficient T in Eq. (10). They must be calculated only once and stored for later use. In a complete plastic-elastic problem, the normal load should be applied stepwise to simulate the accumulation of plastic strain and residual stress [21,24]. To avoid repeated loading processes and conserve computing time, the plastic-elastic contact solution was imported as the initial state in the lubrication problem, which is reasonable because the actual friction test only begins via static loading without shear movement. Notably, the contact solution must be derived from the same applied normal load and surface profile as those in the lubrication problem. Thus, one plastic-elastic contact solution can be used repeatedly as the initial value in the plastic-elastic lubrication problem at different sliding speeds.
The entire numerical procedure comprises two main routines: iterations of water-based lubrication [7] and plastic flow [24]. The equations for lubrication are solved using the Gauss-Seidel method via an under-relaxation strategy. A detailed description of this numerical procedure is available in Refs. [7,28,30]. The rigid-body displacement was adjusted during the iterative procedure to achieve balance between the applied normal load and the summation of the total calculated pressure. The total pressure was utilized to calculate the residual stress and plastic residual displacement via the steps mentioned in the previous section only after pressure convergence was achieved in the lubrication iteration. Once the iterative error of the residual displacement ε vp became sufficiently small, the flag controlling the iteration of the plastic flow was turned off, and the entire process was completed after the final lubrication iteration. This procedure prevents excessive iterations of the plastic flow and ensures that each iteration is based on a stable condition.
The four errors for judging convergence in Figure 1 were calculated using Eq. (12): where w denotes the applied normal load, p sf the surface force in the current iteration, and p sf-cal the precise surface force based on Eq. (3).
The four criteria in Figure 1 can be adjusted in different cases for better precision; however, a lower criterion significantly increases the computing time. Notably, p sf is only determined by the local distance between two surfaces; thus, ε psf can indicate the convergence of the rigidbody displacement iteration.
Because the plastic-elastic model introduces a third dimension (depth direction) compared with the elastic model, the computational effort increases geometrically. Efficiency is an important aspect in the plastic-elastic model. In this study, a simplified calculation procedure was established by canceling the plastic flow iteration shown in Figure 1, which implies that the plastic strain, residual stress, and residual displacement are entirely based on the contact solution and do not evolve after shear movement. This is acceptable because most of the plastic deformation occurred during the static loading stage of the friction test. The explicit error of this simplified model and further discussions regarding the calculation results will be presented in the next section.

Parameters of Surface Force
In a previous water-based lubrication model [7], we used the fitting parameters from a surface force apparatus [10] without considering material selection. Here, an energy conservation method was applied to determine the parameters of specific materials. Assuming that two surfaces approach each other across an aqueous solution, the work induced by overcoming the hydration repulsion force should be equivalent to the energy required to dehydrate all hydrated ions adsorbed on the surfaces [12].
where σ 1 and σ 2 are the charge densities of the two surfaces, ΔG h is the hydration energy of a specific ion, and e is the electron charge. In this study, the ΔG h of Na + was 410 kJ/mol [36], λ 0 was set as 0.2 nm [10], and h 0 was set as 0.27 nm (diameter of water molecule). The surface charge density σ is associated with the surface potential ψ based on the Graham equation [8].
where k B is the Boltzmann constant, and T is the thermodynamic temperature. The surface potential is a theoretical parameter that is difficult to measure directly. It is currently estimated using the zeta potential.  35:117 Considering that the friction pair for the calculation comprises an epoxy hemispherical surface and a sapphire disk, the zeta potentials were −25 mV for epoxy [6] and −20 mV for sapphire [37], both in 1 mM alkali salt solution. Thus, the corresponding surface charge densities were 1.87 mC/m 2 for the epoxy and 1.48 mC/m 2 for the sapphire. Meanwhile, the surface force amplitude p 0 of this friction pair was 71.25 MPa, based on Eq. (13). The zeta potentials in 0.9% NaCl (0.154 M molar concentration) were obtained using Eq. (14) with a fixed surface charge density, which were −2.09 mV for epoxy and −1.65 mV for sapphire. The Hamaker constants were 0.3 × 10 -20 J for the epoxy [38] and 5.32 × 10 −20 J [39] for the sapphire, and the composite constant was their geometric average [8]. Based on the parameters above, the variation in the surface force with the surface spacing is depicted in Figure 2. The surface force decreased significantly as the surface spacing increased and was particularly effective when the surface spacing was less than 1 nm. The decreasing trend of the repulsion surface force at small surface spacings was due to the rapid increase in the van der Waals attraction.

Results of Two Plastic-elastic Lubrication Models
The operating conditions of an actual test [6] were used to validate the feasibility of our new model. Friction occurred between the epoxy spherical pin and sapphire disk across a 0.9% NaCl solution under an applied normal load of 5 N. The elastic moduli of the epoxy and sapphire were 3 and 410 GPa, respectively, and their Poisson's ratios were 0.3 and 0.27, respectively. The curvature radius of the worn hemispherical surface was 46 mm, as deduced from the contact area radius b = 0.375 mm based on the Hertz contact theory. The corresponding maximum Hertz contact pressure was 17 MPa. However, the local pressure can be much higher than this value owing to the high surface roughness of the epoxy surface. The ambient viscosity of the aqueous solution was set as 0.8 MPa•s at room temperature. The yield stress of the epoxy spherical pin was set as 20 MPa, and the material was assumed to be perfectly plastic without work hardening. The roughness of the epoxy surface was represented by a computer-generated discrete Gaussian surface profile comprising 256 × 256 grids and a root mean square roughness Rq of 0.3 μm [30,40].
First, a plastic-elastic contact solution was achieved to obtain the initial state. The normal load was applied in five intervals from 1 to 5 N. Subsequently, the pressure distribution, residual stress, plastic strain, and residual displacement were applied as input data in the lubrication calculation. Five sliding speeds were used: 1.2, 4.7, 23.5, 47.1 and 94.2 mm/s. These speeds are not high and typically result in boundary lubrication or mixed lubrication, which are focused on in this study. The calculated coefficient of friction (COF) curves are shown in Figure 3. The method of calculating the COF is available in Ref. [7]. For comparison, Figure 3 The friction coefficients calculated using the simplified and complete models showed good agreement. Their difference was less significant when the surface force was considered. This is because the friction-force-induced stress was much lower in this case, and the normal pressure distribution was similar to that of the contact solution. This occurred because the surface force separated the surfaces in direct contact when the surface force was not considered. As depicted in Figure 4, the final plastic residual displacement of the complete plastic-elastic model considering the surface force coincides with the residual displacement of the contact solution, which implies that the accumulation of residual stress after the static loading stage is insignificant. Therefore, the friction coefficient does not change significantly.
Meanwhile, significant differences were observed when the surface forces were not considered, particularly in the spike region, where a high local contact pressure occurred with asperities. In cases involving low sliding speeds, the two surfaces were more likely to be in direct contact around the asperities. The friction stress in the direct-contact region was high, which increased the subsurface stress and plastic residual displacement. Therefore, the difference in the friction coefficient between the simplified and complete models without considering the surface force was significant. For the model that considers the surface force, if the surface force is insufficient to separate the asperities, then the deviation of the simplified model will be greater, as shown in Figure 3(b). Because the time consumption of the complete model is typically more than 10 times that of the simplified model, the latter is preferred as it is more efficient in many cases.
The following discussion is based on the results of the complete model.

Effect of Plastic Deformation on Lubrication
As shown in Figure 3, the friction coefficient decreased significantly after the surface force was considered. This is because the surface force undertakes a significant portion of the normal load that the hydrodynamic effect cannot support [7]. However, the effect of plastic flow on lubrication was different in these two cases. Without considering the surface force, the friction coefficient of the plastic-elastic model was slightly higher than that of the purely elastic model. By contrast, the friction coefficient was lower when surface force was considered. To investigate this phenomenon, the surface profiles before and after plastic deformation, as shown in Figure 5, were analyzed. The profile was composed of the general shape of the sphere, surface roughness, and plastic residual displacement after plastic deformation during the static loading stage. Because the maximum Hertz contact pressure was 17 MPa, which is less than the yield stress of 20 MPa, no large-scale plastic deformation occurred in the overall  computational domain. Plastic deformation primarily occurred around sharp asperities where the local contact pressure and subsurface stress were high. The heights of the asperities reduced after plastic deformation. However, as the deformation leveled the surface in the contact area, the possibility of direct contact increased, particularly at low sliding speeds. Therefore, the friction coefficient was slightly higher for the plastic-elastic model that did not consider the surface force. As the sliding speed increased, the difference between the two models decreased gradually (see Figure 3). It can be inferred that when the hydrodynamic effect is stronger at high speeds, the friction coefficient of the plastic-elastic model may be lower than that of the elastic model owing to the decrease in the overall surface roughness, even though the surface force is not considered. However, a different situation is encountered when the surface force is considered. After the plastic deformation, as the overall surface roughness decreased, the local contact pressure decreased as well. Therefore, in the direct-contact region, where the hydrodynamic force was inadequate to accommodate the local pressure, the surface force was sufficient to support the load and separate the two surfaces. The area of the direct-contact region was further reduced, and thus, the friction coefficient decreased. Therefore, plastic deformation improves the function of the surface force and reduces the friction coefficient in water-based lubrication. The effect of plastic deformation can be analyzed comprehensively based on Figure 6, which shows the film thickness and pressure distribution, subsurface von Mises stress field, and residual stress field along the central line y = 0 obtained using the plastic-elastic model.
In Figure 6, the corresponding figures of the stress field are presented in the same color shade, as depicted on the right. Without considering the surface force, the hydrodynamic effect was insufficient to support the normal load, and most of the computational domain belonged to the direct-contact region. Owing to the high shear force in the direct-contact region, the subsurface von Mises stress near the asperities was much higher than the stress for the case without considering the surface force, which is consistent with our previous analysis. The von Mises stresses in the entire domain were limited to less than the yield stress of 20 MPa after plastic flow. The resulting residual stress fields of the two models were similar because they were primarily generated from normal pressure but not shear force, which is consistent with the residual displacement shown in Figure 4. Residual stress primarily occurred at shallow depths. The entire nominal contact area can be partitioned into several small contact regions around each asperity, and the width of the residual-stress region was comparable to their size.
Analogically, the entire friction pair comprising the epoxy spherical pin and sapphire disk can be segmented into small friction pairs between the asperities and disk. Because the surface roughness was relatively high, the hydrodynamic effect at the rough valleys sandwiched between asperities was weak. The surface force separated the region that was originally in contact when it was not considered, resulting in a current film thickness of 2 nm. A hydrodynamic pressure can exist in thin fluid films.
The effect of plastic deformation on water-based lubrication was further investigated. Figure 7 shows the calculation results for different yield stresses. The elastic solution can be regarded as a case involving an infinite yield stress. As the yield stress decreased, the plastic deformation increased, and an increasing number of asperities flattened. Without considering the surface force, when the yield stress was extremely low, the surface roughness decreased owing to plastic deformation, and the hydrodynamic effect was sufficient at a speed of 94.2 mm/s, which reduced the friction coefficient to a value lower than that of the elastic solution, as shown in Figure 7(a).
Based on the surface force in Figure 7(b), the calculated friction coefficient decreased as the yield stress decreased and the plastic deformation increased. The corresponding ratio of the normal load undertaken by the surface force increased, except that at a high speed with a low yield stress, the surface force ratio decreased because the hydrodynamic pressure supported more load. Generally, the friction coefficient should be high at low speeds owing to the weak hydrodynamic effects. In this case, however, the asperities flattened after plastic deformation, and the local contact pressure was extremely low, which implies that the surface force can support most of the normal load. Thus, a lower speed results in a smaller surface spacing and a stronger surface force, as shown in Figure 2. Therefore, the friction coefficient can be reduced by decreasing the sliding speed. The results above suggest that plastic deformation is conducive to friction reduction in water-based lubrication where an aqueous salt solution is used as a lubricant and surface force can be realized.

Effect of Surface Roughness on Lubrication
As discussed, a lower yield stress of the material results in a lower local contact pressure and a lower friction coefficient when a surface force is present. However, plastic materials with low yield stresses exhibit unsatisfactory mechanical properties, thus rendering them difficult to operate under high applied loads. Apart from flattening the asperities via more violent plastic deformations, surfaces with lower roughness levels can be achieved before the friction test via polishing or other methods. To investigate the effect of surface roughness on water-based lubrication between plastic and elastic materials, friction coefficient curves and surface force ratios based on different roughness values Rq were constructed, as shown in Figure 8. Profiles showing different surface roughness levels were prepared by multiplying the original roughness matrix with different coefficients. The friction coefficient decreased over the entire speed range as the surface roughness decreased (see Figure 8(a)), which is a different trend from that shown in Figure 7(a). Compared with only several prominent asperities being flattened by decreasing the yield stress, the number of asperities in the entire domain reduced after the surface roughness decreased, which enhanced the hydrodynamic effect and reduced the friction coefficient. As shown in Figure 8(b), the effect of reducing the surface roughness is similar to the effect of decreasing the yield stress, as shown in Figure 7(b). As the surface roughness decreased, the friction coefficient continuously decreased owing to the simultaneous enhancement in the hydrodynamic and hydration effects. When the surface roughness was 0.2 μm, the surface force ratio decreased from approximately 70% to less than 30% as the speed increased, whereas the friction coefficient remained low, which implies that the hydrodynamic effect was sufficiently strong under the conditions above. Based on the surface force ratio at the lowest speed, the hydration effect improved between surfaces with low roughness levels. The hydration effect was more important under severe conditions, such as low sliding speeds. Using a smoother surface is the preferred method to achieve water-based superlubricity.

Conclusions
A plastic-elastic model of water-based lubrication as well as a simplified model were established in this study. The new models simultaneously consider the surface force and plastic flow in mixed lubrication modeling. Both models can be used to investigate the effect of plastic deformation on water-based lubrication with different precisions. The main conclusions are summarized as follows: (1) When surface force was absent, slight plastic deformation widened some rough asperities and increased the possible area of direct contact where the hydrodynamic effect could not be realized; as such, the friction coefficient increased. (2) When surface force was present, the decrease in the surface roughness and local pressure caused by plastic deformation improved the ratio of the normal load borne by the surface force. In particular, the surface force reduced the friction around the asperities and then decreased the subsurface stress in the corresponding region. (3) As the yield stress decreased, the plastic deformation progressed further and the surface became smoother, thus indicating that plastic deformation is conducive to friction reduction when an aqueous salt solution is used as a lubricant. Utilizing a plastic-elastic material is an effective approach to achieve liquid superlubricity. When designing a practical frictional component, a plastic-elastic material is preferable over an elastic material with similar mechanical properties.