Modeling of static friction in closed-loop kinematic chains—Uniqueness and parametric sensitivity problems

Problems with uniqueness and high parametric sensitivity of the solution of equations of motion, encountered in the static friction regime, are addressed. Friction in joints of a multibody system with closed-loop kinematic chains is discussed. Three different models of friction are studied: the discontinuous Coulomb model with stiction regime represented in terms of additional constraints; the approximate Coulomb model, smoothed in the vicinity of zero relative velocity; and the LuGre model with presliding displacements represented in terms of auxiliary state variables. Firstly, a rigid body model is investigated. It is shown that in the case of constraint addition approach, problems with uniqueness of solution emerge in the static friction regime. In the case of continuous models of friction, the solution in the stiction regime and its vicinity is highly sensitive to some hardly measurable or arbitrarily chosen parameters of the model of friction. Origins of nonuniqueness and high sensitivity are investigated, and the questionable credibility of the stiction regime simulation results is discussed. Secondly, a simplified model of body and joint elasticity is introduced to investigate the impact of flexibility on the mechanism frictional behavior. It is shown that taking the flexibility into consideration may eliminate the uniqueness and sensitivity problems. Moreover, the quantities that represent flexibility may be regarded as the key factors influencing the results of stiction regime simulation. Five examples are provided to illustrate the presented considerations.


Introduction
The rigid body assumption is frequently adopted in multibody system modeling and simulation [1][2][3]. Neglecting body elasticity has many advantages, especially in terms of model complexity and its computational efficiency. There are, however, known limitations of rigid B M. Wojtyra mwojtyra@meil.pw.edu.pl 1 Institute of Aeronautics and Applied Mechanics, Warsaw University of Technology, Nowowiejska 24, 00-665 Warsaw, Poland body approach, associated with redundant constraints existence. They are particularly important when frictional effects are modeled. It has been shown that when a multibody system is overconstrained and Coulomb-like joint friction is present (i.e., friction that depends on normal reactions in joints), the simulated motion of the system may be not unique [4]. As a result, combining joint friction with redundant constraints usually makes the rigid body model useless (the simulated motion is unique only in some special circumstances).
In this study, it is shown that problems with uniqueness of solution may be experienced when closed-loop mechanisms without redundant constraints are modeled and simulated. The investigation is focused on the stiction phase and its vicinity when the relative joint velocities are null or close to zero. The details of problematic behavior depend on methods of friction modeling; however, the results are similar in all cases: neither the stiction forces nor the moment of transition between static and kinetic friction regime can be credibly predicted when the rigid body model is employed.
Three different models of joint friction are investigated. All these models originate from the classical Amontons and Coulomb laws [5] and assume the existence of two friction regimes, static friction (stiction) and kinetic friction (sliding); both static and kinetic friction forces depend on normal reactions in joints. All friction models discussed in this article may be used in multibody dynamics simulations.
The first investigated here technique of modeling joint friction in a multibody system, the method of constraints addition and deletion, was proposed by Haug et al. [6] and was also discussed by Garcia de Jalon and Bayo in their textbook [3]. Due to some numerical difficulties, associated with the necessity of constraints switching and discontinuity of friction laws at zero relative joint velocity, this approach is seldom used in practical calculations. However, examination of this model is beneficial since some issues crucial for this paper are clearly visible. It is shown that in the stiction regime the multibody system becomes overconstrained and calculation of joint reactions is obstructed. As a result, problems with solution uniqueness are encountered.
The second investigated model of joint friction is representative for a group of models with continuous friction force laws. The problems with friction discontinuity at zero relative velocity are often alleviated by smoothing the Coulomb friction law. Models of this kind are frequently used in both theoretical considerations and practical applications; see, for example, [7][8][9][10]. Some drawbacks of smoothed Coulomb models are well known: the set of differential equations of motion is usually stiff, and a drift is observed instead of stiction since the friction force equals zero at zero relative joint velocity. In this paper, it is shown that at the vicinity of zero relative velocity, the multibody model behavior strongly depends on some friction model parameters of secondary importance. The third examined model of frictional effects in joints is the LuGre model [11]. It is representative for a group of models with auxiliary "internal dynamics" state variable that describes the micro-slip during the stiction phase [12]. Apart from the LuGre model, this approach is adopted, for example, in Dankowicz [8,13], Dahl [14], and Leuven [15,16] friction models. In these models, both stiction and sliding regimes are unified and governed by the same equations. At a price of introducing additional differential equations, some drawbacks of smoothed Coulomb models are overcome, and more realistic results are gained [8,17,18]. This article points out that in the stiction regime and its neighborhood, the multibody model behavior is highly sensitive to changes of some less important parameters (often hardly measurable) of the model of joint friction.
It is known that in the case of overconstrained mechanisms, flexibility of bodies is one of the main factors that decide how loads are shared between joints [19,20]. The problems associated with determining friction forces in the stiction regime are closely related to constraints redundancy, and therefore flexibility of mechanism bodies is taken into account in the investigations. The rigid body assumption is released, and the same three models of joint friction are examined once more. It is demonstrated that the secondary parameters of friction models no longer play the essential role in calculation of static friction forces. Instead, the stiffness of mechanism links is crucial for the results, which is observed in the stiction regime for all tested models of friction.
This paper is organized as follows. Firstly, the basics of modeling of joint friction in multibody systems are briefly recalled. Then, the problems emerging when modeling the stiction regime are thoroughly discussed. Next, a simplified flexible body model is investigated in order to analyze the dependencies between joint friction and elasticity. Finally, conclusions summarize the contribution. Every stage of the discussion is accompanied by a relevant example.
2 Modeling friction in joints of a multibody system

