Bifurcation analysis of a rear axle tramp car model

Abstract

Axle tramp is a self-sustaining vibration in the driven axle of a vehicle with a beam axle layout, known to occur under heavy braking or acceleration. A 6DOF mathematical model of this phenomenon is used to identify how the key parameters of driveline stiffness, axle mass and fore/aft stiffness change the system’s dynamics. A bifurcation analysis is performed to study this nonlinear system’s dynamics. Four Hopf bifurcations in the underlying equilibria, along with a fold bifurcation in the outermost limit cycle branch, are shown to bound the parameter space where tramp occurs. The severity of tramp was found to be minimised by increasing drivetrain stiffness, reducing axle mass or increasing fore/aft stiffness: the trade-off for minimising tramp severity is that it may be easier to excite tramp when the drivetrain stiffness is increased, and the speed range over which tramp can occur is increased as fore/aft stiffness is increased. A key outcome from this work is that future electrified powertrains may experience more tramp, albeit at a reduced magnitude, than their combustion-powered counterparts.

1 Introduction

The nonlinear phenomenon known as axle tramp is the undesirable self-sustaining oscillation of a vehicle’s wheel and axle. The phenomena has been studied for almost a century, with early examples in literature highlighting the problem of tramp to be a major issue on vehicles with beam axles [1,2,3,4]. A beam axle, also called a live axle or solid axle, connects the wheels directly via a single beam, meaning any motion acts across the entirety of the assembly. This leads to two types of possible oscillation in the axle and wheels; symmetric oscillations where the whole assembly moves together in the same direction, and asymmetric where one wheel will move in an opposite direction to the other.

Historically, the issue of axle tramp was tackled in both front and rear axles, often alongside shimmy, wheel wobble and front-end shake. Over time, however, enough differentiation between these types of vibration was identified and the literature began to break the problems down into the more specific sub-areas. The differences between related key vibration problems are summarised in Table 1.

Table 1 Model parameters

Conclusions at this time state that the tramping frequency occurred at the same frequency as the rotation of the wheel and therefore the vibration is caused by torsional motions [4] and that tramp oscillations have larger displacements along the longitudinal axis [3]. Vincent [1] stated that tramp occurred in every vehicle produced, and that the problem could only be solved by independently springing the driven wheels. The topic regained traction in the 1960s as the problem was persistent in many vehicles of the time: as an example, see [5]. Ford had many problems with axle tramp on their Ford Cortina Lotus model that reoccurred throughout development. In August 1962, 1 month prior to the planned launch of the Mark 1 model, Lotus completely redesigned the rear axle layout, stating that the current leaf spring layout could not control the axle adequately and that axle tramp would be “inevitable”. The Mark 3 model suffered similar problems and following criticism from public and press, the suspension geometry was revised and redesigned with greater damping added for pitch on the 1974 model. Even with the improvements, at low speeds the rear axle could be described as “lively” [6]. Contemporary studies into the axle tramp phenomena on rear axle-driven vehicles are limited; however, when studied in detail novel approaches and results are still being found [7, 8]. Despite this, most of the recent examples of literature include research on the front axles of vehicles [9,10,11], tractors [12] and how driving conditions effect the system [13].

Although many modern passenger vehicles use independent suspension setups, beam axles are still used in commercial vehicles (such as pickup trucks) and cheaper passenger cars (such as the Ford Fiesta). With the move towards electrification within automotive engineering, there is an increasing industrial interest in how changes in the powertrain will influence the dynamic response of the vehicle. The work presented herein aims to address this interest by studying the dynamics that occur in a mathematical model of axle tramp using numerical continuation methods to identify and trace different types of dynamic features (such as equilibria, periodic orbits and bifurcations) through state–parameter space, under the variation of one or more model parameters. Previous bifurcation studies [14,15,16,17] have shown that numerical continuation is a useful tool for studying nonlinear dynamics in industrially relevant engineering systems, as it allows models developed for engineering purposes to be analysed directly rather than having to simplify a model to make it tractable with analytical nonlinear analysis methods. For the ordinary differential equations presented in this work, the continuation code AUTO, integrated into MATLAB via the Dynamical Systems Toolbox (DST) [18,19,20,21], is used to identify bifurcations in equilibria and periodic orbits, and trace these bifurcations to study how they change in terms of a second parameter of interest. This results in a preliminary analysis of the nonlinear dynamics of axle tramp. For instructions on how to use the toolbox, see [22]. The accuracy of the continuation runs, as defined by convergence criteria, was $1 \times 10^{-6}$.

