Towed wheel shimmy suppression through a nonlinear tuned vibration absorber

The implementation of the nonlinear tuned vibration absorber (NLTVA) for the suppression of shimmy vibration in towed wheels is addressed in this study. We adopt a modified straight tangent tyre model of a single-degree-of-freedom towed wheel system with an attached NLTVA. Stability analysis illustrated that the NLTVA can significantly improve the stability of the equilibrium of the wheel. Bifurcation analysis highlighted the existence of large bistable regions, which undermines the system’s safety. However, numerical continuation analysis, coupled with a dynamical integrity investigation, revealed that the addition of an intentional softening nonlinearity in the absorber restoring force characteristic enables the complete suppression of the bistable regions, also reducing the amplitude of shimmy oscillations in the unstable region. Quasiperiodic motions were also identified; however, their practical relevance seems marginal.


Introduction
Shimmy is a common name for the lateral-torsional vibrations of a towed wheel. It is one of the most common stability problems in vehicle dynamics, which affects several types of vehicles, such as motorcycles and bicycles [1][2][3][4][5][6], trailers [7][8][9][10][11] and aircraft landing gears [12][13][14][15][16]. The vibrations are generated by the dynamically varying forces in the tyre-road contact, and they can arise at various speeds [17]. Although the issue has already been studied for several decades [18], it is still a relevant problem considering the numerous accidents involving trailers undergoing shimmy vibrations [11].
Referring to aircraft, landing gears have several technical requirements, which make them prone to undergo shimmy vibrations. Namely, they must be light, since during the flight landing gear is just a dead weight; they must be long enough to allow a sufficient rotation of the aircraft during take-off; they should be able to absorb the energy of landing impact, which requires shock absorbers with a considerable stroke in order to limit the load occurring during landing [19]. All these requirements result in relatively slender landing gears, which are therefore prone to undergo undesired vibrations, if not properly designed [20][21][22].
The tuned mass damper (TMD), or dynamic vibration absorber, is a device consisting of a small mass attached to a host structure through a spring and a damper [23]. If properly tuned, the TMD can mitigate vibrations in various dynamical systems [24,25]. It is mainly implemented in civil structures, such as tall and slender buildings [26][27][28] and suspended bridges [29][30][31][32], where it suppresses wind-induced vibrations and prevents structural failure in the case of earthquakes. However, it is also adopted in other fields, such as for the mitigation of vibrations in wind turbines [33,34], elimination of the squeal in braking systems [35,36], stability improvement in boring bars for turning machining [37][38][39], train-generated vibrations suppression [40], mitigation of pitch and roll oscillations of floating bodies [41].
The TMD is appealing for shimmy oscillation suppression in landing gear because of its limited weight requirements, particularly if compared to other devices to suppress shimmy, such as shimmy dampers [42], bob masses [12,43] and inerters [1]. Several studies [12,44] tried to estimate the effectiveness of the TMD for suppressing shimmy oscillations in landing gears. These studies focused mainly on the system's stability analysis; indeed, they proved that a properly tuned TMD can significantly reduce the aircraft velocity range presenting shimmy oscillations. This aspect is reexamined in the present study.
Apart from stability, a number of studies [11,14,45,46] illustrated that nonlinear dynamical phenomena are particularly relevant for shimmy dynamics and, if overlooked, they might lead to unexpected behaviours, mining the safe operation of the device at hand.
This study aims to assess the performance of the so-called nonlinear tuned vibration absorber (NLTVA) [47] for suppressing shimmy instabilities. The NLTVA is a specific type of dynamic vibration absorber with linear properties analogous to those of a classical TMD. However, its nonlinear restoring force function also possesses a nonlinear characteristic, which, in general, can be arbitrarily chosen [24]. In this study, the linear and nonlinear behaviours of the system, encompassing a trailing wheel with an attached NLTVA, are studied separately. First, the absorber's linear coefficients are optimized through a stability analysis; then, an analysis of the system's nonlinear behaviour is exploited to optimize the absorber's nonlinear restoring force coef-

