Bifurcation analysis of wheel shimmy with non-smooth effects and time delay in the tyre–ground contact

The nonlinear dynamics of towed wheels is analysed with the help of the brush tyre model. The time delay in the tyre–ground contact and the non-smooth nature of the system caused by contact friction are considered simultaneously. Firstly, the centre manifold reduction is performed on the infinite-dimensional system transforming the governing equations into a normal form containing linear and piecewise-smooth second-order terms. Then, this normal form is used to establish the stability of the non-hyperbolic equilibria of the system and to give an estimation of the limit cycles emerging at the linear stability boundary. This way, it is demonstrated how subcritical Hopf bifurcations in the non-smooth delayed system generate bistable parameter ranges, which are left undetected by standard tyre models.


Introduction
The self-excited vibration of towed wheels is one of the oldest and most thoroughly studied phenomena in vehicle dynamics. This can be explained by the fact that in case of road vehicles, the contact between the vehicle and the road surface is realised by elastic tyres. Consequently, there is a great need to understand the mechanism behind the force generation of the dynamically deformed tyres in vehicle modelling [25]. This is a challenging task due to the fact that the deformed shape of the tyre travelling backwards in the contact region is subject to the effect of partial sticking and sliding governed by friction laws. This is the reason why this task is approached mainly by large-scale numerical studies often using complex tyre models such as the FTire [9], the RMOD-K [22] and finite element-based models [12,14], while limited analytical results are available.
Tyre models can be classified based on their complexity in terms of the physical representation of the tyre-ground contact. The simplest version of these models considers the wheel as a rigid body that has a single contact point at the ground [33,40]. These models are reasonable in case of the rigid wheels of shopping carts, baby strollers or rolling suitcases [8,23], and they provide a relatively simple formulation of the nonlinear governing equations. In these cases, one can consider various friction laws leading to a large variety of behaviours as shown for the example of a rolling rigid disc in [27].
A more sophisticated representation of the contact is possible by considering nonzero contact area between the interacting bodies. In such cases, when the shape of the contact area is well defined (e.g. elliptic [15]), these models were effectively used to study the dynamics of fundamental mechanical examples like the billiard ball or the Celtic stone (see [16,17]).
In case of tyres, the so-called physical tyre models, such as the brush and the stretched string models, are originated in a similar assumption: a contact region of finite area is considered and tyre deformation is represented with the help of spring elements distributed along the contact patch. These models are often simplified by assuming quasi-steady-state deformation. By introducing the side-slip angle calculated from vehicle kinematics, the lateral tyre force and the self-aligning moment characteristics can be derived [25]. This initiated the introduction of semi-empirical tyre models, among others the widely used Magic Formula [26], which capture the tyre force and moment characteristics with shape functions based on extensive experiments. If viscous damping is also considered in the system, these tyre models can provide even quantitatively good results for the linear stability of the rectilinear motion of a caster-wheel system.
However, there are essential discrepancies in the description of the nonlinear behaviour of the tyres in case of the semi-empirical and the physical tyre models. The source of the problem is related to the force characteristics derived by means of a Coulomb-like friction law that leads to non-smooth governing equations.
Note that a Coulomb-like friction law between two rigid bodies usually results in a so-called Filippov system with discontinuous governing equations. These systems may feature bifurcations characteristic only to them (see [30]), and they often exhibit chaos, as shown, for example, for the stick-slip phenomenon in [1,2].
In case of the tyre, however, friction has a direct effect on the lateral deformation, which is integrated along the contact length to calculate the lateral force and the self-aligning moment. Thus, we obtain piecewise-smooth continuous functions in the governing equations where the discontinuity appears in the nonlinear terms only. While this type of discontinuity in a higher derivative is clearly weaker than the discontinuity in case of Filippov systems, it still causes a significant difference compared to smooth systems in the way the limit cycles (self-excited vibrations) emerge at the linear stability boundary [19,29,30]. One difference is that while for a classical Hopf bifurcation in a smooth system, its sub-or supercriticality is decided by the second-and the third-order terms, then in the non-smooth case, the discontinuous second-order terms alone are decisive. This difference is also visible in the corresponding normal forms. If second-order terms are present in a smooth system (i.e. the nonlinearity is non-symmetric), one can eliminate them with a nonlinear near-identity transformation. On the contrary, this is not possible in the non-smooth case: the second-order terms already show odd symmetry and no unique nonlinear transformation exists that is capable to eliminate them in the normal form for both sides of the switching boundary. A further consequence is that with the piecewise-smooth characteristics, one can see a conical structure in the bifurcation diagrams as shown in [7], while using the smooth characteristics of the Magic Formula provides parabolae as expected for classical Hopf bifurcations [11]. Nevertheless, neither model can provide bistable parameter domains in the one degree of freedom model of towed wheels, which were otherwise experimentally presented in [5,24,39]. Additional dry friction at the king pin can be one explanation [24]. However, such parameter regions can also appear as a result of subcritical Hopf bifurcations leading to the coexistence of a stable rectilinear motion and stable self-excited oscillation of the towed wheel. This means that a so-called unsafe parameter region may exist where, depending on the initial conditions, either decaying vibrations or large amplitude oscillations can be experienced for the same set of parameters. Note that using multi-degrees of freedom models considering the compliance of the wheel suspension, subcritical Hopf bifurcations can be detected in narrow parameter regions even with the classical quasi-steady-state tyre models as shown in [41,42] for an aircraft landing gear.
If the assumption of the quasi-steady-state tyre deformation is relaxed and the lateral deformation in the contact zone is described by the governing partial differential equation (PDE), a travelling wave solution can be identified introducing the so-called delayed tyre model [36,38]. In the present paper, we use this dynamically varying tyre deformation in the contact region, which leads to the demonstration of a qualitatively correct nonlinear dynamical behaviour: a subcritical Hopf bifurcation of the straight-line motion is detected, and as a result, bistable parameter domains are explored. A further goal of our analysis is to present how the con-tact delay and dry friction contribute to the complex structure of the stable and unstable rectilinear and periodic motions in the system. With this aim, we perform the bifurcation analysis of the rectilinear motion of the non-smooth system [7] by studying the characteristic normal form [21]. In order to carry out closed-form analytical calculations, the simple brush tyre model is used, where the deformation at the sliding regions can be calculated directly. Although the brush model is less accurate compared to the stretched string models, the analytical results are still relevant either as test examples for advanced numerical codes or as qualitative examples for nonlinear tyre dynamics.
The rest of the paper is organised as follows. In Sect. 2, as a motivation of our study, experimental results are presented to demonstrate the appearance of bistable parameter domains and the strong effect of the tyre-ground contact friction on the dynamics of towed wheels. In Sect. 3, the mechanical model of a towed wheel of elastic tyre and its non-smooth equation of motion are presented. In Sect. 4, we expand the governing equations into a power-series form around the rectilinear motion up to the piecewise-smooth second-order terms. Then, in Sect. 5, critical parameters are presented as a result of linear stability analysis. At these critical parameters, Hopf bifurcations occur, which are analysed by means of centre manifold reduction and normal form transformation. The normal form is studied in Sect. 6 in order to determine the nature of the vibrations that occur around the rectilinear motion close to the critical parameters corresponding to loss of stability. In Sect. 7, the global dynamics of the system is presented by means of both analytical and numerical results, while the conclusions of the paper are summarised in Sect. 8.

