Theoretical and experimental study on the nonlinear dynamics of wheel-shimmy

The dynamics of the one-degree-of-freedom model of a towed wheel is investigated with the help of delayed tyre model. The partial side slip in the contact region is also considered, leading to an infinite-dimensional, piecewise-smooth system of governing equations. The nonlinear behaviour of the system is explored by numerical continuation of the periodic solutions. Thus, we encounter bistable domains in the space of model parameters which are otherwise undetectable by the standard quasi-steady-state tyre models. The presence of these parameter domains is also confirmed by experiments showing qualitatively similar behaviour as predicted by our calculations.

rolling/sliding wheels can occur in various vehicle systems, such as motorcycles [13], trailers [20], cars, aircraft landing gears [17,18,21] or even baby strollers and shopping trolleys. In the past, large effort was made to understand and to avoid these oscillations, since they can be the source of accidents, reduced driving comfort, higher tyre wear and fuel consumption.
As several former studies showed, these oscillations have very rich dynamics and the associated vibrations are essentially nonlinear [17,18]. In experiments, socalled bistable parameters were identified [1] where although the rectilinear motion is stable, it has a finite domain of attraction. This means that depending on the initial conditions, the system can converge either to the straight-line motion or to a large-amplitude periodic solution (shimmy) at the same towing speed. This behaviour can be explained by a subcritical Hopf bifurcation [8] of the rectilinear motion at the critical towing speed, corresponding to the linear stability boundary. Using the simplest, one-degree-of-freedom model of a towed wheel, the standard quasi-steady-state tyre models do not provide such parameter domains. In the meantime, the existence of bistable parameter regions can be shown if the so-called memory effect in the tyreground contact and the contact friction is considered simultaneously [1,4].
In our study, we capture these effects with the help of the brush tyre model (see in [12]) which enables us to describe the dynamically varying lateral tyre deformation with relatively few tyre parameters. The lateral slip of the tyre in the contact region is considered by parabolic deformation limits, resulting a piecewisesmooth continuous system with its specific bifurcation problems (see [9,14]) while the local bifurcation of the rectilinear motion also can be studied by performing the centre manifold reduction and analysing the nonsmooth normal form [3]; in this paper, we concentrate on the numerical calculation of the periodic solutions exploring the global structure of the limit cycles in the system.
The rest of the paper is organised as follows: In Sect. 2, we introduce our mechanical model and the equations of motion. Then, in Sect. 3, we discretise the governing equations using the method of collocation and we present the algorithm calculating the Floquet multipliers to determine the stability of the periodic solutions. In Sect. 4, we give an overview of the analytical study of the non-smooth Hopf bifurcation of the rectilinear motion. In Sect. 5, we present the results of the bifurcation continuation which are compared to measurements in Sect. 6. We summarise our results and give our conclusions in Sect. 7.