Mathematical model
The mechanical system under study consists of a suspension modelled as a one-degree-of-freedom (DoF) system, also known as the "trailing wheel model". As illustrated in Fig. 1, the system comprises a wheel mounted on a trailing arm. The arm can swivel about a vertical rotation axis, which moves at a given velocity V along a straight path. The model includes a yaw damping and stiffness, characterized by the linear coefficients c ψ and k ψ , respectively. Although this is one of the simplest models capable of exhibiting shimmy vibrations, it is sufficient to investigate the general efficiency of the NLTVA in the linear and nonlinear domains, which is the objective of this study.
The NLTVA is attached to the trailing wheel through a viscoelastic element, and it can rotate around the axis of rotation of the suspension. The viscoelastic connection generates a torque acting between the NLTVA and the trail, proportional to their relative angular position and velocity.
Different approaches exist for modelling the tyreground contact, having various levels of complexity, such as the Von Schlippe [48], the straight tangent [18], and the brush [49] tyre models. Intending to provide a general understanding of the performance of an NLTVA, we adopt the rather simple model provided by Pacejka's Magic Formula [46,[50][51][52], which considers the elasticity of the wheel assuming quasi-steady-state tyre deformation and evaluates nonlinear force and self-aligning-torque characteristics in a semiempirical way.
The dynamics of the mechanical system is modelled by the following system of differential equations: where ψ is the yaw angle of the wheel, ϕ is the rotation angle of the NLTVA and α is the tyre deformation angle, i.e. the angle between the velocity direction and the tyre orientation at ground level. F z and M z are the lateral force and moment acting on the wheel, which depend on α according to [50] where B, C, D, E, B t , C t , D t , E t are empirical parameters defined in [50,53]. Namely, B and B t are the stiffness factors, C and C t are the shape factors, D and D t are the peak values, and E and E t are the curvature factors. σ is the relaxation length, which is defined as σ = σ 0 exp −α 2 , as proposed in [53], in order to take into account nonlinear effects of the relaxation length, neglected in [50] due to linearization. t p is the pneumatic trail, which is the offset of the tyre lateral force with respect to the wheel centre in the longitudinal direction. c ψ and k ψ are the yaw damping and stiffness coefficients, respectively. I t is the yaw moment of inertia, calculated at the hinge as I t = I z + e 2 m, where e is the mechanical trail, m is the system mass, and I z is the yaw mass moment of inertia with respect to its centre of mass. a is half of the tyre contact length, V is the towing velocity, c a and k a are the absorber damping and stiffness coefficients, respectively, I a is the absorber mass moment of inertia and k nl is the absorber nonlin-  [12,53], and they are representative of an aircraft landing gear, as better detailed in [12]. The relation is assumed, which is valid for a wide range of twinwheeled landing gears [12]. The absorber mass moment of inertia is assumed to be 5% of I z , acknowledging practical constraints, while c a and k a must be tuned for optimizing absorber performance. k nl is the absorber's nonlinear stiffness, initially set at zero. We remark on some of the assumptions and limitations of the model adopted. Only motions in the horizontal plane and in a straight line are considered; the camber angle, castor angle, turn slip and longitudinal slip are all assumed equal to zero; the vertical load and the towing velocity are constant. Drag forces are neglected. Besides, other vibration modes, apart from yaw motion, which do affect the system stability [11], are also neglected. Regardless of these simplifications, the implemented model is able to capture the essence of shimmy vibrations and enables us to qualitatively and quantitatively evaluate the performance of the TMD, as detailed in the next sections.

Trailing wheel system: dynamical behaviour analysis
At first, we study the stability of the trivial solution of the trailing wheel without the vibration absorber. Accordingly, the system of equations is reduced to

