On the Numerical Modeling of Friction Hysteresis of Conformal Rough Contacts

Abstract

In this work, a numerical model simulating friction hysteresis for lubricated rough and textured surfaces in contact is presented. Friction hysteresis occurs in sliding contacts that are subjected to a non-constant (e.g., sinusoidal) motion. It refers to the phenomenon where the observed friction force during acceleration differs from that during deceleration. Besides the dynamics of the sliding system, a classic mixed friction model is adopted, in which the transient Reynolds equation for the description of the thin lubricant film is combined with a statistical Greenwood–Williamson model for the description of rough surface asperity contacts. The model enables the prediction of the friction hysteresis for predefined contact descriptions (i.e., surface profile and roughness, lubricant, etc.) and allows the study of the physics and parametric influences of dynamically sliding contacts. In this paper, it is shown that (i) friction hysteresis is captured by classic transient models for mixed lubrication; (ii) system parameters, such as roughness, applied load, viscosity and velocity, including the offset, amplitude and motion reversal, influence the shape and area of friction hysteresis; and (iii) the selection of the aforementioned parameters may minimize friction hysteresis.

1. Introduction

When a sliding contact is subjected to a cyclic tangential velocity, the resulting friction force not only depends on the instantaneous value of the sliding velocity but also on its history. Consequently, hysteresis appears in the measured friction force, which depends on the sliding velocity. Such friction hysteresis was first measured by Sampson et al. [1] in 1943. They measured a higher friction force in the acceleration phase in comparison with the deceleration phase during a stick–slip cycle in a lubricated sliding contact. In the literature, different terminology, such as multi-valued friction or friction lag, is used to describe friction hysteresis, depending on the scientific domain wherein this hysteresis behavior is encountered. However, in this work, the authors resort to friction hysteresis to denote the instantaneous friction measured in contact subjected to cyclic motion.

It was previously observed that friction hysteresis phenomena are important in controller design within the control community. In order to conduct correct motion simulations of complex mechanisms, the accurate modeling and simulation of the resulting friction forces is necessary. The modeling of the instantaneous friction force is needed not only to accurately estimate the instantaneous forces but to conduct simulations to allow the accurate estimation of the dynamic behavior of systems [2,3,4,5,6]. Because of the widespread use of fluid film bearings subjected to non-constant (i.e., oscillatory) motion in industrial applications, including pendulum bearings, servo-actuators, cross-head bearing, machine guide ways, etc. [7], this work deals specifically with friction hysteresis.

In 1990, Hess and Soom [8] conducted experiments on a lubricated line contact subjected to a cyclic triangular shaped velocity in observed friction hysteresis. They proposed a simple model starting from a phenomenological Stribeck description, incorporating a time lag to describe the observed friction hysteresis. A time lag between the sliding speed and resulting friction was proposed to describe the observed friction hysteresis. Two different types of lag were considered: (i) a constant time lag independent of the sliding speed and (ii) a time lag inversely proportional to the sliding speed. The constant time delay appeared to match the results more closely. The temporal delay was found to be increasing with load and viscosity but appeared to be independent of the driven oscillation frequency. A second study by Soom in 1995 showed that fluctuations in minimum film thickness (the so-called squeeze effect) are the principal causes for this deviation from the steady friction behavior [9]. In the work of Etsuo et al. [10], only a very small remark was made on the influence of surface roughness on the friction lag, but a systematic study was not undertaken by the authors. In following years, much work was focused on phenomenological modeling for control purposes of the observed friction phenomena under non-steady load and velocity [2,3,11]. These phenomenological models are capable of capturing friction hysteresis behavior but need experiments to determine their model parameters.

A few years later in 1997, Zhai et al.  [12] conducted numerical simulations, using a mixed friction EHL-line model to predict the friction hysteresis. The proposed model was compared to the experimental results of Hess and Zoom [8], and they matched the experimental results quite well. The authors proposed two factors that determine the magnitude of the hysteresis loop: (i) the frequency of velocity variation and (ii) the proportion of the load carried by asperities. They argued that an increase in sliding frequency leads to a higher rate of change in the film thickness, resulting in a higher squeeze force. When a smaller part of the load is supported by the asperities, the squeeze force can contribute more to the total force balance.

Harnoy et al. [13] developed a physics-based model for dynamic friction based upon a closed-form solution of the transient Reynolds model for a short sleeve smooth journal bearing, where presliding effects were included via Dahl’s model. The model was used to predict friction hysteresis. In a follow-up publication, the authors validated the model experimentally and extended the concept to other geometries [14].

Later, Khonsari et al. [15] conducted experiments on an oil-lubricated journal bearing under oscillatory motion. Similar effects to those of Soom et al. [8] were found when multiple loads, oscillation frequencies and lubricant viscosities were tested. It was argued by the authors that the squeeze effect contributes to the friction hysteresis. Continuing, Khonsari [16] developed a model based on decoupling the steady and unsteady terms in the Reynolds equation and Greenwood–Williamson model for rough surface contacts to predict the instantaneous friction in sliding line contacts. This simple semi-analytical model reveals that the squeeze effect plays an important role in contacts, subject to time-dependent sliding velocities.

Lui et al. [17] investigated the film thickness and the rate of change in film thickness of smooth surfaces lubricated under hydrodynamic conditions. A numerical multi-grid model was developed to solve the isothermal Reynolds equation, which was able to predict the observed film thickness hysteresis in hydrodynamic lubrication conditions. The authors showed that the squeeze effect was the dominant mechanism for the time lag occurring in the film formation between smooth surfaces by comparing experimental and numerical results. However, no friction forces were reported, nor was the hysteresis behavior further investigated in boundary or mixed friction lubrication regimes.

Friction hysteresis not only occurs in lubricated contacts but can also be experimentally observed in dry contacts [18], and it is argued by Al-Bender [19] that within dry contacts, friction hysteresis originates from a time dependency in the local adhesion coefficient between interacting asperities.