Equations of motion for systems with friction
Consider a rigid multibody system described by an n-element vector of dependent coordinates q that are subject to m holonomic constraints: It is assumed that constraints Φ are independent, that is, no redundant constraints exist. Note that in the case of a multibody system described by absolute Cartesian coordinates, each joint (kinematic pair) is represented by a subset of Eq. (1), whereas in the case of relative joint coordinates, individual equations comprising Eq. (1) are formulated only for selected, loop-closing joints.
Differentiation of Eq. (1) with respect to time leads to constraints for velocities and accelerations: where Φ q is the constraint Jacobian matrix, and Γ = −(Φ qq ) qq . For the sake of simplicity, let us make two assumptions. First, we assume that the multibody system is described by absolute coordinates [1,2]. As a consequence, each kinematic joint is described by separate constraint equations and all joint reaction forces (i.e., constraint reactions) can be directly obtained. The second assumption is that the considered system has exactly one degree of freedom (m = n − 1) and forms a closed-loop kinematic chain with one or more loops. This guarantees that the relative velocities in all kinematic pairs vanish simultaneously.
When constraints are ideal, that is, the joints are frictionless, the multibody system equations of motion can be derived from Lagrange equations of the first kind [3] or from Newton-Euler equations [1,2] and written in the following form: where M denotes the mass matrix, and Q contains the external forces and all velocitydependent inertial terms. The generalized constraint reactions are represented by the product of the transposed constraint Jacobian Φ T q and the vector of Lagrange multipliers λ. The transposed Jacobian matrix represents the directions of normal reaction forces, whereas the Lagrange multipliers represent the magnitudes of these forces.
Note that Eq. (4) is formulated for the whole system of bodies. Recursive formulations of equations of motion for systems described by absolute coordinates (see, e.g., [21,22]) are computationally more efficient; however, they are not utilized here since closed-form equations are more convenient for discussion.
The mathematical model is more complicated when joint friction is taken into account. Let us first discuss the kinetic friction regime. The generalized forces of sliding friction Q F can be treated as additional external forces that depend on normal reactions [3,4,6], and thus Eq. (4) will be rewritten as The model of friction in joint, used to calculate appropriate elements of vector Q F , establishes the relation between the normal and tangential force. Apart from the frictional properties, the details of calculations strongly depend on the joint type, instantaneous configuration, and geometry (features of design, dimensions, shapes of contacting surfaces, etc.). The details of vector Q F calculations and some examples relevant to various types of joints can be found in Refs. [6] and [4]. Briefly, calculation of the vector Q F (q,q, λ) may be summarized as follows: • The generalized normal reactions magnitudes λ and directions (Φ q (q)) T are used to calculate physical normal reactions N i in each joint i; calculations depend on the pair i type and geometry and often on its instantaneous configuration, and thus N i = N i (q, λ). • The model of friction (e.g., the Coulomb or LuGre model) is utilized to calculate the physical friction force in each joint i : F i = F i (q,q, N i ). Most often, the constructional details of the kinematic pair are taken into account during this calculation. • The physical friction forces acting in joints are transformed into generalized friction forces: There are several approaches to simulate multibody systems described by Eqs. (1)-(3) and Eq. (5). One of the simplest, the so-called index-1 formulation, consists in appending Eqs. (5) and (3) to obtain a system of differential-algebraic equations, and in providing initial conditions consistent with Eqs. (1) and (2). Various numerical methods can be employed to integrate equations of motion. At each integration step, Eq. (6) is solved for accelerationsq and Lagrange multipliers λ. Then, accelerationsq are integrated to obtain velocitiesq and coordinates q. Note that usually some stabilization techniques (see, e.g., Refs. [23][24][25]) are adopted to diminish constraint violation at velocity and position levels.
The generalized friction force Q F is not linear in Lagrange multipliers λ, and thus some iterative methods are utilized to solve Eq. (6) forq and λ (most often, fixed point iterations are applied [3,6]). It is usually assumed that solution exists and is unique. It can be shown, however, that this assumption is false for some systems, especially when friction is high [18,26]. Nevertheless, to focus on problems encountered in the stiction regime, in this paper, it is assumed that exactly one solution to Eq. (6) exists.
In the kinetic friction regime, the differences between various friction models, employed to calculate the generalized friction force Q F , are manifested mainly during calculations of physical friction forces in joints (F 1 , . . . , F k ). In the static friction regime and its vicinity (when relative velocities at joints are close to zero), the differences between friction models are much deeper since in some cases, the structure of equations of motion is affected. Therefore, starting from this point, each friction model considered in this paper must be discussed separately.

Coulomb model with constraints addition and deletion
According to the Coulomb model, the magnitude of kinetic friction force between two contacting surfaces can be calculated as where μ K i is the kinetic friction coefficient, N i is the physical normal force (related to the generalized normal reactions (Φ q (q)) T λ), and F i is the physical friction force. The direction of the friction force vector is such that it opposes the direction of the relative velocity of sliding.
Similar equations may be formulated for a whole joint, with multipoint (or multisurface) contact. For example, in the case of a translational joint, the form of equations remains unchanged, whereas in the case of a revolute joint, the friction torque rather than the force is calculated. In the kinetic friction regime, physical joint friction forces (or torques) F i are used to calculate the generalized friction force Q F that is utilized in Eq. (6).
Equations (6) and (7) are valid as long as the multibody system is in motion. At a certain moment, the relative velocity in a joint can reach the zero value. After that event, motion of the system may continue instantaneously or, when stiction occurs, can be stopped for some nonzero period of time. According to the Coulomb model, no relative joint motion is observed (the joint stays locked) as long as the joint friction force remains within an allowable range. This condition for ith joint can be written as follows [6]: where μ S i is a coefficient of static friction. In multibody modeling, the joint locking due to static friction can be represented by an additional constraint imposed on the system [6] (the constraint is activated when the joint relative velocity is numerically close to zero). This constraint may be regarded as a driving constraint that imposes no relative motion in joint. The reaction associated with the additional driving constraint corresponds to the joint static friction force and can be utilized to calculate the physical joint friction force F i . Condition (8) is repeatedly checked during simulation to determine whether stiction occurs. If condition (8) is not satisfied, then the additional constraint is deactivated, and the joint transits from stiction to sliding.
In the case of a 1-DOF mechanism with one or more kinematic loops, the relative velocities in all joints become equal to zero at the same time. Thus, additional constraints must be simultaneously added to all kinematic pairs with friction. Note that supplementary constraints Ψ consist of k scalar equations, where k is the sum of relative degrees of freedom in all locked joints. During the period of stiction, not only the additional constraints are introduced, but also kinetic friction forces Q F and accelerationsq vanish. As a result, Eq. (6) must be replaced by the following modification of Eq. (5): where q * represents the configuration at which mechanism is stopped (q * = 0), Ψ q denotes the Jacobian matrix of additional constraints, κ is the vector of Lagrange multipliers that represents the reactions of additional constraints, and subscripts denote the sizes of the matrices and vectors.

