On the nonsmooth dynamics of towed wheels

In this paper, a nonsmooth model of towed wheels is analysed; this mechanism can be a part of different kind of vehicles. We focus on the transitions between slipping and rolling in the presence of dry friction. The model leads to a three-dimensional dynamical system with a codimension-2 discontinuity. The systems can be analysed by means of the tools of extended Filippov systems. The essence of the calculation is to find the so-called limit directions, which show the possible directions of slipping-rolling transitions and their properties. By this method, four different scenarios are found. The results are compared to those from the creep models.


Introduction
The towed wheels can be found in vehicles and machines, where the axle of the wheel can freely rotate with respect to the body of the vehicle. Simple everyday examples of towed wheels are the wheels of shopping carts and baby strollers, but the mechanism of towed wheels cover the trailers [23], the aeroplane nose gears [22], and the front wheels of motorbikes [19]. The dynamical properties of these wheels are important from the point of view of comfort and safety.
One of the fundamental phenomena of towed wheels is the self-excited oscillations called shimmy, which has been intensively studied in the literature (see [18,21] and their references). There are many possibilities for modelling the tyre and the vehicle, and the resulting models often lead to complex mathematical models including, for example, time delay [3,4].
In this paper, we try to focus on the simplest relevant models and try to understand the basic dynamical effects coming from the nonsmooth property of the friction forces between the wheel and the road. The resulting nonsmooth dynamical system is analysed by means of the tools developed by the authors [1] recently.
The paper is organized as follows: In Sect. 2, the equations of motion of the system is derived for such coordinates that give the possibility to introduce several different contact models. In Sect. 3, we review and compare a few fundamental models of the tangential contact forces, including smooth and nonsmooth models. In Sect. 4, we analyse the towed wheel model with a linear creep force, which analysis serves as a reference when comparing it to the nonsmooth model. In Sect. 5, we present a throughout analysis of the wheel model with the spatial Coulomb friction, by M. Antali (&) Á G. Stepan Department of Applied Mechanics, Budapest University of Technology and Economics, Budapest 1521, Hungary e-mail: antali@mm.bme.hu using the recently developed calculation methods of this kind of nonsmooth system. Our purpose is to explore the possible scenarios when the system is switching between rolling and slipping states. In Sect. 6, we show that the results from the nonsmooth model can appropriately predict the behaviour of the system with the full nonlinear creep model.

Equations of motion
Consider the model of the towed wheel consisting of the wheel itself and the caster which connects the wheel to the vehicle (see Fig. 1). The wheel-caster mechanism can rotate in the horizontal plane around the king pin. The radius of the wheel is R, the length of the caster is L, and the distance of the centre of gravity B of the caster from the king pin is H. The vehicle is assumed to move with a constant speed v on a straight road. The angle of the caster is measured by w, and the rotation angle of the wheel is measured by #.
This two-degree-of-freedom model is the simplest relevant model of the towed wheel. In this paper, we focus on the nonsmooth dynamical effects caused by the slipping-rolling transitions of the wheel. For this purpose, the differential equations are expressed in such variables that prove to be effective for the subsequent analysis.
The dynamics is described in the inertial reference frame, but the vector quantities are expressed in a rotating coordinate system. The basis vectors i and j are co-rotating with the caster in the horizontal plane (see the Fig. 1). The basis vector i points forward in the longitudinal axis of the caster, and j points to the left. The basis vector k is fixed and points vertically upwards.
In this coordinate system, the geometric relations between the typical points of the mechanism are given by the vectors The velocity of the point A is and the angular velocities of the caster and the wheel are respectively. All along this paper, the dot above a symbol denotes differentiation with respect to time. Then, the kinetic energy of the mechanism can be written into the form where m, S Az and J Az are the total mass, first mass moment, and the (second) moment of inertia of the wheel-caster mechanism, the latter two are calculated with respect to the king pin about the vertical axis. J Cz Fig. 1 The sketch of the mechanical model analysed in the paper. The mechanism consisting of a wheel, a caster and the king pin can be a part of several types of vehicles denotes the moment of inertia of the wheel about its symmetry axis.
The force system acting between the wheel and the road is reduced to a concentrated force F D and a concentrated moment M D acting at the theoretical contact point D. These vectors are denoted by where T x and T y are the friction force components, N is the normal force between the surfaces, and M x ,M y ,M z are the components of the moment. The details of modelling of these quantities are discussed in Sect. 3. Our two generalized coordinates are w and #; the Lagrange equations of the second kind leads to the following equations of motion: Let us transform the equations of motion into a more appropriate form for the subsequent analysis. Consider the velocity v D ¼ u x i þ u y j of the contact point D where the components are the slipping velocites of the wheel on the ground. We describe the position of the bodies by the generalised coordinates w; #. In addition, the velocity state of the bodies is parametrised by the quasivelocities u x ; u y (see [7], p. 217). Then, the equation of motion (7) can be transformed into four first-order differential equations where are the dynamic equations for the evolution of the quasi-velocities and are the kinematic equations for the evolution of the generalised coordinates. In (9)- (10), the definition of the reduced equivalent masses m x and m y are The equations (9)-(12) still do not form a complete system of ordinary differential equations, because we have not defined the quantities T x , T y , M y and M z . In the next subsection, we choose some nonsmooth and smooth models which are used in the analysis of the towed wheel.