Thus far, several studies have revealed the existence of friction hysteresis and have provided some correlation between friction hysteresis and system properties (e.g., load, velocity, viscosity, and inertia); however, (i) the influence of system parameters on friction hysteresis and (ii) the influence of roughness and surface texturing are still open questions. In addition, most of the developed models neglect the inertia of the system. In this work, a numerical model is developed that enables the study of friction hysteresis in rough and textured contacts operating in boundary, mixed and full-film lubrication. The model is used to assess the influence of the aforementioned parameters (e.g., load, velocity, viscosity, and inertia) on the friction hysteresis for (i) generic surface profiles and (ii) generic velocity trajectories.

The paper is organized as follows: Section 2 covers the details of the model developed in this work, whereas Section 2.4 covers the numerical implementation details. In Section 3, numerical friction hysteresis curves obtained with variations in velocity, external applied load, viscosity and roughnessare presented and discussed. Section 5 enlists the major conclusions of this work.

2. Model Description and Implementation

The major objective of the work is to develop and implement a physics-based model that enables the simulation and the study of friction hysteresis, in contrast to the broad class of phenomenological models that are able to describe friction hysteresis.

In mixed lubrication problems, both lubricated and boundary lubrication influences occur simultaneously. This current model proposes to take the squeeze term into account and assumes a constant boundary layer shear strength. This simplification is made due to the limited contribution of the asperity adhesion in a mixed friction contact.

To include the effects of roughness and texture in this model, a multi-scale approach is followed. Three main scales are defined: (i) global, (ii) texture and (iii) roughness scales. This multi-scale approach allows to take the effect of the global and texture scale directly into account in the fluid domain, while the influence of the surface roughness is accounted for by an average flow model as developed by Patir and Cheng [20,21]. The influence of roughness on the contact is included via the statistical Greenwood–Williamson model [22].

Figure 1 schematically shows a slider bearing with a fixed counter surface and a moving sliding bearing shoe. A time-dependent normal force and sliding velocity is applied to the slider bearing, characterized by an inertia value (m). The inertia is assumed to be an equivalent inertia, taking into account the inertia of the bearing and the supported system.

Since the hydrodynamic pressures in such conformal system are typically in the order of 10–200 MPa, rigid surfaces can be safely assumed for most practical bearing materials, with the exception of very soft, rubber-like polymers.

In contrast to other approaches found in the literature ([12,16]), inertia is explicitly included in the determination of friction hysteresis. Hence, this model can be used to determine the contribution of the inertia term to the overall system behavior.

2.1. Rigid Body Motion of the Slider Bearing

2.1.1. Force Equilibrium

Ignoring bearing displacements in the tangential direction (See Figure 1), Newton’s second law is given by

$$m\frac{d^{2}h(x,t)}{d{t}^2}=F_{hyd}(t,h,\dot{h})+F_{asp}(t,h,\dot{h})-F_{ext}\left(t\right),$$

(1)

in this equation, m is the equivalent system inertia, whereas, $F_{hyd},F_{asp}$ denotes, respectively, the hydrodynamic and asperity forces, and $F_{ext}$, the externally applied load. The hydrodynamic and asperity forces are obtained through the integration of the hydrodynamic and asperity pressures over the entire domain:

$$\begin{array}{cc}\hfill F_{hyd}\left(t\right)& =B{\int }_0^L{p}_h(x,t)dx\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill F_{asp}\left(t\right)& =B{\int }_0^L{p}_{asp}(x,t)dx\hfill \end{array}$$

(3)

The total friction force on the slider bearing is determined by summing the viscous shear contribution and the shear stresses at the asperity contact:

$$\begin{array}{cc}\hfill F_f\left(t\right)& =B{\int }_0^L\left[{\tau }_{hyd}\left(t\right)+{\tau }_{asp}\left(t\right)\right]dx\hfill \\ \hfill & =B{\int }_0^L\left[\left(-\frac{1}{2}\frac{\mathsf{\partial }p}{\mathsf{\partial }x}h+\frac{\mu v}{h}\right)\right]dx+\left({\mu }_{bl}{F}_{asp}\right)\hfill \end{array}$$

(4)

A constant boundary lubrication friction coefficient (${\mu }_{bl}$) is assumed, where L and B denote, respectively, the length and width of the contact.

2.1.2. Bearing Clearance

The nominal film thickness is defined as the distance between the nominal surface of the two lubricated surfaces. The total film thickness can be considered the sum of three distinct contributions: the minimum film thickness or clearance, the texture profile, and the nominal bearing geometry. Hence, the film thickness is described as

$$\begin{array}{cc}\hfill h(x,t)& =h_0\left(t\right)+h_g\left(x\right)+h_t\left(x\right)\hfill \end{array}$$

(5)

where $h_0\left(t\right)$ describes the minimum film thickness, $h_g\left(x\right)$ models the film thickness on the global scale, and $h_t\left(x\right)$ describes the texture profile.

2.2. Fluid

A transient isothermal Reynolds equation is implemented to determine the time-dependent fluid pressure within the lubricated contact (see Figure 1). In this section, the models are presented and discussed.

2.2.1. Transient Reynolds Equation

The fluid in the gap is described with the 1D time-dependent average Reynolds equation in the domain $x\in [0,L]$:

$$\frac{\mathsf{\partial }}{\mathsf{\partial }x}\left(\frac{\rho h^3}{12\mu }{\mathsf{\Phi }}_x^P\frac{\mathsf{\partial }p}{\mathsf{\partial }x}\right)=\frac{U}{2}\left(\frac{\mathsf{\partial }\rho h}{\mathsf{\partial }x}+R_q\frac{\mathsf{\partial }{\mathsf{\Phi }}_x^S}{\mathsf{\partial }x}\right)+\frac{\mathsf{\partial }\rho h}{\mathsf{\partial }t}$$

(6)