Motivation
Measurements were carried out on a caster-wheel system running on a treadmill (see Fig. 1). The fork in which the bicycle wheel rotates is fixed to the caster which is mounted to the rack such that it can rotate around the king pin. It can be shown by means of the Galilean transformation that this configuration is identical to the scenario when the wheel is towed in a straight line with a given speed V . During the measurements, a bistable speed domain was explored by applying different force impulses on the caster. These impulses correspond to equivalent initial angular velocities Ω as shown in Fig. 2. Then, observing whether the system converges to stable large amplitude vibrations or to the stable rectilinear motion, we can estimate the angular velocity corresponding to the separatrix between the two stable solutions, which is used as an approximation of the amplitude of an (otherwise undetectable) unstable limit cycle (see red dashed line in Fig. 2). Performing this for different speeds V in the bistable region, the branch of the unstable limit cycles was continued until the point where the mechanical noise from the treadmill and the wheel unbalance resulted large amplitude vibrations without any additional impact hammer excitation. The root mean square (RMS) of the noise characteristic to the experiment is indicated by the dotted horizontal line in Fig. 2. The measurement results are compared with a theoretical bifurcation diagram showing a subcritical non-smooth Hopf bifurcation in the upper panel of Fig. 2. It can be observed that the unstable limit cycle appears to emerge in a conical structure rather than a paraboloid one that is characteristic to Hopf bifurcations in smooth dynamical systems. This indicates that, indeed, the non-smooth nature of the system, caused by dry friction in the tyreground contact, has an essential effect on the way the limit cycles develop.

