A novel dynamic dry friction model for applications in mechanical dynamical systems

Abstract

This work presents a novel continuously differentiable, dynamic dry friction model, with dependence on normal contact force, slip velocity, and static and dynamic coefficients of friction. The state-of-the-art, Brown and McPhee friction model depends on transitional velocity, which is an empirical parameter. A limitation with this model is that there is no generic approach to select the optimal value for the transitional velocity for a certain application. The simulation results presented in this work highlight the sensitivity of friction force with transitional velocity in Brown and McPhee’s model to obtain smooth solutions, which are supportive in control system applications. It is because control systems require jitter-free signals in order to save costs on low-pass filters. The proposed model overcomes such empirical dependence on transition velocity through a combination of an iterative methodology and an empirical parameterization. In this article, a comprehensive analysis of the proposed friction model and the contributing parameters is carried out. The algorithm to implement the proposed friction model is elaborated. The proposed friction model is then simulated and compared against Brown and McPhee’s model. using a stick-slip benchmark problem. Furthermore, the proposed force model is applied to two spatial multibody systems formulated using index 0 tangent space differential-algebraic equations. The results demonstrate how the proposed friction model can be employed in dynamical mechanical systems. The independence on transitional velocity in the proposed friction model is observed to be effective for obtaining smooth solutions that make it suitable for control system applications.

Introduction

Friction is a highly nonlinear and complex phenomenon related to the resistance to the relative motion of bodies in contact. The very first attempts to model friction phenomena include,^1,2 and since then, friction’s properties have been intensively studied, modeled, and used to predict frictional behavior in mechanical systems. The earliest friction force model dates to Coulomb’s friction force model. ³ According to this model, the friction force always opposed the relative motion between two surfaces and it is proportional to the normal force. While it was established that the friction force opposes relative motion, the Stribeck effect laid out that static friction was always higher than kinetic friction^4,5 and showed experimentally that the transition from static to kinetic friction is a continuous process. This means that for low velocities, the friction force decreases with the increase in the relative velocity. This rise of friction force combined with a drop in the relative velocity leads to a phenomenon called “stick-slip” where the contact surfaces stick until the static friction is overcome.

Coulomb friction

Though Coulomb’s model is fundamental, it is quite basic when it comes to differentiating between static and dynamic regions of friction and hence Coulomb’s model is a static model since it disregards the existence of the two distinct types of friction regimes. The Stribeck model, although it considers the transition from static to dynamic friction, it does not do so in a continuous manner. Both, Coulomb’s model and Stribeck model are discontinuous and not continuously differentiable. This makes them further difficult to implement in the case studies where continuous derivatives of friction forces are required, ⁶ for example, in applications related to sensitivity analysis, optimal control, and multibody systems simulations with friction. Furthermore, Haug suggests that introducing friction prevents the formulation of an ordinary differential equation from accurately describing the dynamic characteristics of a system.^7–10 Hence, a problem where friction must be considered needs to be formulated as a system of differential-algebraic equations (DAEs) with Lagrangian multipliers. Since the friction force depends on the normal contact force, the formulation with friction required derivatives of the external forces to satisfy the DAEs. This has been discussed in detail in subsequent sections.

Brown and McPhee’s friction model

As discussed, there is a need for a continuously differentiable friction force model that could also capture the transition of friction forces from a static to a dynamic state. Several continuous models of friction already exist. However, they neglect some of the crucial phenomena associated with friction. Berger’s model ¹¹ neglects time lag in friction. Another important contribution by Armstrong-He’louvry et al. ¹² proposes a continuous velocity-based model which captures the time-lag. However, it takes the absolute magnitudes of static and dynamic friction forces as inputs, instead of evaluating them through the magnitude of normal reaction. Certain models, based on the concept of bristle deflection, propose the usage of a state-variable to capture the micro-slip phenomena.^13,14 However, these models are more complex and often discontinuous because of addition of state variables, which in turn makes them computationally inefficient. Bristle deflection models ¹³ also do not consider the velocity dependence of friction but only displacement dependence. Due to their more complex and less efficient nature, models that include such effects are not considered for this study. Readers are advised to refer to Marques et al.^15,16 for more insight into the comparison of these friction models. Brown and McPhee ¹⁷ proposed a simple velocity-based model where friction has been considered as a continuous function of relative sliding velocity. They also compare their model with certain other contemporary continuous friction models, such as Andersson et al.’s, ¹⁴ Hollars, ¹⁸ and Specker’s model. ¹⁹ To reduce model complexity, only the main velocity-dependent characteristics of friction, that is, the Stribeck effect and viscous friction were included. What makes this friction model especially promising for dynamic systems is its ability to simulate stiction without any discontinuities in the stick-slip transition regime. This model, neglecting the viscous friction, is represented in equation (1)

$$F_f({\mu }_s,{\mu }_d,v_s,v_t)=F_n({\mu }_d\tanh \left(\frac{v_s}{v_t}\right)+\frac{({\mu }_s-{\mu }_d)\left(\frac{v_s}{v_t}\right)}{{({\left(\frac{v_s}{2v_t}\right)}^2+\frac{3}{4})}^2})$$

(1)