$$\begin{array}{cc}\hfill p(x=0)& =p_0\hfill \\ \hfill p(x=L)& =p_L\hfill \end{array}$$

(7)

In the above equations, p is the fluid pressure, h the average lubricant film thickness, U the sliding velocity ($U=\frac{U_1+U_2}{2}$), and $\rho$ and $\mu$ the lubricant density and dynamic viscosity, respectively. We assume one of the surfaces to be stationary, i.e., $U_2=0$ and the other surface to be subjected to a time-dependent changing velocity $U_1=v\left(t\right)$. The lubricant’s density and viscosity are modeled using an equation of state (Section 2.2.2) to account for cavitation. The surface roughness effects are taken into account via a flow-factor method according to the classical Patir and Cheng method (P&C) [20,21].

Two flow factors are introduced: a pressure factor (${\mathsf{\Phi }}_x^P$) accounting for the change in mean flow due to the pressure-induced flow of the surface roughness, and a shear flow factor (${\mathsf{\Phi }}_x^s$) accounting for the fluid transport due to the rough surfaces. The combined root mean square surface roughness is given by $R_q=\sqrt{(R_{q,1}^2+R_{q,2}^2)}$, where $R_{q,1}$ and $R_{q,2}$ are the root mean square surface roughness of the top and bottom surfaces, respectively.

The use of flow factors as proposed by Patir and Cheng [20,21] is only valid for a low load and limited area of contact. Methods have been proposed in the literature to overcome these inherent disadvantages [23,24]; however, these are considered out of the scope of the current work. In this work, the P&C method provides a relatively straightforward and uncomplicated approach for incorporating the impact of roughness on the generated hydrodynamic pressure and film thickness. Furthermore, it is a well-established method for including roughness in a mixed friction model, the advantages, disadvantages, strengths, and weaknesses of which are extensively documented in the literature.

The flow factors are determined as proposed by P&C, using the following relationship (Equation (8)) and parameters for a Gaussian surface (see

Table A1. Flow factor correction coefficients, according to [20,21].

Pressure Flow Factor Coefficients
$\mathit{\gamma }$	C	r			Range
⅑	1.48	0.42			$H>1$
⅙	1.38	0.42			$H>1$
⅓	1.18	0.42			$H>0.75$
1	0.90	0.56			$H>0.50$
3	0.225	1.50			$H>0.50$
6	0.520	1.50			$H>0.50$
9	0.870	1.50			$H>0.50$
Shear Flow Factor Coefficients
$\mathit{\gamma }$	${\mathit{A}}_{\mathbf{1}}$	${\mathit{\alpha }}_{\mathbf{1}}$	${\mathit{\alpha }}_{\mathbf{2}}$	${\mathit{\alpha }}_{\mathbf{3}}$	${\mathit{A}}_{\mathbf{1}}$
⅑	2.046	1.12	0.78	0.03	1.856
⅙	1.962	1.08	0.77	0.03	1.754
⅓	1.858	1.01	0.76	0.03	1.561
1	1.899	0.98	0.92	0.05	1.126
3	1.560	0.85	1.13	0.08	0.556
6	1.290	0.62	1.09	0.08	0.388
9	1.011	0.54	1.07	0.08	0.295

):

