Plastic–elastic Model for Water-based Lubrication Considering Surface Force – Chinese Journal of Mechanical Engineering

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

where A_H is the Hamaker constant; h is the surface spacing or film thickness; ε₀ and ε are the vacuum and relative dielectric constants, respectively; κ is the reciprocal of the Debye length; ψ₁ and ψ₂ are the surface potentials of each surface; p₀ and λ₀ are the amplitude and decay length of the hydration force, respectively; h₀ 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 η₀ 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]:

where ρ₀ 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):

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

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

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