The rest of this paper is arranged as follows: Sect. 2 presents the mathematical model used for this work; Sect. 3 presents a single parameter bifurcation analysis of the model with nominal parameter values; a variety of two-parameter bifurcation diagrams expand on this baseline result in Sects. 4; 5 finishes with some concluding remarks.

2 Live axle suspension model

Fig. 1

Example live axle suspension layout, with key components labelled

Figure 1 presents a schematic view of a live axle setup. The live axles are connected to the vehicle’s chassis via leaf spring suspension elements, which provide an element of both vertical and horizontal compliance. Because leaf springs do not possess any inherent damping, addition damper units are mounted between the axle and the vehicle’s body. Drive from the engine is transferred along a prop shaft, into the differential that is mounted in the centre of the axle: torque from the differential is then transferred to each wheel. The engine itself is mounted to the vehicle body via engine mounts, which provide stiffness and damping to isolate the vehicle body from engine vibrations. The nature of suspending the axle from the chassis means that the live axle is able to pitch relative to the chassis, so the prop shaft has a plunge joint to allow the ends of the shaft to slide relative to one another. This plunge joint has a stiffness associated with it.

Fig. 2

Freebody diagram of the system broken into two parts for readability. In the top-down view a the wheel, with radius R, is connected to axle with mass M. The axle moves longitudinally along X, resisted by stiffness and damping $k_{fa}$ and $c_{fa}$. The axle is connected to the engine mass $M_e$ via $k_p$, which can move along $\eta$. In b the rear view, the axle is shown to move along Z and is connected to the body,$M_b$, via k and c. The body can move along $Z_b$. In c, the axle is shown to pitch $\theta$, with dimensions a, b and e shown

Figure 2 represents the live axle arrangement from Fig. 1 in schematic form, with blocks used to represent the key rigid bodies of the unsprung mass (M—wheels, tyres and axle), sprung mass ($M_b$, the vehicle body) and engine mass ($M_e$). In the mathematical model used for this work, these degrees of freedom are captured via the X and Z displacements of the live axle’s centre of gravity, and the rotation of the axle from the horizontal plane ($\theta$)—see Fig. 2 for an illustration of these degrees of freedom. The axle is connected to the vehicle’s sprung mass via spring elements $k_fa$ and k (Fig 2). The engine is connected via the driveshaft to the differential (housed in the centre of the live axle), which provides additional resistance to the axle’s motion in the longitudinal direction. This resistance is included via a plunge stiffness ($k_p$) between the live axle and engine, and engine mount compliance between the engine and vehicle chassis ($k_e, c_e$). This approach is the same as existing models in literature [5].

From this point, several additions beyond those present in existing literature were made. The sprung mass represents the proportion of the vehicle’s mass carried by the rear axle, and is able to move vertically. Including a sprung mass within the vehicle model allows the vertical load on the tyres to be determined dynamically within the model, rather than set as a parameter within the tyre model at the start of a simulation. The wheel is rotated by the engine, via a gearbox and differential. Initial studies by the authors identified that the “shuffle mode” of the powertrain was a key frequency that drove tramp oscillations, so the torsional elements of the powertrain are represented as a single torsional spring connecting the input speed (as seen leaving the clutch) to the wheel. This spring stiffness was set such that the natural frequency of this simplified system matched that of the whole powertrain. The wheel connects the axle and sprung mass to the road, via an equivalent spring and damper in the vertical direction which represent the tyre’s vertical compliance with the road surface. A nonlinear function was used to relate the friction generated between the tyre and the road to the relative slip of the tyre.

With the above in mind, the equations of motion used to describe the axle’s movement are given as follows:

$$\begin{aligned} {\ddot{Z}}= & {} (c({\dot{Z}}-{\dot{Z}}_b)+k(Z_b-(Z-R)) -k_r(Z-R)\nonumber \\{} & {} -c_r{\dot{Z}}-Mg-c e\theta - k(b-a)\theta )/M \end{aligned}$$

(1)

$$\begin{aligned} {\ddot{X}}= & {} (-c_{fa}{\dot{X}}-k_{fa}X+F \nonumber \\{} & {} +(k_{fa}n+k_ph)\theta +k_ph\eta )/M \end{aligned}$$

(2)

$$\begin{aligned} \ddot{\theta }= & {} (-ce^2{\dot{\theta }} -(k(a^2+b^2)+k_{fa}n^2+k_{p}h^2)\theta -ce{\dot{Z}} \nonumber \\{} & {} -k(b-a)(Z-R)+(kn+k_{p}h)X-k_{p}h\eta )I_y\nonumber \\ \end{aligned}$$