Modelling of the contact forces
The detailed modelling of the contact forces between the wheel and the road leads to complicated models including several phenomena such as dynamical models or time delay (see [3,4] or, for a throughout overview, [18]). In this paper, we focus on the understanding of the fundamental effect of the nonsmooth behaviour of the contact. Thus, we restrict ourselves to the simplest relevant models by considering the following assumptions: • We assume that the ratio of the magnitudes of the moments M y , M z and forces T x , T y is proportional to the typical size q of the contact patch [15]. Then, by considering a small contact patch in the sense of q ( minðR; LÞ, the moments M y and M z can be neglected in (9)- (10), • We consider stationary models, that is, we assume that no additional dynamical variables are considered in the contact. Thus, the contact forces T x , T y are expected to depend on the state variables w; #; u x ; u y only. • We assume that the isotropy of the road and the rotational symmetry of the wheel makes the contact forces independent from the angles w and #.
Finally, we look for the formulae of the tangential forces in the form T x ðu x ; u y Þ and T y ðu x ; u y Þ.

Coulomb friction model
If the stiffness of the wheel-road contact is sufficiently rigid, then the contact forces can be calculated by the simple Coulomb friction model. Then, the slipping friction force is required to be the opposite direction to the slipping velocity and to have a constant magnitude lN where l is the friction coefficient. In planar mechanical problems, it would be described by the formula T C x ðu x Þ ¼ ÀlN Á u x =ju x j. In spatial problems, the formula becomes Note that this model is discontinuous where both slipping velocities u x and u y tend to zero, which issue is the main challenge of the present analysis. The discontinuity at u x ¼ u y ¼ 0 corresponds to the case where the wheel is rolling. The expression of (15) can be extended to the rolling case but then, we require That is, the full dependence of the tangential forces on u x and u y is not a function but a relation (15)-(16), often called a set-valued force law [6]. In the rolling state, the values of T x and T y are not specified explicitly from the contact but they are computed from the dynamics as contact forces with respect to the rolling constraint Note that we implicitly assumed that the static and dynamic friction coefficient coincide. It was shown in [2] that this simplification does not imply loss of generality in the analysis. All things considered, the Coulomb friction model leads to well-defined slipping and rolling states given by (15) and (17), respectively (see the continuous line in Fig. 2). In this paper, we focus on the dynamical consequences of transitions between these two states. However, it is important to compare the results to the models where finite stiffness of the contacting bodies is considered.

Linear creep model
When the stiffness of the contacting bodies is not large enough at the contact point, the local deformation has to be considered, which leads to the phenomenon creep.
If we consider the problem locally, we can experience that the theoretical contact point becomes a contact patch. In a chosen state of the system, slipping and sticking behaviour can occur simultaneously in the different regions of the contact patch [3].
If we consider the motion of the whole wheel as a rigid body, this phenomenon leads to the effect that the pure rolling behaviour does not occur for finite (nonzero) tangential forces. Although the physical background is the co-existence of the local slipping and sticking regions, it seems that the wheel is slipping permanently. This phenomenon is called creep, and in this context, the velocities u x and u y can be called creep velocities.
The essence of this phenomenon can be modelled by a linear creep model (see Fig. 2), where the tangential forces are proportional to the slipping (creep) velocities, T lc x ðu x ; u y Þ ¼ Àc x u x ; T lc y ðu x ; u y Þ ¼ Àc y u y : The creep coefficients c x and c y are influenced by the contact stiffness of the bodies and by the velocity state x ðu x ; 0Þ ¼ Àc x u x is denoted by a dotted line. The Coulomb friction model T x ðu x ; 0Þ ¼ ÀlNu x =ju x j is denoted by a solid line of the contacting bodies. In our model, the velocity of the wheel does not differ much from v for not too large caster angles w, thus, c x and c y can be considered constants for a fixed velocity v. However, the dependence on v can be important when we investigate the effect of different parameters. In the literature, the lateral creep force is often expressed by the creep angle b where tan b ¼ u y =v [3]. Then, a usual empirical linear creep model is expressed as T lc y ¼ Àc y b % Àc y Á u y =v, where the parameterc y is already independent on the velocity v. That is, c y % c y =v is proportional to the reciprocal of v. We expect a similar tendency for c x .

Nonlinear creep model
We can summarize the difference between the previous two models: The Coulomb friction model is appropriate for contact with large contact stiffness (close to a rigid body), and the linear creep model is appropriate for problems with small contact stiffness.
In the literature, these two types of behaviour are usually combined, and nonlinear creep models are obtained (see Fig.2). In two dimensions (u y 0), the nonlinear creep is expected to satisfy ÀlN u x ju x j ju x jc x =ðlNÞ ) 1: The two regions of the function can be connected in the form T nc x ðu x ; 0Þ ¼ ÀlN where sigmoidðxÞ is a monotonous function tending to AE1 at x ! AE1 and having a unit slope at x ¼ 0. For example, such a sigmoid function can be represented by a piecewise linear cut function [8], but piecewise polynomial functions can be also found in the literature [9]. Alternatively, a single smooth sigmoid function can be chosen such as erfð ffiffi , or tanhðxÞ [14]. Other possibilities can be found in [20].
Let us represent the sigmoid function by the hyperbolic tangent function. Then, the three-dimensional (u y 6 ¼ 0) extension of (20) becomes Nonlinear creep models such as (21) are known and widely used in the literature, and they are confirmed by experiments and numerical simulations (see [11][12][13]17]). Thus, although much more complicated models are required in given circumstances, we now consider (21) as a validated model, and (15), (18) as approximating models of (21).