Stability analysis
In order to study the system's trivial solution stability, we linearize the system around the trivial solution, reducing it to the formẋ = Ax, where The stability of the trivial solution can be directly studied by identifying the real part of the eigenvalues of A. If they are all negative, the solution is asymptotically stable; otherwise, it is unstable (excluding the case of non-hyperbolic equilibrium). However, an analytical expression of the stability chart can be obtained by applying the Routh-Hurwitz criterion to the characteristic polynomial of A, that is where The system is stable if a i > 0, for i = 0, . . . , 3, and H 2 = a 2 a 1 −a 3 a 0 > 0. Solving these inequalities with respect to V -and assuming parameter values compatible with a realistic engineering system-shows that the inequalities a i > 0 are always verified for V > 0 and e > 0. While H 2 > 0 provides the following conditions for stability: Stability diagrams in the V, e space (towing velocity and mechanical trail), obtained from these two inequalities, are depicted in Fig. 2. Figure 2a, b illustrate the effect of variations of c ψ and k ψ on the stability, respectively. In general, the trivial solution is unstable for low velocities, although it is stable for V very close to zero. If the mechanical trail e is larger than 1.6, the system is stable for any towing velocity V , regardless of the stiffness k ψ and damping c ψ coefficient values. For c ψ = 0, a critical mechanical trail value can be recognized, for which a larger range of towing velocity generates instability. For the parameter values considered, it is at e ≈ 0.7. As commonly observed in mechanical systems, increasing the damping reduces the unstable region. However, the effect is more significant at the higher boundary of the unstable region, while the lower boundary is only slightly affected by variations of c ψ .
Considering that a realistic landing gear has a nondimensional stiffness coefficient k ψ ≈ 15 [12], Fig. 2b shows that this stiffness provides an unstable area practically identical to a system with quasi-zero stiffness. In other words, the stiffness of the trailing wheel does not contribute to stabilizing the system. In order to improve the stability of the system by stiffening the landing gear, the stiffness coefficient should be increased by one or two orders of magnitudes, which is clearly impractical from an engineering perspective. Nevertheless, increasing the stiffness might reduce the amplitude of limit cycle oscillations (LCOs) in the case of instability, an aspect which is not investigated here.
For a quantitative analysis of the figure, we remind that the distance is normalized with respect to the wheel radius. Therefore, if we consider, for instance, a wheel with a radius of 30 cm, V = 100 corresponds to a velocity of 30 m/s. The analytical expression of the stability chart enables us to perform a much more sophisticated parameter analysis, which, however, is out of the scope of this paper.

Bifurcation analysis and global behaviour
An analysis of the system's eigenvalues reveals that pairs of complex conjugate eigenvalues have their real part becoming positive at the loss of stability. This implies the occurrence of Andranov-Hopf bifurcations. These bifurcations are then analysed utilizing the centre manifold reduction technique, as better explained in [54,55]. The main steps of the procedure consist of the transformation of the linear part of the system into its Jordan normal form, reduction of the dynamics of the system to a two-dimensional space through a centre manifold reduction, elimination of non-resonant terms through a near identity transformation and, finally, a transformation in polar coordinates. The intermediate passages of the procedure are omitted for the sake of brevity. All these steps provide the system's equations of motion in their normal polar form, i.e.
where V cr is the critical velocity at which the system loses stability, λ V is the derivative of the real part of the eigenvalues related to the stability loss, with respect to V (that is, the bifurcation parameter), and r is the amplitude of oscillation in the subspace of the centre manifold. ρ is the so-called Poincaré-Lyapunov constant which characterizes the bifurcation. The bifurcation is supercritical for ρ < 0 and subcritical for ρ > 0; for ρ = 0 the system undergoes a so-called Bautin bifurcation [54]. The analysis was performed along the whole stability boundary for the parameter values indicated in Table 1. Results illustrate that the Poincaré-Lyapunov constant ρ is negative along most of the stability boundary, except in a small region for very small velocities and large mechanical trail e (Fig. 3a). This implies that the bifurcations are mostly supercritical, except for that specific region. From an engineering perspective, this reduces the probability that stable solutions, different from the trivial one, exist in the stable region. This occurrence was verified by adopting the MATLAB toolbox DynIn [56], which aims at iteratively estimate the dynamical integrity [57] of the trivial equilibrium by performing several simulations. However, no stable solution, different from the trivial one, was found in the proximity of the unstable region. The computation revealed only a slight reduction in the dynamical integrity in the proximity of the subcritical bifurcations, which does not seem particularly relevant. Accordingly, the figure illustrating this result is omitted. Figure 3c illustrates the bifurcation diagram for mechanical trail e = 1; the figure was obtained through a numerical continuation technique via the MATLAB toolbox MatCont [58]. Apart from confirming the analytical results, the figure shows that the LCOs have a much larger amplitude for low than for high velocities.