Mechanical model
The nonlinear dynamics around a straight-line rolling of a towed wheel is analysed with the mechanical model shown in Fig. 3. The connection of the casterwheel system to the towing vehicle is realised by a rotational joint at A. The rigid caster is towed along the X -direction by a prescribed constant velocity of V . This can be described by the geometric constraints refer to the X and Y coordinates of the king pin A, respectively. The distance of the wheel centroid T from the joint A is denoted by l, while the distance of the centre of gravity C of the caster and the joint A is l C . The mass of the caster-wheel system is m, while its mass moment of inertia with respect to the centre of gravity C is J C . One can use the angle ψ as a generalised coordinate to determine the orientation of the caster in the (X, Y ) plane. The resultant of the lateral tyre deformation is considered in the form of the lateral force F and the self-aligning moment M, while the effects of The arrows indicate where the solutions converge to. Bottom panel: measurement results. Initial angular velocities (filled markers) and steady-state velocity amplitudes (empty markers) are plotted for different towing speeds. The cases when the system converged to the rectilinear motion are marked by circles, while those cases when persistent large amplitude vibrations occurred are marked with stars. The red dashed line estimates the likely locations of unstable limit cycles longitudinal deformation and longitudinal sliding are assumed to be negligible. This is a standard simplification for the analysis of wheel shimmy when the wheel is towed with constant speed [24]. While it is justified by the good correlation between the theoretical and experimental results, this assumption is not acceptable when the wheel is subject to strong acceleration/deceleration.

Tyre model
The deformation of the elastic tyre is described by means of the brush tyre model [25] shown in Fig. 4. It is assumed that tiny tread elements are distributed along the circumference of the wheel, which can be deformed normal to the wheel centre plane. The corresponding linear elasticity of these elements is charac-

Sticking region
Depending on the vertical load and the deformation of the tyre particles in contact, they can either stick to or slide on the ground. If the elastic elements are considered to stick to the ground, their velocity in the groundfixed (X, Y, Z ) coordinate system is zero, which can be expressed as the kinematic constraint As explained in details in [36,37], this can be used to derive the nonlinear PDĖ which describes the lateral deformation q(x, t) of the centre line. Prime refers to differentiation with respect to the space coordinate x attached to the caster. To this equation, the boundary condition q(a, t) = 0 is attached, where we use the standard assumption for brush models, namely that the leading point of the contact line always sticks to the ground with zero lateral deformation [25].
For this problem, one can construct the travelling wave solution: by introducing the time delay τ (x) distributed along the contact length x, whereas we also obtained a formula: which gives the position along the x axis. Clearly, in linear approximation, the solution of (4) is which is the time needed for a particle to travel from the leading position L to the actual position P at x.

Sliding region
To model the partial sliding of the tyre in the lateral direction, we consider a parabolic vertical force distribution in the contact region, which can be used to define two parabolic deformation boundaries due to the independent deformation of the tyre particles [25]. For a compact formulation, we introduce them as a] where N ± refers to the magnitude of the limiting parabolae, which also depends on the resultant vertical load F z on the tyre and the contact friction coefficient μ. If N ± is considered to be positive, the formula above provides the upper boundary, while a negative value gives the lower boundary.
In what follows, it is assumed that sliding occurs only at the rear end of the contact, which is a reasonable assumption for small vibrations around the rectilinear motion ψ(t) ≡ 0.