Approximate Coulomb model
In practical applications, the constraint addition-deletion approach is seldom used since discontinuity at zero relative joint velocity causes a number of analytical and computational problems [17,18]. More frequently employed models treat both stiction and sliding regimes uniformly, that is, the same friction model is used in both regimes, and no additional constraints are imposed during the stiction phase. In these models, the friction force at zero relative velocity is described by equality rather than inequality. The problems with friction discontinuity at zero relative velocity are often alleviated by smoothing the Coulomb friction law, which is usually accomplished by assuming that the coefficient of friction is a highly nonlinear function of the relative joint velocity. As a consequence, Eq. (7) takes the form where v i is the velocity of sliding at ith joint, and v i depends on generalized coordinates and velocities: v i = v i (q,q).
In practical applications, to approximate the Coulomb model (Fig. 1a), various functions of the velocity of sliding are used. Some of them, for example, saturation [8,17], hyperbolic tangent [7,18], or arctangent [27] do not differ static and kinetic friction coefficients (Fig. 1b). In some other, the gap between static and kinetic friction coefficients is reflected-see, for example, the piecewise linear function [9,28] or its smoothed version [9,29] in Fig. 1c.
In this study, the following piecewise linear function [28] is utilized (see Fig. 1c): where γ i is an arbitrarily chosen (usually very small) characteristic velocity at which the transition from stiction to sliding mode initiates, and the other symbols have been explained earlier.
All variants of approximate Coulomb model suffer from predicting zero friction force at zero relative joint velocity, which results in a drift observed in the static friction regime, especially when the external load is large. Moreover, particularly when the approximating curve is steep in the region of zero velocity of sliding, numerical properties of these models are poor [17,18]. These drawbacks may be diminished by introducing additional state variables, as it is done in the LuGre model of joint friction and other state-variable friction models [12,17].

LuGre model
At the microscopic level, contacting surfaces are irregular and make contact at a number of asperities. These asperities may be visualized as bristles that deflect like springs, when tangential force is applied; this gives rise to the friction force [11]. If the deflection is sufficiently large, then some of the bristles will slip. This process is random, but its results can be treated in an aggregate manner. In the LuGre model of friction [11], an additional state variable z i is introduced; this variable may be interpreted as the average deflection of asperities at contacting surfaces. The average deflection z i is modeled in terms of the first-order differential equationż with function G i that allows the model to accommodate the differing values of static and kinetic coefficients of friction, defined as where σ 0 i is a constant coefficient that represents the aggregate bristle stiffness, and ϑ i is a very small velocity that may be regarded as the border between stiction and sliding. When the LuGre model is applied, Eq. (11) is used to calculate the joint friction force; however, this time the instantaneous friction coefficient μ i is calculated as where the constant σ 1 i is a damping coefficient. Note that in the final mathematical model of the multibody system, equations of motion, Eq. (6), are accompanied by Eqs. (13), formulated for all joints with the LuGre friction. Clearly, the auxiliary state variables z = [z 1 , . . . , z k ] T are utilized to calculate the generalized friction force: Q F = Q F (q,q, z,ż, λ).

Problem statement
Due to its nonlinear nature, friction is always cumbersome in modeling. Encountered difficulties most often pertain to numerical stiffness of equations of motion [12,17]. It can be shown, however, that when friction is present in joints of a closed-loop rigid body mechanism, problems of completely different type, namely problems with uniqueness of solution, may emerge. This type of problems is observed in the stiction regime and its neighborhood.
In the case of friction models without discontinuity at zero slip velocity, these problems (at least formally) do not occur. However, they evolve into robustness problems. High sensitivity of calculated friction forces to some secondary parameters of the model of friction (i.e., parameters difficult to be measured and often arbitrarily chosen) is observed.

Constraints addition-deletion approach
When the Coulomb friction model with constraint addition-deletion is employed, a multibody system in the stiction regime is described by Eq. (10), that is, by n algebraic equations in m + k unknowns (λ and κ, respectively). Recall that in the case of a 1-DOF rigid body mechanism, n = m + 1, where m is the number of permanently imposed constraints (i.e., constraints that represent the kinematic structure of the mechanism, Eq. (1)), and n is the number of coordinates. Hence, when the number of additional constraints k (i.e., the constraints that represent joint locking due to friction, Eq. (9)) is greater than 1, the multibody system becomes redundantly constrained. Mathematically, the system of linear equations (Eq. (10)) is underdetermined. Thus, the Lagrange multipliers cannot be uniquely calculated using Eq. (10). It should be emphasized that, most often, the indeterminacy pertains not only to κ (friction) but also to λ (normal reactions). As a consequence, neither normal reactions nor static friction forces can be uniquely determined in the regime of static friction. Methods described in Refs. [30][31][32] may be employed to obtain precise information on constraint reactions solvability during the stiction occurrence. Note that, since constraints of Eq. (1) are assumed to be independent, when no stiction occurs, normal reactions and kinetic friction forces can be uniquely calculated.
It should be emphasized that when stiction occurs, the motion of the 1-DOF system is unique. Nevertheless, calculation of normal and friction forces is necessary since the conditions for stiction occurrence, Eq. (8), must be repeatedly checked as the simulation proceeds. The central problem is that, in many cases, it is impossible to decide whether or not stiction conditions are fulfilled for given external loads. Two substantially different types of solutions for normal and friction forces-one corresponding to sliding and the other corresponding to stiction-can be found. The rigid body model, combined with constraints addition, offers no means to privilege one solution and discriminate the other, which makes this modeling approach inapplicable. This situation is illustrated in the example that follows. Fig. 2a consists of three moving bodies, two sliders and a connecting rod. Two straight-line sliding rails belong to an unmovable basis of the mechanism. No gravity is considered. At a certain time instant, three external forces are applied to the mechanism: the force S 1 = [−5, 0] T N, perpendicular to the line of sliding, acts on slider 1, the force S 2 = [0, −20] T N, perpendicular to the line of sliding, acts on Absolute coordinates (n = 9) were used to describe the mechanism. Constraint equations that represent two revolute and two translational joints were formulated (m = 8). Next, two additional constraint equations were imposed in order to represent locking of the two translational joints due to friction forces (k = 2). Finally, Eq. (10) was formulated. Since for the discussed mechanism the number of unknown constraint reactions (m + k = 10) exceeds the number of equations (n = 9), the solution is not unique, and infinitely many solutions may be found. The details of calculations are omitted.

Example 1 The mechanism presented in
Let us focus on two solutions, graphically presented in Fig. 2b in the form of a freebody diagram (R ij is the joint reaction force that body i exerts on body j , N k is the normal reaction in joint k, and F k is the static friction force in joint k). The first reaction solution is the following (all forces expressed in newtons): The second solution is as follows: It can be easily verified that the vector sum of forces acting on any of the bodies is zero, for example, F 1 + N 1 + R 31 + S 1 = 0 (slider 1, both solutions). As long as the rigid body assumption holds, both solutions are equally acceptable. There is no justification for rejecting one of them and preserving the other.
In the case of the first solution, condition for stiction occurrence, established in Eq. (8), is not fulfilled for both joints with friction: which suggests (or even requires) that the additional constraints should be released and the friction model should switch to kinetic regime. On the other hand, in the case of the second solution, condition for stiction occurrence is fulfilled for both joints with friction: These results mean that the model of friction should remain in the stiction regime. Apparently, the multibody model gives no clear answer whether the mechanism moves or remains blocked by static friction. This observation concludes the example.
Two more issues concerning the static friction forces should be addressed. Firstly, the directions of friction forces must be coherent, that is, friction in all joints should oppose the upcoming motion (in the case of a 1-DOF mechanism, it is easy to verify this condition). Thus, some of solutions of Eq. (10) must be rejected due to incoherence. Nevertheless, the number of coherent solutions may still be infinite (note that the friction forces analyzed in Example 1 were coherent).
Secondly, the following requirement may be formulated: at the moment of transition from stiction to sliding, friction forces in all joints must reach their maximum allowable values. To justify this requirement, let us consider a joint transiting from stiction to sliding with static friction force less than allowable maximum. If this phenomenon occurs, then the Coulomb static friction law is violated. Thus, the requirement is legitimate as long as the Coulomb friction model is valid. A direct, but not intuitive, consequence is that-during the stiction phase when external forces applied to the mechanism raise-friction forces at a specified joint can increase so that maximum static friction is reached, but cannot go beyond this limit as long as friction forces at all other joints do not reach their respective maxima.
To fulfill the described requirement, the procedure of multibody system simulation must be amended. If, when solving Eq. (10), it is found that the friction in a specific joint is greater than allowable maximum (given by Eq. (8)), the respective constraint representing joint locking should no longer be imposed. Instead, the maximum static friction force (i.e., the force of magnitude μ S i N i (q * , λ) ) should be applied in the form of a generalized force Q F (q * , λ). It is a serious computational complication since it changes the structure of Eq. (10) into the following: where Ψ • represents the modified set of supplementary constraints (without equations for joints with maximum static friction), and κ • is the modified vector of Lagrange multipliers. Note that, unlike Eq. (10), the modified equation is no longer linear in λ.

Handling of redundant constraints in the constraint addition-deletion approach
In the scientific literature, it is difficult to find examples of applications of constraints addition-deletion method to closed-loop mechanisms. In Ref. [3], an open-loop kinematic chain is considered. The authors of Ref. [6] discuss a closed-loop mechanism; however, friction is introduced in just one joint, and thus problems with redundant constraints do not show up. Quite surprisingly, the author of this paper was unable to find a publication addressing, Fig. 3 External loads applied to the mechanism or at least acknowledging, the problem of constraints redundancy caused by introduction of supplementary constraints that represent joint stiction.
There are various methods of handling redundant constraints in multibody modeling. These methods are usually applied to systems that are permanently overconstrained, that is, not only in the stiction phase (see, e.g., [33]). It is known that the calculated joint reactions depend on the method chosen to cope with redundant constraints and, in the case of a rigid body model, no unique reaction solution exists [34]. When joints are frictionless, the simulated motion of the system is unique, despite the nonuniqueness of joint reactions. However, if a joint friction is present in an overconstrained mechanism and if kinetic friction forces depend on normal reactions, then the simulated motion may be not unique [4]. Thus, combining joint friction with redundant constraints usually makes the rigid body model useless (the simulated motion is unique only in some special circumstances).
As it was pointed out, similar problems with redundant constraints causing the nonuniqueness of reaction solution are observed when the constraint addition-deletion method is utilized to analyze the stiction regime. To conduct simulation, diverse methods of handling redundant constraints may be applied. Unfortunately, since the method of redundancy treatment affects the computed constraint reactions, the calculated friction forces depend on the method used. The indeterminacy of Eq. (10) is an inherent feature of the constraint addition-deletion method applied to a rigid body closed-loop kinematic chain, and thus it cannot be resolved by any purely mathematical procedure (i.e., procedure that does not represent physics of the analyzed mechanical system). The next example illustrates the outcome of application of diverse methods of redundant constraints handling. The nonuniqueness of friction and reaction forces and the nonuniqueness of simulated motion is demonstrated.
Example 2 Consider again the mechanism presented in Fig. 2a. The previously given information must be supplemented by the following: characteristic dimension is h = 0.05 m, the coefficients of kinetic friction are μ K 1 = μ K 2 = 0.3, the masses of sliders 1 and 2 are m 1 = m 2 = 20 kg, the mass of the connecting rod is m 3 = 10 kg, and its moment of inertia is J 3 = 2 kg m 2 .
The external forces S 1 , S 2 , and P are piecewise linear functions of time with constant directions and time courses of components presented in Fig. 3. At the beginning of the analyzed period of time, the mechanism is not moving. Note that at time t = 0 s some nonzero loads (S 1 and S 2 ) are applied to sliders 1 and 2 in order to have nonzero normal forces in the translational joints. Then, the force P, applied to the connecting rod and acting along it, rises from zero to reach its final value at t = 1 s. Next, up to t = 3 s, the forces S 1 and S 2 are being changed. Note that at t = 3 s, and after this time instant, the external loads are the loads investigated in Example 1. Two essentially different methods of redundant constraints handling were used during simulation of the mechanism, the constraint elimination method and the method based on pseudoinverse [34]. The results of simulations are presented in Fig. 4 as the calculated friction forces and in terms of the magnitudes of friction forces divided by the magnitudes of respective normal forces.
The constraint elimination method consists in excluding dependent equations from the mathematical model of the multibody system. In the case of the exemplary mechanism, the constraints Φ and Ψ (Eqs. (1) and (9), respectively) constitute a system of 10 algebraic equations. There are several possibilities to eliminate a dependent equation (a dependent equation can never be uniquely indicated). One of the "stiction" equations (Ψ ) could be eliminated; however, the outcome would be too obvious to discuss it. Therefore, some of equations that represent revolute joints were eliminated (when absolute Cartesian coordinates are utilized, each planar revolute joint is represented by two scalar constraint equations). Two different cases were studied.
In the first case, one of the two scalar equations representing the revolute joint on slider 1 was excluded from the mathematical model of the mechanism. As a result, Eq. (10) became a system of nine equations in nine unknowns. During the simulation, the condition for stiction occurrence (Eq. (8)) had been repeatedly checked for both translational joints. At the beginning of the simulation, due to elimination of the redundant constraint, no tangential loads were transmitted to the first translational joint. As a result, the static friction force grew only in the second joint (see Fig. 4a). At time t ≈ 0.8 s, the static friction force in the second joint reached the maximum value allowed by Eq. (8). Since the joint friction force must not go beyond the limit and the mechanism had to remain stopped until the static friction in the other joint reaches maximum, the equations representing the mechanism had to be changed. The previously eliminated equation (representing the revolute joint) was brought back. At the same time, the constraint equation representing the frictional locking of the second joint was eliminated, and since then, the friction force (of magnitude μ S 2 N 2 ) was represented in terms of the generalized force Q F (q, λ). Equation (10), modified in that way, took the form of Eq. (20) and constituted a system of nine equations in nine unknowns (this time, however, nonlinear in λ). The simulation was resumed and continued up to time t ≈ 1.86 s when static friction in the first joint reached its limit. At that moment, the friction regime switched from stiction to sliding. The constraint representing joint locking was deleted, and the equations of motion (Eq. (6)) were integrated during the final stage of simulation.
In the second case, one of the two equations representing the revolute joint on slider 2 was eliminated. In that arrangement of constraints, no tangential loads were transmitted to the second translational joint (see Fig. 4b). The condition for stiction occurrence, Eq. (8), had been repeatedly checked; however, the maximum of static friction force had never been reached. As a result, the investigated mechanism did not move by the end of the simulation.
The pseudoinverse-based method of handling of redundant constraints allows us to preserve all constraints in the mathematical model of the multibody system. The method consists in finding a particular solution to indeterminate Eq. (10): where the symbol [A] + denotes the Moore-Penrose pseudoinverse of a matrix A [35]. In that way, one of infinitely many solutions of Eq. (10) is selected; however, the method of selection is purely mathematical, with no physics of the multibody system involved. Note that the minimum norm solution is obtained. The results calculated using the pseudoinverse method are presented in Fig. 4c. Unlike the previous cases, at the beginning of the simulation, the friction forces grew in both joints simultaneously. Up to time t ≈ 0.19 s, the solution was physically unfeasible since the friction forces had incoherent directions (the minimum norm condition applies to κ and λ jointly, that is why the multipliers κ start from nonzero values). The allowed maximum of static friction force in the first translational joint was achieved at time t ≈ 2.98 s, and thus the constraint equation representing the frictional locking of the first joint was eliminated, and since then, the generalized force Q F (q, λ) represented the joint friction force (of magnitude μ S 1 N 1 ). At time t ≈ 3.59 s, the friction in the second joint reached the limit, and the mechanism started to move.
In conclusions to the example, it should be emphasized that, depending on the utilized method of coping with redundant constraints, the simulated transition from stiction to sliding occurred at time t ≈ 1.86 s, or at time t ≈ 3.59 s, or never. Thus, when the constraint addition-deletion method is used, not only the static friction forces cannot be uniquely calculated, but also motion of the system cannot be realistically predicted.
All models of friction discussed in Section 2 describe, more or less accurately, the same physical phenomena. The Coulomb model with constraints addition-deletion attempts to deal with discontinuous friction law, whereas the other models exhibit no discontinuity at zero slip velocity. Nevertheless, similar results of calculations should be expected. It is reasonable to ask how the problems with solution uniqueness, encountered when the constraint addition-deletion method is used, are reflected in the other models?

Continuous approximations of the Coulomb model
In the approximate Coulomb model, the coefficient of friction depends on the relative velocity of sliding. Regardless of details of the function μ i (v i ), this dependency is particularly strong at the vicinity of zero relative velocity. Let the constant c i (the slope of the function μ i (v i ) at zero relative velocity) be defined as Note that in the case of μ i (v i ) given by Eq. (12), it is simply c i = μ S i /γ i . Generally, for all definitions of the function μ i (v i ) used in practical engineering simulations, the constant c i is a large number.
It is important that particular values of c i usually result from quite arbitrary decisions (made when defining function μ i (v i )) that are not fully justified by the physics of the system. For example, in the case of Eq. (12), the coefficient of static friction μ S i is known from measurements; however, the parameter γ i is subjectively selected to meet some numerical requirements. These arbitrary decisions may influence the properties of the model at the vicinity of zero relative velocity. On the other hand, when the relative velocity is sufficiently large, μ i (v i ) is close to the kinetic coefficient of friction μ K i . Hence, the arbitrarily made decisions do not influence much the results obtained away from zero relative velocity.
When the approximate Coulomb friction model is used, the stiction regime is modeled as sliding with a very small velocity. The magnitude of the physical static friction force in the joint, precisely defined in Eq. (11), can be approximately calculated as where q * represents the configuration at which the mechanism is stopped. The relative joint velocity is a function of generalized coordinates and velocities: v i = v i (q * ,q).
In the case of constraint addition-deletion method, generalized friction forces were represented by the product of transposed constraint Jacobian and Lagrange multipliers (Ψ q (q * ) T κ); now they are represented in terms of an external force Q F (q *  , F 1 , . . . , F k ).
The assumption from Section 2.1 (that exactly one solution to Eq. (6) exists) holds, and thus, at least formally, the problem of solution uniqueness does not emerge.
The crucial questions are: what has happened to nonuniqueness of static friction forces and how the solution will change if the parameters c i are changed? Note that the generalized velocitiesq are dependent since Eq. (2) must be fulfilled. Hence, in our case of a 1-DOF mechanism, if one element of the vectorq (say,q p ) is known, then the remaining elements may be calculated by Eq. (2). Thus, the following relationship between the relative joint velocities holds: where the vector w shows the ratios of joint velocities. These ratios are constant at the given configuration q * . Consequently, the proportions between the coefficients of friction μ i ≈ c i v i are also fixed. Looking at Eq. (23), it is obvious that different selections of constants c i results in different distributions of friction forces among joints (it should be mentioned, however, that the dependency between the joint friction forces and constants c i is nonlinear due to nonlinearity of Eq. (6) with respect to the Lagrange multipliers λ).
There is an analogy between modeling the Coulomb friction by using an approximate model and by constraints addition. The arbitrary selection of c i values may be compared to resolving the indeterminacy problem through imposing proportions between vector κ elements. Both methods unlikely lead to solutions embedded in physics of the mechanical system. To summarize, when the Coulomb model with constraints addition-deletion is used, problems with solution uniqueness are encountered. When the approximate Coulomb model is introduced, some subjective decisions on the slope of functions μ i (v i ) are made (the selection of slopes c i is not unique). These decisions influence the obtained results. Consequently, the nonuniqueness of solutions to Eq. (10) is replaced by the nonuniqueness of parameters of Eq. (6). From the practical point of view, in both types of friction model, the result is quite the same: a number of different solutions can be found, and the model offers no means to select the solution justified by the physics of the system (it is worth noting that this problem may be alleviated if a method of associating the parameters c i with crucial physical properties, for example, stiffness of bodies or joints, is proposed). The discussed problems are illustrated in the example that follows.