System with attached NLTVA
The paper's main objective is to analyse the effectiveness of the NLTVA for suppressing shimmy vibrations in a towed wheel, which is studied in this section. The analysis is divided into two stages. First, we identify the absorber parameter values which minimize the unstable region, considering a given set of parameter values for the trailing wheel. This can be done directly by studying the system's stability in Eq. (1). Then, the system's global behaviour is investigated, utilizing a combination of numerical continuation, direct time integration, and local bifurcation analysis. This second step will enable us to set the nonlinear restoring force coefficient which minimizes the region of existence of LCOs and their amplitude. The same optimization strategy was successfully implemented for tuning the NLTVA in other applications [38,55,59].

Stability analysis
First, we linearise the equations of motion in (1) around the trivial solution, obtaining a system of the formẏ = By, where y = ψψ α ϕφ (11) and Since obtaining an explicit analytical expression of the stability boundary is extremely complicated-if not impossible-the system's stability is studied by numerically evaluating the real part of the eigenvalues of matrix B. If any of them is larger than zero, the trivial position of the system is unstable.
Four parameters characterize the NLTVA: its moment of inertia I a , its damping coefficient c a and its linear and cubic stiffness coefficients k a and k nl , respectively. It is well-known from the literature that the higher the absorber inertia, the better its performance. However, practical constraints usually limit its value, as is clearly the case for aeronautic applications. For this reason, we decided to fix I a at 5 % of the moment of inertia of the trailing wheel, i.e. I a = I z /20. The absorber's cubic stiffness k nl does not affect the stability of the trivial solution but only the nonlinear behaviour. Therefore, it is not involved in this first step of optimization. Accordingly, we are left with only two parameters to be optimized: k a and c a . Having to optimise only two parameter values, brute force optimization is computationally feasible. The optimization functions discussed below are computed for each couple of c a and k a values in a given range.
We optimise the parameters according to three different criteria. Namely, referring to the range of towing velocities which generates instability, we look for the set of absorber's parameter values which minimizes the maximal value of that velocity range (v max ), maximises the minimal value of the range (v min ), or reduces its extent ( v). The optimisation is performed for a fixed configuration of the primary system, corresponding to the values indicated in Table 1. As a reference, we notice that, without the absorber, the system is unstable for V ∈ [0.225, 89.71].
In order to give a better physical insight, we introduce the parameter γ = ω a /ω 1 = k a I t /(I a k ψ ), which represents the ratio of the NLTVA natural frequency divided by the yaw natural frequency of the trailing wheel. Clearly, it is equivalent to optimizing γ or k a . The result of the optimization is presented in Fig. 4. A red dot in the figure marks the optimal configuration for the three cases. Figure 4a shows that, despite the small inertia of the NLTVA, v max can be dramatically reduced from 89.71 to 5.478, with a reduction of almost 94 %. This improvement is obtained for k a = 0.4030 (γ = 1.270) and c a = 0.01393 (ζ a = c a / 2 √ I a k a = 0.06938). Interestingly, this result illustrates that the natural frequency of the NLTVA should be larger than the natural frequency of the trailing wheel. Besides, it is not required for the NLTVA to  Table 1 have a particularly large damping coefficient, simplifying its practical realization. Figure 4b depicts the minimal value of the unstable velocity range, which can be increased from 0.225 to 3.1816. In relative terms, this corresponds to an impressive improvement of 1314 %; in absolute terms, v min is increased by 2.96. The NLTVA parameter values required are k a = 0.3692 (γ = 1.215), c a = 0.01322 (ζ a = 0.06882); comparing these values with the previous case, the damping coefficient is practically identical, while the NLTVA's natural frequency is slightly smaller.
Finally, Fig. 4c depicts the extent of the velocity range generating instability v, which, for the case without absorber, is 89.48. Properly choosing the NLTVA parameters, v can be reduced to 4.74, with an improvement of 94.7 %. The optimal k a and c a parameter values are exactly the same as for the case of Fig. 4a for the considered resolution.
Stability charts corresponding to the optimal absorber configurations are represented in Fig. 5. The dashed line refers to the system without absorber. Despite the NLTVA being tuned for a specific mechanical trail (e = 1), the unstable region is significantly smaller also for other values of e. We also note that the unstable  Table 1 region is divided into two parts for both optimal cases. In Fig. 5a, v max is minimized by pushing the lower unstable region below e = 1, while in Fig. 5b, for maximizing v min , the upper unstable region is pushed above e = 1. In both figures, the direct comparison between the system with and without absorber clearly shows the absorber's effectiveness. Aiming to obtain the same effect by increasing the system's damping instead of using an absorber would require a huge damping coefficient, which is hardly realizable in practice.
In Fig. 6, the sensitivity of the NLTVA performance with respect to variations of its stiffness and damping coefficients is shown. The figures are obtained by selecting the optimal NLTVA parameter values according to the minimal v max and allowing variations of either γ or c a . As visible in Fig. 6a, variations of γ of 10 % increase v max from 5.478 to either 19.27 (increasing γ ) or to 32.5 (decreasing γ ). Conversely, as illustrated in Fig. 6b, increasing c a by 10 %, v max becomes 5.857, while decreasing it by 10 %, v max is 11.45. This result means that the absorber's performance is much more sensitive to variations of its natural frequency rather than its damping. From an engineering perspective, the realisation of the NLTVA is facilitated since accurately tuning the natural frequency is much easier than tuning the damping coefficient [60].
We note that, for the optimal absorber parameter values, the system is stable for e = 1 and V > 5.478; however, it is also very close to the unstable region for V ≈ 9, as visible in Fig. 5a. Accordingly, to have a safer design, it is convenient to select a stiffness value slightly larger than the optimal one to lower that boundary.

