Chapter 8. Models of Friction

The model below was employed by Hess and Soom [75].

where the Coulomb, viscous and static friction is denoted by FC, Fv and Fs respectively. vs is the characteristic velocity of the Stribeck curve. Hess and Soom reported systematic dependence of Fv and vs on lubricant and loading parameters.

Another proposal for modeling the Stribeck effect is a linearized exponential equation of Bo and Pavelescu [21]. The viscous friction term was added by Armstrong-Hélouvry [9].

where δ is an additional empirical parameter. In the literature surveyed by Bo and Pavelescu [21] they find the range of δ from 0.5 to 1. Armstrong-Hélouvry [6], [8] employs δ = 2 and Fuller [57] suggest δ = very large for systems with effective boundary lubricants. Applying δ = 2 the exponential model, Equation (8.20) becomes Gaussian model, which is nearly equivalent of the Lorentzian model of Equation (8.15).

Another simple representation of the Stribeck curve is the Tustin model [144]. The model captures the basic frictional phenomenon and simply superposes them on one another. The three features are the Coulomb friction FC, the static friction Fs with the downward bend of the negative viscous effect determined by α and the viscous friction Fv.

Linearizing the Tustin model we get:

This model although a linear one, leads to a precision of approximated friction only up to 85% in the region where the velocity comparable to its square root. In other cases where the system velocity is not so very low, another model should be used.

A popular approach is simply sewing together the already existing model elements. This is how the seven-parameters model works. Each of the seven parameters of the model represents a different friction phenomenon as presented in Table 8.2. They value principally depends upon the mechanism and the lubrication conditions, but typical values may be offered based on the works of [6], [27], [57], [75], [85], [105].

The friction is given by:

Not sliding (presliding displacement)

Sliding (Coulomb + viscous + Stribeck curve friction with frictional memory)

Rising static friction (frictional level at breakaway):

Where:

Parameter identification. The parameters, which represent a non-linear relationship with the friction force, cannot be identified by using linear methods. Kim et al proposes a possible non-linear identification method for the seven parameter model, the Accelerated Evolutionary Programming (AEP) [88] and gives the description as outlined below. AEP is a stochastic optimization technique for finding the extrema of an arbitrary function. A simple optimization problem can be stated as: determine the values of the ordered set of n parameters

which minimize/maximize the cost function f(z). AEP is a genetic algorithm type method, working with parent/child interactions to converge to the best solution. It uses two variation operators. One of them, the direction operator determines the direction of the search according to the fitness score. The other one is a zero-mean Gaussian operator, used a perturbation added to a parent in order to generate an offspring. AEP also uses a further “age” variable to enhance the diversity of the search and to prevent individuals to remain in the local minima. The Gaussian operator is incorporated with both the direction operator and the variable age. Unlike another genetic methods, AEP generates only one child from its parent. New parents then are selected in one-to-one competition between parent and its child. An offspring is selected if the child wins in the competition with its parent. AEP incorporates the direction of evolution and the age concept into the vector to be optimized in order to increase the convergence speed without loosing the diversity. Thus it brings Equation (8.26) into the form of

where $$dir(z_j^i)\in \left\{-1,1\right\}$$for all$$z_j^i$$$$z_j^i$$is the evolving direction of the jth parameter $$z_j^i$$in $$z_j^i$$the ith vector, and agei denotes the duration of the life in an integer type in the ith vector. The following two rules are applied for perturbing the parents to generate their offspring.

Rule 1

When $$f(z^i[k])<f(z^i[k-1])$$,

then $$dir(z_j^i[k])=\mathrm{sgn}(z_j^i[k]-z_j^i[k-1])$$,

$$ag{e}^i[k]=1$$

or else $$ag{e}^i[k]=ag{e}^i[k-1]+1$$,

Rule 2

When $$ag{e}^i[k]=1$$,

then $${\sigma }^i={\beta }_1\cdot f(z^i[k])$$,

or else $${\sigma }^i={\beta }_2\cdot f(z^i[k])\cdot age$$,

$$z_j^i[k]=z_j^i[k-1]+\left|N(0,{\sigma }^i)\right|$$,

where $$z_j^i[k]$$denotes the jth parameter in the ith vector among Np vectors at the kth generation, $$dir(z_j^i[k])$$denotes the evolving direction of $$z_j^i[k]$$, $$dir(z_j^i[k])\cdot \left|N(0,{\sigma }^i)\right|$$is a realization of a Gaussian-distributed random variable, which is polarized in the direction of $$dir(z_j^i[k])$$, and i, i = 1, 2 are positive constants. Rule 1 indicates that if the performance of a newly generated child is better, then the search will continue towards the promising direction, by setting the direction value and generating a random number in the respective positive or negative region. Otherwise the parent will be kept, by increasing the age value. As Rule 2 “else” part shows, in the latter case the standard deviation of the Gaussian-distributed random variable is getting larger in each step, since the age is included in the product. It also makes more probable that the algorithm doesn’t stuck in the local minima but tends to find the global one.