$$\begin{array}{cc}\hfill {\mathsf{\Phi }}_x(\mathsf{\Lambda },\gamma )& =\left\{\begin{array}{cc}1-C{\mathrm{e}}^{r\mathsf{\Lambda }}\hfill & \mathsf{\Lambda }\le 1\hfill \\ 1+C{\mathrm{e}}^{-r\mathsf{\Lambda }}\hfill & \mathsf{\Lambda }>1\hfill \end{array}\right.\hfill \\ \hfill {\mathsf{\Phi }}_s(\mathsf{\Lambda },\gamma )& =\left\{\begin{array}{cc}{A}_1{\mathsf{\Lambda }}^{{\alpha }_1}{\mathrm{e}}^{-{\alpha }_2+{\alpha }_3{\mathsf{\Lambda }}^2}\hfill & \mathsf{\Lambda }\le 5\hfill \\ A_2{\mathrm{e}}^{-0.25\mathsf{\Lambda }}\hfill & \mathsf{\Lambda }>5\hfill \end{array}\right.\hfill \end{array}$$

(8)

where $\mathsf{\Lambda }$ is the specific film thickness, $\mathsf{\Lambda }={\displaystyle \frac{h}{R_q}}$, and $\gamma$ is the surface lay parameter. In this case, $\gamma =1$, indicating that an isotropic surface is used. The other parameters used in this work are listed in

Table A2. Parameters used for the steady-state simulation.

Parameter	Value	Unit
$F_{ext}$	5	$\mathrm{k}\mathrm{N}$
${\mu }_l$	0.0564	$\mathrm{Pa}\mathrm{s}$
${\rho }_l$	850	$\mathrm{k}\mathrm{g}{\mathrm{m}}^3$
${\mu }_g$	$1.85\times {10}^{-5}$	$\mathrm{Pa}\mathrm{s}$
${\rho }_g$	0.1161	$\mathrm{k}\mathrm{g}{\mathrm{m}}^3$
$p_s$	100,000	$\mathrm{Pa}$
$p_v$	1000	$\mathrm{Pa}$
$C_1$	0	-
$C_2$	1	-
$p_{x=0,x=l}$	100,000	$\mathrm{Pa}$
B	100	$\mathrm{m}\mathrm{m}$
L	$76.2$	$\mathrm{m}\mathrm{m}$
Mesh	$256\times 1$	-
	$v\left(t\right)=V_0$
$V_0$	[−5–5]	$\mathrm{m}{\mathrm{s}}^{-1}$
number of steps	60	-
${\mathsf{\Phi }}_x^P$	1	-
${\mathsf{\Phi }}_x^S$	0	-
$\eta$	$97\times {10}^9$	$\mathrm{m}$
$\sigma$	5	$\mathsf{\mu }\mathrm{m}$
$A_0$	$BL$	$\mathsf{\mu }{\mathrm{m}}^2$
$\beta$	1.5	$\mathsf{\mu }\mathrm{m}$
E	115	$\mathrm{G}\mathrm{Pa}$
${\mu }_{bl}$	0.3	-
${\theta }_1$	2	-
${\theta }_2$	$\frac{4}{5}$	-
${\omega }_{motion}$	0.5	-
Max. Motion Iterations	100	-
Motion Tolerance Error	$1\times {10}^{-4}$	-
${\omega }_{fluid}$	[$5\times {10}^{-4}$–$5\times {10}^{-2}$]	-
Max. Fluid Iterations	10,000	-
Fluid Tolerance Error	$1\times {10}^{-4}$	-

,

Table A3. Operational conditions for the unidirectional sliding simulations.

Parameter	Value	Unit
m	1	$\mathrm{k}\mathrm{g}$
$F_{ext}$	5	$\mathrm{k}\mathrm{N}$
${\mu }_l$	0.0564	$\mathrm{Pa}\mathrm{s}$
${\rho }_l$	850	$\mathrm{k}\mathrm{g}{\mathrm{m}}^3$
${\mu }_g$	$1.85\times {10}^{-5}$	$\mathrm{Pa}\mathrm{s}$
${\rho }_g$	0.1161	$\mathrm{k}\mathrm{g}{\mathrm{m}}^3$
$p_s$	100,000	$\mathrm{Pa}$
$p_v$	1000	$\mathrm{Pa}$
$C_1$	0	-
$C_2$	1	-
$p_{x=0,x=l}$	100,000	$\mathrm{Pa}$
B	100	$\mathrm{m}\mathrm{m}$
L	$76.2$	$\mathrm{m}\mathrm{m}$
Mesh	$256\times 1$	-
	$v\left(t\right)=V_0+V_{1}sin\left(2\pi ft\right)$
$V_0$	$0.25,0.5,1$	$\mathrm{m}{\mathrm{s}}^{-1}$
$V_1$	$\frac{V_0}{2}$	$\mathrm{m}{\mathrm{s}}^{-1}$
f	$\frac{1}{4},1,4$	$\mathrm{Hz}$
${\mathsf{\Phi }}_x^P$	1	-
${\mathsf{\Phi }}_x^S$	0	-
$\eta$	$97\times {10}^9$	$\mathrm{m}$
$\sigma$	5	$\mathsf{\mu }\mathrm{m}$
$A_0$	$BL$	$\mathsf{\mu }{\mathrm{m}}^2$
$\beta$	1.5	$\mathsf{\mu }\mathrm{m}$
E	115	$\mathrm{G}\mathrm{Pa}$
${\mu }_{bl}$	0.3	-
${\theta }_1$	2	-
${\theta }_2$	$\frac{4}{5}$	-
${\omega }_{motion}$	0.5	-
Max. Motion Iterations	100	-
Motion Tolerance Error	$1\times {10}^{-4}$	-
${\omega }_{fluid}$	[$5\times {10}^{-4}$–$5\times {10}^{-2}$]	-
Max. Fluid Iterations	10,000	-
Fluid Tolerance Error	$1\times {10}^{-4}$	-

and

Table A4. Operational conditions for the reciprocating sliding simulations.

Parameter	Value	Unit
m	1	$\mathrm{k}\mathrm{g}$
$F_{ext}$	5	$\mathrm{k}\mathrm{N}$
${\mu }_l$	0.0564	$\mathrm{Pa}\mathrm{s}$
${\rho }_l$	850	$\mathrm{k}\mathrm{g}{\mathrm{m}}^3$
${\mu }_g$	$1.85\times {10}^{-5}$	$\mathrm{Pa}\mathrm{s}$
${\rho }_g$	0.1161	$\mathrm{k}\mathrm{g}{\mathrm{m}}^3$
$p_s$	100,000	$\mathrm{Pa}$
$p_v$	1000	$\mathrm{Pa}$
$C_1$	0	-
$C_2$	1	-
$p_{x=0,x=l}$	100,000	$\mathrm{Pa}$
B	100	$\mathrm{m}\mathrm{m}$
L	$76.2$	$\mathrm{m}\mathrm{m}$
Mesh	$256\times 1$	-
	$v\left(t\right)=V_0+V_{1}sin\left(2\pi ft\right)$
$V_0$	0	$\mathrm{m}{\mathrm{s}}^{-1}$
$V_1$	$0.25,1$	$\mathrm{m}{\mathrm{s}}^{-1}$
f	$\frac{1}{4},1,4$	$\mathrm{Hz}$
${\mathsf{\Phi }}_x^P$	1	-
${\mathsf{\Phi }}_x^S$	0	-
$\eta$	$97\times {10}^9$	$\mathrm{m}$
$\sigma$	5	$\mathsf{\mu }\mathrm{m}$
$A_0$	$BL$	$\mathsf{\mu }{\mathrm{m}}^2$
$\beta$	1.5	$\mathsf{\mu }\mathrm{m}$
E	115	$\mathrm{G}\mathrm{Pa}$
${\mu }_{bl}$	0.3	-
${\theta }_1$	2	-
${\theta }_2$	$\frac{4}{5}$	-
${\omega }_{motion}$	0.5	-
Max. Motion Iterations	100	-
Motion Tolerance Error	$1\times {10}^{-4}$	-
${\omega }_{fluid}$	[$5\times {10}^{-4}$–$5\times {10}^{-2}$]	-
Max. Fluid Iterations	10,000	-
Fluid Tolerance Error	$1\times {10}^{-4}$	-

. We should note that if $\frac{h}{\sigma }\to \infty$, the average Reynolds equation reduces to the classical Reynolds equation [21].

2.2.2. Equation of State

In the case of textured contacts, recent research has shown that taking cavitation should be taken into account to accurately predict the load-bearing capacity as well as the resulting friction [25]. In addition, cavitation should also be taken into account when a contact is subject to transient operation conditions. During transient conditions, changes in load or speed occur, leading to potentially rapid changes in the pressure distribution of the lubricant. Therefore, to accurately predict the actual lubricant pressure distribution, it is essential to take cavitation into account.

In this work, we opted to use a homogeneous mixture approach to model the lubricant when subjected to pressures lower than the saturation pressure:

$$\begin{array}{cc}\hfill \rho \left(\alpha \right)& ={\rho }_l(1-\alpha )+{\rho }_g\alpha \hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill \mu \left(\alpha \right)& ={\mu }_l(1-\alpha )+{\mu }_g\alpha \hfill \end{array}$$

(10)

Therefore, a vapor volume fraction coefficient $\alpha$ is defined (bound in range $0\le \alpha \le 1$). This coefficient reassembles the mixture between the liquid and gaseous–vapor phases present in the lubricated contact. An empirical barotropic equation of state is used, which serves the purpose of a simple phenomenological cavitation model. An analogous formulation was used by [26]:

$$\begin{array}{cc}\hfill \alpha \left(p\right)& =\frac{1}{2}\left[1-tanh\left(\beta \left(p-p0\right)\right)\right]\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \beta & =\frac{1}{2\left(p_s-p_v\right)}ln\left(\frac{{\alpha }_s\left(a_v-1\right)}{a_v\left({\alpha }_s-1\right)}\right)\hfill \end{array}$$

(12)

$$\begin{array}{cc}\hfill p_0& =p_s-\frac{1}{2\beta }ln\left(\frac{-{\alpha }_s}{{\alpha }_s-1}\right)\hfill \end{array}$$

(13)

The coefficients $\beta$ and $p_0$ determine, respectively, the steepness of the transition between vapor/gas and liquid, and the pressure at which the share of vapor and liquid are equal. The $\beta$ coefficient is strictly positive. The parameters follow from the assumption of predefined volume fractions at the saturation and vapor fraction (${\alpha }_s$ and ${\alpha }_v$). The lubricant’s properties, density and viscosity, are subsequently modeled as a homogeneous mixture. Figure 2 shows the variation in the vapor volume fraction with pressure with varying model parameters ($p_0$ and $\beta$).

2.3. Greenwood–Williamson Contact Model

The asperity contact pressure is modeled by means of the statistical Greenwood–Williamson model (GW) [22]. This model leads to a description of the area of contact ($A_r$) and mean contact pressure ($p_c$) as a function of the normalized vertical separation ($\mathsf{\Lambda }=\frac{h}{R_q}$) between the two surfaces:

$$\begin{array}{cc}\hfill A_r\left(\mathsf{\Lambda }\right)& =\pi \eta \beta A_0{F}_1\left(\mathsf{\Lambda }\right)\hfill \end{array}$$

(14)

$$\begin{array}{cc}\hfill p_c\left(\mathsf{\Lambda }\right)& =E^{*}\eta {\beta }^{\frac{1}{2}}{R_q}^{\frac{3}{2}}{A}_0{F}_{\frac{3}{2}}\left(\mathsf{\Lambda }\right)\hfill \end{array}$$

(15)

with

$$\begin{array}{cc}\hfill F_n\left(\mathsf{\Lambda }\right)& =\frac{1}{2\pi }{\int }_{\mathsf{\Lambda }}^{\infty }{\left(u-\mathsf{\Lambda }\right)}^n\varphi \left(u\right)du,\hfill \end{array}$$

(16)

where $A_0$, $\eta$, $\beta$ and $E^{*}$ are, respectively, the real area, asperities density, mean curvature, and the equivalent Young’s modulus. The surface height distribution is given with $\varphi \left(u\right)$.

The GW contact theory was employed in this study as a widely accepted statistical asperity contact model for simulating mixed friction [27]. Statistical multi-asperity contact theories assume that the real contact area is relatively small compared to the nominal or apparent contact area. Since our primary objective is to investigate and understand the friction hysteresis induced by reciprocating and oscillating motions, which primarily arises from the transition between mixed and full-film lubrication, it is not necessary to consider large contact areas relative to the apparent contact area. Therefore, the selection of the commonly used GW model is considered appropriate. In addition, the choice of the GW model is justified due to its ease of implementation compared to more complex and accurate theories, such as the Persson contact model [28], or computationally expensive deterministic contact models like the boundary element method (BEM) and finite element method (FEM) [27]. We do not anticipate significant influences from more complex contact models on the observed trends (e.g., friction hysteresis with viscosity, inertia, or load) and the conclusions drawn in this study. Moreover, the advantages, disadvantages, strengths, and weaknesses of the GW contact theory are extensively documented in the existing literature [27].

2.4. Numerical Implementation

The model developed in the previous section consists of the Reynolds equation, including cavitation (Equation (6)), statistical contact equations (Equation (15)), and the load balance, including inertia (Equation (1)). The resulting system of non-linear equations is solved using the following approach: the equations of motion are solved, and the minimum gap separation is found by solving the load balance using a quasi-Newton secant method. The partitioning between hydrodynamic load and the asperity load is based on the value of the dimensionless film thickness, and follows from the load balance equation.

Algorithm 1 gives an overview of the solution procedure:

Algorithm 1: Main algorithm

2.4.1. Discretization of the Equilibrium Equation

Following Profito et al. [29], the equilibrium equation is discretized using a fourth-order backward difference formula (BDF). The fourth-order BDF scheme is fully implicit, and the stability of it is based on preselected coefficients.

In this study, the same coefficients are adopted as those in the work of Profito et al. [29]; namely, a value of $({\theta }_1,{\theta }_2)=(2,\frac{4}{5})$ is used.

The equilibrium equation (Equation (1)) can be written as an implicit function of the minimum film thickness at time n only. In order to find the equilibrium minimum film thickness at each time instance, a quasi-Newton method (Secant) is used to minimize the residual force ($\mathsf{\Delta }{F}_{res}$) of the following error function:

$$\mathsf{\Delta }{F}_{res}=F_{hyd}(t,h,\dot{h})+F_{asp}(t,h,\dot{h})-F_{ext}\left(t\right)-F_{in}(t,h,\dot{h})$$

(17)

An under-relaxation factor is used to ensure the stability of the secant procedure, and a minimum film height is imposed ($h>\sqrt{(R_{q,1}^2+R_{q,2}^2)}$).

2.4.2. Discretization and Solution of the Reynolds Equation

The transient Reynolds Equation (6) is discretized using a central difference scheme. An iterative procedure with an under-relaxation factor (${\omega }_k^{fluid}$) is used to stabilize the convergence during the solution procedure. Cavitation is taken into account via a specific equation of state (Equation (9)), and negative values of the linearized solution are set to zero. The initial pressure at the n-th time-step ($P_n^0$) is obtained from the last time step $\left(P_{n-1}^{k_{max}}\right)$. The pressure and lubricant properties are updated until convergence is obtained via the solution of the linearized system:

$$\begin{array}{cc}\hfill \widehat{P}& ={\left[LHS\right]}^{-1}\left[RHS\right]\hfill \\ \hfill \delta p& =max(\widehat{P},0)-P_n^k\hfill \\ \hfill P_n^{k+1}& =P_n^k+{\omega }_k^{fluid}\delta p\hfill \end{array}$$

(18)

2.4.3. Discretization and Solution of the Contact Equation

The integrals of the Greenwood–Williamson model (Equations (14)–(16)) are numerically solved for a Gaussian distribution and tabulated in a pre-processing step. A linear 1D interpolation is used to determine the value of the GW integral during the simulation in order to calculate the instantaneous contact area, asperity load and asperity friction.

3. Simulation Results and Discussion

The model is applied to asses friction hysteresis in lubricated contacts subject to a sinusoidal sliding velocity. In this section, the model is validated qualitatively against data in the literature [4,8,15,16]. Due to the lack of complete and detailed experimental data in the literature, a quantitative comparison is not possible. Simulations were carried out for three different (motion) cases: (i) unidirectional, constant velocity sliding; (ii) reciprocating sliding; and (iii) unidirectional, variable velocity sliding.

The simulation case (

Table 1. Velocity profiles used in this work.

Simulation Case	Applied Velocity Profile
Steady State Sliding	$v\left(t\right)=V_0$
Reciprocating Sliding	$v\left(t\right)=V_{1}sin\left(2\pi ft\right)$
Unidirectional Sliding	$v\left(t\right)=V_0+V_{1}sin\left(2\pi ft\right)$ with $V_1\le V_0$

) involves a parabolic slider as described in references [29,30,31], and the parabolic slider is given in Figure 3. The friction hysteresis of the slider is qualitatively compared with the experimental data of Hess et al. and Xiaobin et al. [8,15] and the numerical data of Zhai et al. and Sojoudi et al. [12,16].

3.1. Constant Velocity-Stribeck Curve

First, Stribeck curves are obtained for steady-state sliding motions, which serve as a benchmark solutions for the later unsteady results.

Figure 4 shows six Stribeck curves, corresponding to, respectively, two loads and three lubricants with different viscosities, using Table A2 as the input parameters. The simulation is conducted for two external loads, 5720 N and 20,762 N, and three lubricant viscosities of $0.056\mathrm{Pa}\mathrm{s}$, $0.103\mathrm{Pa}\mathrm{s}$ and $0.143\mathrm{Pa}\mathrm{s}$, which were selected to represent SAE20, SAE30 and SAE40 lubricants at 100 ${}^{\circ }\mathrm{C}$.

The proposed model is able to recreate a typical Stribeck curve with the distinct lubrication regimes of boundary, mixed and full-film lubrication. Starting at a low sliding velocity, a high COF is observed almost equal to the boundary friction coefficient (${\mu }_{bl}$). Increasing the velocity causes a transition to the mixed lubrication regime, which is characterized by a decreasing COF, which eventually reaches a minimum, indicating the demarcation between mixed and full-film lubrication. Increasing the velocity further will result in an increase in the COF due to the increase in viscous losses.

From Figure 4, it is evident that a change in viscosity or applied load causes a shift of the Stribeck curve, leading to different friction values for a specific sliding velocity. This behavior is qualitative in nature and in accordance with other published results in the literature [32].

3.2. Reciprocating Friction Hysteresis

In the following, we show the simulation results of friction hysteresis in the case that a sinusoidal reciprocating motion is applied. The tabulated input parameters (Table A4 in Appendix B) are used unless otherwise noted.

3.2.1. Friction Hysteresis Loops

Friction hysteresis is visualized by plotting the friction force normalized by an external load or film thickness versus the sliding velocity. In the following paragraphs, the influences of frequency, amplitude, load and viscosity are further studied and compared qualitatively with the results of [8,15] and the numerical data of [12,16].

Effect of Velocity Amplitude and Frequency

In order to investigate the influences of the sliding velocity amplitude ($V_1$) of the sliding velocity on the friction hysteresis, steady-state as well as reciprocating sliding conditions (as described in Table 1) were used. One steady-state simulation ($V_0\in [-1,1]$) and nine reciprocating cases with different amplitudes (0.1, 0.5 and 1) m/s for three frequencies of 1/4, 1 and 4 Hz were performed. For each simulation, an external load of 20,762 $\mathrm{N}$ was applied, and a viscosity of the lubricant of $0.056\mathrm{Pa}\mathrm{s}$ was used.

Figure 5 and Figure 6 demonstrate the effect of frequency and amplitude on the frictional and film thickness hysteresis loop of a sliding contact subjected to a sinusoidal sliding velocity. As a reference, the steady-state curves are also added to the plots.

The results presented in Figure 5 demonstrate a clear increase in the area in the friction hysteresis with increasing velocity amplitude, which becomes more pronounced at higher frequencies. This phenomenon is attributed to the squeeze effect: under steady-state conditions, the increasing velocity amplitude leads to an increase in the fluid film thickness and a consequent reduction in the contribution of asperities to the total friction force. However, the build-up of the fluid film is slower than predicted from the steady-state conditions due to the squeeze effect, resulting in a larger friction force than anticipated from the steady-state simulations. The opposite effect occurs for decreasing velocities, resulting in a higher film thickness and a lower friction force. This observation is supported by Figure 6, which shows an increased film height for decreasing sliding velocities and a decreased film height for increasing sliding velocities compared to steady-state conditions. The increase in sliding velocity amplitude further amplifies the squeeze contribution, leading to larger hysteresis loops.

Furthermore, the results indicate the presence of a non-zero friction force at zero sliding velocity, which is believed to be due to the squeeze effect. The squeeze contribution leads to a non-zero film thickness at zero sliding velocity, while the steady-state prediction suggests that the film thickness should be zero at this point. Consequently, the contact operates in a mixed friction regime, resulting in a reduction in the friction force. However, as the film thickness reaches a minimum after a certain time, the contact predominantly operates in a boundary friction regime, leading to the occurrence of the peak friction value at a non-zero velocity as illustrated in Figure 6.

Effect of Externally Applied Load on Friction Hysteresis

To investigate the influence of the externally applied load on the friction hysteresis, nine simulations of a reciprocating motion were performed, combining three external applied loads of $5762\mathrm{N}$, $\mathrm{10,762}\mathrm{N}$ and $\mathrm{20,762}\mathrm{N}$, each at three frequencies of 0.25 Hz, 1 Hz and 4 Hz. The results for the friction hysteresis are displayed in Figure 7.

The obtained friction hysteresis loops are similar to the ones measured by Harnoy and Friedland [13] and Lu and Khonsari [15]. Figure 7 highlights three important observations regarding the effect of an externally applied load on the friction hysteresis. Firstly, the average friction force increases with the increasing load, regardless of the frequency. Secondly, lower frequencies result in smaller loops, and the friction force values near the point of motion reversal are close to each other. Finally, higher frequencies result in larger loops with clearly different values of the friction force near motion reversal, and a maximum COF that shifts towards the maximum sliding velocity. The obtained results indicate that the size of the friction hysteresis loop decreases with the increasing load, while the peak friction force value increases, indicating more asperity contact. This observation is consistent with the fact that an increase in the external load allows a larger part of the load to be carried by the asperities, thereby reducing the friction loops compared to a lower external load. Furthermore, the squeeze contribution for a higher external load has a smaller influence on the change in film thickness, which results in a smaller friction hysteresis loop.

Effect of the Lubricants Viscosity on Friction Hysteresis

In order to investigate the influence of viscosity on friction hysteresis, nine simulations of a reciprocating motion were performed, combining three viscosity values of $0.056\mathrm{Pa}\mathrm{s}$, $0.103\mathrm{Pa}\mathrm{s}$ and $0.143\mathrm{Pa}\mathrm{s}$, each at three frequencies of 0.25 Hz, 1 Hz and 4 Hz. The results for the friction hysteresis are displayed in Figure 8.

Figure 8 reveals that, on average, the friction force decreases as the viscosity increases, regardless of the frequency. Additionally, as the frequency increases, the magnitude of the loops increases for lower viscosity cases but decreases for those with higher viscosities. This phenomenon is accompanied by a change in the shape of the loops, which ultimately results in a shift of the maximum friction force towards the maximum value of the sliding velocity.

Effect of the Slider Inertia on Friction Hysteresis

To illustrate and asses the influence of the slider’s inertia on the resulting friction hysteresis behavior, 12 simulations of a reciprocating motion were performed, combining 4 inertia values of $0.1,1,10$ and $100\mathrm{k}\mathrm{g}$ each at 3 frequencies of 0.25 Hz, 1 Hz and 4 Hz. The results for the friction hysteresis are displayed in Figure 8.

Figure 9 illustrates the negligible differences between the resulting friction hysteresis loops obtained for the four values of the mass. To further investigate the effect of inertia, the maximum relative contribution of the inertial force to the total externally applied force is determined and visualized in Figure 10. From Figure 10, it can be seen that the overall contribution of the inertia term is negligible under the simulated conditions. The inertia contribution ranges from ${10}^{-6}\%$ to ${10}^{-3}\%$:

$$M=max\left(\frac{F_{ini}}{F_{ext}}\right)=max\left(\frac{m}{F_{ext}}\frac{d^{2}Z}{d{t}^2}\right)$$

(19)

By expanding Equation (19), it can be easily shown that the contribution of the inertia scales proportionally to the mass (m) and quadratically to the sliding frequency (f) as shown in Figure 10. Hence, for a system under moderate operational conditions, as simulated in this work, inertia can be safely neglected as was done by other authors ([12,16]).

3.3. Uni-Directional Friction Hysteresis

Next, a sinusoidal velocity is superimposed on the constant-velocity motion of the sliding inertia (called the offset velocity), and the instantaneous friction response is simulated and assessed. Simulations are performed for a unidirectional sliding motion, combining three offset velocities, i.e., 0.05 m s${}^{-1}$, 0.25 m s${}^{-1}$ and 0.50 m s${}^{-1}$ with three oscillation frequencies $0.25$ $\mathrm{Hz}$, 1 $\mathrm{Hz}$ and 4 $\mathrm{Hz}$. The amplitude of the oscillations is assumed to be half of the offset velocity. Table A3 provides an overview of the remaining simulation parameters.

As clearly shown in Figure 11, the friction hysteresis loops are located around the steady-state Stribeck curve. Moreover, the offset velocity plays a critical role in determining the prevailing lubrication regime. Specifically, the regime can vary from boundary to mixed to hydrodynamic lubrication as the offset velocity increases. In addition, the area of the friction hysteresis loops are found to increase with an increase in frequency, with the largest areas being observed for mixed lubrication. Lastly, the study reveals that an increase in frequency near the minimum of the Stribeck curve at the beginning of the hydrodynamic lubrication results in a decrease in the coefficient of friction due to the significant impact of the squeeze effect in this region.

Another interesting observation is that the magnitude of the observed friction hysteresis loop decreases when operating closer to full-film lubrication. The friction hysteresis loop has a rather small magnitude compared to the friction hysteresis loop in the mixed lubrication regime. Variations in the friction force are due to variations in the film thickness, and the observed friction hysteresis loop in uni-directional sliding does not only depend on the amplitude and frequency, as is the case with bi-directional sliding, but also on the offset velocity.

Effect of the Surface Roughness on the Friction Hysteresis

Equivalent conclusions for the influence of load, amplitude and viscosity on the magnitude of the friction hysteresis phenomena can be obtained for the uni-directional sliding case. Hence, this is here omitted for conciseness. To study the influence of the surface roughness, nine simulations are conducted, where the root mean square surface roughness ($R_q$) is varied in steps of $2.5\mathsf{\mu }\mathrm{m}$, $5.0\mathsf{\mu }\mathrm{m}$ and $7.5\mathsf{\mu }\mathrm{m}$, while keeping other roughness parameters constant. Other model parameters are as shown in Table A3 in Appendix B.

Figure 12 presents the influence of surface roughness on the size and shape of the friction hysteresis loop. For an increase in the root mean square surface roughness, while keeping other surface parameters (e.g., average radius $\beta$ and asperity density $\eta$) constant, the friction coefficient will increase due to an increase in the asperity contact [16,22]. This corresponds to a shift of the Stribeck curve to the right. For comparable conditions, the increase in surface roughness leads to a decrease in the size of the hysteresis loop. This is due to the larger contribution of asperities in the load sharing, leading to a lower contribution of the hydrodynamic and, subsequently, squeeze contribution as shown in Figure 13. In Figure 13, we clearly observe that in the case of $\sigma =2.5\mathsf{\mu }\mathrm{m}$, the majority of the external load is supported by the fluid film (>75%), while when the root mean square surface roughness is increased to $\sigma =7.5\mathsf{\mu }\mathrm{m}$, the balance shifts, and the majority of the external load (>75%) is taken by the asperities. Note, however, that for an increase in the sliding velocity, the proportion of asperity force decreases, and the hydrodynamic proportion increases.

4. Normalization—A First Attempt

The use of normalization is often used to condense steady-state results (numerical or experimental) into a characteristic curve, the Stribeck curve. A typically used parameter is the Hersey number, taking into account the influence of applied load ($F_{ext}$), viscosity ($\mu$) and sliding velocity (V) on the friction force. Since, in this work, the sliding velocity is per definition time dependent, the Hersey number is not a single value but a range of Hersey numbers obtained by using Equation (20):

$$H=\frac{\mu V}{F_{ext}}$$

(20)

In Figure 14, eight friction loops for different applied loads, viscosity and frequencies are plotted versus the Hersey number. The goal is to obtain two sets of friction hysteresis curves with a similar range of Hersey numbers by varying the base viscosity and applied external load. From this plot, it can be observed that, for both Hersey ranges investigated, even though the Hersey ranges for different simulations are equal, different shapes and sizes of the COF profiles are predicted because of the unsteady nonlinear squeeze effects that govern the amplitude and frequency of the oscillations. Hence, we conclude that the Hersey number is not sufficient to predict the loop shape and magnitude, necessitating the definition of additional parameters.

5. Conclusions

In this paper, a mixed lubrication model is developed, which involves a transient Reynolds equation, a mass-conserving cavitation model, and a statistical contact equation, combined with a load balance equation. The model is used to study the influence of multiple parameters (e.g., load, viscosity, roughness, amplitude and offset velocity) on the friction hysteresis and other operating conditions. A classic mixed friction model in conjunction with the dynamics of the bearing is sufficient to allow the prediction of the friction hysteresis. The proposed model describes the experimental observed effects for different values of the model parameters (e.g., load, viscosity, velocity, and roughness). Besides performing fundamental studies on friction hysteresis for arbitrary velocity and load parameters, the model can be used to improve control design strategies for transient frictional contacts.

The following general conclusions can be drawn based on the simulation results presented in this work:

• Friction hysteresis loops are caused due to the fluid squeeze effect, leading to fluctuations in the minimum film thickness, and thus a deviation in the instantaneous friction force from the steady-state friction force value.
• With an increasing sliding frequency, the magnitude of the friction hysteresis loops, compared to lower sliding frequencies, increases. However, the exact size and shape of the friction loop depends on the lubrication regime and relative magnitude of the squeeze force. Both terms are influenced by viscosity, surface roughness and external applied load.
• The contribution of the bearing inertia scales proportionally to the bearing mass and quadratically to the frequency.
• System parameters, such as viscosity, load, roughness, etc., influence the friction hysteresis, and can be tailored within a given context and application, leading to better designed tribo-contacts in feedback loops (e.g., precision motion control, and force controlled systems).
• The Hersey number cannot be considered a unique parameter for the scaling of friction loops, as it does not properly account for the magnitude of the squeeze contribution. This will be the object of further research.
• When a uni-directional sliding velocity profile is applied, the friction hysteresis loop is centered around the friction force of the mean sliding velocity as determined from a steady-state Stribeck curve.
• When a bi-directional sliding velocity profile is applied, the instantaneous friction at zero sliding velocity decreases with an increasing frequency. Here, again, the size and shape of the hysteresis loop depends on the instantaneous lubrication regime and relative magnitude of the squeeze force.