where ${\mu }_s$ is the static coefficient of friction, ${\mu }_d$ is the dynamic coefficient of friction, $v_s$ is the slip velocity, and $v_t$ is the transition velocity. The transition velocity $v_t$ is a difficult parameter to estimate in this model. Haug ⁷ recommends $v_t$ to be between 10h and 20h, where h is the time-step size. This model is continuously differentiable and an explicit formulation for evaluating friction force which makes it computationally efficient. This model is continuously differentiable and an explicit formulation for evaluating friction force, which makes it computationally efficient. Such attributes make this model one of the most advanced friction models from a practical perspective. However, a seemingly minor, yet crucial aspect of this model is the selection of $v_t$ which to an extent, is an empirical quantity and can vary depending on the problem being simulated. depending on the problem being simulated. There are no thumb-rules to select the optimal $v_t$ . The importance of selection of $v_t$ is demonstrated in the Results section.

Static friction model

Another friction model has been studied and proposed by Wang and Rui, ²⁰ as provided in equation (2)

$$\mu ={\mu }_s\sin \left(C{\tan }^{-1}(B{v}_s-E(B{v}_s-{\tan }^{-1}B{v}_s))\right)$$

(2)

Here ${\mu }_s$ is the static coefficient of friction and $v_s$ is the slip velocity. The constants C, B, and E are manually adjusted to suit the requirements of the system. The main problem with this model is that it does not include the dynamic coefficient of friction as a functional parameter, and it depends only on the slip velocity. In the Results section, the influence of each parameter on the friction force profile for varying slip velocity is studied in detail.

Limitations of friction models in dynamical systems

Acknowledging the impact of Brown and McPhee’s model and yet considering the drawbacks associated with it, this study focuses on an alternative model which does not include the selection of optimal $v_t$ as an empirical parameter but decides the transition range based on the difference between static and dynamic coefficients of friction and the magnitude of the normal contact force. For the same reasons, the results in this study have been compared to those obtained with Brown and McPhee’s model.

Proposed friction model

Model description

In the previous section, it was outlined why the Brown and McPhee friction model outperforms other dynamic friction models. Also, a static friction model by Wang and Rui ²⁰ was revisited, as per equation (2). It is seen from the Rui and Wang’s model that the friction force is continuously differentiable; however, it is independent of the dynamic terms namely the dynamic coefficient of friction, the normal force, and the transitional velocity, thus rendering the model useful only in static applications. The proposed friction model is motivated by the limitations of these models, and it proposes modifying the Rui and Wang’s model to leverage its benefits in terms of continuous differentiability as well as independence on transitional velocity. The proposed friction model is as follows

$$\begin{array}{rl}{F}_f\left({\mu }_s,{\mu }_d,v_s,F_n\right)=& {\mu }_s{F}_n\sin (C_{cor}{\tan }^{-1}(B_m{v}_s\\ & +2\left(B_m{v}_s-{\tan }^{-1}{B}_m{v}_s\right)))\end{array}$$

(3)

$$B_m=min(100,50e^{({\mu }_s-{\mu }_d)|F_n{|}^{0.2}})$$

(4)

The effect of each parameter in equations (2) to (4) has been studied; from the analysis elaborated in the results section, it is concluded that the parameter C is most sensitive to the friction force and hence needs to be corrected. The parameter B is the next more sensitive parameter, followed by the parameter E, the least sensitive parameter. The procedure to obtain the parameters $C_{cor}$ and $B_m$ is explained in the next section. The value for the parameter E is chosen to be “ $-2$ ” ased on the analysis from the Results section.

Algorithm to obtain corrected model parameters

The value of $B_m$ is obtained using equation (4). The value for $C_{cor}$ , on the other hand, is obtained iteratively as follows:

The function receives the values of coefficients of friction ${\mu }_s$  and ${\mu }_d$ , normal reaction $F_n$ , and relative slip velocity $v_s$ , as inputs.

The factor $B_m$ is calculated using equation (4).

The initial value of friction force $F_f$ at an arbitrary large velocity is calculated with the initial C considered as 1.

Depending on the difference between $F_f$ and the dynamic friction force value ${\mu }_d{F}_n$ , the parameter C is updated iteratively using a constant step size of magnitude 0.001. If the difference is negative, the step size is $-$ 0.001.

The corrected parameter $C_{cor}$ is then used to calculate the final value of friction force $F_f$ at a slip velocity $v_s$ as per equation (3).

Comparison of friction force models

From equations (2) and (3), it is seen that the friction force is a function of the slip velocity. The Brown and McPhee model, on the other hand, also depends on the transitional velocity, as mentioned earlier. The force profiles are plotted in Figure 1. The friction force profile is obtained for a given range of the slip velocity while keeping the normal force constant at $1$ and $100$   N. The transitional velocity in the case of equation (1) is chosen to be $0.015s^{-1}$ . The parameters in equation (2) are obtained using the procedure described in the previous section. The friction force profiles against the slip velocity produced by the two said models are seen to overlap closely in Figure 1. This can be quantified by the obtained root mean squared error (RMSE) (%) between the two models, which is found to be $0.24\mathrm{\%}$ when $F_n=100$  N and $0.88\mathrm{\%}$ when $F_n=1$  N.

Figure 1.

The comparison of friction force profiles against slip velocity.

Hence, it is established, that the proposed friction model behaves similar to the state-of-the-art friction model. This leads to the next section, wherein the proposed model is applied to a benchmark study to demonstrate the applicability of the proposed model.