The procedure of applying the AEP is as follows:

1.) Initialization. Generate an initial population of Np trial solutions with a uniform distribution within the given domain for the ith individual vector.

where $$z_j^i$$és $$dir(z_j^i)$$, $$\forall j\in \left\{1,2,\mathrm{...},n\right\}$$ are randomly initialized and agei is initially set to 1.

2.) Evaluation. Evaluate the fitness of each parent solution by the given cost function f(zi).

3.) Mutation. For each of the Np parents zi, generate a child by using Rules 1 and 2 successively.

4.) Evaluation. Evaluate the fitness of each new offspring.

5.) Selection by one-to-one comparison. If the newly generated child is better than its parent, then select the child and discard the parent, else vice versa. New parents are composed of Np vectors selected among 2Np vectors by one-to-one comparison. The new Np parents are sorted by fitness score.

6.) Termination check. Proceed to step 3 unless available execution time is exhausted or acceptable solution has been discovered.

The most obvious way for parameter identification is to run parallel the experiment and AEP using the same input value. Parameters are identified, when the difference between the real and calculated output values diminish or fall into the previously defined error range. Input value can be position or velocity depending on the used sensor.

For the seven parameter model, first we need to know the input v or (dx/dt) and the state variables t2 and Fs,a. For this let us define a new state variable s, which indicates, whether the system is in stick or slip.

where vcritical is specified by the resolution of the velocity/position sensor built into the real system. Depending on the change of the value of s, t2 and Fs,a are obtained as followed.

When s changes from

where Ts is the sampling time.

Let a vector z denote the set of unknown parameters of Equation (8.24).

where all the parameters are positive. The mismatched displacement error between the plant output and the identification output is

where $$\widehat{z}$$denotes the estimate of z, x(z,ti) is a displacement of the plant, x($$\widehat{z}$$,ti) is a displacement of the identification model at the sample time ti, and Ns is the total number of sample data. Consider the cost function of the squared errors is given by:

Then the parameter identification problem can be posed as an optimization problem as follows:

Research in dynamic friction modeling of rocks in boundary lubrication carried out by geophysicist led to the born of the state variable friction model. [6], [16]. For constant normal stress, the general form including n state variables iis given by

Coming from this representation, a sudden change in velocity does not cause a change in  but it does affect its time derivative. For a single state variable Ruina proposes the following frictional law [16]

Where  is the scalar state variable, L is the characteristic sliding length controlling the evolution of . The pair (v0/f0) corresponds to any convenient point on the steady-state friction-velocity curve. The steady state curve is given by:

and the steady state variable can be related to the mean lifetime of an asperity junction.

Limitations. Unfortunately the model structure itself imposes a limitation on the use of the state variable model. At standstill, the value of v/v0 in Equation (8.40) becomes zero, which gives infinitive friction force due to the result in the logarithmic function (ln 0 = ). In simulations normally it causes a problem in the integrator. A good approach to avoid it is to include a switch, which for zero velocity gives a finite friction value. However since $$\underset{x\mapsto 0}{\lim }\ln (x)\mapsto -\infty$$, it means that in the environment of very small velocities the value of the logarithmic function and thus the friction force changes very fast, resulting a strong dependence of error on the selected environment size. Hence for low velocities this model is not adequate.

Parameter identification. Dupont et al gives a detailed description on the measuring apparatus used for their experiment [49]. In principle any other methods can be used, thus we are not going into details of the original set up. The first step to identify the parameters of Equation (8.38) and Equation (8.39), one must draw the steady state velocity – friction curve (or Stribeck-curve) of the respective system. The change of velocity results a change in the system as well. This can be used to compute the values of A and B on the following way. The system’s response given to a velocity step imposed at time t0, where the pair (v1, f1) corresponds to the initial steady state, the pair (v2, f2) describes to the steady state reached after the step response, while fmax denotes the maximum friction force reached during the transient. The parameters A and B then calculated as:

The exponential decay of the state variable  following a velocity step is given as

The characteristic sliding length L can be computed from a least squares fit of the exponential to the data immediately following the velocity step.

Table 8.4 shows some state variable model parameter values calculated by Dupont et al. The used lubricants are as follows. The used lubricants are paraffin oil with maximum Saybolt viscosity of 158 and a commercial boundary lubricant paste. Test samples were made of heat-treated and oil quenched tool steel, with Rockwell hardness of 59 C.

The magnitude was less then noise level of the measurement.

