Nonsmooth analysis of three-dimensional slipping and rolling in the presence of dry friction

In this paper, the nonsmooth dynamics of two contacting rigid bodies is analysed in the presence of dry friction. In three dimensions, slipping can occur in continuously many directions. Then, the Coulomb friction model leads to a system of differential equations, which has a codimension-2 discontinuity set in the phase space. The new theory of extended Filippov systems is applied to analyse the dynamics of a rigid body moving on a fixed rigid plane to explore the possible transitions between the slipping and rolling behaviour. The paper focuses on finding the so-called limit directions of the slipping equations at the discontinuity. This leads to a complete qualitative description of the possible scenarios of the dynamics in the vicinity of the discontinuity. It is shown that the new approach consistently extends the information provided from the static friction force of the rolling behaviour. The methods are demonstrated on an application example.

discontinuous behaviour. By considering the simple Coulomb model in the two-dimensional (2D) contact problems, the friction force changes sign at zero relative velocity of the surfaces. The situation is similar but more complicated in the three-dimensional (3D) contact dynamics. Then, for infinitesimally small relative velocities, the Coulomb friction model provides continuously many directions of the friction force with a constant finite amplitude.
The direct substitution of the discontinuous friction models into the dynamical equations leads to discontinuous systems of differential equations. In the 2D case, the Coulomb friction leads to Filippov systems (for an overview and examples, see [7]). Considering the friction as a set-valued force law leads to differential inclusions, which is a completely different point of view of modelling (see [13] for an overview). A further approach can be found in [6,16,17].
The generalization of the Filippov systems to codimension-2 discontinuity sets in the phase space leads to the concept of extended Filippov systems (see [1] and [4]). This type of differential equation can be used for modelling and analysis of 3D mechanical systems with dry friction, which was demonstrated in specific mechanical examples in [4]. The early results about two general contacting bodies have been presented by the authors in [3].
In this paper, we analyse the dynamics of a single rigid body in contact with a fixed rigid plane. During the motion of the body, rolling or slipping can occur, and the slipping case reveals to be described by an extended Filippov system. We focus on the transitions between the slipping and rolling dynamics when applying the theory of extended Filippov systems. The so-called limit directions can be determined, which are strongly connected to the possible slipping-rolling transitions. One of our main motivations is to provide a deeper understanding of the qualitative dynamics in the neighbourhood of the discontinuity. But the results makes the possibility for new numerical methods for simulating these mechanical systems, as well.
The paper is organized as follows: In Sect. 2, the dynamic equations of the moving body are derived by appropriate choice of the state variables for the subsequent analysis. In Sect. 3, the basic concepts and definitions of extended Filippov systems are presented. The main part of the paper is Sect. 4, where the theory of extended Filippov systems is applied to the dynamics of the moving body. From the analysis, we get four typical cases of the limit directions. In Sect. 5, the mechanical consequence of the four cases is explained, and the relation with the static friction force is presented. In Sect. 6, the results are demonstrated on a mechanical example. In Sect. 7, an overview can be found about the possible extension of the results to more complicated contact models. This paper is a significantly extended version of the conference paper [3]. The content of Sects. 2-5 has been rearranged and improved, and most importantly, the former conjectures have been developed into a series of proved statements about the possible slippingrolling transitions, as can be found in the current paper. Sections 6 and 7 are completely new.

Dynamics of a rigid body on a flat surface
We consider a rigid body moving in contact with a fixed rigid plane (see Fig. 1). It is assumed that the at any moment, the body and the plane are touching each other in a single contact point, denoted by P. The centre of the gravity of the body is denoted by C. The notation of the important quantities are summarized in Table 1.

Kinematics
Let v C and v P denote the velocities of the points C and P, respectively, and ω denotes the angular velocity of the body. The relation between these quantities of the Fig. 1 Sketch of the analysed mechanical scenario. A rigid body is moving in contact with a fixed plane. In this paper, the transitions between slipping and rolling are investigated Reciprocal of the normal curvature of the body in the direction of motion v C , v P Velocities of C and P u 1 , u 2 Components of the velocity v P (slipping) ω Angular velocity of the body ω 1 , ω 2 , ω 3 Components of ω q Vector of generalized coordinates s Vector of quasi-velocities rigid body is given by where r PC is the position vector of C measured from the contact point P.
During the motion of the body, the contact point P corresponds to different material points of the rigid body. We can take the time derivatives of (1) by using two different approaches: either by following the material point of the body currently located at P, or, following the motion of the instantaneous geometric contact point P.
In the first case of considering P as a material point, the differentiation of (1) gives the usual acceleration formula (2) where a C and a P denote the acceleration of the points C and P, respectively, andω denotes the angular acceleration of the body.
However, when differentiating (1) by considering P as the contact point, we get a C =v P +ω × r PC + ω ×ṙ PC .
(3) The time derivativev P is the rate of change of the velocity of the instantaneous contact point P, which is not equal to the acceleration a P of the material point of P in (2). The time derivative of r PC can be written aṡ (4) In the bracket, the two terms are the velocity v P of the motion of the material point at P and the velocity w P of the translation of the contact point on the surface of the body. This latter quantity depends on the rotation and the local curvature of the body, and it can be calculated by where n is the unit normal vector of the rigid plane and ρ n denotes the reciprocal of the normal curvature of the surface of the body in the plane determined by ω×n. In the general case, the normal curvature is determined by the second fundamental form of the surface of the body (see [12], p. 206). In the case of simple geometries, the quantities ρ or w P can be often found intuitively.
What is the point of writing a C in the form of (3)? In the subsequent calculations, we are using the components of v P as phase variables. As the Coulomb friction model is discontinuous exactly at v P = 0, this choice of variables, the discontinuity set of the resulting differential equation can be treated easily.