Example 3
The mechanism from the previous examples was investigated. The initial conditions and the time history of external loads were the same as in Example 2. This time, however, the approximate Coulomb friction model with function μ i (v i ) defined in Eq. (12) was used. Firstly, the threshold velocities were set to γ 1 = γ 2 = 1 · 10 −5 m/s. It was found that during the simulated period of time, the mechanism did not move (to be exact, did not move significantly). A nonzero friction force is accompanied by a nonzero velocity of sliding, and thus a negligibly small drift, irrelevant to the discussed problems, was observed (in the case of slider 1, it was 0.033 mm at the final time of simulation). The maximum velocities of sliding in translational joints were v 1 ≈ −9.41 · 10 −6 m/s and v 2 ≈ 3.92 · 10 −6 m/s, that is, both velocities were representing the stiction regime. Note that, in the initial configuration of the mechanism, the ratio of velocities of sliding is v 1 : v 2 = −12 : 5 (w ∼ [−12, 5] T ). Since the configuration change was negligible, this ratio of velocities was maintained during the whole simulation.
The calculated friction forces and the magnitudes of friction forces divided by the magnitudes of respective normal forces are presented in Fig. 5a. It is visible that both friction forces remain within the static friction limit.
The second simulation was performed for γ 1 = 2·10 −5 m/s. Note that a customary choice of all γ i equal to each other is the result of a subjective decision (constructional details and operational conditions of joints may be fairly different, and thus in many cases keeping all γ i equal is not sufficiently justified by the physics of the system). Making γ 1 twice greater than γ 2 is just another subjective decision.
The results of the second simulation, presented in Fig. 5b, are not only quantitatively but also qualitatively different from the previous results. At time t ≈ 1.96 s, the static friction force in the first joint reached the maximum and started to decrease since the velocity of sliding continued to increase. Almost simultaneously, the static friction force in the second joint jumped to its maximum allowed value and then dropped to the kinetic friction level. The friction forces were no longer able to prevent gross motion of the mechanism. The process of switching from static to kinetic regime was almost immediate (it took about 0.001 s) since the mechanism started to accelerate significantly.
It is worth noting that the maxima of static friction forces did not occurred simultaneously-the velocities of sliding reached values of γ 1 and γ 2 in different moments. This is an important feature of the approximate Coulomb model-the "full potential" of static friction forces is not exploited.
In conclusions to the example, it should be emphasized that, depending on the choice of parameters γ 1 and γ 2 , the simulated transition from stiction to sliding occurred at time t ≈ 1.96 s or did not occurred at all. Moreover, the distribution of friction and normal forces in joints was strongly influenced by this choice. It should be emphasized that γ 1 and γ 2 are parameters of secondary importance; their values are usually selected due to purely numerical reasons.
In the example, high sensitivity of calculated static friction forces to γ i parameters was illustrated. This unwelcome feature makes the applicability of the approximate Coulomb models questionable. This comment pertains to simulation of the stiction regime only since in the sliding regime, the parameters γ i (or, more generally, c i ) do not affect the calculated fiction importantly.

LuGre model of friction
Equations (13)-(15) form the basis for friction modeling in the LuGre model. The authors of the model in their fundamental paper [11] show that a spring-like behavior is observed in the stiction regime. This property of the LuGre model is crucial for the present discussion.
To get insight into the behavior of the model in the stiction regime, let us start with linearizing Eq. (13) around v i0 = 0 and z i0 = 0. The following is obtained [11]: Let x i denote the relative displacement in the ith joint (x i is the time integral of v i and can be calculated as a function of the multibody system coordinates: x i = x i (q)). From Eq. (25) it may be inferred that when the model state is close to v i0 and z i0 , the change in x i is close to the change in z i : where x i , a small displacement observed during the stiction regime, is often called a microslip or a presliding displacement [12,17]. Substituting Eqs. (25) and (26) into Eq. (15) gives Note that at the given configuration q * , at which the 1-DOF mechanism is locked due to static friction, not only proportions between the joint velocities v i are temporarily fixed, but proportions between the microslips x i are fixed as well. Hence, situation observed in the case of the LuGre model is analogous to the situation discussed in Section 3.4. For a given set of parameters σ 0 i and σ 1 i , at mechanism configuration q * , the ratios of coefficient of friction are settled. Consequently, in the stiction regime, the distribution of friction forces among joints is decided by σ 0 i , σ 1 i , and q * . If the parameters σ 0 i or σ 1 i change, then the solution change accordingly.
It is worth noting that the primary parameters of the model of friction, that is, the coefficients of static and kinetic friction, μ S i and μ K i , respectively, can be relatively easily found in the scientific literature or trustworthily measured. The secondary parameters σ 0 i , σ 1 i , and ϑ i are less frequently measured or published. In many cases, rough estimations of these parameters are used during simulations (especially when a multibody model is not accompanied by a hardware prototype). In that context, high sensitivity of the results to the values of secondary parameters is particularly unwelcome. If test simulations confirm high sensitivity to uncertain parameters, then the results are doubtful. The next example demonstrates the discussed issues.