Governing equations
The equation of motion of the towed wheel can be obtained by means of Lagrange's equation of the second kind as where the mass moment of inertia of the caster-wheel system with respect to the king pin axis at A is J A = J C + ml 2 C , while b t is the torsional viscous damping factor in the joint.
Based on the brush tyre model, the resultant lateral tyre force F and the self-aligning moment M can be calculated by the integral formulae and Thus, substituting into Eq. (7), the governing equations of the mechanical system can be expressed aṡ where Ω is the angular velocity of the caster. The integral formulae that provide the lateral tyre force and the self-aligning moment are divided into two parts as defined in Eq. (11). I 1 corresponds to the sticking part of the contact where the travelling wave solution (3) applies. Thus, instead of coordinate x, the distributed time delay τ can be used for integration by substitution. This means that the limits are changed as x = a ⇔ τ = 0 and x = x * ⇔ τ = τ * , where x * and τ * refer to the actually unknown boundary between the sticking and sliding regions. The derivative dx/dτ used for the integral transformation is based on the travelling wave solution (4) and reads as The integral I 2 includes the lateral force and selfaligning moment generated in the sliding part of the contact. In case of sliding, the deformation is known from the parabolic limits in Eq. (6); this is why we keep the coordinate x as integration variable.
4 Power-series expansion of the tyre force and self-aligning moment

Sticking region
To carry out the local bifurcation analysis of the rectilinear motion, we have to expand the system into a power-series form up to the second-order terms. Note that for Hopf bifurcations in smooth systems, one has to consider the third-order terms as well while the secondorder terms can be eliminated by a near-identity transformation [18]. This cannot be performed in the nonsmooth case due to the piecewise-defined nature of the second-order terms, which are in this sense 'stronger' than the third-order terms as they already determine the vague stability of the non-hyperbolic equilibria of the system [19]. The rectilinear motion can be expressed as ψ(t) ≡ 0, Ω(t) ≡ 0 and q(x, t) ≡ 0. For the coordinate x, we get x = a − V τ since the whole contact region sticks to the ground in this case. This means that the 'boundary' between the sticking and sliding regions is exactly at the rear edge, i.e. x * = −a ⇔ τ * = 2a V =: τ * 0 . In the integral formula I 1 in Eq. (11), both the coordinate x and the lateral deformation q(x, t) are multiplied by the derivative dx/dτ . Consequently, in the non-smooth system, it is sufficient to consider them up to the linear terms: as follows from Eqs. (3), (4) and (12). Using the ansatz ϑ = −τ , the second-order approximation of the resultant moment at joint A from the tyre force F and the self-aligning moment M in the sticking region reads as where the time delay τ * at the boundary between the sticking and sliding regions is state dependent [28]: q(x, t)). After the elimination of q(x, t) by means of the travelling wave solution, this can be transformed into τ represent the infinite-dimensional state-space coordinates of the time-delayed system. Consequently, one can approximate the time delay τ * by the implicit integral formula which can be regarded as a generalised power-series expansion up to the linear terms. The coefficient functions f 1 (ϑ) and f 2 (ϑ) will be determined in the subsequent section.

Sticking-sliding boundary
The boundary of the sticking and sliding regions appears where the travelling wave solution (3) crosses one of the parabolic deformation boundaries in (6).
With the above-explained linear approach, this can be expressed as Substituting the linearised formula for coordinate x into Eq. (14) to obtain Keeping the linear terms only, this equation can be expanded as Comparing the coefficients of ψ(t + ϑ) and Ω(t + ϑ), the functions f 1 (ϑ) and f 2 (ϑ) can be expressed as where δ denotes the Dirac delta function. Thus, for the time delay τ * at the sticking-sliding boundary, we obtain the implicit formula Consequently, the explicit form of this power series can be expressed as where Based on this, also the integral formula I 1 can be expressed in the power-series form where the second subscripts refer to the order of the dependence on the infinite-dimensional state-space variables ψ t and Ω t . This calculation provides and

Sliding region
For the sliding region, the closed-form calculation of the integral I 2 leads to Then, the sticking-sliding boundary x * can be expanded similarly as it was performed for the time delay τ * : Substituting this result into Eq. (31), the integral formula I 2 is considered in the power-series form where where the approximation τ * ≈ τ * 0 = 2a/V is sufficient here for the time delay at the sticking-sliding boundary. This indicates that sliding does not have an effect on the linear stability of the system, which validates the former calculations in [36,37] where pure rolling is considered in the tyre-ground contact during the linear stability analysis.
The power-series expansion of the non-smooth governing Eqs. (2), (10)- (12) can be summarised in the form of the following system of nonlinear delay differential equations (DDEs): where τ * 1 is expressed in (26) and N ± can be found in (6).