Governing equations
The mechanical model of the towed wheel is shown in Fig. 1. We consider an elastic tyre with a rigid wheel rim in a rigid caster. The caster is towed with a constant speed V in the X direction at the rotational joint A. In the joint, a torsional damping b t is assumed. As a result of the tyre deformation, a lateral force F y and a selfaligning moment M z is generated, which are acting on the wheel. Using the yaw angle ψ as generalised coordinate, one can apply Lagrange's equation of the where the yaw mass moment of inertia J A with respect to the joint A can be given as where l C is the distance of the centre of gravity C from the joint, m is the mass of the caster and the wheel while J C is the yaw mass moment of inertia with respect to the centre of gravity. We use the brush tyre model (see Fig. 2) to calculate the lateral tyre force F y and self-aligning moment M z . By means of this model, the elastic tyre is modelled as massless particles with linearly distributed lateral stiffness k. We assume that the effect of longitudinal deformation is negligible. Thus, we can calculate the tyre force and self-aligning moment using the lateral tyre deformation q(x, t). By means of the brush model, it is assumed that the tyre deforms in the x ∈ [−a, a] contact region only. Consequently, introducing the angular velocity the equation of motion (1) can be expressed as While we calculate the lateral deformation, we distinguish between two cases: when the tyre particles stick to or when they slide on the ground. In general, sticking and sliding may occur simultaneously in different parts of the contact region. In the sticking region(s), the lateral tyre deformation is described by the nonlinear partial differential equatioṅ with the boundary condition q(a, t) = 0 [4,15,16]. This condition implies that the tyre particles enter the contact region with zero deformation. For this problem, a travelling wave solution can be provided, introducing time delay τ (x), distributed along the x-axis. The source of the time delay in the contact can be explained by studying the tyre deformation in a caster-fixed frame of reference. While the tyre particles have zero velocity in the ground-fixed frame of reference, they appear to travel 'backwards', i.e. in the direction of the rear edge x = −a with respect to the caster. In this context, the time delay τ (x) gives the time that is needed for a tyre particle to travel to its actual position x from x 0 where it stuck to the ground. Considering these, the lateral displacement in the sticking part of the contact can be given by while we can also derive a formula for the longitudinal position: where q 0 and x 0 denote the lateral deformation and longitudinal position at the instant when the tyre particle in question stuck to the ground. Using the travelling wave solution has a clear advantage over the semidiscretisation of the PDE (5) into a system of ODEs as it does not introduce new state variables. Thus, the dimension of the nonlinear zero problem, which we obtain by applying numerical collocation in time, can be kept relatively low. The effect of partial sliding is considered assuming parabolic vertical force distribution in the x ∈ [−a, a] contact region. Since by means of the brush tyre model it is assumed that the tyre particles can deform independently from each other, we consider parabolic limits for the lateral deformation, too. Depending on the sliding direction, the lateral deformation in the sliding regions can be given as where N ± defines the magnitude of the limiting parabolae. These can be given as a function of the vertical load F z , the lateral tyre stiffness k, the contact-patch half length a and the friction coefficient μ characteristic to the tyre-ground contact: where the positive sign corresponds to the upper and the negative to the lower deformation limit. The derivation of the time-delayed formulae (6) and (7) is given in [15,16], whereas the detailed explanation and the analytical study of the brush tyre model with time delay and side slip can be found in [4]. The main challenge in handling this non-smooth delay differential equation is that unfortunately, there is no closed-form expression available for the position of the sticking/sliding boundary which is otherwise a function of the time-delayed state variables Ψ (t + s), (t + s) for s ∈ [−τ max , 0] with τ max being the largest time delay in the system. This issue is overcome by finding the sticking/sliding boundary numerically in the course of calculating the lateral deformation.

Discretisation and periodic solutions
First, a collocation method is developed to identify the approximate periodic solutions in the discretised mechanical model. Then, the stability of each periodic solutions is determined by means of the calculation of the corresponding Floquet multipliers.