Dynamics
The effect of the rigid plane on the body is modelled by a single force F P acting at P. This force can be separated into a normal force N n and a friction force F f in the form We assume that the normal force N is strictly positive and the body remains in permanent contact with the plane (the scalar product v P · n is zero). That is, the effects of loosing contact, impact without collision, and consistency problems of the Painleve paradox (see [11,15]) are excluded from this analysis. All other forces and torques acting on the body are substituted by a single force F C and a torque T C acting at the centre of gravity C. Let m denote the mass of the body and J is the moment of inertia tensor with respect to the point C. Then, the Newton-Euler equations of the body are In addition, we have to consider supplementary conditions about contact between the body and the plane, according to the rolling or slipping. In the case of rolling, the kinematic constraint is satisfied. In the case of slipping, the friction force F f is modelled by the three-dimensional Coulomb friction law, where μ denotes the friction coefficient. Equations (6)-(9) lead to a system of differential equations in the case of rolling or slipping. In this paper, we focus on the slipping equations with special attention to their behaviour close to the rolling state.

Differential equations for the slipping case
By considering (1)-(5), the Newton-Euler equations (7) can be written into the forṁ We introduce coordinates both on the displacement and the velocity level. Thus, a set of first-order ordinary differential equations (ODEs) is obtained from (10)- (11).
In the vicinity of a chosen initial state, the position and the orientation of the rigid body by five generalized coordinates; the sixth degree of freedom is constrained by the contact between the body and the plane. Let the generalized coordinates be denoted by Note that along the paper, the vectors with a threedimensional physical meaning are denoted by boldface symbols (e.g. F C , v P ), but all other vectors are not distinguished from scalars in notation (e.g. q, s, x). Instead of describing the velocity level by the time derivatives of (12), we choose quasi-velocities (see [14], p. 254) independently from the generalized coordinates (12) to describe the velocity state of the body. By choosing two orthogonal unit vectors t 1 and t 2 parallel to the rigid plane (see Fig. 1), we get an orthonormal basis (t 1 , t 2 , n). In this coordinate system, v P and ω can be written as (13) where the components are chosen as quasi-velocities in the form These five linearly independent variables fully describe the velocity state of the body for any state of general coordinates. That is, the time derivatives of the generalized coordinates can be written aṡ where K (q) is a 5-by-5 matrix depending on the generalized coordinates themselves. By taking the time derivative of (13), we geṫ v P =u 1 t 1 +u 2 t 2 ,ω =ω 1 t 1 +ω 2 t 2 +ω 3 n, which can be substituted into the left-hand sides of (10) and (11). On the right-hand sides of (10)- (11), all quantities can be expressed by q, s, and N in the following way: The geometric quantities r PC (q) and ρ n (q) depend on the generalized coordinates. The moment of inertia tensor J(q) depends on q as well, because of the change of the orientation of the body. With the assumption of no explicit time dependence in the external forces, the resultants F C (q, s) and T C (q, s) are expressed by q and s. In the slipping case, (6), (9), and (13) leads to (17) Consequently, Eqs. (15), (10), and (11) form a set of six differential-algebraic equations in the generalized coordinates (12), the quasi-velocities (14), and the normal force N .
Equations (10)- (11) are linear in the derivatives of the quasi-velocities and N ; it can be solved in the forṁ where f s (q, s) and f N (q, s) denote the formal dependence on the variables. Equations (15) and (18) form a system of ten firstorder ODEs for the variables (12) and (14). Due to the discontinuity of the contact force (17), the system is not defined at u 1 = u 2 = 0, which corresponds to the rolling behaviour. For the rolling states, a different set of differential equations can be derived by excluding the slipping Coulomb law (9) but including the rolling constraint (8). To obtain a deeper insight to the switches between rolling and slipping, this paper focuses on the analysis of slipping system (15) and (18) in the vicinity of the discontinuity u 1 = u 2 = 0.
Note that this discontinuity is located at the states where two variables (u 1 and u 2 ) are zero at the same time. For the analysis of differential equations with such discontinuity, we can use effectively the theory of extended Filippov systems, which is presented briefly in the next section.

Overview of extended Filippov systems
The concept of extended Filippov systems was introduced by the authors in [1] and [4]. Roughly speaking, these dynamical systems are vector fields containing m − 2-dimensional discontinuities in the mdimensional phase space. We will show in Sect. 4 that the contact problem of the rigid body presented in Sect. 2 leads to an extended Filippov system.
In this section, only the most important concepts and definitions of these dynamical systems are presented, which are utilized in the subsequent analysis of the mechanical system. For a more detailed presentation of

F(x)
Vector field of the system Orthogonal space of at x 0 the theory of extended Filippov systems, see [1] and [4]. The notation of the important corresponding quantities can be found in Table 2. Consider a domain D ∈ R m containing an m − 2dimensional smooth manifold ⊂ D. This codimension-2 subset consists of the points where the vector field of the system is discontinuous in the sense of the following definition. At a chosen point x 0 ∈ , let us denote the tangent space by T x 0 , and its orthogonal complement by O x 0 (see Fig. 2). By considering the usual scalar product ., . in R m , the orthogonal complement is defined by Consequently, the direct product of the two-dimensional T x 0 and the m −2-dimensional O x 0 spans the whole vector space R m . Consider a point x 0 ∈ and choose two vectors n 1 (x 0 ) and n 2 (x 0 ) in O x 0 depending smoothly on x 0 with n 1 , n 1 = n 2 , n 2 = 1 and n 1 , n 2 =0. In other words, n 1 (x 0 ) and n 2 (x 0 ) generate an orthonormal basis of O x 0 . Then, let us define in which function n(φ) maps the interval [0, 2π) onto the set unit vectors of O x 0 . The parameter φ ∈ [0, 2π) can be imagined as an angle corresponding to a direction which is orthogonal to at x 0 (see Fig. 2).

Fig. 2
Basic concepts of extended Filippov systems. The codimension-2 discontinuity set is depicted as a curve (1D) in a 3D phase space, but it possesses more dimensions in higher dimensional systems. There are continuously many unit vectors n being perpendicular to the discontinuity set at any point x 0 . The possible normal directions are parameterized by an angle φ. The vector field F is discontinuous at , and it possesses a directional limit F * (φ) for any direction φ. The set of these limits is called the limit vector field Definition 1 Consider the vector fielḋ with the m − 2-dimensional smooth manifold . The system (22) is called an extended Filippov system if the following properties are satisfied: (a) The vector field F is smooth on D\ .
In the sense of Definition 1, is called a codimension-2 discontinuity manifold of F. At a chosen point x 0 ∈ , the function F * (φ) is called the limit vector field of F (see Fig. 2), which contains the directional limits of F from the different directions parameterized by φ.
The three conditions of Definition 1 formally express that there is no discontinuity outside (see (a)), there is indeed discontinuity at any point of (see (c)), and the directional limit does not diverge from any direction (see (b)).
By projecting the limit vector field F * (φ) on O x 0 , we get the components For the subsequent analysis, it is useful to write up the components of F * (φ) considering the component being parallel and perpendicular to the corresponding normal vector n(φ). Let us define The function R(φ) contains the behaviour of vector field in the radially inward or outward direction. The function V (φ) gives the circumferential, rotating behaviour of the vector field around the discontinuity set (see Fig. 3).
Let us now define the concept of limit directions, which are strongly connected to the behaviour of the trajectories at the discontinuity.

Definition 2
Consider an extended Filippov systeṁ x = F(x) and a point x 0 ∈ of the discontinuity manifold. The roots of the equation V (φ) = 0 are called the limit directions of x 0 with respect to F.
This vector can be separated to the components R(φ) and V (φ) by using the direction of the corresponding normal vector n(φ). The function R gives the radial behaviour of the vector field around the discontinuity, and V expresses the circumferential behaviour, which is the key of finding the limit directions It can be proved (see [4]) that if x 0 ∈ possesses at least one limit direction, then all trajectories tending to x 0 (either forward or backward time) approach x 0 along the limit directions.
In this sense, the limit directions are somewhat analogous to the eigenvectors of equilibrium points of smooth systems, but there are fundamental differences. Firstly, an eigenvector of a saddle or node is bi-directional (corresponding to a line), while a limit direction is uni-directional (corresponding to a halfline). Secondly, the eigenvectors of equilibria correspond to infinite-time (exponential) convergence of the solutions in forward or backward time, while trajectories reach x 0 in finite time in forward (attracting) or backward (repelling) direction of time.
By continuing the analogy with the equilibria, we can separate the node-like (sliding) and saddle-like (crossing) behaviour.

Definition 4
Consider a point x 0 ∈ which possesses at least one limit direction. If all the limit directions are either attracting or repelling, then we say that x 0 is located in the sliding region of . If there is at least one attracting and one repelling limit direction, then we say that x 0 is located in the crossing region of .
The terminology of crossing and sliding was introduced in [1] and [4] by generalizing of the crossing and sliding region of classical Filippov systems with codimension-1 discontinuities (see [7]). In the crossing case, there is at least one incoming and one leaving half-trajectory at x 0 , which can be concatenated to a trajectory crossing through at x 0 . In the sliding case, there are either only incoming or only leaving trajectories and the dynamics of F gets stuck into in forward or backward time, respectively. Then, the so-called sliding dynamics generated inside the discontinuity manifold . For the derivations and a more detailed explanation, see [4].
The introduction of the extended Filippov systems was originally motivated by 3D contact problems of rigid bodies. In these mechanical problems, the slipping of the bodies in the presence of Coulomb friction leads to extended Filippov systems, and the rolling or sticking of the bodies corresponds to the sliding dynamics inside the discontinuity manifold. In the following central part of the paper, the analysis of the limit directions is applied to explore the transitions between slipping and rolling between the bodies.

The resulting extended Filippov system
The full phase space of the body consists of the quasivelocities (14) and the generalized coordinates (12). Hence, the state variable vector can be written in the form Consequently, x ∈ D ⊂ R 10 . By composing the vector field from (18) and (15) in the form F = ( f (s, q), K (q) · s), the dynamics of the body can be simply written aṡ The vector field F(x) is discontinuous due to the terms u 1 / u 2 1 + u 2 2 and u 2 / u 2 1 + u 2 2 originating from the contact force (17). Equations (10)- (11) show that the final form of F(x) depends linearly on the contact force F P . Hence, the resulting vector field can be written in the form where A(x), B(x), and C(x) are smooth vector fields.
The system (29) is smooth everywhere except in u 1 = u 2 = 0. That is, the discontinuity set is which is a codimension-2 (8 dimensional) discontinuity.  (21) becomes and direct calculation of (23) leads to These quantities come from the coefficients of the nonsmooth terms in the expression of F P in (17). These coefficients are nonzero because of μ > 0 (there is indeed friction) and N > 0 (required in Sect. 2.2). Therefore, the non-singular linear operations on F P in (10)- (11) show that it is not possible to obtain zero for all components of A(x 0 ) and B(x 0 ) at the same time.
Consequently, all conditions of Definition 1 are satisfied, and thus, (29) is an extended Filippov system.

Analysis of the limit directions of the system
The discontinuity set of (29) is defined by u 1 = u 2 = 0 (see (30)), which corresponds to the rolling of the body on the plane. In this subsection, we categorize the points of according to the number and type of limit directions, which are strongly connected to the transitions between rolling and slipping.
For a chosen point x 0 ∈ , the projection (24) of the limit vector field can be calculated by simply taking first two components of (32). Hence, we get where A 1 . . . C 2 denote the first two components of A(x 0 ), B(x 0 ) and C(x 0 ) from (32) without denoting the (smooth) dependence on x 0 . Then, the functions R and V from (26) become From Definition 2, the limit directions are the zeroes of the function V (φ). In this section, we derive conditions from the coefficients of (33) and (34)-(35), which determine the type and number of the limit directions.

Possible formal simplifications
Proof Let us calculate the coefficients A 2 and B 1 formally from (10)- (11). By using the notations and we get We obtained that B 1 = A 2 .

Proposition 2
The coefficients A 1 and B 2 in (38) are strictly negative for the physically relevant parameters.
Proof By performing similar direct calculation as in the previous proof, we get The moment of inertia tensor J is positive definite. In the inner bracket of the expression of A 1 in (39), the expression can be obtained as the bilinear mapping of the vector r 3 t 2 −r 2 n by the positive definite matrix mJ −1 ; thus, this expression is positive. Similarly, the expression of B 2 contains the bilinear mapping of r 3 t 1 −r 1 n by mJ −1 , which is positive, again. If there is non-vanishing contact force (N > 0) and there is indeed friction (μ > 0), then we obtain A 1 < 0 and B 2 < 0.

Proposition 3
Consider the transformation of the variables u 1 and u 2 defined by u 1 = u 1 cos δ + u 2 sin δ, Then, δ can be chosen such that the coefficients A 2 = B 1 vanish, which denote the transformed coefficients corresponding to A 2 and B 1 in (33).
Proof By performing the transformation (40), the relation of the original and transformed coefficients of (33) is By choosing tan 2δ = 2 A 2 /(A 1 − B 2 ), the coefficients A 2 = B 1 are eliminated.
It follows from Proposition 3 that with an appropriate choice of the basis vectors t 1 and t 2 in the tangent plane of the body, (33) can be written into the form without the loss of generality. In this form, (43)-(44) become

Possible number of limit directions
In the expression of (35) and (44), V (φ) is a truncated Fourier series containing terms up to the second harmonics. According to [10], determining the zeroes of such function leads to the eigenvalue problem of a 4by-4 complex matrix. Alternatively, finding the zeroes of V (φ) is equivalent to solving the following fourthorder polynomial equation.
Proof Equation (45) is a fourth-order polynomial in cos φ, which leads to maximum four different roots for φ on the interval φ ∈ [0, 2π). In the degenerate case when V (φ) is identically zero, all φ ∈ [0, 2π) are limit directions.
Proof In the form (44) of V (φ), the constant term vanishes. Therefore, V (φ) is a periodic continuous function with zero mean value. Thus, it needs to have at least two zeroes on φ ∈ [0, 2π).

Proposition 7 The function V (φ) in (44) has three zeroes if and only if
Proof The periodic differentiable function V (φ) can have odd number of roots only if it has a double root By substituting (44) into the two equations of (47), direct calculation leads to tan φ 1 = −(C 2 /C 1 ) 1/3 and (46).

Theorem 2
Consider a point x 0 ∈ of the system (29). The following cases can occur.

Attracting and repelling limit directions
Theorem 3 Consider a point x 0 ∈ of the system (29). The following cases can occur.
, then x 0 has at least one attracting limit direction and exactly one repelling limit direction.
, then x 0 has attracting limit directions and a limit directions on the boundary of being attracting and repelling.
Proof In Point 3 of the Theorem, we can find the condition of Proposition 8 which ensures the existence of a limit directions between being repelling and attracting (see Definition 3). This condition (48) separates the space of the coefficient A 1 , B 2 , C 1 , C 2 into two regions, and the number of the repelling (or attracting) direction changes by one when crossing the boundary (48).
In the limit case (C 2 1 + C 2 2 )/ max(A 1 , B 2 ) → 0, the terms with sin φ and cos φ are negligible in (44), and we get which is always negative because A 1 , B 2 < 0 (see Proposition 2). That is, there is no possibility for a repelling limit direction φ 1 with R(φ 1 ) > 0, and Point 1 of the Theorem is proved. Proof of Point 2 comes from the fact that exactly one limit direction changes between attracting and repelling on the curve (48).

The four generic cases
The boundary (46)  In the case |A 1 − B 2 | ≤ min(|A 1 |, |B 1 |), the case IV is absent from the plane of C 1 and C 2 (see Fig. 5). Note that the coefficients A 1 , B 2 , C 1 , and C 2 depend smoothly on the state variables ω 1 . . . q 5 from (27). Thus, the boundaries in Figs. 4-5 could be mapped onto seven-dimensional surfaces in the phase space, where they divide the eight-dimensional discontinuity set into regions according to the behaviours I-IV.
The typical structure of the vector field in the four cases can be seen in Figs  The ellipse corresponds to (48) and the star-like boundary corresponds to (46). A similar structure of the figure is obtained for induces that the trajectories approach the discontinuity along the attracting limit directions (denoted by solid lines), and they leave the discontinuity along the repelling limit directions (denoted by dashed lines).

Angularly stable and unstable limit directions
We can see in Figs. 6-9 that the the trajectories are different in the neighbourhood of different attracting directions. For example in Fig. 6, most trajectories seem to be follow the direction φ 1 , and not φ 2 . This difference can be explained by the subsequent analysis.
Let us consider the systeṁ which is an asymptotic approximation of projection of the system (29) into the normal plane of u 1 and u 2 in the limit case u 2 1 + u 2 2 → 0. By using the transformation u := u 2 1 + u 2 2 , φ = arctan(u 2 , u 1 ), and an appropriate transformation of time, the trajectories of (50) are mapped to the trajectories of the system where the dash denoted the differentiation with respect to the new time variable. Note that the solutions of (52) can be determined independently from (51). The Taylor expansion of (52) around an equilibrium φ 1 with V (φ 1 ) = 0 is given by where O 2 denotes the higher order terms. The linear stability of φ 1 of equation (52) is determined by the sign of dV (φ)/dφ at φ = φ 1 . The equilibrium points of (52) corresponds to the limit directions of the original system (see Definition 2). Hence, the terms stable and unstable of the equilibrium point can be transferred to the limit directions.
In the angularly stable case, the limit direction is attracting the adjacent trajectories (see φ 1 in Fig. 6). In the  Case  I  II  III  IV   Total number of limit directions  2  4  2  4 Attracting, angularly stable 1 2 0 1 Attracting, angularly unstable 1 2 1 2 Repelling, angularly stable 0 0 1 1 Repelling, angularly unstable 0 0 0 0 angularly unstable case, the adjacent trajectories get far from the limit directions in the sense of the angle φ (see φ 2 in Fig. 6). The special case dV (φ)/dφ = 0 is a fold bifurcation of (52), which corresponds to the fold of limit directions in (29). This condition was already discussed in Proposition 7 (see (47)); thus, the fold of directions coincides with the condition (46). The next Proposition completes our analysis of the limit directions of (29).

Proposition 9
If a limit direction of (29) is repelling, then it is an angularly stable limit direction.
Proof A limit direction can change from attracting to repelling only on the boundary (48). Let us substitute (48) and the corresponding value of φ 1 from the proof of Proposition 8 into dV (φ)/dφ. Then, we get which is always negative due to Proposition 2. That is, the limit direction is angularly stable.
To summarize our results, the number and properties of the limit trajectories can be found in Table 3.

Mechanical consequence of the limit directions
The system F(x) in (28) was introduced to describe the slipping behaviour of the body. The discontinuity manifold is the set u 1 = u 2 = 0, which coincides with the condition of the rolling constraint (8).
In this subsection, the slipping-rolling transitions are analysed by considering purely the limit directions of the slipping equations determined in Sect. 4. The relation to the dynamical conditions of the rolling equations is presented in the next subsection.

Case I: 2 attracting directions
In this case, all trajectories in the vicinity of x 0 tend to the discontinuity at u 1 = u 2 = 0 (see Fig. 6). That is, the behaviour of the body turns from slipping to rolling. It is proved in [4] that the trajectories reach u 1 = u 2 = 0 in finite time. In some sense, the rolling motion is stable with respect to slipping perturbations, because the effect of a small perturbation in u 1 and u 2 is eliminated by the dynamics in finite time.
Note that almost all solutions reach the rolling state along the angularly stable limit direction (φ 1 in Fig. 6). The angle φ can be imagined not only as an angle in the phase space but as an angle of the slipping velocity v P , as well. Therefore, the dominant behaviour of the limit direction φ 1 causes that the slipping velocity points typically into the direction φ 1 when the motion changes from slipping into rolling. There is only a single trajectory which approaches the state u 1 = u 2 = 0 from the direction φ 2 . The trajectories close to φ 2 contain a high-curvature turning when reaching x 0 . Hence, the direction of the friction force changes rapidly just before the transition from slipping to rolling.

Case II: 4 attracting directions
This case has a behaviour similar to Case I: all surrounding trajectories tend to the discontinuity (u 1 = u 2 = 0) in finite time (see Fig. 7. From mechanical point of view, this means that the rolling motion is realizable, because small perturbations causing slipping are eliminated by the dynamics and the body starts rolling again. But in contrast to Case I, here, we have four limit directions and there are two angularly stable limit directions (φ 1 and φ 3 in Fig. 7). That is, there are two typical directions of the slipping velocity when the body is in transition from slipping to rolling. The angularly unstable directions φ 2 and φ 4 behave as separatrices. In the regions φ 4 < φ < φ 1 and φ 1 < φ < φ 2 , all trajectories approach the rolling state along the direction φ 1 . In the regions φ 2 < φ < φ 3 and φ 3 < φ < φ 4 , the trajectories tend to the limit direction φ 3 .

Case III: 1 repelling and 1 attracting direction
When an attracting limit direction turns into repelling, the structure of the phase changes significantly. In Case III, we can find an attracting and a repelling limit direction (see Fig. 8). The attracting direction is angularly unstable and the repelling direction is angularly stable (according to Proposition 9). That is, the typical behaviour of the system is slipping, and almost all trajectories avoid the discontinuity at u 1 = u 2 = 0.
In the vicinity of the discontinuity set (rolling behaviour), the trajectories tend to the repelling limit direction φ 1 and they diverge from the rolling state. That is, a slipping motion is generated with a typical direction φ 1 of the slipping velocity. There exists one single trajectory which reaches the discontinuity, and this happens along the limit direction φ 2 . But the system reaches the rolling state just for a moment, and it starts slipping immediately in the direction of φ 2 .

Case IV: 1 repelling and 3 attracting directions
In Case IV, there are a repelling and three attracting limit directions, which lead to the most complicated scenario of the four cases (see Fig. 9. The two angularly unstable attracting directions φ 2 and φ 4 are the separatrices between two typical regions of the phase plane. Between φ 4 < φ < φ 1 and φ 1 < φ < φ 2 , all trajectories avoid the discontinuity and the trajectories tend to the repelling limit direction φ 1 . In this sense, this case is similar to Case III, and the typical behaviour of the system is slipping.
However, in the regions φ 2 < φ < φ 3 and φ 3 < φ < φ 4 , the trajectories tend to the angularly stable attracting direction φ 3 , and they all reach the discontinuity at u 1 = u 2 = 0. There is rolling only for a moment, and the system starts slipping with a slipping velocity described by the direction φ 1 . In contrast to Case III, not just a single trajectory is connected to the discontinuity, but a large portion of the phase plane tends to u 1 = u 2 = 0. That is, the typical long-time behaviour is slipping, but for many initial conditions, rolling can occur for a moment.

Summary of the typical types of behaviour
After the detailed survey of the possible types of solution, let us summarize the typical four cases of behaviour from the mechanical point of view.

Corollary 1 Consider a rolling state of the body and let us perturb the motion by a small amount of slipping.
According to the chosen state of rolling, the typical behaviour of the body is the following: -In Cases I and II, the perturbed body returns to rolling in finite time and then it maintains the rolling state. -In Case I, the slipping velocity vanishes from a certain direction for almost all perturbations (see Fig. 6). -In Case II, the slipping velocity vanishes from two certain directions for almost all perturbations (see Fig. 7). -In Cases III and IV, the perturbed body is unable to maintain a lasting rolling state and it continues slipping. -In Case III, the slipping velocity remains finite for almost all perturbations, thus, pure slipping continues (see Fig. 8).

-In Case IV, two types of behaviour occur accord-
ing to the direction of the perturbation. Either, the body continues pure slipping like in Case III. Or, the slipping velocity vanishes in finite time, rolling motion occurs for a single moment, and then, the body continues slipping again (see Fig. 9).

Comparison with the rolling condition of the friction law
Up to this point, we analysed the rolling-slipping transitions based on purely the phase space of the slipping system (29). It can be seen that a detailed, consistent structure of the behaviour can be obtained from this analysis. But what is the relation between these results and the ones from the rolling condition of the friction law with the static friction force?
The equations of the rolling vector field can be derived either from the Newton-Euler equations (7) with the rolling constraint (8), or directly from the limit vector field (32) of the slipping case.
In the latter case, we consider the dynamics on generated by F * (x) by a convex combination, which is called sliding dynamics in the terminology of Filippov systems and extended Filippov systems (see [1,7] and [4]). In mechanical problems, we have to be careful with this terminology because sliding dynamics correspond to the mechanical rolling and not to the mechanical slipping.
At a point x 0 ∈ , we search for the sliding vector F (x 0 ) ∈ T x 0 in the form where α is a [0, 2π) → [0, 1] function which satisfies 2π 0 α(φ)dφ = 1 (a convex combination). By direct calculation from (32)-(33), we get that the only such sliding vector is By using the reduced form (42) in the appropriately chosen coordinates, the formula (56) simplifies to If the body is rolling then u 1 = u 2 = 0, and the dynamics of the other variables is described by the system (56). Then, the formula (9) is not valid, and the friction force F f can be obtained as a constraint force. The dynamic condition of the rolling is usually expressed by the maximal admissible friction force in the form where μ 0 is the static friction coefficient. When the static and dynamic friction coefficient is equal (μ 0 = μ) then we can state the following theorem.
Theorem 4 Consider a state x 0 ∈ of rolling, when the static friction force is known from the constraints and μ 0 = μ.
1. If the rolling is strictly admitted by (58), that is, |F f | < μ 0 N , then x 0 possesses no repelling limit directions. 2. If the rolling is not admitted by (58), that is, |F f | > μ 0 N , then x 0 possesses a repelling limit direction. 3. In the special case |F f | = μ 0 N , x 0 possesses attracting limit directions and a limit direction on the boundary of being attracting and repelling.
Proof The rolling dynamics ensures u 1 = u 2 = 0 permanently, that is, the derivativesu 1 andu 2 have to be zero, as well. If μ 0 = μ then third statement of the theorem with |F f | = μ 0 N gives back the condition (9) of the dynamic friction. That is, the rolling and slipping dynamics is valid at the same time. In the slipping dynamics, the conditionu 1 =u 2 = 0 is equivalent to the condition (48) (see the proof of Proposition 8), which decides whether there exists a repelling limit direction or not (see Theorem 3). The magnitude |F f | of the friction force tends to zero when A 1 → 0 and B 2 → 0 [compare (6), (11), (17) and (29)]. Consequently, the three cases of Theorems 3 and 4 are pairwise equivalent. The theorem is valid for the simple Coulomb model with a uniform friction coefficient μ 0 = μ (see the topleft panel of Fig. 10). However, it is not valid for the stiction model where there are two distinct values of μ and μ 0 in (9) and (58), respectively (see the top-right panel of Fig. 10). In some sense, this stiction model provides inconsistent friction forces in the rolling and slipping cases. In this model, there are three discontinuities at |v P | = 0: the change of the sign of the friction force and the change between the static and dynamic friction in both directions. This degeneracy can be avoided by replacing the constant dynamic coefficient μ by a functioñ and then, the slipping Coulomb model (9) becomes The model (59)-(60) provides a Stribeck friction model without the viscous effect (see the bottom panel of Fig. 10). The parameter γ can be estimated empirically, and the limit case γ → ∞ gives back the stiction model. In fact, (59) can be replaced by any smooth functionμ(|v P |) satisfying  Proof By the first limit of (61), the dynamic friction coefficient tends to μ 0 when |v P | → 0. Then, μ can be replaced by μ 0 in F * and the related quantities all along the analysis of the paper. Then, the proof of Theorem 4 can be repeated.
Theorem 4 and Proposition 10 state that the analysis of the limit directions in Sect. 4 is consistent with the checking of the maximal admissible friction force. In Cases I and II with only attracting limit directions (Figs. 6-7), the condition (58) is satisfied (rolling is realizable). In Cases III and IV with a repelling limit direction (Figs. [8][9], the condition (58) is violated (rolling is not realizable).
That is, Theorem 3 decides whether rolling is possible or not based on the slipping equations, and we do not need to calculate the rolling dynamics (56) and check the condition (58) of rolling. This property can be useful especially in those systems where it is not possible to calculate the static friction force (see [1,2]). Moreover, the analysis of limit directions gives more detailed information than just deciding between rolling or slipping. The number and location of limit directions characterize the possible transitions between the slipping-rolling states and the direction of the slipping velocities at the transition.

Application example
Consider a wheel moving on a plane (see Fig. 11), where the symmetry axis of the wheel remains parallel with the plane (no tilting of the wheel). The wheel is modelled by a rigid disc with a radius ρ and a negligible thickness. The external forces acting on the wheel are the gravity force (mg), the steering moment (M s ), the driving moment (M d ), and the balancing moment (M b ), ensuring the horizontal orientation of the symmetry axis. All the other notations are the same as in Sects. 2-4 (see Tables 1). It is shown in [9] that the rolling motion of this model is equivalent to the motion of the Chaplygin-sleigh, which is an important benchmark problem of nonholonomic mechanics.
It is useful to fix the basis vectors t 1 and t 2 to the wheel such that t 2 is parallel to the symmetry axis of the wheel. Then, the Newton-Euler equations of the wheels in the form (10)- (11) becomė The location and orientation of the wheel on the plane do not appear in (62) due to the symmetry properties of the problem. By eliminating the trivial coordinate ω 1 ≡ 0, as well, the state vector (27) of the system can be written into the reduced form Then, the vector field F(x) becomes The discontinuity set defined by u 1 = u 2 = 0 is a plane of the variables ω 2 and ω 3 , and its selected point is denoted by x 0 = (0, 0, ω 2 , ω 3 ). The limit vector field (23) becomes In the form (32), the vectors A, B, andC at x 0 are given by and the coefficients (33) are As A 2 = B 1 = 0, the system is already in the form (42), and it is not necessary to transform the variables according to (40). By substituting these into the boundary curve (48) between slipping and rolling, we get The boundary curve (46) between the 2 and 4 limit directions becomes

M d mρ
These curves are visualized in Fig. 12 in a similar diagram similar to Fig. 4. There are three special values of |M d |, which is the absolute value of the driving moment: |M d | = 3μmgρ, |M d | = 2μmgρ, and |M d | ≈ 0.459μmgρ (see Fig. 12). By selecting typical values of M d in between these special values, the sketch of the discontinuity set can be found in Fig. 13. The discontinuity set is the plane of the variables ω 2 and ω 3 . In this plane, the boundary curves are transformed into hyperbolas due to the product ω 2 ω 3 in the expressions. It is not surprising that the rolling is realizable (Cases I and II) when the product of the angular  Fig. 12), the different regions correspond to the four cases I-IV of behaviour presented above (see Figs. [6][7][8][9] velocities are not larger than a critical value. Note that the 4 limit directions (Cases II and IV) appear at low values at the angular velocities.

Overview of more complex contact models
The results of the paper are based on the following modelling assumptions of the contact between the bodies: 1. the assumption of planar geometry of the fixed body, 2. the assumption of rigid bodies, 3. assuming that the contact state is either slipping or rolling, 4. the assumption of the Coulomb friction model, 5. the assumption of a well-defined, single contact point of the unloaded bodies.
In this section, we overview the possibilities of the extension of the analysis in the case of any of these restrictions are released.
If we replace the fixed plane by a rigid body with an arbitrary curved fixed surface, then the normal plane of the contacting surfaces is changing during the motion. It would make the dynamic equations more complicated, but it would not be a structural modification of the model. Thus, the different scenarios of the limit directions are expected to be preserved in this case.
If we do not neglect the deformation of the bodies, then we have to consider the formation of the contact area around the theoretical contact point. In the literature, it is usual to separate the motion into two parts: the rigid body motion of the whole body and the local deformations in the vicinity of the contact area. The Hertz theory (see, e.g. [25], p. 55 or [19] p. 84) assumes elliptical contact area and a parabolic distribution of the normal pressure between the bodies in the frictionless case. Similar but higher order theories exist for the normal pressures (see, e.g. [20]). Combining these models with friction leads to theoretical and computational challenges [19,20,25]. For the purpose of dynamical applications, analytical and semianalytical models of the contact forces can be derived from these theories. When we consider the combined effect of the slipping and drilling motion, the contact laws are determined by the Coulomb-Contensou friction model (see, e.g. [21]). The combined effect of slipping and rolling motion leads to the contact laws of creep models (see, e.g. [18]). By improving our analysis by these models, higher-codimensional (3)(4)(5) discontinuities are expected to appear. This can be the topic of the further research work. The concept of limit directions probably remain important in these cases, as well, to find the possible transitions between rolling and slipping.
We assumed that the normal contact force N is strictly positive, and the surfaces remain in permanent contact at the contact point. Then, depending on the state of the system, the friction model decides whether slipping or rolling behaviour occurs. However, when the contact force N decreases to zero, the bodies can separate from each other (lift-off ), and the nonsmooth behaviour of the dynamics with the discontinuity set vanishes. However, the switching of the contact and the no-contact states introduces a further discontinuity, containing the impacts of the bodies, as well. The generalization of the results of the paper to these cases would need additional extensive research work. The Painlevé paradox of the contact states [11,15] causes further complications.
In the analysis, we first considered the simple Coulomb law for modelling dry friction. Then, in Proposition 10, the results are generalized for a class of friction models similar to the Stribeck friction law. However, several different friction models can be found in the literature (see [22] and [23] for an overview). It is an open question how further effects like hysteresis (e.g. the Karnopp model) or internal variables (e.g. the Dahl model) modify the qualitative structure of the dynamics at the discontinuity.
A further complication can be the coexistence of multiple contact points between the contacting bodies. The results in [1] show that the concept of limit directions is applicable to two contact points, but still, a throughout analysis would be necessary. An even further case is the contact of conforming bodies, where there is a finite contact area even with the rigid body assumption. (The simplest example is a block moving on a plane.) Then, information is needed about the normal pressure distribution, and the integration on the contact area is expected to lead to highercodimensional discontinuities as we expected in the deformable models, too.

Conclusion
The dynamical equations of a rigid body were derived, in which body is in 3D slipping or rolling contact with a rigid plane in the presence of dry friction. It was shown that by assuming Coulomb friction model, the differential equations of the dynamics of the body lead to an extended Filippov system. That is, the phase space of the system contains a codimension-2 discontinuity set.
The nonsmooth differential equations of slipping were analysed by the methods of new theory of extended Filippov systems. The possible number and type of limit directions of the points of the discontinuity were determined, where the transitions between slipping and rolling occur. We got four structurally different cases of limit directions; there can be two or four limit directions, from which maximum one can be repelling. The effect of these scenarios of the mechanical behaviour was discussed in detail. It was shown that in case of simple Coulomb model and the Stribeck model, the limit directions lead to such conditions of rolling, which are consistent with the condition of maximum admissible friction force. Furthermore, the results of the new approach provides more information about the qualitative behaviour of these mechanical systems near the discontinuity. The result was demonstrated on an example of a wheel. A part of the further work would be to apply the results to other well-known systems such as the rolling disc (see [5,8]) and the classical skate problem (see, e.g. [9]). The considered contact model is clearly the simplest mechanical model which is capable to describe the problem of the different possible directions of transitions between rolling and slipping in three dimensions. However, several possibilities were presented in the last section to improve the contact model in different ways.
One further important direction of the subsequent research would be to utilize these results to develop effective and reliable numerical methods for simulation of these systems. The information of the structure of the trajectories and limit directions could help to find appropriate event-driven strategies similar to those of classical Filippov systems [24].