Validation of proposed model through Rabinowicz’s stick-slip experiment

Figure 2 shows the stick-slip experimental set-up. This benchmark test was used by Brown and McPhee to compare their model against Andersson’s and Hollars’ models. Hence, this simulation study was carried out to compare results of the proposed friction model against Brown and McPhee’s model. A rigid body of mass $M$ , is connected to a spring of stiffness $K$ over a belt moving with a constant speed $v_0$ . The static friction coefficient between the rigid body and the belt is given by ${\mu }_s$ and the dynamic friction coefficient is given by ${\mu }_d$ . The spring-mass system is a second-order system that can be modeled as in equation (5) and the system parameters are tabulated in Table 1

$$M\ddot{x}+K(x-x_0)+F_f=0$$

(5)

Figure 2.

Rabinowicz’s stick-slip experimental setup.

Table 1.

Stick-slip experiment data.

Parameter	Value	Unit
$M$	500	kg
$K$	40	${\mathrm{Nm}}^{-1}$
$v_0$	5	${\mathrm{ms}}^{-1}$
${\mu }_d$	0.2	No unit
${\mu }_s$	0.3	No unit

Figure 3 compares the results obtained with the proposed friction model and the Brown and McPhee model on the key dynamic components of the system: position, velocity, and acceleration of the mass. It is concluded from the sudden drops in the acceleration curve that stiction was captured by both models. The transition velocity for the Brown and McPhee model was kept at 15 $h$ , where $h$ is the time step. Figure 4 compares the resisting friction force obtained with both models. It can be observed that the friction force profiles are equivalent for both models, which can be seen from the RMSE values included in Table 2.

Figure 3.

The comparison of position, velocity, and acceleration.

Figure 4.

The comparison of friction force.

Table 2.

Root mean squared errors (RMSEs) of dynamic parameters.

Parameter	RMSE (%)
$x(t)$	0.018
$\dot{x}(t)$	6.25
$\ddot{x}(t)$	2.27
$F_f(t)$	0.047

Since the system’s dynamic characteristics (position, velocity, and acceleration) almost overlap each other respectively, it can be concluded that both models perform equally well for this case study. However, such a comparison is incomplete, as in this experiment the normal contact force between the block and the belt does not change with time. Hence, case studies have been performed by applying the proposed friction model to multibody systems where the contact normal forces are time dependent, as discussed in the following section.

Multibody dynamic system formulation

The multibody dynamic systems considered in this study are made of rigid links constrained by joints. Such constrained dynamic systems are usually formulated using index 3, index 2, or index 1 DAEs, however, these DAE formulations consider only one of the position, velocity, or acceleration constraints, respectively. ⁹ Haug ¹⁰ proved that the index 0 tangent space approach captures the highly nonlinear nature of constrained multibody dynamic systems and exhibits better error control with all the dynamic parameters after integration. This is because it includes constraints on position, velocity, as well as acceleration. Hence, the case studies were modeled using this approach.

Formulation of the tangent space index 0 DAE

Consider an orthonormal inertial reference frame $\mathbb{I}\triangleq \{O;X,Y,Z\}$ centered at the origin, $O$ with axes $X,Y,Z\in {\mathbb{R}}^3$ . The body reference frame is given by $\mathbb{B}\triangleq \{O^{\prime };X^{\prime },Y^{\prime },Z^{\prime }\}$ centered at $O^{\prime }$ with axes $X^{\prime },Y^{\prime },X^{\prime }$ . The generalized coordinates are given by $\mathbf{q}\triangleq \mathbf{[}\mathbf{r}\mathbf{p}{\mathbf{]}}^{\mathbf{T}}$ , where $\mathbf{r}\triangleq \mathbf{[}\mathbf{x}\mathbf{y}\mathbf{z}{\mathbf{]}}^{\mathbf{T}}$ and the Euler parameters, $\mathbf{p}\triangleq \mathbf{[}{\mathbf{e}}_{\mathbf{0}}{\mathbf{e}}_{\mathbf{1}}{\mathbf{e}}_{\mathbf{2}}{\mathbf{e}}_{\mathbf{3}}{\mathbf{]}}^{\mathbf{T}}$ . The constraint manifold is given by $\mathrm{\Phi }(\mathbf{q}\mathbf{,}\mathbf{t}\mathbf{)}\triangleq \mathbf{0}$ which captures the constraint mappings of the connected links in the multibody system.

For any generalized coordinates, $\mathbf{q}\in {\mathbb{R}}^{\mathbf{n}}$ , the equations of motion for multibody systems with $m$ holonomic constraints can be written in Langrange multiplier form as a DAE as follows^7–9

$$\mathbf{M}(\mathbf{q})\ddot{\mathbf{q}}+{\mathrm{\Phi }}_{\mathbf{q}}^T(\mathbf{q},t)\lambda -\mathbf{S}(\mathbf{q},\dot{\mathbf{q}},t)-{\mathbf{Q}}^A(\lambda ,\mathbf{q},\dot{\mathbf{q}},t)=0$$

(6)