Numerical collocation
We calculate numerically the periodic solutions of the dynamical system Eqs. (3) and (4) as a two-point boundary value problem over a time period with the help of the method of collocation [19]. Accordingly, the solutions are approximated by m piecewise Lagrange polynomials of degree n [10]. Thus, the time t ∈ [0, T ] is discretised by using nm + 1 equidistant representation points: Note that alternatively, one could apply an equidistant, periodic grid by using a finite segment of Fourier series for the discretisation [2]. In this case, the periodicity of the solutions would be assumed by definition; hence, no additional periodicity condition is needed. However, if a periodic grid is used for the discretisation, the numerical stability analysis of the periodic orbits cannot be performed.
We introduce the dimensionless time ϑ := t/T to project the time period onto the ϑ ∈ [0, 1] interval. The collocation points in the dimensionless time are introduced as with the coefficients p i, j defined by Based on this, the time derivativeu(t) can be calculated: where prime denotes derivation with respect to the dimensionless time ϑ. Evaluating this equation in the representation points given in (10), one can introduce the differentiation matrix D L m,n corresponding to the piecewise discretisation relying on m Lagrange polynomials of degree n. This procedure provides a regular Fig. 3 The structure of the differentiation matrix with piecewise Lagrange polynomials. (In this example 4 sub-intervals are assumed). The shaded segments belong to differentiation matrices for Lagrange polynomials of degree n. In the overlapping parts the matrix elements are added together (mn+1)×(mn+1) matrix which has a block-diagonallike structure (see Fig. 3) where the nonzero blocks are obtained from Lagrange differential matrices D L n relying on n intervals. Introducing the solution vector u contains the function values u i = u(t i ) at the representation points t i for i = 0, ..., mn + 1). The vector of the time derivatives can be given aṡ where the elements of the derivative vectoru are defined asu i =u(t i ). Based on this, the angular velocity function (t) and the yaw angle function ψ(t) are represented in the solution vectors and ψ and their derivatives are calculated aṡ Assuming that the function values i and ψ i for i = 0, ..., nm + 1 are available, we can provide the corresponding lateral deformation in the contact region using the travelling wave solution (6), (7) and the limiting parabolae (8). This procedure is performed starting from the leading edge x = a because q(a, t i ) = 0 is given by the boundary condition attached to the PDE (5). Then, we follow the path of each tyre particles in the wheel-fixed frame of reference (see Fig. 4) in time steps defined by the discretisation of the time period. In each step, we assume rolling to calculate  (3) and (4) with (6) and (7), i.e. when wheel-shimmy occurs. The green curves belong to sticking while the red ones to sliding parts of the tyre-ground contact. The thick blue curve shows how a given point travels towards the rear edge x = −a. The application of Lagrange polynomials provides an equidistant mesh in the time domain while the spatial discretisation along the x ∈ [−a, a] contact patch is a result of the travelling wave solution (7). (Colour figure online) the longitudinal position with the travelling wave solution and where τ i, j is the time delay needed for the tyre particle in question after it sticks to the ground. The longitudinal position and lateral deformation corresponding to this delay are denoted by x 0 i, j and q 0 i, j , respectively. The time-delayed yaw angle can be given as Using the periodicity of the solution, this can be obtained by the modulo operation as Then, the calculated deformation is compared with the parabolic limits (8). Namely, if the deformation is between these parabolic boundaries, then the travelling wave solution is accepted. Otherwise, considering equal friction coefficients for sticking and sliding, the deformation is equal to the corresponding sliding limit, that is, . (21) This procedure is repeated until all the tyre particles leave the x ∈ [−a, a] contact region at the rear end x = −a of the contact line. As a result of this calculation, the deformed shape is available at each time instant t i . We use these values to calculate the integrals I i := I (t i ) at the right-hand side in (3): We approximate this integral by means of the trapezoid rule: where the weights are calculated as with ν i denoting the number of discretisation points in the contact region x ∈ [−a, a] corresponding to the representation point t i . Thus, the discretised system of the equation of motions (3), (4) can be given as To ensure the periodicity of the solution, we attach the following conditions to the discretised system: while in order to fix the phase of the periodic solutions, we also add the condition Using equations (25)-(29), we define the nonlinear zero problem where the solution vector x can be expressed as while the vector λ contains the bifurcation parameters. We carry out the continuation of the periodic solutions by means of two bifurcation parameters: the towing speed V and the caster length l. This two-dimensional continuation is realised by the Continuation Core and Toolboxes (COCO) in MATLAB relying on the 2D atlas algorithm [6].