In some model representations friction may be a discontinuous function of the steady-state velocity. In case of integrating discontinuous ordinary differential equations it may cause a problem, though with embedded switching functions this problem is still manageable. In addition, the integrator also has to bear with the ability to locate the point of discontinuity within the given interval. Variable-step-size variable-order methods are appropriate for this purpose and they are built in elements in all commercially available simulation software. However in the physical representation of a controller in the majority of cases we have to use fixed-step-size methods as the sampling time is fixed.

The Karnopp [84] model solves this problem by introducing a pseudo velocity and a small neighborhood of zero velocity ±DV in which the speed is set to zero. The pseudo velocity, p(t) which in the original material is referred as momentum, is integrated on the standard way.

where Fa is the force acting on the system, while Fm is the modeled friction force. The velocity of the system is set according to Equation (8.45), where M is the system mass.

The friction law is given as:

This includes the Coulomb, static and viscous terms (FC, Fs, Fv respectively) but can be upgraded with any other convenient curves (e.g. Stribeck curve).

Unlike the seven parameter model, the LuGre model is an integrated dynamic friction model. In some sources it is also referred as “integrated dynamic model”. The name comes from the abbreviation of the Lund Institute of Technology and INPG Grenoble, the two universities hosting the cooperating scientists. The LuGre model is based on the elasticity in the contact [33]. The variable z(t) represents the friction state and can be interpreted as the mean deflection of the junctions between the two sliding surfaces. The instantaneous friction Ff is given by

where σ0 describes the stiffness for a spring-like behavior for small displacements, and σ1 gives the negative viscous friction term. The state variable is updated according to:

where

Fc, Fs and Fv represents the Coulomb, Stribeck and static friction respectively. For some special cases, an additional element Fg can be added to the sum, representing the effect of gravity (e.g. in case of dealing with a vertical machining axis) or amplifier offsets. The steady-state friction is as follows:

The advantage of this model, is that besides it is valid in both rotational and linear coordinates, it describes the phenomena for normally flat surfaces as well as for rolling bearings element

Parameter identification. An identification method is presented in [2] based on the general methods used for linear dynamic systems. Note that the method includes Fg contrary to Equation (8.47)-(8.49). The parameters Fs and σ0 are easily deduced from the Dahl’s curve.

The assumption $$\dot{v}=0$$ implies Fs=F. Furthermore, if $$v<<v_s$$, then $$F\approx {\sigma }_{0}z+F_g$$and $${\sigma }_{0}g(v)\approx F_s$$. Combining (8.50) with these assumptions leads to

Using Euler’s approximation for time derivatives, a linear regression results

The authors of [1] concluded that these algorithms turned to be sensitive to position measurement noise, thus they propose to use interferometers instead of inductive sensors to eliminate the bias on $${\widehat{\sigma }}_0$$.

The original Dahl model [40], [20] predicts a frictional lag between the velocity reversals and leads to hysteresis loops. It does not contain however the Stribeck effect. To overcome on this deficiency Canudas et al proposed to transform the Dahl model into a second order linear space invariant model, which emulates the Stribeck friction by producing an overshoot in the response of friction forces. However, doing so, they also eliminated the dissipative properties of the first order model. The modified Dahl model is a first order model redesigned by Canudas et al and this time combining the Stribeck phenomenon without increasing the system state dimension.

The friction force represented by Ff is described by

where ε is a time constant proportional to 1/vsimilarly to the original Dahl model, z is a state variable should be interpreted as an “internal” friction state, vs is the characteristic velocity of the Stribeck curve, while α1, α2 and α3 are positive constants representing the Coulomb, Stribeck and viscous friction terms respectively.

The steady state friction is given by:

The M2 model was developed by the same researcher group, than that of the modified Dahl model. Yet it has not got any catchy name, even in the original source it simply referred as M2. [33]. The model is based on the separate representation of stiction Fs and dynamic friction Fd. The two states are connected through a dynamic weighting factor s, which depends on the relative velocity. Thus the total friction Ff is given as

The weighting factor s is defined by solving the equation

where τs is a time constant representing the frictional lag. s = 1 means sliding whereas s = 0 implies sticking. Hence, s can be interpreted as the “measure” of sliding, having a uniform centered deterministic distribution with respect to relative displacement [33].

The so-called Dahl effect models the behavior of sticking by employing a stiff spring ks accompanied by a damper ds. Stiction however is limited by the break-away force, thus a saturation function is included in the form of

The saturation is described by

where

The variable  is an internal variable representing the relative displacement from the position where stiction occurred. The second term serves to reset  to zero when sliding started

The Stribeck effect is a result of the combination of the saturation and the weighting function, but can be traced in the steady state characteristics as well:

where: