The Nonsmooth Dynamics of Combined Slip and Spin Motion Under Dry Friction

We investigate the motion of rigid bodies subject to combined slipping and spinning over a rough flat surface in the presence of dry friction. Integration of Coulomb friction forces over the contact area gives rise to a dynamical system with an isolated discontinuity of codimension 3. Recent results about such vector fields are applied to the motion of flat bodies under the assumption of known, time-independent distributions of normal contact forces and to general bodies where kinematic constraints enforce a state-dependent normal contact force distribution with a discontinuity at the sticking state. In both cases, the equations of motion are transformed into a smooth slow–fast dynamical system. The fixed points of the fast flow indicate the possible directions of combined slip–spin motion immediately before a body stops. We also introduce an approximation of the frictional forces and moments by the leading-order term of a spherical harmonic expansion, which allows for an explicit formulation of the equations of motion. The approximate model captures important empirically observed features of the motion. It is proven analytically and illustrated by examples that the number of fixed points in the approximate model is 2, 4, or 6.


Introduction
The presence of dry friction in mechanics often leads to dynamic (slipping) and static (sticking) contact states. In dynamic problems, these mechanical states can switch into each other. Slip-stick transitions bring about nonsmooth behaviour of the system. The conditions of transitions between stick and slip can be implemented in the models in different ways. If we focus on the sticking state, then the (static) friction forces are determined by the dynamical equations as constraint forces, and the transitions can be modelled to occur at the boundary of the set of admissible contact forces. Alternatively, if we focus on the slipping state at transitions, then the slipping dynamics can be analysed in the vicinity of the sticking states to determine the presence and type of transitions between stick and slip. The two approaches often complement each other, and they can be implemented in different ways in friction models.
Here, we focus on the latter, dynamic approach by analysing the state space of slipping dynamics to explore the stick-slip transitions. A two-dimensional model of two rigid bodies with a single contact point and Coulomb friction gives rise to a piecewise smooth dynamical system (often called Filippov system) containing a codimension-1 discontinuity set in the phase space. In this case, the stick-slip transitions can occur in two directions ("left" and "right"), and the corresponding mathematical theories are well established in the literature (see Leine and Nijmeijer 2004;di Bernardo et al. 2008;Jeffrey 2018). In three-dimensional problems with a single contact point, the state space contains a codimension-2 discontinuity. Then, from the continuously many possible slipping directions, one can find the so-called limit directions where the slip-stick transitions are realized Stepan 2018, 2019;Batlle 1996;Cheesman et al. 2021), and this analysis can be naturally extended for multiple contact points (Antali and Stepan 2021). In this paper, we consider finite contact regions with a combined slip and spin motion, where the integration of locally prescribed Coulomb friction forces leads to a codimension-3 discontinuity in the state space (Marques et al. 2016). Vector fields with higher-codimensional discontinuities have been recently analysed by the authors (Varkonyi and Antali 2021); those mathematical tools can be applied to the slip-stick transitions of bodies with a finite contact region.
As proposed by Varkonyi and Antali (2021), we transform the equations of slip motion to spherical coordinates, whereby the discontinuity is replaced by smooth dynamics with two characteristic time scales. The positions and the attractivity of fixed points of the fast dynamics-or limit directions in our terminology-determine the direction of terminal motion at slip-stick transitions, including the direction of slip and the ratio of velocity and angular velocity at the time of the transition. The limit directions of some simple mechanical systems have already been discovered (Goyal et al. 1991b); however, our work is the first one to describe their role in a systematic and general way.
In some simple cases, the resultant of frictional forces over the contact region can be given in terms of elliptic integral functions (Leine and Glocker 2003;Farkas et al. 2003;Malhotra 2005, 2007;Wolf et al. 2007;Magyari and Weidman 2001). Another class of problem that lends itself for analytical investigation is the case when the contact forces are concentrated at points or along lines introducing new forms of nonsmoothness (Goyal et al. 1991a;Hammond 2010;Borisov et al. 2016). However in general, the resultant can only be evaluated numerically. In order to enable analytical investigation, several closed-form approximations of frictional forces and torques have been proposed. These include Padé approximations (Kireenkov 2002a(Kireenkov , b, 2017Zhuvarlev and Kireenkov 2005), hyperbolic tangent functions Awrejcewicz 2011, 2012), concatenated Taylor series (Leine and Glocker 2003) and implicit quadratic equations (Goyal et al. 1991a;Howe and Cutkosky 1996). In addition, some works have proposed a priori closed-form friction models, neither derived by integration of local Coulomb forces nor derived by approximations of such integrals (Zhou et al. 2018;Walker and Leine 2019). In this paper, we propose a new and natural approximation, which is found to be the leading-order term in a spherical harmonic expansion of the exact model. Our approximation is represented by a linear function in some appropriately chosen variables. The simple form of the approximate model allows us to find explicit algorithms for identifying limit directions, and we also establish upper and lower bounds for the number of limit directions.
The structure of the paper is the following: In Sect. 2, the kinematic and dynamic description of a finite contact region is presented, and we identify the properties of the discontinuity, which originates from the resultant of the distributed dry friction forces. In Sect. 3, we review some recent mathematical tools for analysing codimension-3 discontinuities of vector fields (Varkonyi and Antali 2021). In Sect. 4, we introduce an approximate ellipsoidal model for the mapping between the variables of the velocity state and the friction forces. Section 5 contains the detailed analysis of the nonsmooth dynamics of a single rigid body under the approximate friction model, focusing on the sticking-slipping transitions. The new methods and results are demonstrated through mechanical examples in Sect. 6.

Problem Statement
We assume that the motion of the bodies can be decomposed into rigid body motion of the bodies and small local elastic deformations in the contact region. Therefore, the global elastic deformations of the bodies are neglected, and the effect of the local deformations is included in the friction forces (and torques) applied to the dynamical model of rigid bodies. This approach is common in the literature (see, e.g. Leine and Glocker 2003).
Friction phenomena greatly depend on the geometry of the contacting surfaces of the bodies. From the several possible situations, we focus on two typical cases: • (regular) point contact In this case, we consider a single theoretical contact point, and at this point, the curvatures of the contacting bodies are substantially different (nonconforming surfaces). Then, by considering the elastic deformation, a finitesized contact region is initiated around the theoretical contact point. We assume that the contact region is a plane shape, and its size is considered to be much less than the typical size of the bodies. A prototype example of this case is a ball moving on a plane.
• plane contact In this case, we consider two bodies with contacting planar surfaces. Then, the contact region is the overlapping part of the planar surfaces of the two bodies. The prototype example of this case is a block moving on a plane.
There are other possible situations beyond the scope of the present analysis, such as: • (general) conforming contact If the surfaces have approximately the same curvatures at the theoretical contact point, then the contact region has comparable size to the typical size of the bodies, and it can be a significantly curved surface. (Plane contact is the special case of conforming contact where the curvatures of both surfaces vanish.) • (regular) line contact The contacting surfaces can be conforming in one direction and nonconforming in another direction at the theoretical contact point. Then, the initiating contact region is a thin stripe. A typical example of this type can be the contact of two cylinders with parallel symmetry axes. • degenerate point and line contact If the edges or vertices of the bodies are involved in the contact, then the deformations and the geometry of the contact region can be irregular compared to smooth surfaces.
Although the cases of point contact and the plane contact are substantially different from each other in several aspects, they have some similarities, which make it possible to analyse them in a common calculation framework. In both cases, we can find a well-defined contact plane tangential to both surfaces, and the contact region can be considered as a subset of this plane. Throughout this section, these two cases are considered together unless otherwise stated. When there are differences between the point contact and plane contact case, the analysis is shown parallel.

Resultant of the Friction Forces
Without loss of generality, we describe the relative motion of the bodies. That is, one of the bodies is considered as a fixed body, and the other one is referred to as a moving body.
In the contact plane (tangential plane) between the bodies, a planar contact region A is initiated. The motion of the moving body is described in a steady reference frame attached to the fixed body but in a rotating coordinate system attached to the moving body. The origin of the coordinate system is located at a chosen reference point O ∈ A. The coordinate axes are denoted by x, y, z where x, y lie in the contact plane and z points in the normal direction towards the moving body (see Fig. 1). The corresponding unit basis vectors are denoted by i, j, k, respectively. Then, the contact region is parametrized as a subset A = {(x, y)} ⊂ R 2 . During the motion, the coordinate system, the contact region A and all the contact forces may change in time. Time-dependence is not considered in this section (see, however, Sect. 5). Instead, we analyse a single time instant when the contact forces are calculated.
We assume that after the initiation of the contact region, the tangential motion of the points of the moving body inside the contact region is planar rigid body motion. That is, from the velocity v O = u x i + u y j of the reference point O and the normal component Left panel: Sketch of the contact region of two contacting bodies. One of the bodies is considered to be fixed. The coordinate direction z points into the normal direction towards the moving body. Right panel: The velocity components u x , u y of the reference point O and the angular velocity component ω z in the normal direction. With the assumption of planar rigid body motion in the contact plane, the velocity field can be calculated from these variables ω z of the angular velocity of the moving body, we can express the tangential velocity Equation (1) includes the assumption that either the whole contact region is sticking (when u x = u y = ω z = 0), or almost all points are slipping (when any components from u x , u y , ω z are nonzero). This assumption is widely used in the literature (see, e.g. Leine and Glocker 2003). To describe sticking and slipping zones in the contact region, a more detailed analysis of the deformations is necessary (see the details in, for example, Johnson 1985).
In this section, we assume that within the contact region, the normal pressure distribution is given, and it is denoted by p(x, y). Then, the total normal force P z k between the bodies is determined by (2) Moreover, the pressure distribution creates a resultant torque T x i + T y j about the point O, where the components are given by: We assume spatial Coulomb friction locally at each point of the contact region, and the friction coefficient μ is assumed to be the same for dynamic and static friction. Then, in the slipping case of the contact region, the distribution q xy (x, y) = q x (x, y)i + q y (x, y)j of the tangential contact forces is given by: Then, by substituting (1) into (4), the components of the tangential pressures become: Let us calculate the moment of the tangential force field with respect to the reference point O, where r = xi + yj and Finally, the components of the resultant force and torque of the tangential force distribution are given by: and the resultant of all contact forces with respect to the reference point O becomes: It can be seen from formulae (5) and (7), that, in general, the integrals in (8) cannot be computed analytically. Analytical solutions are still possible in some special cases; a few of them are shown in Sect. 4.5.

Analysis of the Discontinuity
Equation (1) shows that the velocity distribution in the contact region is determined by the variables u x , u y and ω z . These three variables are considered in the vector form where u z = ρω z is the normalized angular velocity. The multiplier ρ is defined by: i.e. ρ represents average distance from the reference point, or equivalently the characteristic length of the contact region. The rescaling of ω z by ρ makes the dimension of each component of (10) identical; furthermore, they also have a similar order of magnitude whenever u x , u y , and ω z are of the same order of magnitude [see the velocity reduction formula (1)].
In this paper, we use the boldface letters both for physical vectors such as P and for general vectors and matrices such as x o . For the shorthand notation, the n-tuple (x, y, z) is considered to be equivalent to a column matrix x y z . The notation . refers to the usual Euclidean 2-norm of vectors.
By the subscript 'o' in the notation x o , we anticipate the fact that in the full state space of a dynamical system, the three variables in (10) correspond to the orthogonal subspace of the discontinuity set (for details, see Sect. 3).
The norm of x o is: which expresses the magnitude of slipping. In the limit case r = 0, the velocity of all points in the contact region becomes zero, and the bodies are sticking to each other. Let us transform the variables (10) into spherical coordinates in the form where The vector w expresses the direction of slipping including information about the angular velocity, as well. Note that the variables w 1 , w 2 , w 3 are not independent because they have to satisfy the equation w = w 2 1 + w 2 2 + w 2 3 = 1. Thereby, x o is decomposed into the magnitude r ≥ 0 of slipping and the direction w of slipping in the space (u x , u y , ρω z ) (see the left panel of Fig. 2).
The direction vector w = (w 1 , w 2 , w 3 ) on the unit sphere S 2 ⊂ R 3 can be parametrized by two angles in the form: In (15), φ ∈ [0, 2π) denotes the "longitude" on the unit sphere, expressing the direction of the velocity of the reference point, and θ ∈ [−π/2, π/2] denotes the "latitude", expressing the ratio between ω z and the velocity magnitude of the reference point (see the right panel of Fig. 2). The parametrization is singular at the "poles" w = (0, 0, ±1), which states correspond to the pure rotating contact state with respect to the reference point (v O = 0). In most part of the paper, we use the parameter set w = (w 1 , w 2 , w 3 ), which is redundant due to w = w 2 1 + w 2 2 + w 2 3 = 1. Sometimes, it is useful to transform to the usual longitude φ ∈ [0, 2π) (expressing the direction of the translational motion) and latitude θ ∈ [−π/2, π/2] (expressing the ratio between the rotational and translational motion). Then, the "equator" of the sphere corresponds to the pure translational motion, while the north and south "poles" of the sphere corresponds to the pure rotational motion with respect to the reference point It can be seen from (5) and (7) that the tangential force distribution and its resultant become discontinuous at the sticking state x o = 0 (equivalent to r = 0). From (13), let us express u x = r w 1 , u y = r w 2 , ω z = r w 3 /ρ. By substituting these into (5)-(8), we get expressions, which depend formally on the direction w only, and they are independent of the magnitude r of slipping velocity: By using the vectorial form Q = (P x , P y , T z /ρ) ∈ R 3 , formulae (16) of the tangential forces can be written as Equation (17) maps the unit sphere S 2 onto R 3 , and the image E = Q(S 2 ) is the set of possible friction force resultants while slipping. Now, let us highlight some properties of (16)-(17), which are crucial for our analysis: • The integrand of (17) is equivalent of the local form of the Coulomb friction law (5), which has codimension-2 discontinuities at each point (x, y) ∈ A. The integration in (17) usually eliminates these (zero-measure) discontinuities of the integrand, and hence, (17) becomes smooth. (A similar phenomenon occurs in the following elementary integral: 1 0 (x − ξ)/|x − ξ |dξ = 2x for 0 ≤ x ≤ 1 where the integrand is discontinuous but the result is smooth.) Thus, by requiring some regularity of the pressure distribution p(x, y) and the geometry A of the contact region, the image E = Q(S 2 ) forms a closed, smooth 2-surface in R 3 .
• As a notable exception, pressure distributions concentrated along lines (see Sect. 4.5.4) may result in nonsmooth function Q(w), and pressure distributions concentrated in isolated points (see Goyal et al. 1991a;Hammond 2010;Borisov et al. 2016) may generate discontinuous Q(w). • Despite smoothness of Q(w), the friction force as a funtion of the velocity variables x o = (u x , u y , u z ) has a codimension-3 discontinuity at x o = 0. This is caused by simultaneous discontinuity of the local friction force at all points in the trivial (motionless) state. Thus, mapping (17) itself is not discontinuous, but the discontinuity is included into the variable w [see (14)]. • The state x o = 0 cannot be directly substituted into (16). The limit lim r 0 (r w) = x o = 0 gives the same (sticking) state independently from w; however, substituting this limit into (16) gives continuously many directional limits depending on w. The dynamics of vector fields in the vicinity of such discontinuities can be explored effectively by the mathematical tools developed recently in Varkonyi and Antali (2021), which we summarize briefly in Sect. 3.

Codimension-Discontinuities of Vector Fields
In a recent paper (Varkonyi and Antali 2021), the authors of the present work investigated the general properties of dynamical systems with an isolated codimension-n discontinuity. In this section, we review the main results related to the n = 3 case, which results are applied to slipping and spinning bodies in Sect. 5.

Codimension-3 Discontinuities
Consider the space R m where the vectors are denoted by x = (x 1 , x 2 . . . x m ) ∈ R m , and the following codimension-3 subspace of R m : Consider the differential equationẋ where F is an R m \ → R m vector field with components F = (F 1 , F 2 , . . . F m ). Dot means differentiation with respect to time t; the explicit time dependence x(t) is not denoted except when it is necessary. We assume that the vector field F is undefined in , and it has the following properties: 1 F is smooth everywhere in R m \ . 2 In all pointsx ∈ , the limit exists for any vector v ∈ R m orthogonal to . 3 The limit (20) depends smoothly on v andx.
In this context, we call a codimension-3 discontinuity set with respect to F. For a selected pointx ∈ , we call F (x, v) the limit vector field of F atx. Note that the linear spaces R m and can be generalized to smooth manifolds if necessary.
In the case of contacting bodies, F describes slipping motion, and the discontinuity set corresponds to rolling or sticking states. The dynamics of F can be extended to the discontinuity set , either by using the dynamical equation of sticking (or rolling) or by calculating the sliding dynamics (in the sense of Filippov theory) generated by F on (see Varkonyi and Antali 2021).
Our goal is to investigate the behaviour of the differential equation (19) in a small neighbourhood of with a focus on trajectories starting or ending at a pointx ∈ . As we will see, such trajectories correspond to slip-stick and stick-slip transitions of the bodies in the mechanical model.

Reformulation of the Vector Field
To achieve an appropriate form of the equations, we reformulate the vector field in four steps.

Decomposition into Orthogonal and Tangential Parts
Let us separate the state variable x into parts orthogonal and tangential to , where x o ∈ R 3 is the orthogonal part and x t ∈ R m−3 the tangential part. Then, the vector field can be written in the form

Spherical Decomposition of the Orthogonal Part
Next, a spherical decomposition of the orthogonal variable x o is considered: where r = x o represents the distance of x from , and is the unit vector showing the direction of x around . This set of variables is redundant and subject to the constraint w = 1. That is, w is located on a unit sphere. Note that in (10)- (14), the same quantities were defined from mechanical point of view.

Asymptotic Dynamics
Note that due to the requirement of the smooth dependence of (20) on the variables, the functions R, V and F t depend smoothly on the variables r , w, x t . Furthermore, they have well-defined limits which are smooth functions in w and x t . Then, we can define the asymptotic dynam- It was demonstrated in Varkonyi and Antali (2021) that in some important cases, the local dynamics (23) in an infinitesimal neighbourhood of a pointx t ∈ can be described via investigation of the asymptotic dynamics (25).

Time Rescaling at the Singularity
The singularity associated with r = 0 is still present in the second equation of (25). This singularity can be removed by a singular rescaling of time, rendering (25) into where dash means derivation with respect to the transformed time variable τ . The new set of equations is a multiple time scale dynamical system for small r , where w is a fast variable, and r and x t are slow variables. Note that (27) does not fit into the standard form (see Kuehn 2015, p. 8), of multiple time scale systems. Standard slow-fast systems are written as: where y contains the fast variables, z contains the slow variables, and 0 < 1 is the so-called time-scale parameter. In (27), x t is a usual slow variable, but the slow variables r behave differently: We cannot see a time-scale parameter in (27), but the separation of the time scales is caused by variable r itself when 0 < r 1. Equation (27) appears to be outside the scope of (28); nevertheless, it has been found in Varkonyi and Antali (2021) that the overall behaviour of the system can be investigated via separation of fast and slow variables, similarly to the standard case.

Analysis of the Asymptotic Dynamics
In (27), the fast dynamics of w fully decouples from the slow variables and the equation w = V (w,x t ) can be solved independently. According to the Poincaré-Bendixson theorem, trajectories of smooth dynamical systems over S 2 converge to fixed points, limit cycles, or polycycles (Lopez and Llibre 2007). By investigating a large number of cases, we found that trajectories corresponding to dynamics of slipping bodies always converge to fixed points, albeit this property is not proven here. In what follows, we only review results related to fast dynamics converging to fixed points. We begin by introducing the definition of limit directions: Definition 1 Consider a fixed pointŵ of the circumferential dynamics satisfying V (ŵ,x t ) = 0. Then, we callŵ a limit direction of the system (19) atx t . In the cases R (ŵ,x t ) < 0 or R (ŵ,x t ) > 0,ŵ is called an attracting or repelling limit direction, respectively.
The names defined above are inspired by the following pair of statements:

Theorem 1 (from Varkonyi and Antali 2021) Consider an attracting [repelling] limit
Theorem 1 states the existence of a single trivial trajectory at the fixed pointŵ of the circumferential dynamics. The limit directions can be further classified based on the number and type of trajectories satisfying (29).
For this, let us linearize the circumferential dynamics at the limit directionŵ: where The second term of (31) makes the linear term in (30) a projection on the unit sphere of w; we can check w V w w = 0. (Note that alternatively for this analysis, the unit sphere w = 1 could be considered as a smooth two-manifold.) Thus, in general, the 3 × 3 matrix V w has two nontrivial eigenvalues λ 1 and λ 2 and a trivial eigenvalue λ 3 = 0, the latter correspond to the radial dynamics parallel to w. We callŵ a hyperbolic fixed point if Reλ 1 = 0 and Reλ 2 = 0. In this case, we can define the following types of limit directions: Definition 2 Consider a limit directionŵ atx = (0, 0,x t ) ∈ . Assume thatŵ is a hyperbolic fixed point of (25) with corresponding nontrivial eigenvalues λ 1 and λ 2 . The limit direction is called In other words, a stable node or focus of w = V (w,x t ) corresponds to an attracting-dominant or a repelling-isolated limit direction; an unstable node or focus corresponds to an attracting-isolated or repelling-dominant limit direction, and a saddle corresponds to a saddle-type limit direction. This categorization shows whether a limit direction attracts the nearby trajectories in the time direction of approaching the discontinuity set. At a dominant limit direction, the trivial trajectory from Theorem 1 has a neighbourhood where all trajectories are connected to along the limit direction. At an isolated limit direction, there is a single isolated trajectory that is connected to . At a saddle-type limit direction, there is a mixed behaviour according to the stable and unstable directions of the saddle.

Ellipsoidal Approximation of the Slipping Friction Force Resultants
It can be seen from the integrals (16) that generally the expressions of the resultants of the friction forces cannot be computed analytically. Analytical solutions are still possible in some special cases; a few of them are shown as reference results in Sect. 4.5. For many purposes of dynamical analysis, it could be appropriate to calculate (8) numerically for the desired accuracy. Instead, in this section, we present a natural leading-order approximate model of (16), which makes it possible to apply the nonsmooth mathematical tools of Sect. 3 analytically in order to explore the different dynamical scenarios of nonsmooth behaviour of the contacting bodies (see Sect. 5).