The Jacobian of the constraint manifold, ${\mathrm{\Phi }}_{\mathbf{q}}\in {\mathbb{R}}^{\mathbf{m}\times \mathbf{n}}$ of the three times continuously differentiable function $\mathrm{\Phi }(\mathbf{q}\mathbf{,}\mathbf{t}\mathbf{)}$ must have full row rank. Equation (6), which is in the generalized coordinate configuration space, can be solved by mapping the equation to a tangent space configuration space as equation (7)

$$\mathbf{M}\left(\mathbf{q}\right)\ddot{\mathbf{v}}+{\mathrm{\Phi }}_{\mathbf{q}}^T(\mathbf{q},t)\lambda -\mathbf{M}\left(\mathbf{q}\right)\mathbf{U}\mathbf{B}\gamma -{\mathbf{Q}}^{\mathbf{A}}(\mu ,\lambda ,\mathbf{q},\dot{\mathbf{q}},t)=0$$

(7)

where

$$\gamma =-{\left({\mathrm{\Phi }}_{\mathbf{q}}\left(\mathbf{q}\right)\dot{\mathbf{q}}\right)}_{\mathbf{q}}\dot{\mathbf{q}}-{\mathrm{\Phi }}_{tt}$$

(8)

$$\mathbf{U}={\mathrm{\Phi }{\left({\mathbf{q}}_{\mathbf{0}}\right)}_{\mathbf{q}}}^T$$

(9)

$$\mathbf{B}=\mathrm{\Phi }{\left(\mathbf{q}\right)}_{\mathbf{q}}{\mathbf{U}}^{-1}$$

(10)

$$\mathbf{D}=\mathbf{I}-\mathbf{U}\mathbf{B}\mathrm{\Phi }{\left(\mathbf{q}\right)}_{\mathbf{q}}\mathbf{V}$$

(11)

Here $\mathbf{M}\mathbf{(}\mathbf{q}\mathbf{)}$ is the combined mass-inertia matrix of the system, $\mathbf{V}$ is the null space of the constraint Jacobian $\mathrm{\Phi }(\mathbf{q}{\mathbf{)}}_{\mathbf{q}}$ , and ${\mathbf{q}}_{\mathbf{0}}$ is the orientation of the system at any time step in ${\mathbb{R}}^3$ , where reparameterization happens, ⁷ and $\mathbf{U}$ and $\mathbf{V}$ together form the tangent space planes at point ${\mathbf{q}}_{\mathbf{0}}$ . Here $\mu$ represents the vector of static and dynamic coeffiecients of friction. In equation (5), it may be observed that the Lagrangian multiplier $\lambda$ is an input to the external force function ${\mathbf{Q}}^{\mathbf{A}}$ . This means that there is no explicit method to evaluate $\lambda$ . Hence, at every timestep, equation (5) is satisfied by Newton-Raphson iterations that require Jacobian matrix of the complete equation (5) with respect to acceleration of independent coordinates $\ddot{\mathbf{v}}$ in the tangent plane and $\lambda$ . Interested readers are advised to refer to Haug^7,9 for a detailed description of Index-0 tagent space DAE formulation. For this study, the above presented equations of motion were numerically integrated using explicit Runge-Kutta Fehlberg-Nystrom 45 integrator with adaptive time-step scheme. Throughout the integration, it was ensured that the norm of the constraint manifold was not violated

Multibody system case studies

In this section, two multibody dynamic systems that were formulated using index 0 tangent space DAEs are described. The constraint equations of relevant joints are presented along with the geometric and inertial parameters of both systems.

Spatial four mass system

As shown in Figure 5, the spatial four mass system consists of four point masses represented by $m_1,m_2,m_3$ , and $m_4$ . These masses form translational joints with the ground. The mass pairs $(m_1$ , $m_2)$ and $(m_2,m_3)$ are connected by fixed length massless rods, the masses $m_1,m_2$ , and $m_3$ are connected to the origin via springs of stiffness $k_1,k_2$ , and $k_3$ , respectively.

Figure 5.

Spatial four-mass system setup.

The masses $m_1$ and $m_4$ are connected by another spring of stiffness $k_4$ . The system comes into motion when mass $m_3$ starts falling under the action of gravity. Each point mass has seven generalized coordinates given by $\mathbf{q}\in {\mathbb{R}}^{\mathbf{7}}$ . Hence the configuration space of the system, $\mathit{\text{C}}\in {\mathbb{R}}^{28}$ . The system is constrained by four translational joints and two distance constraints. The distance constraint equations between two bodies with generalized coordinates ${\mathbf{q}}_{\mathbf{i}}$ and ${\mathbf{q}}_{\mathbf{j}}$ are shown in equations (11) and (12) ⁸

$${\mathrm{\Phi }}^{dist}({\mathbf{q}}_i,{\mathbf{q}}_j,d)=\frac{({\mathbf{d}}_{ij}^T{\mathbf{d}}_{ij}-d^2)}{2}=0$$

(12)

$${\mathbf{d}}_{ij}={\mathbf{r}}_j+{\mathbf{A}}_j{\mathbf{s}}_j^{\prime }-{\mathbf{r}}_i-{\mathbf{A}}_i{\mathbf{s}}_i^{\prime }$$

(13)

