Basic Equations

The elastic deformation v is calculated via the Boussinesq integration in the entire computational domain, as Eq. (1) [16]:

$$ v_{x,y} = frac{2}{{uppi E^{prime}}}iintlimits_{Omega } {frac{{p_{xi ,eta }^{{{text{tot}}}} }}{{sqrt {left( {x – xi } right)^{2} + left( {y – eta } right)^{2} } }}{text{d}}xi {text{d}}eta }, $$


where ptot is the total pressure, which comprises the surface force pressure psf 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):

$$ V_{alpha ,beta } = sumlimits_{xi } {sumlimits_{zeta } {C_{alpha – zeta ,beta – xi } P_{zeta ,xi }^{{{text{tot}}}} } } , $$


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

$$ begin{aligned} p_{{{text{sf}}}} & = p_{{{text{vdW}}}} + p_{{{text{dl}}}} + p_{{{text{hf}}}} \ & = – {{A_{{text{H}}} } mathord{left/ {vphantom {{A_{{text{H}}} } {6uppi h^{3} }}} right. kern-nulldelimiterspace} {6uppi h^{3} }} + 2varepsilon varepsilon_{0} kappa^{2} psi_{1} psi_{2} e^{ – kappa h} + p_{0} e^{{ – {{left( {h – h_{0} } right)} mathord{left/ {vphantom {{left( {h – h_{0} } right)} {lambda_{0} }}} right. kern-nulldelimiterspace} {lambda_{0} }}}} , \ end{aligned} $$


where AH 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; p0 and λ0 are the amplitude and decay length of the hydration force, respectively; h0 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.

$$ frac{partial }{partial x}left( {frac{{rho h^{3} }}{12eta }frac{partial p}{{partial x}}} right) + frac{partial }{partial y}left( {frac{{rho h^{3} }}{12eta }frac{partial p}{{partial y}}} right) = Ufrac{{partial left( {rho h} right)}}{partial x}, $$


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

$$ eta = eta_{0} e^{alpha p} , $$


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 water-based lubrication with pressure obeys the following relationship proposed by Dowson and Higginson [33]:

$$ rho = rho_{0} left( {1 + frac{{0.6 times 10^{ – 9} p}}{{1 + 1.7 times 10^{ – 9} p}}} right),, $$


where ρ0 is the ambient density of water.

In the direct-contact region, the direct-contact pressure is calculated using Eq. (7) [30]:

$$ P_{i,j} = {{left( {0 – H_{i,j} + C_{0,0} P_{i,j}^{{{text{old}}}} } right)} mathord{left/ {vphantom {{left( {0 – H_{i,j} + C_{0,0} P_{i,j}^{{{text{old}}}} } right)} {C_{0,0} }}} right. kern-nulldelimiterspace} {C_{0,0} }}, $$


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

$$ h = frac{1}{2R}left( {x^{2} + y^{2} } right) + delta left( {x,y} right) + h_{{text{r}}} + vleft( {x,y} right) + v^{{text{p}}} left( {x,y} right), $$


where the first term represents a general spherical profile, whereas the second to fifth terms represent the surface roughness profile δ, rigid-body displacement hr, elastic deformation v, and plastic residual displacement vp, 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].

$$ sigma_{alpha ,beta ,gamma }^{ij} = sumlimits_{xi } {sumlimits_{zeta } {left( {D_{alpha – zeta ,beta – xi ,gamma }^{{{text{N}}ij}} p_{zeta ,xi }^{{{text{tot}}}} + D_{alpha – zeta ,beta – xi ,gamma }^{{{text{S}}ij}} s_{zeta ,xi } } right)} } , $$


where γ corresponds to the coordinate along the depth direction; ij denotes the stress direction; DN and DS 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].

$$begin{aligned} sigma_{alpha ,beta ,gamma }^{*ij} & = frac{ – mu }{{4uppi (1 – nu )}}sumlimits_{vartheta } {sumlimits_{xi } {sumlimits_{zeta } {T_{alpha – zeta ,beta – xi ,gamma – vartheta }^{ijkl} varepsilon_{zeta ,xi ,vartheta }^{*kl} } } } \ & quad + frac{ – mu }{{4uppi (1 – nu )}}sumlimits_{vartheta } {sumlimits_{xi } {sumlimits_{zeta } {T_{alpha – zeta ,beta – xi ,gamma + vartheta }^{ijkl} {{varepsilon^{prime}}_{zeta ,xi ,vartheta }^{*kl}} } } } \ & quad + 2sumlimits_{xi } {sumlimits_{zeta } {D_{alpha – zeta ,beta – xi ,gamma }^{{_{Nij} }} tau_{zeta ,xi }^{zz} } } , \ end{aligned}$$


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.

$$ V_{alpha ,beta }^{{text{p}}} = 2sumlimits_{xi } {sumlimits_{zeta } {C_{alpha – zeta ,beta – xi } tau_{zeta ,xi }^{{{text{zz}}}} } } . $$


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 stress-pressure coefficients DN and DS 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.

Figure 1
figure 1

Flow chart of numerical procedure

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

$$ left{ begin{gathered} varepsilon_{{text{p}}} = frac{{sumlimits_{y} {sumlimits_{x} {left| {p_{x,y} – p_{x,y}^{{{text{old}}}} } right|} } }}{{sumlimits_{y} {sumlimits_{x} {p_{x,y} } } }}, hfill \ varepsilon_{{text{w}}} = frac{{left| {w – sumlimits_{y} {sumlimits_{x} {p_{x,y} } } } right|}}{w}, hfill \ varepsilon_{{{text{psf}}}} = frac{{sumlimits_{y} {sumlimits_{x} {left| {p_{x,y}^{{{text{sf}}}} – p_{x,y}^{{text{sf – cal}}} } right|} } }}{{sumlimits_{y} {sumlimits_{x} {p_{x,y}^{{text{sf – cal}}} } } }}, hfill \ varepsilon_{{{text{vp}}}} = frac{{sumlimits_{y} {sumlimits_{x} {left| {v_{x,y}^{{text{p}}} – v_{x,y}^{{text{p – old}}} } right|} } }}{{sumlimits_{y} {sumlimits_{x} {left| {v_{x,y}^{{text{p}}} } right|} } }}, hfill \ end{gathered} right. $$


where w denotes the applied normal load, psf the surface force in the current iteration, and psf-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, psf is only determined by the local distance between two surfaces; thus, εpsf can indicate the convergence of the rigid-body 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.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit


This article is autogenerated using RSS feeds and has not been created or edited by OA JF.

Click here for Source link (