Form of the Approximation
Let us approximate the closed surface E = Q(S 2 ) to replace the exact integration of (16). A natural approximation of a closed surface is an ellipsoid. IfẼ =Q(S 2 ) ≈ Q(S 2 ) is an ellipsoid, thenQ must be simply a linear function in w in the form whereˆ is a 3 × 3 matrix. The ellipsoidal approximation (35) can be considered as the leading term of a harmonic expansion of the function Q(w) on the unit sphere, which is the higherdimensional extension of the harmonic Fourier expansion on a unit circle. Such expansion would be achieved by using the Laplace spherical harmonics (see Brechbühler et al. 1995 andArfken et al. 2013), which leads to higher-order terms of w in (35). (Note that this is not a Taylor expansion because w is located on the unit sphere, and it cannot have the value 0.) In this paper, we focus on the first-order approximation (35) by an ellipsoid.
Note than in Howe and Cutkosky (1996) and (Arkadani et al. 2020), a similar ellipsoidal approximation of the surface E is proposed. Those authors take a different approach: E is generated as the boundary of the set of admissible static contact forces, and it is called limit surface (Goyal et al. 1991a). The relation between Q and the corresponding slipping direction w is established via an orthogonality assumption Q (w)w = 0. Those authors also approximate E by an ellipsoid; nevertheless, their orthogonality condition is not consistent with (35), and hence, our model and the model of Howe and Cutkosky (1996) and Arkadani et al. (2020) are not equivalent.

Fitting the Ellipsoid to Special Motions
There are several possibilities to fit the ellipsoidal approximation (35) onto the exact results (16)-(17). One possibility is to perform a proper harmonic expansion on the sphere and take the leading term. A second possibility would be to calculate (16) numerically at a sufficient number of points w ∈ S 2 , and the matrixˆ would be determined by, for example, the method of least squares. We follow a third way: we select those three special motions for which only one component of w is nonzero and chooseˆ to give exact results for these motions.
Let us first define the following integrals, which are used in the expressions below: Note the connection with (3). The elements ofˆ can be found by calculating the images of three linearly independent elements in R 3 corresponding to the special motions. In these cases, (16) leads to simple integrals expressed in terms of the quantities (36): • When u y = ω z = 0, the reference point is moving into the x direction with a velocity u x and the body is not rotating in the normal direction z. Then, we get By substituting these into (16), we get • When u x = ω z = 0, the reference point is moving into the y direction with a velocity u y and the body is not rotating in the normal direction z. Then, we get x o = (0, u y , 0), which leads to x o = |u y |, w 2 = u y /|u y |, and w = (0, w 2 , 0). By substituting these into (16), we get • When u x = u y = 0, the reference point is not moving but the body is rotating with an angular velocity ω z . Then, we get x o = (0, 0, ρω z ), which leads to x o = ρ|ω z |, w 3 = ω z /|ω z |, and w = (0, 0, w 3 ). By substituting these into (16), we get By comparing (35) with (37)-(39), the matrixˆ becomes: where =ˆ /(−μP z ). That is, we determined the matrixˆ in such a way that the ellipsoidal approximate model (35) provides an accurate result for the integrals (16) in the case of special motions described above.

The Proposed General Ellipsoidal Model
From (35), (40) and (36), the ellipsoidal approximation of the set of friction forces leads to the model The components of (41) are given by: If the region A and the pressure distribution p(x, y) do not change in time with respect to the co-moving coordinate system, then the elements of are parameters that can be calculated from different moments and averages (36) of the normal pressure distribution p(x, y). When the pressure distribution p(x, y) depends on other state variables of the system, then the elements of will change in time (see Sect. 5).
The resultant effect of the friction forces (16)  The components (42) are discontinuous at u x = u y = ω z . That is, the application of this model in mechanical systems creates a codimension-3 discontinuity in the vector field (see Sect. 3).

Special Cases
Let us consider some important special cases of the ellipsoidal model: • If the region A is symmetric with respect to the reference point and p(−x, −y) = p(x, y), then, according to (36), π x = π y = π x0 = π y0 = 0. Thus, becomes diagonal and (42) leads to • The single-point contact is obtained as a limit case when A is very small; then, ρ → 0 and π xy → 0. As the tiny contact region can be considered symmetric, (43) simplifies to This contact model excludes the effect of drilling friction, and it leads to a codimension-2 discontinuity analysed in Stepan (2018, 2019). • When the pressure distribution p(x, y) does not change in time, it is useful to choose the reference point O and the coordinate directions x and y appropriately to make the form of as simple as possible. Assume now an initial choice of O for which has the general form (40). By a translation of O in the contact region, we can find a unique point when the resultant torques (3) vanish, and then, we get π x = π y = 0 [see (36)]. Then, by an appropriate rotation of the coordinate axes x and y, we can eliminate π x0 (or, alternatively π y0 ). Finally, we get the model in the form

Case Studies
In this subsection, we calculate the ellipsoidal approximate model for some simple geometries and pressure distributions, where we can find analytical reference solutions for comparison.