Stability of the periodic solutions
The stability of the numerically calculated periodic solutions can be decided with the help of the corresponding Floquet multipliers that are the eigenvalues of the monodromy matrix [10]. Let us consider our governing equations (3) and (4) in a delay differential equation form aṡ with a periodic solution where τ max stands for the largest time delay, while the vector of state variables is defined as Let us consider a solution y as a sum of the periodic solution y p and a perturbation v: Substituting the expression above in the delay differential equation (32), we obtaiṅ Linearisation of the system around the periodic solution provideṡ where Df(.) is the Fréchet derivative of the right-hand side. Since the periodic solution y p satisfies the equation of motion (32) by definition, this equation can be simplified aṡ Discretising this equation with the help of piecewise Lagrange polynomials as described in Sect. 3.1, the following equation is obtained: where the vector F contains the right-hand-side f evaluated at the representation points given in (10), while w contains the state variables at the representation points . . .
where the time instant t −r corresponds to the maximal delay τ max and the differentiation matrix D can be given in a block-diagonal form The discretised variational equations (39) can be expressed also in the form The coefficient matrix J = D−grad w F can be obtained from the Jacobian of the zero problem (30) without considering the phase and periodicity conditions and the modulo operation in the time-delayed variables [10]. The Jacobian J mod of the zero problem (30) has a block matrix structure as shown in Fig. 5 where we already omitted the periodicity and the phase conditions, as well as the derivations with respect to the period T . Based on equations (25) and (26), three of the four blocks can be given analytically. However, we cannot calculate the top-right block (shaded in light grey in Fig. 5) directly since this corresponds to the derivatives of (25) with respect to the yaw angle ψ. Accordingly, these derivatives depend on the lateral deformation and the consideration of sliding makes it impossible to provide a formula that links the deformation to the timedelayed state variables in a straightforward way. This is why we calculate this block by finite differences. Nevertheless, we still know that due to the distributed contact delay, matrix J mod is expected to have a band-like The structure of the Jacobian matrix with and without the modulo operation for the time-delayed terms. The light grey shading indicates the part of the matrix which has to be calculated by finite differences. The dark grey bands show the structure of the nonzero derivatives from the terms with distributed delay structure, which, as a result of the modulo operation, appears also in the upper-right triangle of the block.
In order to omit the modulo operation, one can reorganise this matrix resulting the extended Jacobian J in case the largest delay is smaller than the period: τ max < T . However, this cannot be performed directly for τ max > T since in this case the periodic band in the top-right block overlaps itself. Thus, one has to calculate the top-right block of the extended Jacobian by finite differences evaluating the zero problem (30) without the modulo operation.
The monodromy matrix can be calculated by defining a map between the initial and the final solution segments v(t 0 + s) and v(t 0 + T + s), s ∈ [−τ max , 0], respectively. In the discretised system (39), this can be provided by the extended Jacobian J. Reorganising the columns of matrix J (see Fig. 6), Eq. (42) can be expressed in a block matrix form This equation is underdetermined, with a 2(r + 1)dimensional solution space. Thus, the final solution segment w 1 can be given as a function of the initial solution segment w 0 where and the coefficient matrix M is the monodromy matrix. Note that if τ max > T , the solution segments w 1 and w 0 overlap and one part of the monodromy matrix is the identity. The stability of the periodic solutions can be decided with the help of the eigenvalues of the monodromy matrix M. A periodic solution is stable, if the spectral radius of the mondoromy matrix satisfies ρ(M) ≤ 1, where ρ(M) = max abs(λ) , for any Ms = λs. Note that for a periodic solution, 1 is always a trivial eigenvalue of the monodromy matrix. Asymptotic stability is ensured if the eigenvalue λ = 1 is multiplicity one while for the non-trivial eigenvaluesρ(M) < 1 holds, whereρ(M) = max abs(λ) , for any Ms =λs,λ = 1.
Theoretically, the boundary of the bistable domain(s) could be determined by the continuation of the saddlenode bifurcation (fold) of the periodic solutions. How-ever, the calculation of the monodromy matrix is computationally expensive while determining its critical eigenvalue is also an ill-conditioned problem due to the non-smooth nature of the system. Therefore, in our analysis, we determined the stability of the periodic solutions after the solution manifold is calculated.