(3)

Equations (1, 2 and 3) capture the vertical (Z), longitudinal (X) and rotational ($\theta$) motion of the axle. In these equations: k and C are vertical suspension stiffness and damping; $k_\textrm{fa}$ and $c_\textrm{fa}$ are the longitudinal stiffness and damping; $k_r$ and $c_r$ are vertical tyre stiffness and damping; M and $I_y$ are the unsprung mass and axle pitch inertia respectively; $k_p$ is the pinion stiffness; a, b, e, n and h are moment arms for terms $\frac{k}{2}$ (a and b, shown in Fig. 2), c, $k_\textrm{fa}$ and $k_p$ (for e,n and h respectively, not shown in figure to maintain clarity); g is acceleration due to gravity.

As described at the start of the section, the axle’s motion is influenced by the longitudinal position of the engine, $\eta$. To account for this influence, the following equation is added to describe the motion of the engine with respect to the vehicle’s body:

$$\begin{aligned} \ddot{\eta }=(-c_e{\dot{\eta }}-k_p\eta -k_e\eta +k_{p}X-k_ph\theta )/M_e \end{aligned}$$

(4)

Here, $M_e$ is the mass of the engine and $c_e$ and $k_e$ are the longitudinal engine stiffness and damping.

The influence of the vehicle’s sprung mass is captured as follows:

$$\begin{aligned} \ddot{Z_B}= & {} (-c({\dot{Z}}_b-{\dot{Z}}) \nonumber \\{} & {} -k(Z_b-(Z-R)-M_bg)/M_b \end{aligned}$$

(5)

where $M_b$ is the body mass.

The final equation of motion uses a single DOF to model the wheel’s torsional displacement, $\delta _x$, as a function of input speed, $\varOmega$. This captures the change in angular displacement of the wheel relative to the driveshaft’s input position, rather than an absolute angular rotation with respect to a stationary reference frame, to enable a steady-state to exist even when the wheel is rotating.

$$\begin{aligned} \delta \ddot{x}=(-c_1\delta {\dot{x}}-k_1 \delta x-F\ \cdot R)/i_y \end{aligned}$$

(6)

where $i_y$ is the referred wheel inertia and $k_1$ and $c_1$ are referred torsional stiffness and damping.

Fig. 3

Slip ratio, $\sigma$ vs normalised longitudinal force, $F^*$

Nonlinearities enter the system through the tyre’s ability to generate force. The force from the tyre (F), which acts at a distance R from the wheel’s centre, is obtained via a nonlinear function that captures typical tyre behaviour [23,24,25]. The equation that describes the force function is:

$$\begin{aligned} F=\dfrac{2\cdot \mu \cdot F_N \cdot \sigma _\textrm{peak}\cdot \sigma }{ \sigma ^2_\textrm{peak} + \sigma ^2} \end{aligned}$$

(7)

where $\mu$ is the coefficient of friction, $F_N$ is the vertical normal force. $\sigma$ is the slip ratio, with $\sigma _\textrm{peak}$ representing the value where the slip ratio is at its maximum.

Figure 3 shows how the normalised tyre force $F^*$ varies as a function of normalised tyre slip ratio. The slip ratio is the ratio of tyre rotational speed to longitudinal vehicle speed, defined by the following equation:

$$\begin{aligned} \sigma =(\delta {\dot{x}}+\varOmega )R-{\dot{X}} \end{aligned}$$

(8)

Here, R is the wheel radius, $\varOmega$ is the system’s input speed (equivalent to the speed of the driveshaft after the clutch), ${\dot{x}}$ is the wheel’s rotational speed deviation from the system’s input speed $\varOmega$ and ${\dot{X}}$ is the fore/aft speed of the axle and wheel assembly. The resulting slip value is expressed as a multiple of vehicle speed: in the subsequent analysis, the vehicle is assumed to be moving at a constant speed throughout the tramp response. Equation (8) therefore couples Eqs. (1, 2 , 3) and (6).

The normal force on the tyre is applied using the following equation.

$$\begin{aligned} F_N=-k_r(\delta x-R) \end{aligned}$$

(9)

To reflect the change in vertical force on the tyre when it leaves the ground, the tyre’s vertical stiffness $k_r$ becomes zero when the wheel leaves contact with the ground via the following piecewise function. Otherwise, stiffness value is taken from [5].

$$\begin{aligned} k_{r}= {\left\{ \begin{array}{ll} 0 &{} Z > R \\ 166371.6 &{} Z \le R \end{array}\right. } \end{aligned}$$

(10)

Table 2 Model parameters

Equations (1)–(9) describe the symmetric motion of a beam axle system.

In matrix form, the system can be expressed as;

$$\begin{aligned} {\dot{x}}=Ax+b \end{aligned}$$

(11)

Where

$$\begin{aligned} x= & {} [\delta x_1 \, \delta {\dot{x}}_1 \, \delta x_2 \, \delta {\dot{x}}_2 \,\ldots \, \nonumber \\{} & {} \ldots \, \delta x_{(n-1)} \, \delta {\dot{x}}_{(n-1)} \, \delta x_{(n)} \, \delta {\dot{x}}_{(n)}] \end{aligned}$$

(12)

where A is the stiffness and inertia system matrix, b is the nonlinear matrix, and I is the identity matrix.

$$\begin{aligned} A= \left[ \begin{matrix} 0 &{} 1 &{} 0 &{} 0&{} \ldots &{} \ldots &{} \ldots &{} \ldots &{} \ldots &{} 0 \\ -k_1/I_1 &{} 0 &{} k_2/I_2 &{} 0&{} \ldots &{} \ddots &{} \ddots &{} \ddots &{} \ddots &{} \vdots \\ 0 &{} 0 &{} 0 &{} 1 &{} 0 &{} \ldots &{} \ddots &{} \ddots &{} \ddots &{} \vdots \\ k_1/I_2 &{} 0 &{} -(k_2+k_1)/I_2 &{} 0 &{} k_2/I_2 &{} 0 &{} \ldots &{} \ddots &{} \ddots &{} \vdots \\ 0 &{} 0 &{} 0 &{} \ldots &{} \ddots &{} \ddots &{} \ddots &{} \ddots &{} \ddots &{} \vdots \\ \vdots &{} \ddots &{} \ddots &{} \ddots &{} \ddots &{} 0 &{} 0 &{} 1 &{} 0 &{} 0 \\ \vdots &{} \ddots &{} \ddots &{} \ddots &{} \ddots &{} k_{(n-2)}/I_{(n-1)} &{} 0 &{} -(k_{(n-1)}+k_{(n-2)})/I_{(n-1)} &{} 0 &{} k_{(n-1)}/I_{(n-1)} \\ \vdots &{} \ddots &{} \ddots &{} \ddots &{} \ddots &{} \ldots &{} 0 &{} 0&{} 0 &{} 1 \\ 0 &{} \ldots &{} \ldots &{} \ldots &{} \ldots &{} \ldots &{} 0 &{} k_{(n-1)}I_{n} &{} 0 &{} -(k_{n}+k_{(n-1)})I_n \end{matrix} \right] \end{aligned}$$

(13)

$$\begin{aligned} b=\begin{matrix} [0 \ 0 \ 0 \ F/M \ 0 \ 0 \ 0 \ 0 \ 0 \ 0 \ 0 \ -FR/i_y]' \end{matrix} \end{aligned}$$

(14)

The model is parameterised using values from previous work [5], as described in Table 2. The final model therefore takes the symmetric DOF from Sharp [5] with the following changes; the wheel is driven by input speed to ensure an equilibria can be obtained requiring a new torsional system, the wheel dynamics are updated with modern equations for longitudinal slip [11,12,13], the longitudinal engine dynamics from Sharp are combined in Eq. (4), the mass of the body is added as a DOF to account for any vertical motion in Eq. (5). The following section studies the bifurcations that are present in this model, with a view to explaining some of the physical mechanisms that drive tramp-like oscillations in automotive systems.

3 Bifurcation analysis of beam axle model

In this section, initial bifurcation results are presented to show the bifurcations that occur in the nominal parameterisation of the tramp model. Input speed, $\varOmega$, is used as the primary continuation parameter with relevant output states being plotted as a function of $\varOmega$. The bifurcation diagrams are shown in terms of X, Z and $\delta x$ as these are the prominent motions for axle tramp. Showing both X and Z is required to observe the wheel orbit. The graphical conventions used for the bifurcations are summarised below in Table 3.

Table 3 Graphical conventions for bifurcation diagrams

Fig. 4

One-parameter bifurcation diagrams for a Z, b X and c $\delta x$ as a function of input parameter clutch speed, $\varOmega$ with zoomed view provided

Figure 4 presents the bifurcation analysis for a longitudinal position Z, b vertical position X and c torsional displacement $\delta x$, all as functions of input speed, $\varOmega$. There are several features to note in the analysis. There is a single band of equilibria that is located at the datum position. This band is stable (blue) for higher speeds, and unstable (red dashed) at low speeds. The change in stability is identified with a Hopf bifurcation at $\varOmega$=16.6 rad/s. This subcritical Hopf bifurcation has a limit cycle associated with it that is either stable (green) or unstable (dashed cyan): the stability change occurs at a fold bifurcation. This set of limit cycles has a frequency of 1.86Hz (the same as the torsional eigenfrequency), and is associated with large amplitudes in the torsional displacement. There is a second Hopf bifurcation near zero which bounds these limit cycles. As the system approaches this second Hopf bifurcation point, nonlinear resonances are observed (see Fig. 4$a_{ii}$–$c_{ii}$). These peaks are most visible in the translational degrees of freedom of the model, suggesting that the torsional motion is undergoing some form of resonance with the wheel’s vertical and longitudinal motion. In this nonlinear resonance region, there is a second pair of Hopf bifurcations that occur at $\varOmega$=1.16 & $\varOmega$=4.85. These have a limit cycle frequency of 10.77 Hz, the same as the longitudinal eigenfrequency, and the amplitude of the limit cycles is the largest of all limit cycles when viewed in terms of longitudinal and vertical wheel displacements.

Regions of behaviour that are bound by key bifurcations emerge from the initial bifurcation diagrams. Region 1 is bound by the Hopf bifurcations near zero and $\varOmega$=16.6 rad/s, and contains unstable equilibria and a combination of stable and unstable limit cycles. The frequency of these oscillations can either be 1.86Hz or 10.77 Hz, with the 1.86 Hz oscillation occurring across this entire speed range and the 10.77 Hz oscillations being a subset of this region bound by the locations of the additional Hopf bifurcations. The magnitude of the oscillations can be determined from the size of the minimum and maximum amplitude shown on the diagram. Region 2 is bound by the Hopf bifurcation at $\varOmega$=16.6 rad/s and the fold bifurcation at $\varOmega$=39.2 rad/s. This region contains both dynamically stable equilibria and a combination of stable and unstable limit cycles (with a frequency of 1.86 Hz). Therefore, an axle that operates in this region may or may not tramp, depending on the initial conditions. Finally, region 3, which exists for any input speed higher than the corresponding speed of the fold bifurcation, only contains a dynamically stable equilibrium meaning any speed input would result in no limit cycle oscillations regardless of perturbation. The existence of these regions is similar to results obtained by the authors in previous work on brake vibrations [17], which contained a similar nonlinear friction force.

Fig. 5

A series of simulations showing the behaviour in each region

These regions can be verified by conducting simulations in each region. Figure 5 contains four different simulations showing the behaviour of X over time. In Fig. 5a, the input speed is set to $\varOmega$=10 m/s. This simulation is in region 1, and shows a self-sustaining oscillation as predicted. Figure 5b shows a simulation for $\varOmega$=50 m/s, which is in region 3 and shows the system decaying to rest at the equilibrium observed in the bifurcation diagram. Simulations in Fig. 5c and d are both set to $\varOmega$=25 m/s which is in region two where stable equilibria coexist with stable limit cycle. In Fig. 5c, the system observes no self-sustaining vibration. However, when the system is perturbed in Fig. 5d, a self-sustaining oscillations can be observed.

Fig. 6

Phase plots in the angular and longituinal DOF

Given the complicated limit cycle behaviour at low speeds, observing the limit cycles in phase space with the surrounding vector field can provide additional insight. Figure X contains a phase portrait at $\varOmega$=1.5 rad/s viewed in: (a) the $\delta X$ and $\delta {\dot{x}}$ plane; (b) the X and ${\dot{X}}$ plane. In Fig. 6a, the outer limit cycle corresponds to the torsional mode (at 1.86Hz). In this plane, the vector field appears to be influenced by this mode more than the longitudinal response (smaller limit cycle, bottom left of figure). In Figure (b), the vector field tends to head past the larger limit cycle (the longitudinal mode) towards the smaller limit cycle (which is the 1.86Hz torsional mode). This suggests that the torsional oscillation is more likely to be observed at low speeds than the longitudinal oscillation.

Fig. 7

Two-parameter bifurcation diagram showing change in location of bifurcations as input speed, $\varOmega$, and torsional stiffness, $k_1$ are varied simultaneously. Panel a shows Hopf and fold bifurcations in parameter space and panel b shows the $\delta x$ value of the fold bifurcation as a function of $k_1$. The background hatching shows the frequency observed

Fig. 8

One-parameter bifurcation diagrams showing the behaviour a Z, b X and c $\delta x$ as a function of input parameter clutch speed, $\varOmega$, in second gear, with zoom viewed added. For reference, the results for first gear are greyed out in the background

Fig. 9

Two-parameter bifurcation diagram showing change in location of bifurcations as input speed, $\varOmega$, and axle mass, M are varied simultaneously. Panel a shows Hopf and fold bifurcations in parameter space and panel b shows the $\delta x$ value of the fold bifurcation as a function of M. The background hatching shows the frequency observed

Fig. 10

One-parameter continuation in Z and X for a reduced axle mass of 15Kg (a$_i$,b$_i$) and an increased axle mass of 60Kg (c$_i$,d$_i$), with zoomed counterpart for each shown below in figures $(a-d)_{ii}$

Fig. 11

The initial two-parameter bifurcation diagram showing change in location of bifurcations as input speed, $\varOmega$, and longitudinal stiffness, $k_\textrm{fa}$ are varied simultaneously. Panel a shows Hopf and fold bifurcations in parameter space and panel b shows the $\delta x$ value of the fold bifurcation as a function of M. The background hatching shows the frequency observed

Fig. 12

The one-parameter continuation for (a$_i$) Z and (b$_i$) X for a reduced longitudinal stiffness value of $k_\textrm{fa}=1\times 10^5$ N/m, with zoom view a $_{ii}$ and b $_{ii}$ below

Fig. 13

One-parameter continuation for an increased longitudinal stiffness value of $k_\textrm{fa}=8\times 10^5$ N/m in a Z and b X with zoomed view provided

Fig. 14

The final two-parameter bifurcation diagram showing change in location of bifurcations as input speed, $\varOmega$, and longitudinal stiffness, $k_\textrm{fa}$ are varied simultaneously. Panel a shows Hopf and fold bifurcations in parameter space and panel b shows the $\delta x$ value of the fold bifurcation as a function of M. The background hatching shows the frequency observed

In light of the observations made above, the following section will explore how this model’s tramp response changes as key parameters are altered. The aim of this subsequent investigation is to identify how live axle designs could be improved to mitigate tramp from persisting, or minimise its effects if it were to occur.

4 Two-parameter bifurcation analysis

With a one-parameter study presented and the key bifurcations detected that indicate the regions of behaviour, several additional two-parameter studies can be conducted to determine how alterations in the model may affect these bifurcations and therefore the nature of axle tramp. Three scenarios are considered: changing the torsional stiffness (which captures the effects of changing gear); changing the mass of the beam axle; changing the stiffness between the beam axle and chassis. Three parameters are chosen: the first change can be instigated by the driver, as in the torsional stiffness study it is shown how the system will behave when the system undergoes a gear change; the changes to longitudinal stiffness and axle mass could be altered in the design phase and are chosen for study here as the system contains a large amount of longitudinal motion, directly related to the longitudinal stiffness; the axle mass is directly linked to the force applied on the wheel which likely influences the stick-slip behaviour involved in the axle tramp phenomena.

4.1 The influence of torsional stiffness

Figure 7 shows the bifurcation diagram in the parameter space a and state parameter space b in terms of the torsional stiffness, for the key bifurcations in this system. Following on from the analysis presented in the previous section, different regions are observable: two Hopf bifurcations are traced to show the origin of tramp oscillation at the longitudinal eigenfrequency; two additional Hopf bifurcations are traced to show the origin of the tramp that occurs at the torsional eigenfrequency; the key fold bifurcation, which occurs at high speeds, is traced to determine the speed boundary on tramp.

The fold bifurcation at higher speeds indicates the peak amplitude that occurs in torsional displacement during a tramp event. This bifurcation can therefore be shown in the state parameter space as a marker for the peak amplitude in $\delta x$. As shown in Fig. 7b, the maximum amplitude of the change in wheel position decreases as torsional stiffness increases, although it appears to be tending towards a limit. This result suggests that a higher torsional stiffness will reduce the size of the torsional component of tramp, but it does not eliminate the oscillations or change the speed range over which they occur.

Whilst the driveline’s stiffness could be set at the design stage, its effective stiffness will change during operation: changing gears alters the ratio between the wheels’ and engine’s displacement. As the torsional system uses a referred stiffness parameter (to capture the shuffle mode), a change in gear can be represented by a change in torsional stiffness. An eigenvalue analysis of a multi-degree-of-freedom torsional system created by the authors identified that the driveline shuffle mode increases to 2.37 Hz from 1.87 Hz as the system gores from first to second gear for the vehicle under consideration. This new torsional system can be reduced to a 1DOF system and applied in the same way outlined above in the modelling section.

Figure 8 shows a one-parameter bifurcation diagram of the system in 2nd gear. For comparison, the previous diagram in 1st gear is shown in grey. The bifurcations that bound regions of the phase space behaviour outlined earlier do not change, as the gearing alteration does not move the Hopf or fold bifurcation in relation to the x-axis. The frequency of the limit cycle oscillations is at the longitudinal and new torsional eigenfrequency, 10.77 Hz and 2.37 Hz respectively.

There are significant changes in the amplitude of the limit cycles in all DOF. In X and Z, there was a local peak in amplitude for $\varOmega \in$[10,20] rad/s (bounded by a pair of fold bifurcations). The gear change shifts this region to occur at the lower speed range, $\varOmega \in$[5,10] rad/s. Elsewhere, the resonance spikes still occur but at lower speeds, and at higher speeds similar limit cycle behaviour is observed but at a reduced amplitude. The amplitude of the oscillations reduces in $\delta x$ across the speed range, due to the change in torsional system resistance.

This result suggests that tramp could be more severe if provoked in 1st gear than 2nd, given the amplitude across the entire speed range is generally lower as the gear is increased from first to second; however, the change in the unstable limit cycle emanating from the highest speed Hopf bifurcation suggests that tramp could be easier to provoke in second gear: the unstable limit cycle branch will bound the stable limit cycle’s basin of attraction; as it occurs at a lower $\delta x$ value for a given $\varOmega$ (and similarly in other states too), smaller perturbations could tend towards the stable limit cycle.

4.2 The influence of axle mass

Figure 9 presents a two-parameter bifurcation diagram for axle mass, with panel (a) showing the two-parameter bifurcation runs in parameter space. The trend in the fold bifurcation indicates that, for a lower beam axle mass, the regions in which tramp will occur is reduced across the speed range. Furthermore, plotting the fold bifurcation (which represents peak torsional vibration) against axle mass in state–parameter space shows a small decrease in the torsional displacement when the system with a low axle mass tramps. This result has potential consequences for heavy vehicles, buses and trucks which are likely to have a beam axle setup involving heavy masses, suggesting that tramp in those vehicles could be exacerbated by their large axle mass.

Two discrete choices of axle mass can provide further insight into the tramp amplitude in Z and X. Reducing the axle mass to 15 Kg, shown in Fig. 10a,b demonstrates that a smaller axle mass reduces the tramp amplitude at low speeds (associated with the longitudinal mode), but increases the amplitude at higher speeds (associated with the torsional mode). Alongside this quantitative change in tramp, there is a qualitative change in the limit cycle branch, as a torus bifurcation present in the nominal mass case no longer exists. This qualitative change from the nominal case also occurs when the axle mass is increases: increasing the axle mass to 60kg (shown in Fig. 10c,d) also removes the torus bifurcation from the limit cycle branch. The change in amplitude of the limit cycles follows the observation made for the low mass to nominal mass comparison: at high speeds (associated with the torsional mode) tramp reduces; at low speeds (associated with the longitudinal mode) tramp increases.

4.3 The influence of longitudinal stiffness

Figure 11 presents an initial two parameter bifurcation analysis for longitudinal stiffness. In Fig. 11a, the two-parameter study traces the evolution the key fold and Hopf bifurcations and indicates regions of tramp and the potential frequency. In Fig. 11b, the state parameter space is shown for the torsional wheel displacement, where the fold bifurcation indicates the peak amplitude for this state. Regions are carved in the parameter space to show where the limit cycle occurs and at what frequency. The majority of the space contains limit cycles at 1.86 Hz (the same frequency as the torsional mode) but with a subsection that also contains an oscillation at the longitudinal mode, which changes as a function of $k_\textrm{fa}$. This region of coexistence is bound by two Hopf bifurcations, which now meet near $k_\textrm{fa}=2\times 10^5$ N/m. This result, unlike the other parameter variations, includes a parameter range where only a single tramp frequency is present, warranting further investigation.

Figure 12 contains the bifurcation diagram for a low longitudinal stiffness of $k=1 \times 10^5$ N/m. This places it in the region where the previous result suggested there was only a single tramp frequency, yet there still appears to be multiple limit cycle branches at low speed. There is a general growth in limit cycle amplitude in both Z and X, resulting in larger axle tramp oscillations in this system. The large cluster of bifurcations found that the default system also remains and the resonance spikes still occur; however, the increased amplitude magnifies the effect felt by the vehicle. The continued presence of large-amplitude limit cycles at low speed suggests that the previous result in Fig. 11 needs to be reconsidered at low stiffness values.

Figure 13 shows the one-parameter continuation for (a) Z and (b) X with zooms for each below, for the case of longitudinal stiffness increased to $k=8\times 10^5$ N/m. For reference, the original result is shown in grey. There are three significant changes to the limit cycle behaviour. First, the tramp oscillations that are at the torsional eigenfrequency in Z and X are smaller by an order of magnitude (compared to the nominal case). Second, the tramp at the longitudinal eigenfrequency is larger and the fold bifurcation associated with this branch is now located at a higher speed value. With the new fold location exceeding the location of the highest speed fold bifurcation, it is insufficient to use a single fold bifurcation as a boundary for tramp when considering longitudinal stiffness as a secondary continuation parameter. Finally, the increase in longitudinal stiffness has removed the fold bifurcations at low speeds and smoothed this branch, resulting in no resonance spike and a more comfortable transition when operating at low speeds.

In light of the results from the one-parameter continuation runs, two additional bifurcations need to be considered to correctly capture the range of tramp behaviour in a single diagram. The result obtained by tracing these additional bifurcations is shown in Fig. 14, which provides an overview of all the behaviour discovered. Notably, this increases the size of the parameter space in which oscillation at the longitudinal eigenfrequency can be observed: the second fold bifurcation acts a boundary on this behaviour for higher stiffness values, and the high amplitude oscillations at the longitudinal eigenfrequency are bound by new Hopf bifurcations at low stiffness values and the previous ones when the stiffness is higher. This result has implications for conventional tramp suppression approaches, which propose increasing the stiffness between the axle and chassis as a solution: at some point, increasing stiffness may broaden the range over which tramp can occur.

5 Concluding remarks

This paper presents a bifurcation analysis of a mathematical model of axle tramp. Nonlinearities enter the model through the tyre, where the friction force varies nonlinearly with longitudinal slip. Five key bifurcations are identified: two pairs of Hopf bifurcations lead to tramp oscillations associated with two different frequencies; one fold bifurcation bounds the maximum speed at which tramp can be observed. Using these key bifurcations, the influence of drivetrain stiffness, axle mass and axle fore/aft stiffness are considered. It is shown that the severity of tramp can be minimised by increasing drivetrain stiffness, reducing axle mass or increasing fore/aft stiffness. The trade-off for minimising tramp severity is that it may be easier to excite tramp when the drivetrain stiffness is increased, and the speed range over which tramp can occur is increased as fore/aft stiffness is increased.

The broad motivation for this work stems from the current trend in the automotive sector to electrify vehicles. Such a change tends to significantly increase the stiffness of the driveline, as gearing and driveshafts are reduced/removed. The outcome from this initial study suggests that such a change in stiffness could make tramp easier to provoke, although the resulting oscillations would be of a lower amplitude than in an equivalent combustion-driven vehicle.

The breadth of parameters considered in this work has discovered plenty of scope for future work to investigate certain aspects in more detail. With the current model, the existence and disappearance of a torus bifurcation under the action of axle mass warrants further study to identify the region over which this bifurcation exists, the dynamic behaviour that arises because of its existence, and the reason why it can be brought into/out of existence by changing the axle mass. The nonlinear resonances at low speed could also be studied to identify parameters that move the fold bifurcations in these limit cycles: knowledge of how these resonances “stiffen”, and potentially switch in nature, would be of academic interest. The influence of tyre model parameterisation on the nonlinear response could also be studied.

A semi-analytical study may provide different insights into the influence of model parameters on key bifurcation points. Such a study could adjust the model to capture tramp with a minimum number of states and parameters, and then explore the influence of this simplified model’s parameters on the location of the Hopf bifurcations in the model. This could be achieved either by linearising the tyre model and introducing its gradient and intercept as parameters that drive stability, or by analytically differentiating the tyre function and using a graphical approach to study the conditions that lead to zero crossings of the eigenvalues of the nonlinear system’s Jacobian matrix.

The model presented in this work assumes that the tyre’s response is linear up to the point it leaves the ground. In reality, the tyre’s stiffness will change with vertical displacement: this would be influenced by tyre construction, size and pressure. Future work could identify parameter sets that lead to oscillations where the wheel can leave the ground, and study the influence of the model used to capture this critical change in system dynamics. Finally, future experimental work could look to confirm the existence of some of the key behaviours identified in this model (such as the presence of a subcritical Hopf bifurcation, or two different tramp frequencies).