Centre manifold reduction
In order to analyse the dynamics of the above-described non-smooth system, first the stability analysis of the rectilinear motion is carried out, which is followed by a detailed and rigorous Hopf bifurcation calculation.
Since the system of DDEs (36), (37) is defined in the infinite-dimensional state space of the continuous functions ψ t , Ω t , this system is transformed to the operator differential equation form (see [31]) where y t is defined by the shift y t = y(t + ϑ) for ϑ ∈ [−2a/V, 0], while y represents the vector of the caster angle and angular velocity: The operator A stands for the linear part, while F corresponds to the nonlinear part of the system. The linear operator A is defined as (see [10]) for any χ : [−2a/V, 0] → C 2 where the coefficient matrices A 0 and A τ are calculated form the linear part of the DDE (37): The nonlinear operator F in (38) is approximated till second order by where the vector function f 2 is calculated from the nonlinear part of the DDE (37) [see also (30), (35)] as . (44)

Stability analysis
The eigenvalue-eigenfunction problem for the linear operator A can be formulated as where λ ∈ C is the eigenvalue and φ is the right eigenfunction of the operator defined above [−2a/V, 0] ϑ. Based on (40), this equation can be expressed as the ordinary differential equation (ODE) with the boundary condition The general solution of ODE (46) has the exponential form where c ∈ C 2 is a constant vector. Substituting this into the boundary condition (47), we obtain This leads to the characteristic function D(λ) and the characteristic equation where I refers to the identity matrix while A 0 and A τ are given in (41) and (42), respectively. The characteristic Eq. (50) is identical to the one that is obtained after substituting the exponential trial solution into the system of DDEs (36), (37). After performing the integration and expanding the determinant, the characteristic function in (50) can be expressed as The characteristic equation can be used to study the linear stability of the rectilinear motion: all the real parts of the infinitely many characteristic exponents are negative in case of exponential stability. Based on this, a necessary stability condition l > −a can be derived. This critical caster length corresponds to saddle-node bifurcation, i.e. a real characteristic exponent crosses the imaginary axis at the stability boundary (λ = 0). In practice, however, the stability boundaries related to oscillatory loss of stability are the essential ones, that is, when a pure imaginary pair of characteristic exponents crosses the imaginary axis (λ = ±iω). These boundaries can be found by the D-subdivision method [38]. Calculating the real and imaginary parts of the characteristic function D(iω) provides the equations respectively. For the undamped system (b t = 0), the stability boundaries can be calculated analytically as shown also in [36,38]. On the one hand, l = a satisfies both (53) and (54) with the frequency of On the other hand, one can derive formulae for the critical towing speed V cr and caster length l cr by introducing the dimensionless frequency α := 2aω/V cr as a parameter in the formulae V cr = 4a 2 ak 3J A × 2α 2 + 2α(α 2 + 6) cos α + α(α 2 − 12) cos(2α) − 6(α 2 + 2) sin α + 6(α 2 − 1) sin(2α) Equations (55) and (56) provide the stability boundaries that divide the parameter space into linearly stable and unstable domains. The sense of the different parameter regions can be determined by methods presented in [32]. Stability charts constructed this way in the plane of the towing velocity V and the caster length l are shown in Figs. 6 and 9 in Sect. 7. In Fig. 6, the critical vibration frequencies ω in (57) at the stability boundaries are also presented.

Tangent space of centre manifold
The Hopf bifurcation calculation requires the determination of the tangent space of the centre manifold at the critical parameters corresponding to the oscillatory loss of stability. Accordingly, consider the case of the pure imaginary pair of characteristic exponents λ 1,2 = ±iω corresponding to the Hopf bifurcation of the rectilinear motion, where ω denotes the critical angular frequency of the system. Then, the characteristic exponent λ 1 = iω satisfies Eq. (45) with the eigenfunction φ 1 : Substituting it back into Eq. (50), we obtain from where the two eigenvectors in the ODE (46) can be expressed as and the corresponding eigenfunction of the linear part of the operator differential Eq. (38) is Note that φ 2 = φ 1 , where the overbar refers to complex-conjugate. The real and imaginary parts of the eigenfunction φ 1 provide a basis of the centre eigensubspace of the system. We represent them in the matrix function In order to reduce the system to the centre manifold, one has to calculate the adjoint basis as well using the left eigenfunctions of the linear operator A . The adequate eigenvalue problem can be expressed as where ψ is the left eigenfunction. Alternatively, we can use the adjoint operator A * as where A * is defined as for any ρ : [0, 2a/V ] → C 2 with A * 0 and A * τ as transposed matrices. Substituting into the eigenvalue problem in Eq. (64) at the critical characteristic exponent with boundary condition The general solution of this ODE for the left critical eigenfunction ψ 1 of the operator A can be expressed in the exponential form where d 1 ∈ C 2 is a constant vector. The boundary condition (67) leads to which has the solution and similarly, d 2 = d 1 . With the help of the real and imaginary parts of the eigenfunction ψ 1 , the adjoint basis is formed, which is represented by the matrix