where for bodies $i$ and $j$ , respectively, ${\mathbf{d}}_{\mathbf{i}\mathbf{j}}\in {\mathbb{R}}^{\mathbf{3}}$ is the distance constraint, ${\mathbf{s}}_{\mathbf{i}}\in {\mathbb{R}}^{\mathbf{3}}$ and ${\mathbf{s}}_{\mathbf{j}}\in {\mathbb{R}}^{\mathbf{3}}$ are the vectors to the centroids in body fixed reference frame, ${\mathbf{A}}_{\mathbf{i}}\in {\mathbb{R}}^{\mathbf{3}\times \mathbf{3}}\mathbf{,}{\mathbf{A}}_{\mathbf{j}}\in {\mathbb{R}}^{\mathbf{3}\times \mathbf{3}}$ are the rotation matrices and ${\mathbf{r}}_{\mathbf{i}}$ , ${\mathbf{r}}_{\mathbf{j}}$ are the vectors to the centroids in the inertial reference frame. $d$ is the length of the connecting rod, ${\mathrm{\Phi }}^{dist}\in {\mathbb{R}}^1$ is the spatial distance constraint manifold. The translational joint constraints are given in equations (14) to (16) ⁸

$${\mathrm{\Phi }}^{trans}({\mathbf{q}}_i,{\mathbf{q}}_j)=\left[\begin{array}{c}{\mathrm{\Phi }}^{dot2}({\mathbf{u}}_j,{\mathbf{d}}_{ij})\\ {\mathrm{\Phi }}^{dot2}({\mathbf{v}}_j,{\mathbf{d}}_{ij})\\ {\mathrm{\Phi }}^{dot1}({\mathbf{w}}_j,{\mathbf{u}}_j)\\ {\mathrm{\Phi }}^{dot1}({\mathbf{w}}_j,{\mathbf{v}}_j)\\ {\mathrm{\Phi }}^{dot1}({\mathbf{v}}_i,{\mathbf{u}}_j)\end{array}\right]=0$$

(14)

$${\mathrm{\Phi }}^{dot1}({\mathbf{a}}_i,{\mathbf{a}}_j)={\mathbf{a}}_i^{\prime T}{\mathbf{A}}_i^T{\mathbf{A}}_j{\mathbf{a}}_j^{\prime }=0$$

(15)

$${\mathrm{\Phi }}^{dot2}({\mathbf{a}}_j,{\mathbf{d}}_{ij})={\mathbf{a}}_j^{\prime T}{\mathbf{A}}_j^{\prime T}{\mathbf{d}}_{ij}=0$$

(16)

where for bodies $i$ and $j$ , $\mathbf{u}\mathbf{,}\mathbf{v}\mathbf{,}\mathbf{w}$ represent unit vectors in body fixed reference frames and $\mathbf{a}$ is any non-zero body fixed vector orthogonal to ${\mathbf{d}}_{\mathbf{i}\mathbf{j}}$ . Hence, ${\mathrm{\Phi }}^{trans}\in {\mathbb{R}}^5$ , ${\mathrm{\Phi }}^{dot1}\in {\mathbb{R}}^1$ and ${\mathrm{\Phi }}^{dot2}\in {\mathbb{R}}^1$ .

Thus, the system has 22 constraint equations, which results in the system to have only two degrees of freedom. The constraint manifold of the system, $\mathrm{\Phi }(\mathbf{q}\mathbf{)}\mathbf{=}\mathbf{0}$ , therefore, is in ${\mathbb{R}}^{22}$ with $\mathbf{q}\in {\mathbb{R}}^{\mathbf{28}}$ . Hence the dimension of the constraint Jacobian, $\mathrm{\Phi }(\mathbf{q}{\mathbf{)}}_{\mathbf{q}}\in {\mathbb{R}}^{\mathbf{22}\times \mathbf{28}}$ . The geometric and inertial parameters of the four-mass spatial system is tabulated in Table 3. The coefficient of friction ${\mu }_s$ was taken to be 0.3 and ${\mu }_d$ was taken to be 0.1.

Table 3.

Geometric and inertial data for four-mass system.

Parameter	Value	Unit
$m_1$	2	kg
$m_2$	2	kg
$m_3$	6	kg
$m_4$	6	kg
$k_1$	2	${\mathrm{Nm}}^{-1}$
$k_2$	2	${\mathrm{Nm}}^{-1}$
$k_3$	10	${\mathrm{Nm}}^{-1}$
$k_4$	10	${\mathrm{Nm}}^{-1}$
${\ell }_1$	5	m
${\ell }_2$	7	m

Spatial slider-crank mechanism

As shown in Figure 6, the three-dimensional slider-crank mechanism is composed of a crank, massless connecting rod and slider. The mass and the inertia of the crank about its rotational axis are given by $m_c$ and $I_{rr}$ , respectively. The length of the connecting rod is given by $d$ and mass of slider is given by $m_s$ . The crank length is given by ${\ell }_c$ . The slider is connected to the ground through a translational joint. The crank rod is constrained to ground through a revolute joint. The crank and slider are constrained by a distance constraint. The crank is allowed to fall under the action of gravity, in the absence of any driving constraint.

Figure 6.

Spatial slider–crank mechanism.