Role of the approximating models
The parameters express the dominant character of the nonlinear creep model (21), whereũ is the typical amplitude of ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi u 2 x þ u 2 y q in the given problem. If C x ; C y ( 1 then the nonlinear creep model (21) tends to the linear creep model (18). However, if C x ; C y ) 1 then the nonlinear creep model (21) tends to the Coulomb model (15). That is, the linear creep and the Coulomb model can be considered opposite limit cases of the more realistic nonlinear creep model (compare the curves in Fig. 2). This fact can be shown more clearly from the expansion of the nonlinear characteristic curve. Consider the hyperbolic tangent function, which was chosen as the sigmoid function when constructing (21). At x ¼ 0, the Taylor expansion of tanhðxÞ is where the coefficients of the terms can be expressed by the Bernoulli numbers ( [16], p. 291). Note that the convergence of this series is extremely slow. Even for considering the first hundred terms from the series, the truncated series is impossible to make an appropriate approximation of the saturation of the curve at AE1.
(see the left panel in Fig. 3). This phenomenon seems to be typical for other sigmoid functions, as well.
Therefore, if an approximate analysis is needed, including the saturation effect, other methods need to be found. It has been shown (see p. 18 of [10]) that the sigmoid function can be effectively approximated by asymptotic series where the leading term is the sign function. This asymptotic expansion of the hyperbolic tangent function is [16] tanhðxÞ ¼ sgnðxÞ Á 1 À 2 expðÀ2jxjÞ þ 2 expðÀ4jxjÞ ð À 2 expðÀ6jxjÞ þ 2 expðÀ8jxjÞ À . . .Þ: This series initially gives a good approximation for large absolute values of x, and by increasing the terms considered in the series, the convergence is relatively fast towards the small absolute values of x (see the right panel of Fig. 3). The first approximation is the sign function itself, which immediately contains the saturation of the original curve at x ¼ AE1.
In this sense, the linear term x and the sign function sgnðxÞ are the leading terms of the two expansions (23)-(24) of the sigmoid function tanhðxÞ. By comparing this to the formulae (15), (18), (21) and Fig. 2, we can get a similar conclusion for the contact models. The linear creep model and the Coulomb friction model can be considered the leading approximations of two different expansions of the nonlinear creep model.
The advantage of these two approximate models is the availability of analytic methods to explore the behaviour of the system. The linear creep model can be investigated by usual linear stability analysis and multi-scale considerations. The main new result of the paper is the analysis of the system with the Coulomb model by using the nonsmooth methods introduced recently by the authors [1].

General set of ODEs
If we formally insert any of the models (15), (18) and (21) into (9)- (12) in the form T x ðu x ; u y Þ and T y ðu x ; u y Þ, we get a system of four ODEs. Note that the rotation angle # does not appear on the right-hand side of these equations. Thus, the ODE (12) of this cyclic variable can be omitted, and the dynamics of the other three variables is determined by the system In vector form, (25) can be written as The following part of the paper contains the analysis of the system (25) by assuming the models (18) and (15). Section 4 contains the linear stability analysis of the system with the linear creep model. In Sect. 5, a throughout analysis is carried out on the system with the Coulomb model by using the nonsmooth analytical tools of extended Filippov systems.