Bilinear form
In order to set an orthonormal basis and adjoint basis of the tangent space of the centre manifold, we need the scalar product defined by the bilinear form [10] , where the new basis of the left eigenfunctions is introduced with the constant multiplier matrix K ∈ R 2×2 . This coefficient matrix is to be determined from the orthonormality condition Based on this, the inverse transpose of the matrix K can be given as The elements of this matrix are given as κ jk j, k ∈ {1, 2}, where With the help of Eq. (54), one can show that at the stability boundaries κ 22 = −κ 11 . Introducing the notation the constant matrix K in (73) that normalises the tangent basis of the centre manifold can be given as

Decomposition and reduction to centre manifold
Using the bases and , an arbitrary solution y t is decomposed as where u is the vector of the state variables in the centre manifold: while w t refers to the rest of the state variables in the infinite-dimensional complementer space. As a result, the functions of the deflection angle ψ t and the angular velocity Ω t are approximated in the centre manifold by harmonic functions: (ϑ)u(t) = cos ωϑ sin ωϑ −ω sin ωϑ ω cos ωϑ Thus, the reduced system can be expressed in the forṁ where w t depends on the second-order terms of u 1 and u 2 in the centre manifold. Since the operator F 2 includes second-degree terms only, the approximation w t ≈ 0 can be considered, which leads to with coefficients (89)