The system has two bodies whose generalized coordinates are given by $\mathbf{q}\in {\mathbb{R}}^{\mathbf{7}}$ . Hence the configuration space of the system, $\mathit{\text{C}}\in {\mathbb{R}}^{14}$ . The system is constrained by one translational joint, one revolute joint, and one distance constraint. The equations for distance constraint and translational joint constraint are exlpained the previous subsection. The revolute joint constraint equations between any two bodies $i$ and $j$ having generalized coordinates, ${\mathbf{q}}_{\mathbf{i}}$ and ${\mathbf{q}}_{\mathbf{j}}$ are shown in equation (17) ⁸

$${\mathrm{\Phi }}^{rev}({\mathbf{q}}_i,{\mathbf{q}}_j)=\left[\begin{array}{c}{\mathrm{\Phi }}^{dot2}({\mathbf{u}}_j,{\mathbf{d}}_{ij})\\ {\mathrm{\Phi }}^{dot2}({\mathbf{v}}_j,{\mathbf{d}}_{ij})\\ {\mathrm{\Phi }}^{dot1}({\mathbf{w}}_j,{\mathbf{u}}_j)\\ {\mathrm{\Phi }}^{dot1}({\mathbf{w}}_j,{\mathbf{v}}_j)\\ {\mathrm{\Phi }}^{dot2}({\mathbf{w}}_j,{\mathbf{d}}_{ij})\end{array}\right]=0$$

(17)

where ${\mathrm{\Phi }}^{rev}({\mathbf{q}}_i,{\mathbf{q}}_j)\in {\mathbb{R}}^5$ is the revolute constraint manifold.

Thus, the system has 11 constraint equations, which results in the system to have only one degree of freedom. The constraint manifold of the system, $\mathrm{\Phi }(\mathbf{q}\mathbf{)}\mathbf{=}\mathbf{0}$ , therefore, is in ${\mathbb{R}}^{11}$ with $\mathbf{q}\in {\mathbb{R}}^{\mathbf{14}}$ . Hence the dimension of the constraint Jacobian, $\mathrm{\Phi }(\mathbf{q}{\mathbf{)}}_{\mathbf{q}}\in {\mathbb{R}}^{\mathbf{11}\times \mathbf{14}}$ . The geometric and inertial parameters of the spatial slider crank system are tabulated in Table 4. The coefficient of friction ${\mu }_s$ was taken to be 0.5 and ${\mu }_d$ was taken to be 0.2.

Table 4.

Geometric and inertial data for slider–crank mechanism.

Parameter	Value	Unit
$I_{rr}$	0.2	$\mathrm{kg}$ - ${\mathrm{m}}^2$
$m_c$	5	kg
$m_s$	10	kg
${\ell }_c$	0.08	m
$d$	0.23	m

The two case studies discussed in this section were simulated using numerical integration. The observations are discussed in detail in the next section.

Analysis and inferences

Anaslysis of the proposed friction model

The Coulomb’s friction model is a nonlinear function with discontinuities. In reality, the highly nonlinear nature of friction arises due to complex phenomena observed at a molecular level, ²¹ such as surface asperities, material properties, and temperature. However, modeling such phenomena for multibody dynamics and control applications is exceptionally challenging. The question then arises, how to use available information to model friction as close to accurate as possible while meeting requirements for dynamic simulations and controls? The work reviewed, as described in the Introduction section delineated approaches that attempted to remove the discontinuity, yet capturing the nonlinearity using independent quantities, such as slip velocity and coefficient of friction and introducing empirical or manually tunable parameters like transition velocity.

This leads to the next question, which is how to model friction with tunable parameters that do not have to be guessed but that can be obtained by using available information? The proposed model is able to circumvent the issue of having to guess such tunable or empirical parameters by using available information.

This is achieved by computing the tunable parameters using all available information, that is the magnitude of the normal reaction force and the surface properties, such as the static and the dynamic friction coefficients along with the slip velocity. The sine and inverse tangent function in equation (3) capture the nonlinearity while also maintaining $C^n$ continuity and differentiability. The stick-slip transition range is determined by evaluating the difference between static and dynamic coefficients of friction and the magnitude of normal reaction force, instead of manually estimating it. This is achieved by calculating the term $B_m$ as described in equation (4), and iteratively determining $C_{cor}$ as described in one of the previous sections. The following subsection describes the influence of the tunable parameters in equation (3) in detail.

Investigation into the proposed friction model parameters

To understand the effect of each parameter on the friction force model, simulations were performed by varying the query parameter while maintaining other parameters constant.

The terms are analogous to Pacejka’s tire model. ²² The parameter C is the “shape factor.” It can be seen in Figure 7 that when C is set to 0, the curve remains a flat horizontal line, indicating there is no friction force. Hence, there is no curvature to the profile, however, as the value of C is increased from 0 to 2, it is seen that curve resembles a negative hyperbolic tangent function until C = 1 and then begins to show a profile as one would expect in a friction force curve against slip velocity. Further, it is seen that, when C is increased beyond 1.4, the peaks become sharper. Since, C controls the shape of the friction force plot, this is a key parameter.

Figure 7.

Influence of $C$ on $F_f$ .

Next, we study the influence of parameter B. The values of B are changed from 0 to 1 as shown in Figure 8. It can be seen that, at B = 0, a flat horizontal line is obtained for the friction force. At B = 0.1 the slope changes but with minimum curvature. Beyond 0.3, the profile remains almost the same however, the negative gradient at zero slip velocity begins to increase. This influences the stiffness of the curve and hence it is termed as the “stiffness factor.”