Analysis with the linear creep model
Let us substitute the formulae of the linear creep model (18) into the set (25) of general equations. Then, we get The trivial equilibrium of this smooth system is located at ðu x ; u y ; wÞ ¼ ð0; 0; 0Þ. The linearised differential equation around the trivial equilibrium can be written in the form The eigenvalues of the coefficient matrix are determined by the characteristic equation and the eigenvalues are The real part of all eigenvalues is negative, and the trivial equilibrium is asymptotically stable.
The structure of the dynamics is significantly affected by the parameter Note, that v/L is the coefficient of the decay of the caster angle at the pure rolling motion (see (36)), which we call the kinematic effect. At v ¼ 0, the lateral slipping velocity (creep velocity) u y decays with a coefficient c y =m y (see (27)), which we call the lateral creep effect. Thus, the parameter S expresses the ratio of the time scale of these two effects. As it was explained at the end of Sect. 3.2, the parameters c x and c y are supposed to be approximately proportional to 1/ v. If S is relatively large (high velocity or low wheel stiffness) in the sense S [ 1=4 then (30) shows that the eigenvalues k 2;3 become complex, which causes oscillatory decay of the variables. When S\1=4 (low velocity or high wheel stiffness), all eigenvalues are real, and there are no oscillations during the exponential decay.
If S is very small in the sense of S ( 1 then the eigenvalues tend to k 2 ! Àc y =m y and k 3 ! Àv=L. In this case, the creep effect and the kinematic effect separates to different time scales. As c x =m x and c y =m y are supposed to have the same order of magnitude, we can say that we have a multi-scale system with the fast eigenvalues k 1 and k 2 and the slow eigenvalue k 3 ( k 1 ; k 2 . In this case, the eigenvector of k 3 tends to the w axis, and the corresponding stable manifold lays close to the line u x ¼ u y ¼ 0, as well (see Fig. 4).
How do the trajectories of the system behave in the vicinity of the equilibrium at the origin? Due to the multi-scale dynamics with k 1 ; k 2 ) k 3 , the creep velocities u x and u y decay much faster than the caster angle w. (The eigenvectors are approximately parallel to the coordinate directions.) Thus, the trajectories approach the line u x ¼ u y ¼ 0 very fast, and then, they converge to the origin along this line governed by the slow eigenvalue k 3 ¼ Àv=L. Although pure rolling does not occur, the trajectories are quantitatively very close to that.
It can be shown (see later in Sect. 5.4.2) that m x \m y and m y is typically much larger than m x . If the values of c x and c y are close to each other then c x =m x is typically larger than c y =m y . That is, during the fast dynamics in the plane u x À u y , most trajectories approach the w axis along the eigendirection of the eigenvalue k 2 % Àc y =m y . Later, this property of the eigenvectors can be compared to the results from the nonsmooth analysis with the Coulomb model.

Nonsmooth vector field from the Coulomb model
By substituting the expression of Coulomb model (15) into (25), we get This set of differential equations is valid for the case of slipping. In the 3D phase space of the system (32), the vector field has a discontinuity on the line u x ¼ u y ¼ 0, which line is called discontinuity set and denoted by R (see the left panel of Fig. 5). We can say that R is a codimension-2 discontinuity set because it has two dimensions less than the full phase space. Thus, the discontinuous dynamical system (32) can be categorized as an extended Filippov system. The mathematical background of these systems, more general definitions and solution methods can be found in [1]. Now, we restrict ourselves to the introduction of the concepts and formalisms which are necessary for the analysis of our mechanical problem.
In the discontinuity set R, and the limit of the vector field (32) does not exist. (Note the singularity at u x ¼ u y ¼ 0 at the denominator of the friction terms.) However, we can define the directional limit of the vector field in the phase space.
Denote the polar angle in the plane u x À u y by /. (See the right panel of Fig. 5.) Then, the slipping velocities can be written in polar coordinates as is the magnitude of slipping and the physical meaning of / is the direction of slipping.
Let us select a value w and a point x R ¼ ð0; 0; wÞ in the discontinuity set R. Then, the function nð/Þ ¼ ðcos /; sin /; 0Þ provides the possible unit vectors perpendicular to R. Consider the directional limit which we call limit vector field at x R , depending on the angle /. By using the notation f Ã ¼ ð _ u Ã x ; _ u Ã y ; _ w Ã Þ, the components of the limit vector field of (32) become This limit vector field can be used for two purposes. On the one hand, it provides the rolling dynamics, which process is almost trivial in our example. On the other hand, the limit vector field (35) provides the possibility of a detailed analysis of the rolling-slipping transitions.

Rolling dynamics
The rolling constraint (17) shows that the rolling state of the system corresponds to u x ¼ u y ¼ 0, and the In the vicinity of the origin, the trajectories converge fast to the line u x ¼ u y ¼ 0 and then, they approach the origin being asymptotically close to this line dynamics in the rolling state is described purely by the caster angle w.
If the third component of the limit vector field (35) depended on the direction /, then the rolling dynamics could be calculated by a convex combination in an integral form (see [1] for details). However, in our system, the third component _ w Ã is independent of /, and thus, it directly gives the rolling dynamics in the form (Note that an alternative way to derive (36) is the pure kinematic description of the system in the presence of the rolling constraint (17).) The solution of (36) can be given analytically, the solution of the linearised equation _ w ¼ Àv=L Á w is a very good approximation of the full solution wðtÞ ¼ arccos where w 0 ¼ wð0Þ is the initial condition. In the case of pure rolling, the caster exponentially tends to the centre position w ¼ 0.
In the next subsection, we analyse the solutions of the system when transitions occur between rolling and slipping.

Transitions between slipping and rolling
Let us express the first two components of the limit vector field (35) by using polar coordinates according to (33). Then, we get where According to (39)-(40), the function Rð/Þ shows the rate of change of the magnitude u of the slipping velocity, and Vð/Þ is related to the rate of change of the angle of the slipping velocity. By direct calculation, we get Vð This polar coordinate form is useful to find and analyse the possible directions (angles) where the transitions between rolling and slipping occur.
The system in the original form (32) is not defined at the discontinuity at u x ¼ u y ¼ 0. That is the reason why the limit (40) is not defined except for some special values of /. This limit exists at the roots of the equation Vð/Þ ¼ 0. These roots are called limit directions.
Definition 1 (Limit direction) An angle / 1 2 ½0; 2p and the corresponding direction is called a limit direction of the dynamics at the point x R 2 R if Vð/ 1 Þ ¼ 0.
By finding the limit directions, we can determine the possible directions of transition between slipping and rolling. When the slipping wheel starts rolling, the slipping velocity vanishes along the corresponding limit direction. When the rolling wheel starts slipping, the slipping velocity is initiated along the corresponding limit direction.
Our analysis has three main steps to achieve a qualitative overview of the slipping-rolling transitions: 1. Determine the number and location of limit directions, which show the direction of the slipping-rolling transitions. 2. Determine the attracting or repelling property of the limit directions, which shows whether slipping turns into rolling or vice versa. 3. Determine the dominant or isolated property of the limit directions, which show whether the limit direction corresponds to a typical or a special solution.
These calculations can be performed analytically if the linearised system of (32) is considered where the caster angle is assumed to be small in the sense of w ( 1. Then, we get the system where the higher-order terms in w are neglected. Thus, the functions R and V from (43)-(44) become

Number and location of the limit directions
According to Definition 1, the limit directions are determined by Vð/Þ ¼ 0. The roots of (47) can be written in the form The limit directions / 1 and / 2 always exist, but the limit directions / 3 and / 4 exist only for a certain range of the caster angle w. It will be shown in Sect. 5.4.2 that for the physically relevant mass distribution of the wheel, we have 1=m x [ 1=m y , and thus, the parameter m is always positive. Then, we can state the following proposition: Proposition 1 A point x R ¼ ð0; 0; wÞ of the system (45) possesses four limit directions given by (48)-(49) if If jwj [ lNL=ðv 2m Þ then the system has two limit directions given by (48).

Attracting and repelling limit directions
In the attracting case, the trajectories approach the discontinuity set along the limit direction, which corresponds to the slipping-to-rolling transition. In the repelling case, the trajectories leave the discontinuity set, which corresponds to the rolling-to-slipping transition. By substituting the values from (48)-(49), we get Rð/ 3;4 Þ ¼ ÀlN Thus, we can state the following proposition: The limit directions / 3 and / 4 in (49) are always attracting. The limit direction / 1 is attracting for w\lNL=ðv 2 m y Þ, and / 2 is attracting for w [ À lNL=ðv 2 m y Þ.
If all limit directions are attracting (jwj\lNL=v 2 m y ), then a transition from rolling to slipping is not possible, and from a rolling initial condition, the wheel sustains the rolling state. Moreover, if a perturbation initiates a small amplitude of slipping, the system robustly comes back to the rolling state along one of the attracting limit directions. However, if a repelling limit direction appears, then rolling-slipping transition becomes possible. We can say that the rolling behaviour suffers from some kind of instability because after a small perturbation, an increasing amount of slipping initiated along the repelling limit direction.
As an interest, one can check that the boundary of these two scenarios coincides to the maximum admissible friction force calculated from the Coulomb model when the static friction coefficient is assumed to be l 0 ¼ l. That is, a part of the analysis of the limit directions provides the alternative method to determine the possibility of rolling purely from the slipping differential equations.

Dominant and isolated limit directions
In the dominant case, the limit direction attracts the surrounding trajectories when approaching the discontinuity in the positive or negative direction of time. Then, the limit direction is an asymptote of continuously many trajectories. In the isolated case, the limit direction is repelling the surrounding trajectories at the discontinuity. Thus, only a single trajectory is connected to the discontinuity along the limit direction.
From the mechanical point of view, the dominant attracting limit directions are the typical directions of transitions from slipping to rolling, and at the isolated attracting limit directions, the transition occurs in a special case only. In case of repelling directions, the authors have found only isolated ones in mechanical systems with Coulomb friction, and this seems to be a generic property (see a partial proof in [2]).
By substituting the values from (48), (49) into the conditions of the definition, we get The consequences can be summarized in the following proposition: The limit directions / 3 and / 4 in (49) are always isolated. The limit directions / 1 and / 2 are dominant for sufficiently small absolute values of w, and this property can be changed at the values jwj ¼ lNL=ðv 2 m y Þ and jwj ¼ lNL=ðv 2m Þ.
That is, for very small caster angles w, the transition from slipping to rolling typically occurs along the limit directions / 1 and / 2 . By increasing the absolute value of w, the situation gets more complicated. The different scenarios are summarized in the next subsection through a bifurcation analysis.

Limit directions and bifurcations
During the previous subsection, we analysed the dynamics of the slipping velocities u x and u y only, and we considered the caster angle as a parameter. This approach and the subsequent bifurcation analysis needs to be confirmed from the point of view of the strength of the leading terms of the dynamics.
When we calculate the Taylor expansion of _ w from (45) at a x R ¼ ð0; 0; wÞ, the leading term is linear. When this linear term is non-zero, the dynamics near this point is dominated by the linear part, the higherorder terms are less significant, and the tendency of the solution is exponential.
In the case of the time derivatives _ u x and _ u y in (45), the dominant term is the nonsmooth term containing the square root. It can be showed by the transformation of the time (see [1]), that in the vicinity of the discontinuity set R, the nonsmooth terms generate a convergence that is faster than exponential, and the trajectories reach the discontinuity in finite time in the forward or backward direction of the time.
Compared to these nonsmooth terms, the effect of the linear and higher-order terms can be neglected when we are close enough to the discontinuity. Then, as the dynamics of the caster angle w does not contain nonsmooth terms, its change can be neglected, and w can be fixed as a parameter, and the system can be reduced to the plane of u x and u y . (Of course, when we are further from the discontinuity, the dynamics of w becomes significant.) From this consideration, we can treat the state w as a parameter, and this can be used as a bifurcation parameter for categorizing the typical behaviours of the system. A value w is considered as a bifurcation point from the point of view of limit directions then either the number of the limit directions or their properties (attracting-repelling, dominant-isolated) change. From summarizing the previous subsection, we can state that bifurcations occur at (The subscripts refer to the type of the bifurcation, which is explained later.) To compare the location of these bifurcation points, we have to determine the relationship between the reduced parameters m y and m.

On the inertial parameters of the model
We have not specified the modelled system; the towed wheel can be a part of several wheeled machines from a shopping trolley to an aeroplane. Thus, the values of the parameters can be significantly different in the applications. Still, we can obtain information about the relation between the equivalent masses m x , m y (see (13)) andm (see (50)) which appear in the formulae (56), (57) of the bifurcation points.
To minimize the number of the parameters and make the formulae easier, we used the moments of inertia J Cy and J Az . If we want to express them with the mass m c of the caster and the mass m w of the wheel, we get J Cy ¼ j wy m w R 2 ; The dimensionless moments of inertia of the wheel with respect to its centre of gravity are denoted by 0\j wy \1 and 0\j wz \1=2. The dimensionless moment of inertia of the caster with respect to the king pin A is denoted by 0\j cz \1. Then, from (13), By direct calculation from (50), it can be shown that m\m y is satisfied when This expression would be negative if compared to the caster, the wheel is very small and very heavy at the same time, which is physically unrealistic. That is, for physically relevant parameters,m\m y , Thus, by increasing the amplitude of the caster angle in the linearised model (45), the system first goes through the bifurcation point w T and after that, the bifurcation point w P .

The bifurcations of the limit directions
Now, we have all the information to analyse the bifurcations of the limit directions in the system (45), which is summarized in Fig. 6. The region jwj\w T is denoted by number 1 in Fig. 6. Here, there are four limit directions, all of them are attracting, two of them are dominant (/ 1 and / 2 ) and the other two is isolated (/ 3 and / 4 ).
At jwj ¼ w T , one of the dominant directions turn around and becomes an isolated repelling direction. According to the similar name from the literature (see [10] or [5]), we call this bifurcation a tangency bifurcation, because, at the bifurcation point, one trajectory is exactly tangent to the discontinuity set.
The region w T \jwj\w P is denoted by number 2 in Fig. 6. Here, there are still four limit directions, but one of them is repelling.
At jwj ¼ w P , the remaining three attracting directions collide at one point. After that, a single attracting, isolated limit direction remains. We call this bifurcation a pitchfork bifurcation of limit directions.
The region w P \jwj is denoted by number 3 in Fig. 6. Here, there are only two limit directions; both are isolated, one attracting and one repelling.
After the overview of the bifurcations, let us explore what happens in the phase plane of the regions 1-3 and what are the mechanical consequences.

The three typical scenarios of the system
If we sketch the parameter map of the caster angle w and the velocity v of the vehicle, we get the graph of the left side of Fig. 7. Here we can detect the three regions and the two bifurcations from Fig. 6. In each scenario, the limit directions help us to understand what happens with the system from the point of slipping-rolling transitions. Both the figure and the explanation corresponds to the positive caster angle w, but the opposite direction can be treated in a similar way.
In region 1, all the four limit directions are attracting, and thus, there is no possibility to escape from the discontinuity. Therefore, the rolling motion of the wheel is realizable. If a perturbation initiates a small amount of slipping, the system returns back to rolling in finite time. Then, the transition from slipping to rolling occurs along one of the dominant directions. The isolated directions behave as separatrices between the areas denoted by 'a' and 'b' in Fig. 7. From the area 'a', the slipping-rolling transition occurs along the Fig. 6 Bifurcation diagram showing the location and type of the limit directions / i against the caster angle w. The solid line denotes the dominant attracting limit directions. The dashed line denotes the isolated attracting limit directions. The double line denotes the isolated repelling limit directions. The tangency bifurcation is denoted by T, and the pitchfork bifurcation of limit directions is denoted by P limit direction / 1 , that is, the slipping velocity vanishes from the positive u y direction. From the area 'b', the slipping velocity vanishes from the negative u y direction, along the limit direction / 2 .
In region 2, one of the limit directions becomes repelling, and thus, the continuous rolling becomes impossible. When a perturbation pushes the system into the area denoted by 'c', the solution escapes along the limit direction / 2 , and slipping initiated in the negative u y direction. In this region, the dominant attracting direction / 1 still exists, which has interesting consequences. If a perturbation pushes the system into the area 'd', the solution first converges to the origin along the attracting direction, and then, it escapes along the repelling one. That is, there is a transition from slipping to rolling, but at the same moment, the system switches from rolling to slipping, again.
This latter behaviour vanishes in region 3, where only two limit directions remain. The continuous rolling is not realizable in this region either. For almost any small perturbations, the solutions escape along the limit direction / 2 . That is, the system just starts slipping without any complications.

Results for the full nonlinear system
In Sects. 5.3, 5.4, the linearised model (45) was analysed assuming small caster angles. For the nonlinear system (32) with the Coulomb model, the procedure cannot be repeated analytically but just numerically. The qualitative result of this process can be seen in Fig. 8. Two significant differences appear compared to the linear model.
On the one hand, a new, fourth region appears for large caster angles. In this region, there are two attracting limit direction, one dominant and one isolated. Unlike in region 1, there is only one dominant direction for the slipping-rolling transition. For almost all perturbations, the slipping velocity vanishes along the limit direction / 2 .
On the other hand, for large caster angles, the limit direction / 1 and / 2 do not remain in AEp=2, but they can have a significant deviation from that (see Fig. 8). ) For the first sight, we can see some analogy between the eigendirections of the smooth system (with creep) and the limit directions of the nonsmooth system (with the Coulomb friction model). They both act as organizing directions of the trajectories towards a critical set-the smooth, stable manifold or the discontinuity set. However, there are two fundamental differences between eigendirections and limit directions.
1. The eigendirections are two-directional (lines), but the limit directions are unidirectional (halflines). The two sides of an eigendirection correspond to the same eigenvalue and have similar behaviour. But each limit direction can have independent dynamical properties.
2. The eigendirections correspond to an exponential decay characterised by the corresponding eigenvalues, where the full convergence to the critical set requires infinite time. In the case of the limit directions, the convergence is faster than exponential and solutions converge in finite time.
These properties explain the main qualitative differences and the quantitative similarities between the linear creep model and the corresponding similar region of the Coulomb model. In the system with the Coulomb model, the trajectories reach the discontinuity set in finite time. Then, pure rolling begins, and the dynamics is determined by purely the caster angle w. In the creep model with large creep coefficients, the slipping velocities decay much faster than the caster angle, but with a similar exponential tendency. Thus, the system gets very close to the rolling state but reaches the pure rolling only when the caster angle decays to zero.
The main difference is that the Coulomb friction model contains the saturation of the contact force which leads to a new phenomenon in the scenarios 2-4 including the vanishing of the continuous rolling state and the different kinds of slipping-rolling transitions.

Comparison of the Coulomb model and the nonlinear creep model
The system with the nonlinear creep model is analysed by means of numerical simulations. As we did in Fig. 7, the vector field is projected to the plane u x À u y , and the resulting phase portrait shows the behaviour of the system in the vicinity of a chosen state x R ¼ ð0; 0; wÞ. In Figs. 9, 10 and 11, the system parameters are set to L ¼ 0:4 m, N ¼ 18:75 N, m x ¼ 0:5 kg, m y ¼ 1:458 kg, l ¼ 0:8, c x ¼ c y ¼ 100 Ns/m, v ¼ 10 m/s. The state w in Figs. 9, 10 and 11 is chosen to w ¼ 0:024 ¼ 1:38 , w ¼ 0:048 ¼ 2:75 and w ¼ 0:096 ¼ 5:50 , respectively. These values are related to scenarios 1, 2, and 3, respectively (see Fig. 7).
In each diagram, the dashed circle denotes the region where the relative difference between the forces from the nonsmooth (Coulomb) and nonlinear creep models is more than 10%. Inside the circle, the creep effect is significant, but outside this circle, the phase portraits are very similar to those obtained from  Fig. 7 when the original nonlinear model (32) is considered. The continuous lines denote the bifurcations of the nonlinear system; the dashed lines correspond to previous results of the linearised system. A new, fourth region appears with two attracting directions the Coulomb model (see Fig. 7). In each diagram, the trajectories with the thick lines are related to the limit directions of the Coulomb model (see Fig. 7).
In Fig. 9, the structure of the phase space is similar to that of the Coulomb model. The small circle shows the nullcline of the dynamics, which is an equilibrium point of the projected dynamics. In the full system, this state is very close to the pure rolling motion, which would be located in the origin of the figure. The two dominant attracting limit directions of the nonsmooth model coincides with one of the eigendirections of the nonlinear creep model (the vertical line in the figure). The continuation of the other eigendirection becomes curved (denoted by a dashed line), and its asymptote is parallel to the isolated limit directions of the nonsmooth model.
In Fig. 10, we can see the effect of the tangency bifurcation of the nonsmooth system, which is related to a structural change in the system with the nonlinear creep model. The nullcline vanishes, and a special trajectory is created pointing downwards from the origin, which attracts the trajectories. This behaviour cannot be shown directly from the nonlinear creep model, but it can be predicted from the repelling limit direction of the nonsmooth system (see Fig. 7). Physically, it means that the u y component of the slipping velocity vanishes, and the slipping typically becomes parallel to the motion of the vehicle. Figure 11 is not much different from Fig. 10; thus, the pitchfork bifurcation of the nonsmooth system does not create a structural change in the nonlinear creep system. However, there is a geometric change: In Fig. 11, all trajectories possess an inflection point. In Fig. 10, there are trajectories with and without an inflection points, corresponding to regions 'c' and 'd' in Fig. 7. That is, the isolated limit directions of the nonsmooth systems separate the regions with different tendencies of the vanishing of the lateral slipping velocity u y that can be seen in the nonlinear creep model.
In all three cases in Figs. 9, 10 and 11, the analytical results from the nonsmooth Coulomb model could be used to predict the qualitative behaviour of the nonlinear creep model.

Conclusion
In this paper, we investigated the simplest relevant model of the towed wheel while focusing on the dynamical effects of different contact models. The model leads to a system of three first-order differential equations for the angle of the caster and the components of the slipping velocity of the wheel. It is shown that, at asymptotic expansion, the Coulomb model can be considered as the leading approximation of nonlinear creep models. In another way, the linear creep model can also be used as an approximation. The analytical investigation of the system was carried out, considering these two models. With the linear creep model, the linear stability analysis and multi-scale considerations are used to obtain qualitative information about the system. With the Coulomb friction model, a special nonsmooth system was created, which was analysed by the recent methods of the extended Filippov systems. By using the concept of limit directions, we determined four typical scenarios regarding the types of possible slipping-rolling transitions.
It was shown that the analytical results from the nonsmooth model could be used to get a throughout overview of the qualitative behaviour of the full nonlinear system close to the rolling behaviour. These analytical results can be used to validate multibody software that use more complex contact models.
In future research, the global dynamics of the system should be investigated. It is also important to introduce stiffness to the king pin to create the possibility of the self-excited oscillations [21]. The presented nonsmooth methods hopefully help to analyse and understand such situations when the self-excited instability is coupled with the dissipative effects from the friction models. From the general analysis of the towed wheel, more specific results can be obtained if the parameters are gained from the selected application.