Bifurcation analysis
The value of the absorber nonlinear stiffness coefficient k nl is completely irrelevant with respect to the system's stability. Therefore, the parameter was not included in the optimization performed in the previous section. However, it can have an important effect on the system's global behaviour and on the bifurcations occurring at the loss of stability. Accordingly, a bifurcation analysis is performed to study the system's local nonlinear behaviour and optimize k nl .
The bifurcation analysis of the system with the attached NLTVA is limited to the optimal absorber parameter values for minimizing v max , i.e. k a = 0.4030 (γ = 1.270) and c a = 0.01393 (ζ a = c a / 2 √ I a k a = 0.06938). The mechanical trail e is fixed to 1 since the optimization was obtained for that value. The analysis was performed as for the case of the system without absorber in Sect. 3.2. All parameters are treated as constant numerical values, except for the velocity V , that is the bifurcation parameter, and the NLTVA cubic stiffness coefficient k nl .
The analysis provided the following Poincaré-Lyapunov coefficients for the lower and the upper stability boundaries:  Table 1 ρ lower = 0.0007302 + 0.02810 k nl (13) ρ upper = 0.01142 + 0.2583 k nl .
For a linear absorber (k nl = 0), ρ is positive both at the lower and upper limits, implying that the Andranov-Hopf bifurcations are subcritical. However, for both stability boundaries, a sufficiently large softening characteristic of the absorber's restoring force (k nl < 0) will make both bifurcations supercritical. In particular, for the stability boundary at low speed, it is required to have k nl < −0.02599, while for the one at high speed k nl < −0.04421. The accuracy of the analytical bifurcation analysis is verified by direct comparison with results obtained through numerical continuation, as illustrated in Fig. 7. Numerical results (thick black lines) were obtained via the MATLAB toolbox MatCont [58]. The numerical and analytical results display an excellent matching in the vicinity of the bifurcations, confirming both approaches' validity. In particular, we note that for k nl = −0.04421, the Andranov-Hopf bifurcation at higher speed has a transition between subcritical and supercritical characteristic-which corresponds to a Bautin bifurcation [54]-and the emerging branch of limit cycles is almost vertical (Fig. 7b). The corresponding curve computed numerically bends to the right, meaning that high-order terms lead to a subcritical behaviour for this limit case. A similar phenomenon, not illustrated here, is experienced for k nl < −0.02599 at the lower stability boundary.
From an engineering perspective, this result hints that, in the vicinity of the bifurcation, but in the stable region, the dynamical integrity of the trivial solution is most probably bounded. This might cause unexpected shimmy oscillations in parameter regions judged safe from the stability analysis, making it practically unsafe. The transition of the Andranov-Hopf bifurcation from subcritical to supercritical suggests that a softening characteristic of the absorber most probably increases the dynamical integrity of the trivial solution. However, this conjecture should be verified by analysing the system's global dynamics.