Analytical study of the non-smooth Hopf bifurcation
One can also get an image of the stability of the periodic solutions by topologically extrapolating the properties of the local non-smooth Hopf bifurcation of the rectilinear motion ψ(t) ≡ 0, (t) ≡ 0 at the linear stability boundaries. This bifurcation can be studied analytically as described in detail in [3,4]. For this purpose, we consider only one sliding region in the x ∈ [−a, a] contact range at the vicinity of the rear edge x = −a. Thus, it is assumed that in the interval x ∈ [−x * , a] all particles stick to the ground, while the in x ∈ [−a, x * ] all points are sliding. It can be shown that this assumption is sufficient for studying small oscillations around the rectilinear motion. Consequently, the equation of motion (3) can be expressed in the following form: where the integral I 1 refers to the sticking region, where the travelling wave solution (6), (7) applies, while I 2 gives the effect of the sliding region, where the deformation is given by the limiting parabolae (8). Note that the boundary x * and the corresponding time delay τ * are state dependent, i.e. they are varying as the tyre oscillates. It is shown in [4] that, by performing centre manifold reduction in the non-smooth delay differential equation, (46) can be transformed to a normal form with smooth linear and piecewise-smooth second-order terms where μ and ω are the real and imaginary parts of the critical characteristic exponent, while h 2 (.) contains the non-smooth second-order terms. Analysing the normal form, one can determine the sub-or supercritical nature of the bifurcation, while it also provides the tangent lines of the arising branches of limit cycles at the non-hyperbolic equilibrium. The amplitude of the limit cycles around the critical equilibrium can be given as a function of the towing speed V by the formula where the parameter δ is equivalent to the Poincare-Lyapunov constant in smooth systems. That is, if δ < 0 the non-hyperbolic equilibrium is asymptotically stable due to the nonlinear terms, and the non-smooth Hopf bifurcation is supercritical and the emerging periodic solutions are orbitally asymptotically stable, and similarly, if δ > 0 the non-hyperbolic equilibrium is unstable, the bifurcation is subcritical and the limit cycles are unstable. This formula also shows that the periodic orbits emerge in a conical structure rather than in a paraboloid one that is characteristic in Hopf bifurcations of smooth systems [9].

Results
Before the nonlinear analysis, we divided the (V, l) parameter plane into linearly stable and unstable domains. Studying Eqs. (3) and (4), one can establish that lateral sliding does not contribute to the linear part of the equation. Thus, the linear stability of the rectilinear motion can be decided by semi-discretisation [7] of the delay differential equations and by the eigenvalues of the resultant coefficient matrix. We carried out the two parameter continuation of the periodic solutions, emerging from the linear stability boundaries by means of the towing speed V and the caster length l. For this analysis, we used the parameters given in Table 1. Figure 7 shows the vibration amplitudes ψ max of the periodic solutions in the undamped system (b t = 0). We also plotted the linear stability boundaries of the rectilinear motion in the (V, l) plane. One can show analytically that, as common in time-delay systems, an infinite number of stable and unstable domains exist, appearing in a chessboardlike structure. In our analysis, however, we only focus on the periodic solutions above the rightmost unstable domains. The nonlinear analysis reveals that the subsequent critical speeds correspond to sub-and supercritical bifurcations of the rectilinear motion in an alternating fashion. The critical nature of the stability boundary at l = a also becomes apparent. This line is found to be 'neutral', that is, neither sub-or supercritical by analysing the non-smooth second-order terms. In a sense, this boundary behaves similarly to a coincident Hopf bifurcation of the rectilinear motion and the fold bifurcation of the periodic solutions since it separates both stable and unstable equilibria and limit cycles. It can be seen that each surface of periodic solutions has a generic 'fold' line as well, corresponding to the fold of the periodic orbits. Accordingly, bistable parameter domains appear which correspond to the subcritical bifurcation of the rectilinear motion.
Thus, we can distinguish four types of parameter domains as shown if Fig. 8. In the white 'stable' regions, the rectilinear motion is linearly stable and has an infinitely large domain of attraction. In the 'unstable' regions (shaded in medium red), the rectilinear motion is unstable, whereas in the 'bistable' domains (shaded in light red) the rectilinear motion is stable but has a finite domain of attraction. That is, the system converges to a stable limit cycle for large enough perturbation. This structure is also visible in the oneparameter bifurcation diagrams (panels (a) and (b) of Fig. 8) where only the towing speed V is varied while the caster length l is kept fixed. In panel (b), it is also apparent that there exist a parameter range where the rectilinear motion is unstable, but the system has two stable limit cycles. This type of parameter domain is marked by dark red shading in panel (c). Figure 9 shows that if damping is considered in the rotational joint A, the chessboard-like structure of the linearly stable and unstable domains collapses into one D-shaped unstable parameter region. Our analysis revealed that the structure of periodic orbits is largely preserved even in this case; namely, a bistable domain occurs beside the stability boundary corresponding to subcritical bifurcation of the rectilinear motion. It is also visible that all the calculated branches of limit cycles emerge in a conical structure from the linear stability as predicted by the analytical calculations.
Unfortunately, when compared to numerical simulations, the calculation of the Floquet multipliers turned out to be unreliable as indicators of stability in certain cases. While they turned out to be robust at the solution surface at the highest speed, they tend to provide contradictory results at lower speeds, especially at the solutions close to linear stability boundaries. The suspected reason of this behaviour is that due to the non-smooth nature of the system, the continuation of the periodic solutions can only be carried out with relaxed error tolerances of the solution algorithm while, among others, the boundary between the sticking and sliding regions is found only approximately. Although these approximations are tolerable for calculating the periodic solutions, the Jacobian matrix will be effected by these inaccuracies as it is partly calculated by finite differences, involving the evaluation of the zero problem (30) with the variation of each state variable. As a conclusion, the periodic solutions are marked as stable/unstable based on the topology of the manifold using the fold line.