Example 4
The same mechanism was investigated once more. The initial conditions and the time history of external loads were the same as in Examples 2 and 3. The LuGre model of joint friction was used. The first simulation was performed with the following set of parameters: σ 0 1 = σ 0 2 = 1 · 10 5 1/m, σ 1 1 = σ 1 2 = 300 s/m, ϑ 1 = ϑ 2 = 1 · 10 −3 m/s. At time t ≈ 1.94 s, the static friction force in the first joint reached the maximum and started to decrease. The static friction force in the second joint reached its maximum a moment later, at time t ≈ 1.96 s. Both friction forces dropped to the kinetic friction level at time t ≈ 1.99 s.
The results of the first simulation are presented in Fig. 6a. As in the previous examples, the calculated friction forces and the magnitudes of friction forces divided by the magnitudes of respective normal forces are presented.
The second simulation was performed with doubled value of the parameter σ 0 1 (i.e., σ 0 1 = 2 · 10 5 1/m). The other parameters of the model of friction remained unchanged. It was found that during the simulated period of time, the mechanism did not move (at least, did not move significantly). Both friction forces remained in the stiction regime. The results of simulation are depicted in Fig. 6b.
Conclusions to this example are similar to conclusions drawn in the previous cases. In the stiction regime, the distribution of friction and normal forces in joints is strongly influenced by the choice of some parameters of the model of friction. These parameters are usually considered as parameters of secondary importance. The problem of solution uniqueness, experienced in the case of the Coulomb model with constraints addition, is replaced by the problem of high sensitivity of the solution to some parameters with hardly measurable values.
The discussion in Section 3 shows that as long as relative joint velocities of a rigid body model are in the vicinity of zero, that is, when the stiction regime is simulated, the obtained results are doubtful. Depending on the friction model applied, solutions are either nonunique or sensitive to some arbitrarily chosen or uncertain parameters. Hence, the rigid body models of multibody systems with joint friction have inherent limitations of great importance.

The role of flexibility in joint forces distribution
It is known that flexibility of bodies or joints is one of the main factors deciding on the distribution of reactions in overconstrained mechanisms [19,36]. Similarly, an important role of flexibility may be expected when redundant constraints are introduced to represent locking of the mechanism due to friction. It is not important whether it is "true" locking, represented by additional constraints, or "approximate" locking, represented by sufficiently large friction forces accompanied by microslips at joints. The point is that flexibility may strongly influence (or even determine) normal and frictional reactions in the stiction regime.
It should be emphasized that the flexibility introduced to the multibody model should properly reflect the real system properties. Subjective selections of bodies to be modeled as the flexible ones or arbitrary decisions on stiffness properties usually lead to unrealistic reaction solutions [19] and thus cannot be accepted in modeling of friction.
When flexibility of bodies or joints is taken into account, the number of system degrees of freedom increases. Thus, the impact of flexibility introduction on the modeling of friction with the use of the constraint addition-deletion method is easily predictable. Dependent equations do not appear in the mathematical model of a multibody system in the stiction regime, and therefore the problems with solution uniqueness do not emerge (usually, coordinates far outnumber the constraint equations). Normal reactions and friction forces in joints are unique and correspond to physics of the flexible multibody system. Joint locking due to friction is modeled in exactly the same way as in the case of rigid-body models. This time, however, the joints are not expected to transit form stiction to sliding simultaneously.
The impact of flexibility introduction on the modeling of friction with the use of the approximate Coulomb or the LuGre model needs to be separately commented. The deflections of flexible bodies are usually small; however, microslips or drifts due to stiction regime velocities at joints are small as well. Thus, deflections of flexible parts may be large enough to eliminate the coupling between joint displacements and the coupling between joint relative velocities. As it was pointed out in Sections 3.4 and 3.5, coupling boosts the role of secondary parameters of friction models and lets them to importantly influence the distribution of friction among joints. When modeling of joint friction is to be accompanied by modeling of flexibility of the system elements, it is very important to decide which bodies are to be modeled as flexible ones. The most reliable results are expected when elasticity of all parts is taken into account (in the case of the exemplary mechanism, it would be the connecting rod, two sliders, and two siding rails). In many situations, however, especially when some bodies are more flexible than the other, only selected parts may be treated as flexible ones, whereas the other may remain rigid (note that caution is recommended when dealing with joint loads in mixed rigid-flexible systems [19]).
Another important issue is how to model the flexibility of parts. There are various approaches to modeling of flexibility of bodies [37,38]. There are practically endless possibilities (e.g., when the finite element method is used, the same rigid body may be in many ways divided into finite elements, and different types of elements may be used). That is why it is difficult to discuss or investigate a generic flexible multibody system with joint friction.
In this paper, a very simple model of flexibility is used, just to show how flexibility influences the distribution of joint loads (with attention focused on friction forces in the stiction regime). For the exemplary mechanism, it is assumed that flexibility of the connecting rod dominates and only this part deflects importantly.

Simplified model of a flexible link
To model flexibility of a beam, advanced techniques can be used (e.g., the model discussed in Refs. [39,40] fits well absolute coordinates approach adopted here); however, to keep the model as easy as possible, much simpler approach is proposed here.
Consider a flexible link connected via revolute joints at its ends to other parts of a multibody system and assume that the link is carrying mainly axial loads. Let us build a simplified model of the link in which only axial deflections are considered and other elastic effects are neglected. To account for both elasticity of link and flexibility of joints at link ends, the simple model presented in Fig. 7 is developed.
In the simplified model, the link is split into three segments. The masses and lengths of the extreme segments constitute small fractions of mass and length of the middle segment. The segments are connected by translational joints and massless spring-damper elements (small viscous damping is added to avoid high-frequency vibrations). Coefficients of stiffness k 1 and k 2 may have different values since revolute joints at link ends may differ in construction details and the link itself may have a nonuniform cross sections along its length. The same applies to damping coefficients d 1 and d 2 .
The model of elasticity is extremely simple; nevertheless, it sufficiently captures properties essential for considerations presented in the next section. Moreover, flexibility of the link is modeled without resigning from the rigid body approach.