Global behaviour
The system's global dynamics is investigated by combining a numerical continuation analysis with direct estimation of the trivial solution's dynamical integrity.
The numerical continuation of the periodic solutions emerging at the Andranov-Hopf bifurcations was performed through the MatCont toolbox [58] by advancing the curves already shown in Fig. 7 for higher amplitudes. The result is depicted in Fig. 8. The figure illustrates how the branches of periodic solutions arising at the bifurcations merge, forming a unique branch. For k nl = 0, the periodic solution branch experiences two folds, one on the right of the unstable region and one on its left. This scenario, quite common for subcritical Andranov-Hopf bifurcations, implies that the dynamical integrity of the system's trivial solution is limited for velocity values between the folds and the Andranov-Hopf bifurcations, while between the two Andranov-Hopf bifurcations it is unstable. On the left-hand side of the unstable region, the fold is very close to the stability boundary (fold at V = 0.68, stability boundary at V = 0.75); conversely, the right fold is relatively far from the stability boundary (fold at V = 13.3, stability boundary at V = 5.48), making the unsafe region decidedly relevant from a practical perspective. Indeed, the unsafe region (between the stability boundaries and the fold bifurcations) is larger than the unstable region itself.
Imposing k nl < 0, the two folds move towards the centre, reducing the extent of the unsafe region. The left-hand side unsafe region disappears when the corresponding Andranov-Hopf bifurcation turns supercritical. This is not the case for the right-hand side unsafe region. The transition between subcritical and supercritical occurs for k nl = −0.04421, and for that value, the fold bifurcation is at V = 8.81. A continuation of the fold bifurcations (green line in Fig. 8) illustrates that the unsafe region disappears for k nl = −0.0574, i.e.
when it occurs at the same velocity of the stability loss. The fold completely disappears for k nl < −0.05926.
The numerical continuation also showed that, for k nl sufficiently small, the periodic solutions lose stability through a couple of Neimark-Sacker bifurcations, generating quasiperiodic solutions. The nature of the quasiperiodic solutions is verified in Fig. 9. The Poincaré sections of the steady-state solution display closed curves in the phase space, which is a clear mark of a quasiperiodic solution (red line in Fig. 9). The blue lines in Fig. 8a, representing the maximal ψ amplitude of the quasiperiodic solutions, were obtained through direct numerical simulations. Although this kind of motion is detrimental, the quasiperiodic motion's amplitude is not significantly larger than the one of the unstable periodic motions. Besides, quasiperiodic motions exist only for velocities within the unstable region, therefore already marked as unsafe from the stability analysis. Accordingly, these quasiperiodic motions do not seem particularly relevant from an engineering perspective. The red curves in Fig. 8 are the locus of the Neimark-Sacker bifurcations, also obtained through the MatCont toolbox. According to this curve, Neimark-Sacker bifurcations exist for k nl ≤ −0.0306.
We also remark that introducing a softening nonlinearity reduces the amplitude of the stable periodic solutions within the unstable region, which is practically desirable. However, softening nonlinearity generally has a destabilizing effect, as also suggested in this system by the appearance of Neimark-Sacker bifurcations. Additionally, nonlinearity can always generate unexpected dynamical phenomena, a condition that rarely can be ruled out a priori. Therefore, the absolute value of k nl should be kept as small as possible, providing that the unsafe region is eliminated or minimized.
In order to obtain a more comprehensive picture of the system's global dynamics, the dynamical integrity of the stable trivial solution was computed [57], this enables us to quantify the robustness of the equilibrium against external perturbations, and verify if other undetected solutions exist. For this purpose, the local integrity measure (LIM) was implemented, which corresponds to the radius of the largest hypersphere entirely included in the basin of attraction of the solution and centred in it [61]. The LIM was computed through the algorithm proposed in [56] and the relative MATLAB toolbox DynIn. Results are depicted in Fig. 10. Figure 10 confirms the results obtained with the bifurcation analysis and the numerical continuation, i.e. hardening nonlinearity (k nl > 0) is detrimental concerning the system dynamical integrity while softening nonlinearity (k nl < 0) is beneficial. Indeed, a relatively small softening nonlinearity can make the trivial solution practically globally stable (within the limit of validity of the mechanical model considered). Figure 10a, b show a similar trend of the limit of the unsafe regions on the left and the right of the unstable area.
According to our computation, the LIM slightly decreases for excessively large k nl absolute values. This effect is not related to the existence of other attractors but to the fact that, for large values of the differential rotation of the NLTVA |ψ − φ|, the restoring force becomes negative if the cubic term is negative, causing static losses of stability. This is a limitation of the math-ematical model adopted since, in general, it is unreasonable to devise an NLTVA with a negative restoring force. Accordingly, except for the dark regions, the system can be considered globally stable within the validity of the adopted mathematical model.