Measurements
Measurements were carried out on a towed bicycle wheel running on a treadmill. The experimental rig is shown in Fig. 10. The wheel is rotating in the fork which is mounted to the caster, such that the caster length l can be modified by a slider. The caster is mounted to the rack by a rotational joint. The deflection angle ψ(t) is measured by an angle sensor built on the joint. By means of the Galilean transformation, it can be showed that the scenario when the conveyor belt has a speed V is identical to the configuration when the caster wheel system is towed with speed V with respect to the ground. By varying the speed of the conveyor belt, we explore the linear stability boundary where the system spontaneously diverges from the rectilinear motion. Then, after the system converged to large-amplitude oscillation, we decrease the treadmill speed until the oscillations spontaneously decay. This speed corresponds to the fold of the periodic solutions which defines the boundary of a bistable parameter domain. This procedure is repeated for different caster lengths. Thus, the linear stability boundary and the bistable parameter range can be explored in the (V, l) parameter plane. In Fig. 11, the measurement results are compared to our calculations. It can be seen that our mechanical model and the measurements provided qualitatively identical diagrams and that the bistable parameter domains appear in a similar structure in the two nonlinear stability charts. Nevertheless, there are quantitative differences related to a known error of the brush tyre model: since this model neglects tyre deformation outside the contact region, the unstable and bistable domains appear at a somewhat lower speed domain than in the experiment. This error could be adjusted via the stretched string tyre model [12]. As shown in [1], this model provides qualitatively similar, but quantitatively better results for the critical towing speed. However, by means of the stretched string model, sliding can be considered only with differential equations instead of pre-defined sliding limits. As a result, one may find it challenging to perform a similar theoretical study with the stretched string model.

Conclusions
In our study, we investigated the nonlinear dynamics of wheel-shimmy. By performing the continuation of the periodic solutions in two bifurcation parameters, the towing speed and the caster length, beside the linearly stable and unstable regions, we were able to identify bistable parameter domains which are undetectable by using the standard quasi-steady-state tyre models with the one-degree-of freedom model of the caster wheel system. The local bifurcation of the rectilinear motion can also be investigated analytically. This reveals that, at the stability boundaries corresponding to the bistable domains, the straight-line motion undergoes a subcritical Hopf bifurcation. Additionally, it is also proven that the limit cycles emerge in a conical structure due to the non-smooth nature of the system. These results were also corroborated by numerical analysis.
By conducting experiments on a treadmill, we showed that the structure of the periodic solutions is qualitatively similar in the real-life system as it is predicted by the theoretical results. Nevertheless, there are quantitative differences between the theoretical and experimental results due to the neglected tyre relaxation (deformation outside the contact region). Moreover, obtaining the periodic solutions numerically has its own limitations, too. The continuation of the solution manifold, as well as determining their stability, requires the Jacobian matrix of the system equations. However, the calculation of the Jacobian is ill-conditioned due to the non-smooth nature of the system. Therefore, besides the Floquet multipliers, we also relied on the topology of the periodic solution surfaces and results of analytical calculations when we determined the stability of the solutions.
The part of the parameter domain where two nearby stability boundaries with different critical frequencies intersect each other is still relatively rarely identified. One can encounter quasi-periodic vibrations in those parameter domains where two complex conjugate pairs of roots of the characteristic equation have a positive real part. The analysis corresponding Hopf-Hopf bifurcations and the continuation of the quasi-periodic solutions is subject of future research.