Figure 8.

Influence of $B$ on $F_f$ .

The effect of parameter E on the friction force is presented in Figure 9. It can be seen that for values of E > 0.5, the friction force profile begins to show a sinusoidal behavior, but for values of E < 0.5, the friction force profile is what is as expected. This parameter thus controls the curvature of the friction curve and hence it is known as the “curvature factor.”

Figure 9.

Influence of $E$ on $F_f$ .

From the analyses of the three contributing parameters, it can be seen that the choice of C has the strongest influence on the magnitude of the friction. As a consequence, the closer $C$ is to its “optimal” value, the more accurate the magnitude of the friction force will be, specifically in the dynamic range. The value of B on the other hand, does not affect the friction force magnitude significantly, however, the gradient of the curve at zero slip velocity changes, so $B$ is ranked as the second significant parameter. Considering this, the model in equation (2) is modified as presented in equations (3) and (4) for the proposed friction model.

An investigation into the effect of $B_m$

The parameter $B_m$ affects the stiffness of the system. It was observed that as $B_m$ was increased from $0$ to $1$ , the gradient became steeper. The limit values of $B_m$ , $(50-100)$ was obtained after a separate analysis on the sensitivity of the friction force on $B_m$ alone. The sensitivity of $F_f$ with respect to $B_m$ is shown in equation (18)

$$\frac{d{F}_f}{d{B}_m}=AB$$

(18)

where

$$\begin{array}{rl}A& =\left(\frac{C{F}_n{\mu }_s\cos (C{\tan }^{-1}\left(2{\tan }^{-1}\left(B_m{v}_s\right)\right)-3B_m{v}_s)}{{(2{\tan }^{-1}\left(B_m{v}_s\right)-3B_m{v}_s)}^2+1}\right)\\ B& =(3B_m-\frac{2B_m}{B_m^2{v}_s^2+1})\end{array}$$

From the sensitivity plot in Figure 10, it can be seen that there is a dramatic change in the gradient when $B_m$ is below 50, however, beyond 100, the standard deviation of the sensitivity is below 2%. Hence, the maximum possible value for $B_m$ is set to 100. Further, it was seen that, by increasing the value of $B_m$ over 100, the friction force vs. slip velocity gradient is extremely steep, which renders the problem stiff. Hence, with values of $B_m$ larger than 100, the stiffness of the problem made the integrator to not converge. On the other hand, when the lower limit for $B_m$ was chosen to be below 50, the gradient was smaller in the friction force–slip velocity curve, and the model was not able to accurately capture the friction phenomena.

Figure 10.

Sensitivity of friction force to $B_m$ .

Alongside the information we gather from bristle deflection-based models, ¹³ the contact surfaces A and B sliding with relative slip velocity $\dot{\theta }$ , can be visualized as bristles as shown in Figure 11. A higher $B_m$ corresponds to a higher shear resistance, hence the bristles between the body and the surface are expected to show a wider transition range. On the other hand, lower $B_m$ corresponds to lower shear resistance, which indicates a faster transition from static to the dynamic regime. It can also be inferred that a higher $B_m$ corresponds to tighter interlocking between the bristles in contact which manifests as higher shear resistance. This means, that the term $B_m$ should be directly proportional to the magnitude of normal reaction force.

Figure 11.

Bristle deflection model.

Secondly, the difference between the static and dynamic coefficients of friction reflects the structural property of individual bristles. A large difference in the coefficients of friction assumes a brittle structure of the bristle. This means, that a large difference in friction coefficients makes the transition faster from static to dynamic regimes as brittle bristles are more prone to fail in shear. However, a small difference in coefficients of friction corresponds to a more ductile structure of bristles and it allows for a slower transition from static to the dynamic regime as bristles go through elastic bending by displacement $z$ . Therefore, the stiffness factor $B_m$ can be expected to be directly proportional to the difference between the static and the dynamic coefficients of friction.

Based on the aforementioned inferences, it was deemed necessary to develop a function that evaluates $B_m$ as a function of the contact force and the difference between static and dynamic friction coefficients. The model proposed is given in equation (4). Figure 12 further visualizes the above-explained variation of $B_m$ with respect to change in the magnitude of normal reaction force. It can be seen that the equation can be used to calculate an optimal choice for $B_m$ rather than having to guess an optimal value; however, the methodology developed to determine $B_m$ can be further improved.

Figure 12.

Variation of $B_m$ with magnitude of normal force.

Case study simulation results

Spatial four-mass system

This section illustrates the simulation results of the spatial four-mass system case study. For the entirety of the simulation, the joint constraints were not violated and is quantified by the norm of the constraint, $\Vert \mathrm{\Phi }(\mathbf{q}){\Vert }_2=3.05\times {10}^{-13}$ which was below the set tolerance of ${10}^{-6}$ . For the Brown and McPhee friction model, the recommended value of transition velocity which is, $10-20$  h ¹⁰ gave a jittery behavior, as shown in Figure 13, in the acceleration plot after $1.5$ s, even though the system comes to a halt, as can be inferred from Figure 14.

Figure 13.

Acceleration versus time: Brown and McPhee model, $v_t=20h$ .

Figure 14.

Position versus time: four-mass system.