Second-order non-smooth normal form
Equation (85) obtained from the centre manifold reduction is the specific case of the general form where the real part μ of the critical eigenvalues is zero for the non-hyperbolic equilibrium [as given in (85)], and h 2 contains the non-smooth second-order terms. These can be given as where the coefficients c jk (0) and d jk (0) for j, k ∈ {1, 2} are given as The function H stands for the switching manifold, which is defined with the help of the limiting parabolae (6). Since we assumed that sliding occurs only at the rear section of the contact region, it follows that a transition from the upper parabola to the lower one can take place only at the rear edge R (at x = −a) as the sliding region shrinks to that single point. This means that for the time instantt of switching, the whole contact region is sticking to the ground, that is, x * = −a and q(x * ,t) = 0. Consequently, the travelling wave solution (3) related to the sticking region can be used at x = −a to express the switching condition. Since in (91) it is satisfactory to have its linear approximation only, the switching condition can be expressed with the linearised travelling wave solution (13) by projecting the state variables ψ and Ω =ψ to the centre manifold with u 1 and u 2 as given in Eq. (82): This expression is calculated at the stability boundary for the critical towing speed V cr when μ = 0 in (90). Equations (90) and (91) define the normal form of the local bifurcation of the rectilinear motion. This reduced system is thoroughly examined in [7], to which we rely on in the present study. The introduction of the polar coordinates r, ϕ with u 1 = r cos ϕ and u 2 = r sin ϕ in (90), (91) leads to the non-smooth ODĖ where the coefficient of the second-order terms can be expressed as f (ϕ; μ) := c 11 (μ) cos 3 ϕ + (c 12 (μ) while ϕ 0 refers to the orientation angle of the switching line H (see Fig. 5): For the angle ϕ, the polar coordinate transformation in (90), (91) leads to the approximationφ ≈ −ω.
Using the odd symmetry of the flow, one can compose a map of consecutive intersections of the switching line and a trajectory corresponding to an initial condition r 0 as shown in Fig. 5. The non-hyperbolic equilibrium at μ = 0 is stable if r k+1 < r k for any k ∈ N. It can be shown [7] that its stability is decided by the integral Namely, the equilibrium is stable (unstable) if δ < 0 (δ > 0). The substitution of (96) into the integral (98) leads to the Poincarè-Lyapunov constant δ of the nonsmooth Hopf bifurcation in the form where, by abuse of notation, c jk = c jk (0) and d jk = d jk (0). As in case of the classical Hopf bifurcation theorem [18], the stability of the non-hyperbolic equilibrium is topologically extrapolated to the arising limit cycles. The approximate amplitude r 0 for μ = 0 can be derived as the nonzero trivial (constant in time) solution of (90), which is approximated with the help of formula (98): This indicates that if the non-hyperbolic equilibrium is asymptotically stable due to the nonlinear terms (δ < 0), then the non-smooth Hopf bifurcation is supercritical and the emerging periodic solutions are orbitally asymptotically stable, and similarly, if the non-hyperbolic equilibrium is unstable (δ > 0), then the bifurcation is subcritical and the limit cycles are unstable. Formula (100) of the vibration amplitude also represents a common feature of piecewise-smooth systems [19]: the limit cycles emerge in a conical structure at non-smooth Hopf bifurcations instead of the paraboloid structure at classical Hopf bifurcations of smooth dynamical systems [18].
In practical examples, bifurcation parameters are selected, which help to visualise the structure of the emerging periodic solutions as a function of some relevant physical parameters. In case of the DDE model (36) and (37) of the shimmying wheel, a natural choice for this bifurcation parameter is the towing speed V since this parameter can be varied in the simplest way in experiments. The real part of the critical characteristic exponents λ 1,2 (V ) = μ(V ) ± iω(V ) depends on the towing speed V and μ(V cr ) = 0 at the stability boundary. The root tendency can be calculated by means of the implicit derivation of the characteristic function D(μ + iω) in (52). The bifurcation diagrams are constructed by transforming (100) with respect to the bifurcation parameter in the form

Results
Before investigating the nonlinear behaviour of the system, the linear stability analysis of the rectilinear motion has to be performed. Stability charts are created in the plane of the towing speed V and the caster length l. Figure 6 shows the stability chart and the critical frequencies (57) for the case when no damping is considered at the king pin (b t = 0). The rest of the system parameters are listed in Table 1.
The stability boundaries (55), (56), all corresponding to Hopf bifurcation, divide the (V, l) parame- ter plane into linearly stable and unstable parameter domains. One can use semi-discretisation [13] to establish the stability of the rectilinear motion for each domain. Thus, linearly stable and unstable parameter domains are obtained alternating in a chessboard-like structure. For large caster lengths, the unstable parameter boundaries shrink to asymptotes corresponding to l → ∞ which can be given asṼ In the undamped case, for V → 0, an infinite number of linear stability boundaries exists which behaviour is typical in delay differential equations. After establishing linear stability, the stability of the non-hyperbolic equilibria is investigated with the help of the non-smooth second-order terms as explained in Sect. 6. We found that by varying the towing speed  (10), (11). The thick curves are results of numerical collocation, while the thin straight lines represent the results of the analytical bifurcation calculation (102). Red and blue colours refer to unstable and stable branches, respectively  (10), (11). The thick curves are results of numerical collocation, while the thin straight lines represent the results of the analytical bifurcation calculation (102). Red and blue colours refer to unstable and stable branches, respectively V , the nonlinearly stable (shown in blue) and unstable (shown in red) non-hyperbolic equilibria also appear in an alternating structure in Fig. 6. The line l = a was found to be 'neutral' (δ = 0) for the second-order terms; therefore, its stability could be decided by the higher-order terms. Note that in previous studies [34], using a quasi-steady-state tyre model, this line proved to be unstable due to the third-order terms. There are also points in the (V, l) parameter plane where two linear stability boundaries intersect each other. At such points, two complex-conjugate pairs of characteristic exponents cross the imaginary axis resulting possible quasi-periodic in oscillations having two frequencies.
For an analytical investigation of these 'double-Hopf' points (see [20,35]) an extension of the presented algorithm is needed, which is beyond the scope of this study together with the numerical analysis of the further possible bifurcations of the periodic motions emerging from the presented primary bifurcations. As one can observe, for each linearly unstable domain, the left boundary relates to subcritical bifurcation, while the right boundary corresponds to supercritical bifurcation. Thus, at each unstable domain, both branches of the arising limit cycles are tilted to the left (i.e. to smaller velocities). Figures 7 and 8 show two examples of the structure of the periodic solutions for small (l = 0.02 m) and large (l = 0.2 m) caster length. In these bifurcation diagrams, the branches of limit cycles are calculated using numerical collocation [4,6] while their tangent lines at the bifurcation points are determined by considering the non-smooth second-order terms according to expression (102). It can be seen that the numerical collocation and the semi-analytical approximation are in good agreement with each other for small amplitudes in the vicinity of the bifurcation points. For both small and large caster lengths, bistable parameter domains occur where stable rectilinear motion and periodic oscillation coexist with domains of attraction separated by unstable limit cycles. This behaviour is already predicted by the second-order terms, though our semi-analytical calculation does not capture the numerically observed Fig. 9 Stability chart of the rectilinear motion with viscous damping considered at the king pin. The white domains are linearly stable, while the light red domains are unstable. The linear stability boundaries are coloured by means of their nonlinear stability considering the second-order terms. The blue curves are stable, the red ones are unstable, while the black one is neutral in this sense. The thin curves correspond to the linear stability boundaries of the undamped system saddle-node bifurcation of the periodic orbits (see the folds in Figs. 7 and 8 ).
As it is denoted in Fig. 7 for short caster (l = 0.02 m), the vibration amplitudes may reach π/2 in the unstable parameter region of V ≈ 0.35 m/s. At this point, the mathematical model has a singularity and it cannot describe those motions where the wheel is rotating in the opposite direction. In other cases, the amplitudes remain below 5 × 10 −3 rad, e.g. in the long caster case. While these so-called 'micro-shimmy' vibrations [39] are unlikely to cause structural vibrations and accidents, they are still undesirable due to the related energy loss, high tyre wear and noise generation.
If viscous damping b t is considered at the king pin, the stability boundaries can only be found numerically, for example, with the multidimensional bisection method [3]. In Fig. 9, we present how the stability boundaries change if an increasing value of damping is added to the system. As one would expect, the unstable parameter domains shrink and for large enough damping only the rightmost unstable domains survive. In the meantime, the nonlinear analysis revealed that the damped system largely preserves the structure of the nonlinearly stable and unstable boundaries. In Fig. 10, the bifurcation diagram is shown for l = 0.02 m and b t = 0.61 Nms where the system also features a bistable parameter domain.

Conclusions
In this paper, local bifurcation analysis was performed with respect to the rectilinear motion of a towed wheel. It was confirmed analytically that the simultaneous presence of time delay and partial sliding in the tyre-ground contact leads to an alternating structure of supercritical and subcritical bifurcations in the undamped system. If viscous damping at the king pin is considered, this structure breaks down to unstable islands in the parameter plane. Another, yet theoretically relevant, result is that discarding the side-slip in linear stability analysis is correct even if time delay is considered in the brush tyre model.
The experimental, analytical and numerical results all confirmed the existence of subcritical bifurcations of the rectilinear motion, which presents bistable parameter domains for certain towing speeds. These are undetected by simple quasi-steady-state tyre models. The results follow the rule of thumb that time delay tends to destabilise dynamical systems not only in terms of linear stability but also in the nonlinear domain by inducing subcritical bifurcations. Moreover, the results also underline why the consideration of the contact memory effect could give us a better understanding of the tyre deformation dynamics and help to develop more accurate tyre models.
The unique analytical calculations in this infinitedimensional and non-smooth problem were possible because it was shown that the evaluation of the secondorder terms only is satisfactory for the Hopf bifurcation analysis. In future work, this could be followed by considering the third-order terms, which would give an estimation for the saddle-node bifurcation of the periodic orbits; this defines the fold boundary of the bistable parameter domain. Another topic of further studies should be the application of the presented technique for the stretched string tyre model. As presented in [36] for linear stability analysis, the brush and the stretched string tyre models are qualitatively identical and it is expected that their behaviour in the nonlinear domain is similar, too. Although these future extensions cause further difficulties in both the analytical and the numerical calculations, the results are expected to be close to the experimental ones not just quantitatively but also qualitatively.