Conclusions
This study addressed the implementation of a dynamic vibration absorber for suppressing shimmy vibrations in a towed wheel, specifically focusing on a simplified model of an aircraft landing gear. Nevertheless, the model is highly general and can be applied to other applications as well [46]. The vibration absorber considered is the nonlinear tuned vibration absorber (NLTVA) [47], which is a nonlinear extension of the widely used tuned mass damper.
The analysis showed that the NLTVA can significantly improve the system's stability. However, bifurcation and global analysis also revealed the existence of unsafe regions where a stable trivial solution coexists with a stable periodic solution, and their extent is significant. This phenomenon mines the functional safety of the system and highlights that a stability analysis is insufficient to assess it.
An analysis of the system's dynamical integrity illustrated that, in these unsafe regions, the trivial solution's dynamical integrity is quite small. However, the analysis also revealed that introducing a softening nonlinearity in the absorber stiffness can reduce and eventually eliminate these unsafe regions, while also reducing the amplitude of the periodic solutions in the unstable region. Overall, the results obtained are promising towards the practical implementation of the NLTVA for shimmy suppression.
As a limitation of the approach used, it should be noted that the adopted tyre model may not always correctly capture the system's nonlinear dynamics, as illustrated in [45]. Therefore, it is necessary to validate the results through more sophisticated tyre models, such as the brush tyre model [52]. Additionally, analysing the NLTVA for shimmy suppression in rigid wheel models [62] would further enhance our understanding of the real potential of this technical solution.
Practical limitations of the system, such as the fabrication of a vibration absorber with a specific nonlinear characteristic [63][64][65], or the space and weight requirement of the NLTVA in comparison to other shimmy  Table 1; k a = 0.4030 and c a = 0.01393. The trivial solution is stable in all the range considered in the figure suppression solutions [12], were not analysed in this study. These aspects will be the subject of future studies. a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.