To fathom further into the sensitivity of the transitional velocity in the Brown and McPhee model, the value for $v_t$ was taken to be 50 h and then 100 h which is beyond the recommended range. Figures 14 to 16, respectively, represent the position, velocity, and acceleration curves of the system, respectively.

Figure 15.

Velocity versus time: four-mass system.

Figure 16.

Acceleration versus time: four-mass system.

At $v_t=50$  h, the acceleration plots still show jittery behavior which reduces when $v_t$ was increased to 100 h. A large $v_t$ value means lower sensitivity to change in sliding velocity. This shows that the selection of appropriate $v_t$ value is not only crucial for obtaining smooth solutions, but also dependent on the system’s characteristics and is subject to change with every mechanism. This case study thus re-emphasizes the need for an algorithm to choose a model dependent empirical parameter which significantly influences the dynamic simulations. On the other hand, the jittery behavior was not observed, when using the proposed model, as can be observed in the position, velocity, acceleration plots. By utilizing an iterative regulation of shape and stiffness factors (B and C) which adapt for slight changes in normal reaction and relative velocities, respectively, the simulations produce smooth profiles in acceleration, velocity, and position.

As can be observed from the acceleration plots, the proposed model inherently possesses the capacity to tune the key parameters automatically by taking into account all the information about the forces acting on the body. Because of this inherent characteristic, it provides smooth signals irrespective of the choice of input parameters. This ability could be useful for implementing control systems into the systems with dry friction. The modern state-space control systems employ a number of mechanisms to generate smooth signals, such as anti-aliasing filters and low-pass filters. ²³ Especially, in the multi-input and multi-output (MIMO) control systems applications, the number of signal processing devices to employ is directly correlated to the number of signals. It could be a cost effective measure in such cases if the dynamic signals can be smoothened out naturally.

Spatial slider–crank mechanism

For crank and slider, respectively, Figures 17 and 18 show the accelerations, and Figures 19 and 20 show the velocities. Figure 21 shows the position of slider against time.

Figure 17.

Crank acceleration against time.

Figure 18.

Slider acceleration against time.

Figure 19.

Slider velocity against time.

Figure 20.

Crank velocity against time.

Figure 21.

Slider position against time.

As in the previous case study, the joint constraints were not violated and is quantified by the norm of the constraint, $\Vert \mathrm{\Phi }(\mathbf{q}){\Vert }_2=2.86\times {10}^{-12}$ which was below the set tolerance of ${10}^{-6}$ . When using Brown and McPhee’s friction model, $v_t$ was kept within the recommended range $15$  h, ¹⁰ and the solutions with both friction models were observed to agree with each other.

Unlike the spatial four-mass system, no jittery behavior was observed in the acceleration plot. The acceleration and velocity profiles are smooth, and it can be seen that they are very similar, overlapping almost completely. It is, however, observed that after $0.5$  s into the simulation, the position of the slider crank has a slight offset. This can be attributed to the truncation error that builds up during the numerical integration of acceleration and velocity.

Conclusions

To summarize, a novel friction model which possesses continuity and differentiability. The proposed friction model has been studied critically for the effect of each of its parameters. The proposed model brings novelty in the terms that, unlike its contemporary friction models, the proposed model does not depend on any empirical parameter such as transition velocity. This makes this model very useful for simulating large-scale multibody systems while considering all the major aspects of friction phenomena.

Further, it has the inherent property to tune its parameters such that the produced signals are smooth in nature. This characteristic of the proposed model could prove vital for control system applications where smooth signals are preferred to save costs on anti-aliasing and low-pass filters.

The theoretical and conceptual reasoning for the criteria and algorithms to select the model-dependent parameters have been described in detail. The proposed friction model was simulated for the Benchmark Rabinowicz’s stick-slip experiment, which validated the potential for the proposed friction model to be used in the simulation of dynamical systems.

The performance of the proposed model was also compared to the state-of-the-art Brown and McPhee’s friction model and it was shown that both models are in agreement. The applicability of the proposed friction model was then demonstrated by simulating two multibody dynamic systems.

The index-0 tangent space DAE formulation method was adopted for developing equations of motion of the spatial four-mass system and of the spatial slider–crank system which incorporated into the friction force models. It was concluded that the choice of the transitional velocity $v_t$ can impact the dynamic simulations, lower values of transitional velocity induced jitter in the acceleration profile for the spatial four-mass case study.

Acceleration, could either be an output from a sensor or just a state-space variable in the modern control systems. To deal with multiple jittery signals, the modern control systems often use low-pass filters and anti-aliasing devices. In such control systems where every signal needs to be filtered with a device of appropriate band-width and Nyquist frequency, a lot of cost could be saved by having jitter-free signals.

From the simulation results incorporating the proposed friction model, it was seen that similar results can be indeed achieved, while no jittery behavior can be seen. This was accomplished by eliminating the need for manually tunable parameters and instead using parameters that use available information from the system.

The algorithm has been explained as to how it inherently accounts for the shape and stiffness factors thus adapting to slight changes in normal force and slip velocity. Hence, the proposed friction model can be used for multibody dynamic system applications.

The future scope of this work includes implementing the model into designing control systems and testing this model against one that includes bristle deflection theory, such as the LuGre and Gonthier’s model.

Further improvements that can be carried out on the proposed model is to include the viscous friction regime.