The General Axisymmetric Case
Consider a circular contact region with a radius R with an axisymmetric pressure distribution p(x, y) =p( x 2 + y 2 ). Assume that the reference point O is located in the centre of the circle. It can be shown from axisymmetric properties of the contact region and pressure distribution that in this special case, (16) can be written in the form where P and T are 2π periodic functions; P w and T w are scalar functions over the respective intervals (0, 1) and (−1, 1). From (11) and (36), the parameters ρ and π xy are given by Note that from transformation (15), the functions defined in (46) satisfy In the literature, a further parametrization can be found (see, e.g., Leine and Glocker 2003;Farkas et al. 2003;Weidman and Malhotra 2005). Let be a dimensionless ratio between the variables defined by and consider the re-parametrized function P ( ) and T ( ) defined by
In Farkas et al. (2003), Weidman and Malhotra (2005), we can find the analytical formulae of (46) for the uniform pressure distribution. According to Weidman and Malhotra (2005), the functions P ( ) are given by: where E(m) and K (m) are the incomplete elliptic integrals expressed by the elliptic parameter m (see Abramowitz andStegun 1972, p. 589 or Oldham et al. 2009, p. 637). (Note that the elliptic integrals can be also defined by elliptic modulus √ m in the formĒ( √ m) = E(m) Abramowitz and Stegun 1972;Oldham et al. 2009, which parametrization in used, e.g., in Farkas et al. (2003.) The sections of the piecewise defined functions (52)-(53) are smoothly connected at = 1, and each can be expressed by a single formula containing incomplete elliptic integrals (Weidman and Malhotra 2005).
Having this reference analytical solution, let us apply the approximate model in (41). From the symmetry of the problem with respect to the reference point, the conditions of (43) are satisfied. By comparing (46), (43), and (51), we get that the ellipsoidal approximation leads to or, with the two types of parametrization from (46), As ρ = π xy , the matrix from (40) becomes a unit matrix.
The exact solution (52)-(53) and its ellipsoidal approximation (54)-(56) can be seen in Figs. 3 and 4. On the diagrams, we observe that the ellipsoidal model provides an acceptable approximation with a very simple algebraic formula, which is useful in dynamical calculations, for example, presented in Sect. 5. Fig. 3 The dimensionless friction force (P) and friction torque (T ) of a circular contact. On the three diagrams, the same quantities are shown by using different parametrizations (see (49) and (51) Fig. 3. Due to the fitting method, the ellipsoidal approximation is accurate at pure translational and rotational motions with respect to the reference point. Left panel: the dimensionless friction force and torque values. We can see that these graphs for the constant and parabolic distributions are almost identical. Right panel: graph the friction force and torque components of the constant distribution, which we can call friction ball (Leine and Glocker 2003). The friction balls (exact and approximated) for the parabolic distribution have a very similar shape, but the axis of the ellipsoid in the direction of T z is somewhat smaller (the value of π xy is different)

Circular Contact with Parabolic Pressure Distribution
Consider the parabolic pressure distributioñ for which (48) leads to π xy = 3πr /16. Based on Leine and Glocker (2003), the analytical solution to the problem is When calculating the approximate solution (41), we get the same results as (54)-(56). It is not surprising as the exact dimensionless functions for constant and parabolic pressure distributions are close to each other (see Fig. 3).

Long, Thin Rectangular Contact
Consider a thin rectangular contact region A = [−L/2, L/2] × [−δ, δ] with δ/L 1 with a constant pressure distribution p(x, y) = P z /(δL) (see Fig. 6). In the limit case δ/L → 0, the integration in (16) leads to a single-variable integral. Then, the friction forces in (16) become The integrals (60) have lengthy closed forms containing square root and area hyperbolic sine functions. As we did in (46), let us define the dimensionless friction force P and friction torque T , which quantities now depend on the angle φ due to the anisotropic geometry of the contact region [see (14)-(15)]. In the two special cases when the velocity of the reference point is parallel (φ = 0) or perpendicular (φ = π/2) to the x axis, we get simple analytical results from (60)-(61): where γ = 2 tan θ . Let us calculate the ellipsoidal approximation according to (41). From (11), we get and because of the uniform pressure distribution, (36) gives the same value for π xy = L/4. Then, matrix in (40) becomes a unit matrix, and the approximate model is given by (43) with the substitution ρ = π xy = L/4. The results can be seen in Fig. 5. The exact friction law is anisotropic, and thus, the functions P and T depend on the angle φ, as well. However, the ellipsoidal model behaves isotropically here and gives the same approximation for each direction. The rectangle is assumed to be very thin (δ L), which tends to the case of line contact. Even in this singular case, the approximation is acceptable at φ = 0; the error becomes significant at φ = π/2.

Long, Thin Rectangular Contact-Change of the Reference Point
Let us check the effect of the choice of the reference point O on the previous example. Let us consider the new reference point O at x = L/8 (see Fig. 6). The corresponding coordinate axes are x and y , and the quantities with respect to O are denoted by Q , w and . Then, the form of the ellipsoidal model becomes By using the reduction formula of velocities and moments of rigid bodies, the relationship between the two sets of variables becomes Thus, although we transformed back the ellipsoidal model to the reference point O, the transformed model W W still slightly deviates from the unit matrix obtained in the previous case. This example demonstrates that the choice of the reference point modifies the fitting of the ellipsoidal approximate model.

Circular Contact with Linear Pressure Distribution
Consider now a circular contact region of radius R with a linear pressure distribution in the form p(x, y) = p 0 + c x x + c y y.
Then, from (2), we get Compared to the previous examples, this pressure distribution has fewer symmetries, making both the exact and the approximate solution more irregular.
From the resulting friction force components (16), let us define the dimensionless friction components P x , P y , T z by which are now available only from numerical computation. In the ellipsoidal model, the integrals (36) can be calculated by using polar coordinates, and the matrix from (40) becomes Now, the matrix becomes nondiagonal; that is, the friction force characteristic is anisotropic.
The comparison of the ellipsoidal model and the numerically computed exact solution can be seen in Fig. 7. We can see that in several different cases, the ellipsoidal model provides an approximation that captures the qualitative behaviour of the friction forces.

Dynamic Equations of a Rigid Body in Conforming Contact with a Plane
Now, we demonstrate the application of the mathematical tools of Sect. 3 to the approximate model presented in Sect. 4. In this section, we consider the case of plane contact from Sect. 2 when one of the objects is a fixed plane, and a flat face of a moving body is in conformal contact with this plane. In this section, the word body refers to the moving body.

Notation
Recall that the angular velocity of the body is ω = ω z k and the velocity of the reference point is v o = u x i + u y j. Thus, the velocity state of the body is fully determined by the three variables u x , u y and ω z .
Assume that the planar frictional contact is the only nonsmooth effect in the dynamics of the rigid body. Then, we expect that the state space of the body has the form (21), where x o = (u x , u y , ρω z ) [see (10)], and x t contains the other state variables, which has smooth effect on the dynamics. In the simplest case, x t contains the three generalized coordinates describing the configuration of the body.
Our goal is to obtain a set of first-order differential equations for x o in the forṁ which is the basis of the analysis presented in Sect. 3. Then, it is useful to transform kinematic quantities with respect to the centre of gravity C of the body. Let the location of point C be r c = x C i + y C j + z C k, where x C , y C , z C are parameters because O, C, and the coordinate system are all fixed to the body. The velocity of C is given by: Then, the velocity state of the rigid body can also be determined by From (75), the relation between (76) and (10) becomes where the transformation matrix W is given by It is assumed that an external force system is acting on the body whose resultant in O is described by the concentrated force and torque P e = (P ex , P ey , P ez ), respectively, which are assumed to depend smoothly on the variables x o and x t . In addition to that, the body is subject to normal and frictional forces as described in Sect. 2. The resultant force and torque of the friction forces at the reference point O are given by P = (P x , P y , P z ) and T = (T x , T y , T z ), respectively.

Newton-Euler Equations of Motion
The angular momentum of the body with respect to the centre of gravity C can be written in the form: where J xz , J yz and J z are the corresponding components of the mass moment inertia matrix of the body with respect to the point C. Note that the axes x, y, z are fixed to the body. The Newton-Euler equations of the body become mv c = P + P e ,L = T + T e − r c × (P + P e ), where m is the mass of the body. The six scalar equations corresponding to (81) are where the terms ω z v cy , ω z v cy , J yz ω 2 z , J xz ω 2 z account for the rotating coordinate axes i, j having the time derivativesi = ω × i andj = ω × j.

Transformation of Variables
In matrix form, (82), (83), and (87) can be written as: where and W is given by (78), and Q = (P x , P y , T z /ρ) was defined in (17). Now, we express the equation (88) Let κ denote the radius of gyration defined by J z = mκ 2 . Then, we can define the mass matrix and the gyroscopic matrix whereby (84) becomes However, (93) is still not fully defined because we have to calculate the vector Q of friction force components.

Calculation of the Contact Forces
In the ellipsoidal model (41), we proposed the approximation Q = −μP z w where w = x o / x o . For the calculation of the elements of the matrix , we used the integrals (36) containing the pressure distribution p(x, y). In Sect. 2, the pressure distribution p(x, y) was declared to be known. However, in the case of conformal surface contact between quasi-rigid bodies, p(x, y) is determined by microscopic deformations of the objects involved in the contact; moreover, it depends extremely sensitively on geometric imperfections of the surfaces (Huang et al. 2017). Hence, we need to find p(x, y) and the acceleration of the moving body simultaneously.
Obviously, this problem is not tractable within the context of rigid body dynamics. One possible approach is considering detailed models of the material behaviour and geometric imperfection. To avoid the complexity of this approach, we make additional assumptions about p.
According to (2)-(3), p(x, y) has to satisfy the conditions where P z , T x , T y are determined by constraints of the motion. In particular, P z is given by (84) as and the torque components can be obtained from combining (85)- (87): where and j yz = J yz /J z , j xz = J xz /J z are the dimensionless nondiagonal components of the inertia tensor. For now, we restrict the analysis to the limit case N = 0 of flat bodies when the point C lies in the contact region (z c = 0) and the nondiagonal components of the inertia tensor vanish ( j xy = j yz = 0).
If N = 0, then T x , T y can be directly expressed from (96), and thus, the conditions (94) allow to identify three parameters of p(x, y). For example, one may assume that p is linear as in (69). Then, the parameters p 0 , c x , c y are given by (94), and thus, p(x, y) becomes completely determined. Finally, we can calculate the elements of from (36), and model (41) provides the friction forces for the differential equation for (93).
In the more general case with N = 0, further complications occur because (96) already contains and w. This situation is analysed later in Sect. 5.6.

Analysis of the Discontinuous Dynamics
By considering the friction model Q = −μP w from (41), the time derivativeẋ o can be expressed from (93) in the form: where and Based on the previous assumptions, B and c depend smoothly on x = (x o , x t ). Based on the decomposition (21), Eq. (98) is complemented with the dynamics of the smooth variables in the formẋ t = F t (x). The resulting vector field F = (F o , F t ) gives a differential equation of form (19), and thus, the analysis and results from Sect. 3 are applicable to the system.

Radial and Circumferential Dynamics
The discontinuity of the system is located at x o = 0, that is, x o is the orthogonal subspace of variables according to the separation in (21). Consider the spherical decomposition x = (r , w, x t ) from (22) with x o = r w. By using (23), the orthogonal dynamics (98) can be decomposed to the radial dynamicṡ and the circumferential dynamics where I is the identity matrix.

Limit Directions
Let us choose a pointx = (0, w, x t ) of the discontinuity set , and letB = B(x), c = c(x). Then, according to (25), the local approximation of the dynamics at x =x becomesṙ It can be seen from (104) that w evolves on a fast time scale in the vicinity r = 0 of the discontinuity.
At r = 0, the fixed points of (104) are determined by V (w,x t ) = 0, which are the limit directions of the dynamics at x t (see Definition 1). To obtain information about these fixed points, considerẋ o =(r w) =ṙ w + rẇ. (105) By comparing (98) and (105) in the case x =x andẇ = 0, we get The solution of equation (106) are pairs (w i ,ṙ i ), where w i . . . w k are the limit directions andṙ i . . .ṙ k are the rate of evolution of the radial dynamics along the limit directions. From (106), we can draw consequences about the number of the limit directions.

Limit Directions Without External Accelerating Forces
Consider first the special case ofc = 0. As the gyroscopic matrix G vanishes at r = 0 (ω z = 0), the casec = 0 is equivalent to Q e = 0, that is, there are no "accelerating components" of the external forces. The other external force components P ez , T ex , T ey can be nonzero, which affect the dynamics through modifying the pressure distribution. Ifc = 0, then (106) becomes which is the eigenvalue-eigenvector problem ofB leading to the following statement: Proposition 1 Let λ ∈ R be a real eigenvalue ofB and let γ be an associated eigenvector satisfying λI −B γ = 0. Then, the eigenvector determines w = ± γ −1 γ , which are two fixed points of (104) in the special casec = 0. That is, in the number of limit directions is twice the number of the real eigenvectors ofB.
Although the other two eigenvalues ofB can be expressed explicitly, they are lengthy algebraic expressions, and it cannot be checked directly whether they are real or not. However, from the physical behaviour of the friction law, we expect that all three eigenvalues ofB are real. Therefore, our conjecture is that in general, the body without external forces Q e has six limit directions.
Let us consider a singular case that was excluded from Proposition 2.

Limit Directions with External Forces
In the general case when c = 0, (106) is a so-called inhomogeneous eigenvalue problem (Mattheij and Soderlind 1987) in the form Assume thatṙ is not an eigenvalue ofB. Then, from (108), we can express and we get the equation From (109)- (110), we can reach some important conclusions:

Proposition 4
The system (104) has a maximum of six fixed points.
Proof Equation (110) be rewritten into a 6 th degree polynomial inṙ . (Note that the inverse of a matrix can be expressed by the adjugate matrix and the determinant.) Consequently, it has maximum 6 real roots forṙ , and then, and the related fixed points w are determined by (109).

Proposition 5
The system (104) has minimum two fixed points.
Note that in (109)-(110), we explicitly assumed thatṙ cannot be an eigenvalue ofB. We can use the following result from (Mattheij and Soderlind 1987), which includes these degenerate cases, as well.
Proposition 6 Consider the inhomogeneous eigenvalue problem (108) whereB has at least one real eigenvalue. For the number k of the real solutionsṙ ∈ R and w ∈ R 3 (multiplicity counted), we have k = 2, k = 4 or k = 6.
It can also be found in Mattheij and Soderlind (1987) that for the real solutions, (108) can be transformed into which is a usual eigenvalue-eigenvector problem in R 6 for the eigenvaluesṙ and eigenvectors (w,ṙ w).

Case of Bodies with Significant Height
In Sect. 5.5, we assumed that the matrix N vanishes in (96), which occurs when the height of the centre of gravity is negligible (z c = 0), and z is a principal axis of the inertia tensor ( j xz = j yz = 0). If these conditions are not satisfied, then N = 0. Then, Eq. (96) cannot be used directly to express T x and T y because the matrix depends on the pressure distribution p(x, y) , which should be calculated from T x and T y . Therefore, there is feedback from the dynamical equations, which makes the previous direct calculation impossible.
However, the problem can be resolved if we consider the assumed form (69) of the pressure distribution containing the parameters. In this case or for any other pressure distribution containing three parameters as linear coefficients, we can choose these three parameters as unknowns. Let us substitute (69) into the integrals in (36) and (94). Then, the parts of the integrals P z , T x , T z , and can be computed purely from the geometry of the contact region A, and the parameters p 0 , c x , c y appear as unknown coefficients. Then, by substituting P z , T x , T z , and into (95)-(96), we get a set of three linear equations, which can be solved for the parameters p 0 , c x , c y .
Moreover, this solution depends on the slipping direction w, as well. Thus, we get solutions in the form p 0 (w), c x (w), c y (w). Consequently, this case leads to a friction model in the form This expression introduces nonlinearity in w into the dynamics (98). Then, the results in Sect. 5.5 cannot be applied directly. However, we demonstrate in the example in Sect. 6.2 that several details can be explored about the limit directions even in this case.

Examples
In the previous section, we surveyed several cases of finding the limit directions of a moving body in planar contact with a plane. These limit directions express the specific directions where transitions can occur between the slipping and sticking states of the moving body. In what follows, we demonstrate the results on mechanical examples.

Thin Disk on a Flat Surface with an In-plane Loading Force
Consider a flat, circular disk of radius R with mass m and homogeneous mass density, implying radius of gyration κ = R/ √ 2, and J xz = J yz = 0. The disk moves on a horizontal plane, and the (vertical) gravitational acceleration is g (see the left panel of Fig. 8). The height of the disk is assumed to be much smaller than the radius, and the reference point is the midpoint of the base facet, that is, r c = 0. An in-plane constant loading force is acting at the centre of gravity C of the disk; this force has a constant value F and a constant direction of the fixed plane. Then, the resultant of the external forces is P e = (F cos ψ, F sin ψ, −mg), T e = 0, where ψ denotes the orientation of the disk on the plane withψ = ω z . (Note that the coordinate system is fixed to the disk.) Then, (89)-(92) imply From (84)-(86), the resultant of the normal contact forces becomes P z = mg, T x = T y = 0. We expect that the pressure distribution is constant with p(x, y) = p 0 = mg/(R 2 π). [This can be obtained by direct calculation of (69), as well.] From Sect. 4.5.2, the uniform pressure distribution implies that is a unit matrix.
In addition to the orthogonal variables x o = (u x , u y , ρω z ), the only other state variable in the expressions is the orientation ψ of the disk, that is, x t = (ψ). As the problem is symmetric with respect to the rotation about the z axis, the dynamics is practically identical to any ψ ∈ [0, 2π). Without loss of generality, we can take ψ = 0 andx t = (0). Then, direct calculation of (98)-(105) yields: By substituting (114) into (103)-(104), the radial and circumferential dynamics becomes Fig. 9 Sketch of the fast dynamics and limit directions on a unit sphere in the case of the thin disk with an in-plane loading force. Encircled numbers correspond to the three cases of behaviour identified in Sect. 6.1. The colour codes are the following: red and blue dots correspond to repelling and attracting behaviour, respectively. Red and blue arrows correspond to unstable and stable behaviour on the unit sphere, respectively. A limit direction is dominant if the colours of the dot and the corresponding arrows match, and otherwise, it is isolated (Color figure online) respectively. The condition (106) of limit directions is recast as three scalar equations (ṙ + μg)w 1 = F/m, (ṙ + μg)w 2 = 0, (ṙ + 8/9 · μg)w 3 = 0.

Fig. 10
Sketch of the fast dynamics and limit directions on the invariant circle in the case of the thin disk with an in-plane loading force. The invariant circle depicted to these sub-figures are denoted by green colour in Fig. 9. The circled numbers correspond to the three cases of behaviour identified in Sect. 6.1 The trivial limit directionsŵ 1 andŵ 2 always exist, while the nontrivialŵ 3 and w 4 exist for F < μmg/9 only. Note that all limit directions are located on the circle w 2 = 0 on the unit sphere (see Figs. 9,10). Moreover, it can be seen from (116) that the circle w 2 = 0 is an invariant set of the fast dynamics. Physically,ŵ 1 and w 2 correspond to pure translational motions in the same direction (case ofŵ 1 ) or the opposite direction of F (case ofŵ 2 ). In contrast,ŵ 3 andŵ 4 are related to motions with a certain fixed ratio between the velocity u x and the angular velocity ω z . Consider the loading force F ∈ R as a bifurcation parameter of the system. Then, according to the number and type of the limit directions, we can identify three structurally different regions: Case 1 If |F| > μmg, then there are two limit directions:ŵ 1 is repelling-isolated, w 2 is attracting-isolated. As the only attracting limit directionŵ 2 is isolated, there is a single trajectory characterized by u x < 0 and u y = ω z = 0, which reaches the discontinuity (the sticking state). Even in this case, the disk sticks for a moment and starts slipping immediately, prescribed by the repelling directionŵ 1 . All the other trajectories avoid the discontinuity, and thus, the sticking of the disk is practically unrealizable. The condition |F| > μmg of impossible sticking coincides with the trivial condition from elementary dynamics where only the velocity component u x is considered without the transverse velocity or rotation. This case is depicted in the left sub-figures of Figs. 9 and 10.
Case 2 If μmg > |F| > μmg/9, then there are two limit directions,ŵ 1 attractingdominant, andŵ 2 is attracting-isolated. As both limit directions are attracting, all nearby trajectories end at the discontinuity, and thus, sticking is a realizable longterm behaviour of the disk. A single trajectory exists-the same trajectory mentioned in the previous cases-which reaches the sticking state from the isolated direction w 2 . However, asŵ 1 is a dominant limit direction, almost all trajectories approach the sticking state from this direction. Physically, it means that for almost all initial conditions of slipping, the disk reaches the sticking state in finite time, and just before sticking, the disk has a pure translational slipping with a positive velocity u x . This case is depicted in the middle sub-figures of Figs. 9 and 10.
Case 3 If μmg/9 > |F| > 0, then there are four limit directions:ŵ 1 is attracting saddle type,ŵ 2 is attracting-isolated andŵ 3 andŵ 4 are both attracting-dominant. Fig. 11 Bifurcation diagram of the limit directions when the loading force F is varied. The blue and red lines denote the attracting and repelling limit directions, respectively. The solid lines denote the dominant limit directions, while the dashed lines correspond to the isolated or saddle-type limit directions. The circled numbers of the regions correspond to Cases 1-3 identified in Sect. 6.3, and the corresponding phase portraits can be found in Figs. 9 and 10 (Color figure online) The main difference from the previous case is that now, almost all trajectories tend to the discontinuity alongŵ 3 andŵ 4 . On the unit sphere of the fast dynamics, the basins of attraction of these fixed points are separated by the "equator" w 3 = 0 of the unit sphere, which circle is an invariant set of the fast dynamics [see (116)]. Trajectories starting from this plane-without rotation-typically tend to the discontinuity along the saddle-type directionŵ 1 , except for the single trajectory ofŵ 2 mentioned before. However, the typical behaviour of the disk is reaching the sticking state along the nontrivial equilibriaŵ 3 andŵ 4 . That is, just before coming to rest, the motion of the disk is a certain combination of simultaneous translational and rotational motion, where the ratio of the velocity u x and the angular velocity ω z is prescribed by the limit directions in (119). This case is depicted in the right sub-figures of Figs. 9 and 10.
In Magyari and Weidman (2001), a thin disk on a slope was analysed, which is analogous to the present example. The authors of Magyari and Weidman (2001) used a different approach analysing a reduced phase space implicitly assuming a zero transversal velocity u y , and the exact solution of the tangential force system was considered (see also Farkas et al. 2003;Weidman and Malhotra 2005). Our results about Cases 2 and 3 are consistent with the results of Magyari and Weidman (2001).
Let us take a look at the bifurcations occurring between the cases presented above (see Fig. 11). At |F| = μmg/9, the trivial limit directionŵ 1 interacts with the nontrivial onesŵ 3 andŵ 4 in a way that we can call a pitchfork bifurcation of limit directions. At |F| = μmg, neither the number, nor the fast eigenvalues of the limit direction change, butŵ 1 turns from attracting into repelling.
A special degeneracy occurs at F = 0, which is the force-free motion of the disk on the plane. Then, (117) is satisfied by a continuous set of limit directionŝ w 0 (ζ ) = (ζ, ± 1 − ζ 2 , 0) for any −1 ≤ ζ ≤ 1 as well as the pairŵ 3,4 = (0, 0, ±1). The limit directionsŵ 3,4 = (0, 0, ±1) are located at the "north and south poles" of the unit sphere, corresponding to a pure rotational motion. These discrete limit directions are attracting-dominant, while the directions of the continuous setŵ 0 (ζ ) are on the limit state between attracting-saddle and attracting-isolated. It can be shown that for initial conditions without rotation (ω z = 0), the disk reaches the sticking state along one of these degenerate directions. Therefore, the model predicts that the typical behaviour of the force-free disc just before stopping is a pure rotational motion. However, this degenerate case is analysed in Farkas et al. (2003) with the exact contact forces, and it was found that the typical behaviour should be a combined motion with translational and rotational slipping (there, theŵ 3,4 = (0, 0, ±1) would not reach the "poles" at F = 0). This shows that for low values of F, higher-order terms of w, neglected in the ellipsoidal model (41), have a significant effect in this example.

Thick Disk on a Flat Surface with an In-plane Loading Force
Consider now an extension of the previous example, where the disk has a height 2h comparable to its radius R (see the middle panel of Fig. 8). Assume that the loading force acting at the centre of gravity, which is at the height z c = h from the plane (uniform mass distribution is assumed). Then, r c = (0, 0, h), and the external force system becomes P e = (F cos ψ, F sin ψ, −mg), T e = Fh(− sin ψ, cos ψ, 0), where ψ is the orientation angle of the disk, again.
Expression (126) is a third-degree polynomial in w 2 3 , and we search for the nontrivial roots 0 < w 2 3 < 1. We can investigate the roots by using the discriminant and Descartes' Rule of Signs. Straightforward but lengthy calculations lead to the following result: Consider the critical value of the loading force. In the case |F| > F 2/4 crit , the polynomial (126) does not have a positive root. For 0 < |F| < F 2/4 crit (126) has a single positive root in the interval 0 < w 2 3 < 1. Each root corresponds to two limit directions in the form ± w 2 3 . Finally, the two types of limit directions are two trivial limit directions w =ŵ 1,2 = (±1, 0, 0) and two nontrivial limit directions nontrivial limit directions w =ŵ 3,4 determined by (125). After checking the eigenvalues of the fast dynamics, we get similar cases of behaviour to those of the previous example: Case 1 If |F| > μmg, then there are two limit directions:ŵ 1 is repelling-isolated, w 2 is attracting-isolated.
Thus, the dynamics of a thick disk is structurally similar to that of the thin disk. However, there are significant differences in Case 3 with relevant physical consequences: • In the limit case h/R → 0, formula (127) of the critical loading force gives the value previously determined to the thin disk. However, (127) shows a significant increase in the parameter value when we increase the height h. Thus, the thicker the disk is, the more typical it is that the motion will stop at the limit directionŝ w 3,4 corresponding to a combined translational-rotational motion.

Fig. 12
Left panel: Sketch of the fast dynamics and limit directions on a unit sphere in the case of the thick disk on a flat surface. Cases 1 and 2 are the same as those in the previous example (see Fig. 9). In Case 3, the attracting-dominant limit directions leave the circle w 2 = 0, and the former invariant circle becomes a twisted invariant curve between the limit directions (denoted by the dashed green line). Physically this effect causes the rotating thick disk to stop slipping after turning away from the direction of the loading force by a certain angle (see the right panel). This angle can be read from the unit sphere of slipping directions (see the projected unit sphere in the middle panel) (Color figure online) • In the case of the thin disk, the limit directionsŵ 3,4 are located on the invariant circle w 2 = 0 (see the right panel of Fig. 9). Now, (125) shows that in general, these nontrivial limit directions move away from this circle, and the component w 2 = 0 grows with the increasing η = μh/R. Therefore, at the stopping of the motion of the thick disk along the nontrivial limit directionsŵ 3,4 , the velocity of the reference point is not parallel to the direction of the loading force, and the trajectory of the centre of the disk turns away from the direction of the loading force. This result is consistent with the heuristic estimation of Farkas et al. (2003) and with numerical results of Penner (2001).
The results are summarized in Fig. 12. We close the analysis of this example by checking whether the pressure distribution (69) gives a nonnegative pressure value at all points of the contact region. For the circular contact with this linear pressure distribution, the pressure is always nonnegative if p 2 0 − R 2 (c 2 x + c 2 y ) ≥ 0. By direct calculation from (122), this is equivalent to 160η 2 w 2 3 − 144η 2 + 9 ≥ 0. If it is satisfied for all 0 ≤ w 3 ≤ 1, we should prescribe η ≤ 1/4. Therefore, the analysis is valid below a certain thickness given by μh/R < 1/4; otherwise, the assumption of the linear pressure distribution could give negative normal pressure at some points.
• The condition F o = 0 shows the boundary case between the attracting and repelling limit directions (see Def. 1). By direct substitution, we get the critical force F 2/2 crit = 1 − 2ζ 2 . It can be checked that at |F| = F 2/2 crit , the limit directionŵ 1 turns from attracting into repelling.
After summarizing the results, we get the following structurally different regions: Fig. 13 Left panel: sketch of the fast dynamics and limit directions on a unit sphere in the case of the thin disk with an in-plane loading force and an in-plane torque. Compared to the thin disk without the loading torque (see Fig. 9, a new, fourth scenario appears with six limit directions. Right panel: The four different dynamical cases depicted on the parameter space of the loading force and torque. On the coloured lines, bifurcations occur involving the limit directions (Color figure online) Case 1 If |F| > F 2/2 crit , we get a similar behaviour that of Case 1 in Sect. 6.1 without the loading torque.
Case 2  crit > |F| > 0, we obtain a new type of behaviour, which was not present in the example without the loading torque. It can be checked that at |F| = F 4/6 crit , the limit directionsŵ 5,6 are born fromŵ 1 through a pitchfork bifurcation. Thus, we get six limit directions:ŵ 1 andŵ 2 are attracting-isolated,ŵ 3 andŵ 4 are attracting-dominant, w 5 andŵ 6 are attracting saddle-type. The fundamental qualitative behaviour does not change much compared to Case 3; the body typically stops slipping along the limit directionsŵ 3 andŵ 4 . However, the birth of the new two limit directions changes the structure of the fast dynamics around these special motions. The different regions in the parameter space can be seen in Fig. 13.
Finally, let us check the validity of the linear pressure distribution. As we showed in the previous example in Sect. 6.2, we can use the condition p 2 0 − R 2 (c 2 x + c 2 y ) ≥ 0 to check whether the normal pressure is nonnegative everywhere in the contact region. In this example, this condition is equivalent to |ζ | = |T |/(mgr) ≤ 1/4.

Conclusion
We considered a finite planar contact region between two contacting bodies, where the (relative) velocity distribution in the contact region is assumed to be described by planar rigid body motion. By integrating the Coulomb friction model applied locally to the normal and tangential pressure distributions, the resultant force and torque components depend on three kinematic variables of the slip and spin motion in the contact region. When applied to problems with rigid body dynamics, the resultant friction force and torque characteristics lead to a codimension-3 discontinuity in the vector field of slipping, which discontinuity set corresponds to the case of static (sticking or rolling) contact. By utilizing recently developed mathematical tools, the dynamics of this vector field at the discontinuity was investigated by spherical decomposition of the variables and by using the blow-up method. It was shown that the resulting multiple time scale system can be decomposed into fast and slow dynamics at the discontinuity, which revealed the structure of transitions between slipping and sticking states. Although kinematically, the stick-slip transitions could occur in any direction, the fixed points of the fast subsystem generate specific limit directions, where the transitions are dynamically realized. By finding the number and type of limit directions, we can explore the typical types of behaviour of the contacting bodies.
We created an ellipsoidal approximate model in the space of frictional force and torque components. This leading-order approximation is realized as linear formulae in the fast variables containing the discontinuity. The model was compared with several analytical reference solutions for simple geometries. By applying the ellipsoidal friction model for rigid body dynamics of a flat body moving on a rough plane, we could prove that the number of limit directions is 2, 4 or 6 in the generic cases. The results were demonstrated through mechanical examples where the different types of behaviour were explored.
For further research work, it would be useful to select some special benchmark problems and create detailed numeric and experimental studies to check the analytical results. A further major goal would be to include the coupled discontinuous effect from the rolling resistance at nonconforming bodies.
Finally, the ellipsoidal model could be developed by adding nonlinear terms in order to improve the accuracy. Some examples showed that nonlinear terms can have a significant effect on the bifurcation scenarios. Therefore, appropriate forms of nonlinear extensions are indeed useful.