Joint friction in flexible body models
The limitations of rigid body models that manifest themselves when stiction regime is simulated may be alleviated or overcome when flexibility of bodies is taken into account. In this paper, the importance of flexibility and its influence of joint friction distribution, briefly discussed in Section 4.1, is investigated on the basis of a case study. The example that follows presents the results of flexibility introduction.
Example 5 The multibody model from previous examples was adjusted to account for elasticity of the link connecting the two sliders. The simplified model of link flexibility, described in Section 4.2, was used. The simulations from Examples 2-4 were repeated two times for two different sets of the flexible link stiffness parameters.
During the first series of simulations, the following set of flexibility parameters was utilized: k 1 = k 2 = 1 · 10 5 N/m, d 1 = d 2 = 2 · 10 3 Ns/m. Firstly, the constraint additiondeletion method (as in Example 2) was applied. It was found that the constraint equations were independent during the whole simulation, even when the supplementary constraints were added to represent joints locking due to stiction. Hence, there was no need to eliminate redundant constraints or to calculate the Moore-Penrose pseudoinverse. Consequently, just one variant of simulation was performed. At time t ≈ 1.90 s, friction in the first translational joint switched from stiction to sliding regime (the additional constraint was deactivated); a moment later, at time t ≈ 1.92 s, the same happened with friction in the second joint. The exemplary results of simulation are presented in Fig. 8a. It is apparent that the results do not resemble any of the results depicted in Fig. 4.
Secondly, the approximate Coulomb model of friction was utilized. Both simulations from Example 3 were repeated with exactly the same settings of friction model parameters. The results of simulations are shown in Fig. 8b. In the scale of Fig. 8b, the results obtained for different values of parameter γ 1 are indistinguishable (the differences between calculated friction forces are less than 0.5%). The curves from Fig. 8b are closely similar to those from Fig. 8a but fairly different from the results presented in Fig. 5. Thirdly, the LuGre model was employed. Two simulations from Example 4 were repeated (the parameters of model of friction remained unchanged). The results are presented in Fig. 8c. Minor differences between the two simulations are visible just after switching from stiction to sliding; nevertheless, the role of the parameter σ 0 1 is far less important than in the simulations of Example 4. Clearly, the results presented in Fig. 8c are similar to those of Fig. 8a and Fig. 8b and completely dissimilar to results shown in Fig. 6.
During the second series of simulations, the stiffness k 1 was tripled (k 1 = 3k 2 = 3 · 10 5 N/m), and the remaining flexibility parameters were left unchanged. Simulations with three different models of friction were repeated. Of course, all parameters of friction models used in particular simulations were kept unchanged. The results are presented in Fig. 9.
There were substantial consequences of change in the stiffness coefficient k 1 . In all conducted simulations, irrespectively of the model of friction used and irrespectively of values of secondary parameters, the joint friction remained in the stiction regime.
It is noticeable that almost the same results were obtained with the use of all tested models (only in the case of the approximate Coulomb model, small effects of drift are visible). On the contrary, comparison of curves in Fig. 9a with those in Fig. 4 reveals huge differences. Similarly, Fig. 9b can be compared with Fig. 5, and Fig. 9c with Fig. 6, to draw the same conclusions. Completely dissimilar results of the previous examples are contrasted with results very close to each other, obtained when flexibility of the multibody system is taken into account.
All results obtained in this example suggest that flexibility of the multibody system plays a crucial role in calculations of friction forces distribution among joints. Unlike the case of rigid-body mechanism, the results obtained for different models of joint friction are close to each other. Moreover, in the case of the Coulomb model with constraint addition-deletion, no problems with uniqueness of solution emerge, and in the case of the other models, sensitivity of the results to the values of secondary parameters of friction models is very weak.
In comments to the discussed example, it should be emphasized that the results obtained with the use of a flexible model can be recognized as unique and credible only when the flexibility of the real system is properly reflected by the model. The physics of the system must be sufficiently captured. Otherwise, flexibility is just another parameter (uncertain or sometimes quite arbitrarily decided) that influences the distribution of loads among the joints of the simulated multibody system.

Conclusions
Attempts to mathematically describe the phenomenon of friction originate from at least 15th century [5,12]. Nowadays, friction is well recognized: frictional effects in multibody systems are commonly taken into consideration in engineering practice, and some aspects of modeling of friction are taught during multibody dynamics courses at various universities [3,41]. On the other hand, friction in multibody systems is extensively investigated by 21st century researchers (e.g., [7,8,[16][17][18]28]). This paper was aimed at pointing the attention to some important properties of various approaches to friction modeling that manifest themselves in the stiction regime.
In this study, it is shown that application of well-known models of friction to a rigid multibody system with closed-loop kinematic chains causes some unwelcome effects. Three different models of friction, representing three typical approaches frequently used in multibody dynamics, were studied: a discontinuous model (Coulomb), a continuous model (approximate Coulomb), and a model with internal dynamics (LuGre). In each case, it was shown that rigid-body models of multibody systems offer no means to credibly predict friction forces in the stiction regime and its vicinity. Note that most recent studies show that similar effects are observed also for the Dahl model [42].
The Coulomb model of stiction regime, applied to rigid-body mechanisms by means of constraints addition-deletion, leads to a system of equations with no unique solution. Consequently, neither friction forces nor the moment of stiction to sliding transition can be uniquely calculated. In the case of an approximate Coulomb model, although formally equations have the unique solution, a high sensitivity of the results to the friction model parameters of secondary importance (whose values are often quite arbitrarily chosen, predominantly to meet some numerical requirements) is observed. Similarly, in the case of LuGre model, a high sensitivity of results to changes of values of some hardly measurable parameters of the friction model is observed. In the paper, origins of this high sensitivity are investigated and explained.
The first essential conclusion is that results of stiction regime simulations, obtained using purely rigid-body models, may be considered as doubtful.
Since rigid-body models provide no clear information about the multibody system state in the stiction regime, a simplified flexible body model was studied. It was demonstrated that when the flexibility of mechanism bodies is taken into account, the selection of values of "secondary" parameters of models of friction only negligibly influences the results of simulations. On the contrary, the quantities that represent the flexibility of bodies and joints are crucial for friction forces distribution among joints and consequently for the simulated moment of stiction to sliding transition.
The second essential conclusion is that the flexibility of bodies may be the key factor deciding on friction forces distribution during the stiction phase.
Note that the coefficients of static and kinetic friction were recognized as primary factors characterizing frictional properties of the system. Their influence on the system behavior was regarded as obvious, and thus sensitivity of the simulated results to changes of the coefficients of friction was not investigated.
In future works, closed-loop systems with more than one degree of freedom should be investigated. Problems with high parametric sensitivity may be expected when the number of joints simultaneously locked by frictional forces is greater than the number of degrees of freedom and the whole system does not move. However, when some joints are in static and some other in kinetic regime, a generalization of results obtained here is not straightforward. Some additional problems, possibly of combinatorial nature may be foreseen. Another interesting question is whether other models of friction (e.g., those discussed in the recent review paper [43]) exhibit similar problems of high parametric sensitivity in the vicinity of stiction regime.