Using continuation analysis to identify shimmy-suppression devices for an aircraft main landing

This paper considers several passive shimmy-suppression devices for a dual-wheel main landing gear (MLG) and proposes a method of selecting the device parameter values for which no shimmy occurs. Two of these devices include an inerter, a novel mechanical element with the property that the applied force is proportional to the relative acceleration between its terminals. A nonlinear mathematical model is developed to represent the MLG dynamics. A bifurcation study is then carried out to investigate the eﬀects of the shimmy-suppression devices on the gear steady-state response. The aircraft forward speed and the device damping are chosen as the continuation parameters. A range of device parameter values that ensure the aircraft is free from shimmy instability for any forward speed within its operating region are identiﬁed. It is shown that the use of a proposed spring-damper conﬁguration can result in a more robust device in terms of the device damping over that of a conventional shimmy damper. Two inerter-based shimmy-suppression devices are then considered and yield further beneﬁts on expanding the zero-shimmy regions in the two-parameter bifurcation diagrams.


Introduction
The landing gear of any civil aircraft is required to be free from excessive vibrations and any dynamic instabilities over a conservative range of operating conditions [1]. A key source of such vibrations or instabilities is the phenomenon called shimmy. In the design process, the demand for suppressing shimmy instability may impose several design constraints on the structural stiffness 5 and geometry of landing gear [2]. However, if modifications to geometry, stiffness or weight are infeasible or undesirable, a shimmy damper is often introduced to alter the response [3]. Normally passive dampers are used to suppress shimmy oscillations. However alternatives are available. For example, it has been proposed that that the control orifice, present in some nose-gear hydraulic steering actuators, can be used to suppress shimmy [4]. Typically the shimmy damper is modelled 10 as a damping coefficient in parallel with the gear torsional stiffness [5]. Recently, studies have proposed semi-active or active shimmy-suppression strategies, such as using fuzzy adaptive control [6], magnetorheological damping [7], and sliding mode control [8]. However, while such semi-active or active controllers outperform passive ones, passive devices do have some advantages. They are typically simpler, requiring no power source, and are unconditionally stable. For example, under 15 some circumstance when the electrical power is lost, the powered active damping systems may fail to function. So current shimmy-suppression methods are typically still passive shimmy dampers [8].
Concentrating on passive solutions, studies have been reported in which the performance of various passive devices have been assessed in conjunction with linear gear model, see for example 20 [1,9,10]. In this paper, a method of selecting the network layout and parameter regions for shimmysuppression devices is proposed. This method ensures no sustained shimmy oscillations will occur over the aircraft operating velocity range. This approach is different from previous design methods in the literature and is applicable to different mechanical structures, for example, the lag damper for helicopters. 25 Firstly a passive device consisting of a linear spring and damper in parallel (which we term the shimmy damper ) is considered. The response of a dual-wheel MLG equipped with this shimmy damper is assessed. Then we investigate the effects of a proposed layout which adds a linear spring in series with the shimmy damper. In addition, configurations which include an inerter will be considered. The inerter is a commercially-available component, first proposed by Smith [11], that 30 generates a reaction force proportional to the relative acceleration between its two terminals. It completes the analogy between mechanical and electrical systems, allowing a wide range of passive absorber structures to be realised by mechanical networks. Performance advantages of suppression devices that include inertance have been identified for various systems, including vehicle suspensions [12,13], motorcycle steering systems [14], building vibration-suppression systems [15,16] and 35 railway suspensions [17,18]. The effects of the inerter on landing gear shimmy behaviour have been reported in [19,20], with [20] discussing the advantages of inerter-based shimmy-suppression configurations in terms of landing gear transient response.
The MLG steady-state response will be analysed using a nonlinear low-order model in this work. Here the tyre will be modelled using the exact stretched-string formulation [21], an extension of 40 the model proposed by Von Schlippe and Dietrich in [22]. For the representation of the landing gear structure, the torsional motion is a vital consideration in capturing the shimmy mechanism, see [23] for example. Since real landing gear systems exhibit various nonlinearities, the nonlinear dynamics of landing gear have attracted significant research interest. One approach to studying the gear's nonlinear response is to use a bifurcation analysis [24]. In [25,26], Thota et al. performed 45 a bifurcation study to investigate the effects of the geometric nonlinearity raised by a non-zero rake angle. They found that a lateral bending motion becomes coupled with the torsional motion. Further examples of bifurcation analysis applied in nonlinear systems can be found in [27][28][29][30] where [28-30] focused on aircraft shimmy analysis in particular. Other analytical analyses have been conducted in the literature to investigate the nonlinear aircraft shimmy problem, based on 50 the techniques of perturbation analysis [31] and the incremental harmonic balance method [32].
In this paper, we investigate the influences of the passive shimmy-suppression devices on the MLG steady-state response via continuation analysis. In Section 2, a nonlinear mathematical model of a typical MLG configuration reported in [28] is discussed. In Section 3, a passive device consisting of a linear spring and damper in parallel (the shimmy damper) is considered. A bifurcation 55 study is carried out using the continuation software AUTO [33], which is integrated into a Matlab environment via the Dynamical Systems Toolbox [34]. This allows us to identify the device parameter region in which no sustained shimmy oscillations occur over the entire operating speed range. This region in parameter space will be defined as the zero-shimmy region. A beneficial shimmy-suppression device with spring-damper layout is introduced, and its ability in expanding the zero-shimmy region are assessed. Based on this layout, two inerter-combined devices are proposed in Section 4. Their effects on the bifurcation diagrams are then studied and the performance advantages discussed. Finally, in Section 5 we draw some conclusions.

Main landing gear shimmy model
In this section, a typical dual-wheel MLG system reported in [28] is considered and the formulation of a low-order mathematical model of the MLG is presented. The MLG motion is modelled in terms of two degrees of freedom (DOFs); the gear torsional rotation and the lateral bending. The dynamics of the shimmy-suppression device and elastic tyres are also considered.

Dynamics of a MLG system
A sketch of the dual-wheel MLG is shown in Fig. 1  attached to the fuselage at the point F, as illustrated in Fig. 1(a). The main strut which is inclined to the Z axis by a non-zero rake angle φ is constructed of two cylinders or tubes. To keep the alignment of the wheels, a pair of torque links is employed, with the upper link attached to the upper strut cylinder and the lower one to the lower cylinder (piston), see Fig. 1(b). The end point of the piston is labelled C and the wheel axle is offset from the point C via a caster of length e.  As can be seen from Fig. 2(a), a body frame Bξηζ is used in order to describe the dynamics of the MLG system in the disturbed state. The axis ζ is rotated from the Z axis along the Y axis by the rake angle φ, and is aligned with the strut axis. The axis ξ is parallel with the caster while the axis η completes the right-handed coordinate system. With a radius l δ , the lower gear is allowed to bend laterally about point B along the ξ axis by the angle δ to represent the lateral compliance 85 of the gear. We consider k δ and c δ to capture the structural stiffness and damping of such lateral motion. The gear below point B has a center of gravity located at the point D. The wheel and axle assembly may rotate about the ζ axis with the torsion angle ψ, as shown in Fig. 2(b). This represents the rotational compliance of the torque links that span the upper and lower parts of the strut and are present to provide a rotational stiffness k ψ . A torsional damping c ψ is also introduced 90 to capture the MLG rotational damping. The MLG is allowed to move vertically and the vertical displacement of point B is z B . A constraint that the tyres must always contact the ground is assumed. Note that as in [28] we do not consider any axial deflection here. Hence, δ and ψ are the two MLG DOFs. Further, we consider a shimmy-suppression device fitted in the apex location between the upper and lower torque links, as shown in Fig. 2(a). Note that the conventional shimmy 95 damper is a translational device which generates force to control the relative movement of the two torque links. Such characteristics of the shimmy damper can be converted into equivalent torsional characteristics about the strut center line. All the following shimmy-suppression devices discussed have equivalent torsional properties. The DOF 'ε' across the device is introduced to represent the equivalent torsional motion of the shimmy-suppression device. The overall equivalent torque T ψ 100 generated by torsional damping, torque link stiffnesses and shimmy-suppression device is illustrated in Fig. 2(b). As presented in [1,9,20], the torsional stiffness k ψ is treated as a series connection between the suppression device and the MLG torsional motion. Such MLG torsional mechanism is illustrated in Fig. 3(a).
(a) (b) Figure 2: Schematic of (a) the MLG system in the disturbed state and the location of shimmy-suppression device, (b) the ψ degree of freedom in the Cξη plane. The toque T ψ represents the overall equivalent torque generated by torsional damping, torque link stiffness and shimmy-suppression device.
In this work the MLG mathematical model is established using the Lagrangian method. For each DOF, Lagrange's equation holds, where L is Lagrangian and L = T − U , T and U represent the kinetic and potential energy of the 105 MLG system, respectively. D is Rayleigh's dissipative function, Q i is the generalised force applied to the MLG system and q i is the generalised coordinate. When deriving the equations of motion for the system, z B is temporarily treated as a MLG DOF before being eliminated using a compatibility equation of the tyre-ground contact constraint.
The MLG kinetic energy is where m D is the mass of the lower gear (below B), J D is the inertia matrix at point D in the global frame, v D is the velocity vector of the point D, ω δ and ω ψ are the angular velocity vectors of the points D (lateral bending motion) and C (torsional motion), respectively. The potential energy is The ε DOF, representing the motion across the shimmy-suppression device, can be eliminated using a compatibility equation. The torque generated by the shimmy-suppression device T d is balanced by the torsional spring torque k ψ (ψ − ε), as illustrated in Fig. 3(a). The expression for T d is dependent on the device layout. For the shimmy damper layout (labelled L1) shown in Fig. 3(b), the compatibility equation is given by where k d and c d denote the stiffness and damping of the device. It should be noted that k d is included and fixed at a default value for all the shimmy-suppression device layouts considered. This is to ensure the gear has sufficient rotational stiffness for centering the gear. The default value of k d was determined via a scaling calculation based on [9], which is 1.09 × 10 5 Nm/rad. The Rayleigh's dissipative function is given by where c δ and c ψ are the damping coefficients of the respective DOF. Now we need to determine the exact expression of T based on Eq. (2). We consider all the terms in the global frame, so it is desirable to define the rotation matrix transforming the vectors from the body frame to the global one. As the global frame can be transformed to the body frame via three rotations φ-δ-ψ in sequence, this matrix is given by The position and velocity vectors of point B in the global frame are given by where V is the aircraft forward speed, t is the time, L A is the MLG height, and l B is the distance from points A to B. The position vector of point D in the body frame is given by where l D is the distance from B to D. Then, we obtain the position and velocity vectors of point D in the global frame by Due to the sequenced rotations applied from global to body frames, the two angular velocity vectors may be written as The inertia matrix tensor is given by where J b DW W represents the moment of inertia at point D with respect to W axis of the body frame. Hence, the inertia matrix tensor in the global frame can be given by Fig. 4 shows the external forces applied to the MLG system. There is a vertical force F N acting on the point B and the gravitational force at point D, giving 6 where F N = M g, M is the mass of the fuselage and the upper part of the MLG (that above point B), and g is the gravitational constant. At the ground-tyre contact points, there are vertical forces F zi and lateral forces F yi , where i = L or R represents the force applied to the left or right tyre. These are given by where F zi is the magnitude of F zi , θ = ψ cos φ cos δ is the angle between the wheel actual travelling direction and the X axis. The expression of the coefficient Λ i is dependent on the tyre model and will be given in Section 2.2. The vertical forces are given by where k t is the tyre vertical stiffness, F z is the total vertical force acting on the MLG from the ground, and a is the track width, as illustrated in Fig. 4. The angle µ is given by Apart from the forces, the tyres also experience self-aligning moments M ki , as given by where C ki will be given in Section 2.2. The velocity vectors at ground-tyre contact points E i are also needed for the calculation of the generalised forces. The positions of E i in the body frame are given by where R i is loaded radius of the left/right tyre. R i can be expressed by where R is the tyre unloaded radius, d i = ∓ 1 2 aµ, which represents the tyre deflections (see Daugherty [35]). Then, the position and velocity vectors can be calculated by Thus, the generalised force for each coordinate is calculated by This completes the set of terms needed for Eq. (1).

115
As previously mentioned, the ground-contact constraint ensures the MLG always contacts the ground. Therefore, the velocities of the contact points in the vertical direction should equal 0, where v Ei.Z andṙ BEi.Z are the Z components of v Ei andṙ G BEi (r G BEi = Hr b BEi ), respectively. Therefore, to satisfy this constraint for the two contact points, we havė (a) (b) Figure 4: View of (a) the general forces applied to the MLG, (b) the vertical forces acting at the ground contact points of the two wheels.

Tyre model
The exact stretched-string model detailed in [21] is used here to capture the tyre dynamics, and will be discussed in two parts. The first part involves the forces applied to the MLG, i.e., the lateral forces F yi and the self-aligning moments M ki . According to the tyre model, to capture these forces, the definitions of coefficients Λ i and C ki need to be provided, written as Λ i = k λ tan −1 (7 tan α i cos(0.95 tan −1 (7 tan α i ))), (25) and where k λ is the tyre restoring coefficient, α i = tan −1 ( λi l ) is the tyre slip angle and l is the tyre relaxation length. In Eq. (26), k α is the tyre self-aligning coefficient and α m is the maximum slip angle. These definitions are taken from [26].
The second part relates to the tyre dynamics at the ground-tyre contact plane. As shown in Fig. 5, the local tyre frame E i xy is used and the actual shape of the tyre contact line is defined by a function P i (x, t), x ∈ [−h, h], where 2h is the contact patch length. Consider a point P i between the leading and trailing points (P 1 and P 2 , respectively) of the contact patch, it has a 8 lateral deformation y i (x, t) along y axis. Then the absolute position of the point P i is given by where r Ei.X , r Ei.Y and r Ei.Z denote the X, Y and Z components of the vector r G Ei . Then its 120 velocity vector is given by where v Ei.X , v Ei.Y and v Ei.Z denote the X, Y and Z components of the vector v G Ei . Assuming the tyres are fully adhered to the ground leads to a zero velocity condition at the point P i , whereby each component of the vector v pi equals 0. Further algebraic manipulation giveṡ The slope at P 1i is related to the lateral deformation of the leading point λ i by the approximation Substitution of Eq. (30) into Eq. (29) giveṡ Figure 5: The stretched-string tyre model and the tyre deformation λ i .

Equations of motion
The equations of motion of the MLG system for q i = δ, ψ and z B are now written in the compact form Via necessary substitutions, we obtain wherez Bi is obtained via differentiation of the two constraints Eqs. (23) and (24). The expression of the shimmy-suppression dynamics, given by Eq. (4), is specific to the type of the device considered.

125
Hence, the three second-order differential equations Eqs. (32)-(34), the first order differential equation Eq. (31) for tyres, along with Eq. (4), complete the equations of motion for the MLG system. The coefficients N b aa are detailed in the Appendix. The parameter values shown in Table 1 are used for the MLG system and are mostly taken from [28]. The stiffness for the shimmy-suppression device, k d , is calculated via the scaling analysis based on [9]. Note that all the parameters in Table 1 130 are fixed in the following discussions.

Analysis of non-inerter shimmy-suppression devices
Using the model presented in Section 2, two shimmy-suppression devices are investigated in this section, the first of which is the traditional shimmy-suppression device, the shimmy damper, shown in Fig. 3(b). Its effects on the MLG steady-state response are studied via continuation.

135
The continuation approach allows the region in parameter space in which no steady-state shimmy oscillation solutions exist at any velocity within the operating range, the zero-shimmy region, to be identified. This is achieved via two-parameter bifurcation analysis in terms of forward velocity V and shimmy-suppression device damping c d . The second proposed device comprises a springdamper layout and the influence of this configuration on expanding the zero-shimmy region is 140 explored. This zero-shimmy region, in which no sustained oscillations exist, has previously been studied in terms of operating conditions, see for example [36,37]. Here it is considered also in terms of device parameter values. We note that establishing such zero-shimmy parameter region is the first design stage for shimmy-suppression device. Once it has been identified that no steady-state shimmy solutions exist, the emphasis of the device design shifts to its second stage, which is to 145 minimise the maximum transient shimmy deflections due to perturbation, see for example [8,20]. However, it is the first phase that is the focus of the work presented here.

The default shimmy damper
In this subsection, the effects of a shimmy damper are analysed. The device stiffness k d is fixed to 1.09 × 10 5 Nm/rad and the influence of the device damping c d on the MLG steady-state dynamic behaviour is investigated. First, choosing the aircraft forward speed as the continuation parameter while setting c d to three example values, one-parameter bifurcation diagrams are constructed showing both the MLG steady-state zero-amplitude and periodic solutions within a certain range of speeds. After that, we treat c d as a second continuation parameter and the resulting two-parameter bifurcation diagram is given in the parameter plane of the MLG forward speed and c d . Regions 155 in which different types of shimmy oscillations occur as well as stable zero-shimmy behaviour are identified in this parameter space.

One-parameter continuation
A bifurcation study is conducted here to investigate the MLG steay-state solutions; see [38] as an example for details of bifurcation theory and [39] as a review of applications of bifurcation analysis for aircraft dynamics. The continuation software AUTO [33], which has been integrated into Matlab via the Dynamical Systems Toolbox [34], is used in the analysis. A one-parameter bifurcation analysis is now performed with the aircraft forward speed V as the continuation parameter. Fig. 6(a) shows with respect to V the maximum amplitudes of MLG steady-state solutions in terms of torsional motion ψ and lateral bending motion δ * when c d = 3.6 × 10 4 Nms/rad. Here, δ * represents the resulting deflection of δ DOF at the point C. The stable solutions are represented by blue solid lines and unstable by red dashed ones. At low forward speeds, the straight-rolling solution is observed, which represents zero shimmy. When V is 13.4 m/s, the zero-shimmy solution loses stability and shimmy occurs. This qualitative change in behaviour is due to a Hopf bifurcation (labelled H) corresponding to a pair of complex conjugate eigenvalues 170 crossing the imaginary axis in the complex plane; this gives birth to a periodic solution. This periodic branch ends with another Hopf point at V = 33.8 m/s. Along this branch, stability is changed at V = 30.8 m/s, with the existence of a torus bifurcation point (labelled T), resulting from a pair of complex eigenvalues with unit modulus. Note that for V = 20 m/s (point A in Fig. 6(a)) the maximum deflection angle of ψ is 3.7 degrees while the maximum amplitude of δ * is only of the order of 10 −4 m. We define this solution as a kind of torsional-dominated shimmy (we refer to this as ψ-shimmy). It can be observed that along this branch the maximum deflection angle of ψ is always more significant than that of δ * . Hence this branch represents ψ-shimmy. Moreover, there is another branch of periodic solutions which is connected by a further pair of Hopf bifurcation points. Following this branch, initially the solution is unstable but regains stability when a torus point is  Fig. 6(a)) as an example, and hence is referred to as δ-shimmy. Note that there is a velocity region bounded by the two torus bifurcations, where both branches are stable; within this bistable region the initial perturbation determines which solution branch will be observed. Fig. 6(b) and (c) shows one-parameter bifurcation diagrams for 185 smaller device damping: (b1) and (b2) for c d = 5.0 × 10 3 Nms/rad; (c1) and (c2) for c d = 2.0 × 10 3 Nms/rad. In contrast to the results shown in Fig. 6(a), for both cases, only one branch is observed and along this branch the solutions are stable. Since the δ-component is more significant in Fig. 6(b), the branch is δ-shimmy while that shown in (c) is a ψ-shimmy branch.

190
As observed in Fig. 6, the device damping c d plays a significant role in the MLG dynamic behaviour. Hence, in order to investigate the influence of the device damping c d , we choose c d as the second continuation parameter and construct a two-dimensional bifurcation diagram in the (V , c d )-plane. This is illustrated in Fig. 7. Fig. 6 represents the three slices c d = 3.6 × 10 4 , 5.0 × 10 3 and 2.0 × 10 3 Nms/rad in the two-dimensional plane. The red curves in Fig. 7 are Hopf curves which are 195 formed via the continuation of the Hopf point H in the one-dimensional plane. In the (V , c d )-plane, we observe two Hopf curves intersected with a pair of Hopf-Hopf points (labelled HH). Two curves of torus bifurcations emerge from the two HH points. The third Hopf curve is observed in the lower area of the (V , c d )-plane. As discussed in the one-parameter bifurcation analysis, the existence of Hopf points represents the onset of shimmy oscillations. Similarly, in Fig. 7, these Hopf curves 200 bound regions of different shimmy behaviour. In particular, we find two ψ-shimmy (left-shaded) and one δ-shimmy (right-shaded) region. Grey-shaded areas indicate (V , c d ) regions in which no shimmy solutions exist. Consequently, a region c d ∈ (2754.1, 3019.4) bounded by the two blue lines, see insert plot in Fig. 7, is also observed. In this region the landing gear will not have sustained shimmy oscillations at any forward speed over the operating range 0-200 m/s. Specifically, this 205 indicates that if the damping of the shimmy damper is designed in the range (2754.1, 3019.4), the MLG will be free of shimmy for any value of V -this is the zero-shimmy region.

A beneficial spring-damper layout
From Fig. 7, we found that if the damping of the shimmy damper is designed in the zero-shimmy region, the MLG will be free from instabilities for any forward velocity. However, this zero-shimmy 210 region (2754.1, 3019.4) is quite narrow, so placing tight requirements for the design of the damper, such that it remains in this region for all operating conditions while also accounting for other damping contributions such as joint friction. Hence, it is desirable to expand this zero-shimmy region to provide a more flexible design requirement. To achieve this, we now investigate how varying the stiffness properties will affect the zero-shimmy region.

215
Note that the focus of this work is the shimmy-suppression device while the structural stiffness as listed in Table 1 will not be changed. The suppression device comprises a torsional stiffness k s connected in series to a spring-damper pair. This configuration is denoted L2 and is depicted in Fig. 8. Panel (a) illustrates that, when equipped with the shimmy damper, k ψ and k s contribute to give the effective stiffness k ψe . These stiffnesses are related by the expression We note that if k s = ∞, k ψe would equal k ψ , otherwise k ψe < k ψ . Now we can investigate the effects of varying k ψe via varying k s of L2 on the two-parameter bifurcation diagrams. Fig. 9 shows a series of continuations in the (V , c d )-plane when varying the values of k s . The black lines in Fig. 9(a) which represents k s = ∞ case is equivalent to that in Fig. 7. Note that the torus bifurcations curves are omitted for Fig. 9 and the following two-parameter bifurcation diagrams. The dominant types shrinks and the two ψ-shimmy regions become more significant. This indicates that a smaller k s , and hence a smaller torsional stiffness, would stabilise the MLG δ-shimmy while destabilising the ψ-shimmy. Fig. 9(b) shows that if k s is decreased to 5.9 × 10 5 Nm/rad, the δ-shimmy region 225 disappears and a wide zero-shimmy region is found. However, if k s is reduced further, the two ψ-shimmy regions intersect each other eliminating the zero-shimmy region. To identify the sensitivity of the system to the variations of k s in L2 qualitatively, four boundary points of the three shimmy regions in the (V , c d )-plane are defined, as shown in Fig. 10(a). BP1 and BP3 are the boundary points of the upper and lower ψ-shimmy regions, respectively. BP2 u 230 and BP2 l represent the upper and lower boundaries of the δ-shimmy region. As k s is varied the c d values of these four boundary points form four curves in Fig. 10(b), and the grey area represents the zero-shimmy region. It can be observed that the red (BP3) and blue (BP1) lines intersect at k s = 4.73 × 10 5 Nms/rad. When k s is increased from this value, the zero-shimmy region is growing to the maximum until the green lines, i.e. the δ-shimmy region, arise at k s = 7.86 × 10 5 Nms/rad. The grey area then splits into two parts, one is bounded by BP3 and BP2 u , the other by BP1 and BP2 l . Fig. 10(c) shows the width of the grey area with respect to the value of k s and we define this width as the width in damping of the zero-shimmy region. It can be seen that the zero-shimmy region is largest in terms of damping values (c d ∈ (3813.7, 14032.7)) at k s = 7.86 × 10 5 Nms/rad; at this point the width in damping of this region is 1.02 × 10 4 Nms/rad (14032.7 − 3813.7), which 240 is approximately 38 times what is obtainable with the shimmy damper alone (k s = ∞). After the maximum value, the width in damping of the zero-shimmy region reduces sharply, as shown by the two green lines of Fig. 10(c). Therefore, comparing with the original shimmy damper case, this beneficial spring-damper layout L2 is seen to provide a greater allowable damping range over which the MLG will be free from shimmy for all operating velocities.

Effects of inerter-based shimmy-suppression devices
The performance improvements obtained by the inerter-based vibration-suppression devices have been identified for many industrial applications, such as [12-18, 20, 40-43]. In this section, two variations on layout L2 that include inerters are considered. The MLG's sensitivity to the parameters of these two inerter-based layouts are investigated and their performance advantages 250 on expanding the zero-shimmy regions are then discussed. Fig. 11 shows two inerter-based shimmy-suppression layouts proposed in this work, labelled LI1 and LI2. The layouts are inspired by work on the vibration suppression of buildings. LI1 is referred as a TVMD-type layout, which includes the inerter in parallel with c d and k d of L2. This layout 255 was first proposed by Ikago et al. in [44] and was called the tuned-viscous-mass-damper (TVMD). The original layout excludes k d ; we propose k d here to function as a centering spring. LI2 is a TID-type layout; the original TID (tuned-inerter-damper) layout was introduced by Lazar et al. in [15]. The spring k s was not introduced in the original layout but is considered here, again acting as a centering stiffness. Note that for both layouts, we focus on the effects of k s and b while fixing 260 k d to the value used for the default shimmy damper, 1.09 × 10 5 Nm/rad. Note also that the two layouts are equivalent to layout L2 if b is set to zero.  Fig. 12 shows two-parameter bifurcation diagrams for the layout LI1 for several representative values of k s when fixing b to different values in the range of b ∈ [20, 150] Nms 2 /rad. Note that 265 b = 0 result is equivalent to the L2 bifurcation plots of Fig. 9. Several qualitative changes in the appearance of the δ-shimmy regions are observed here. When the inertance b is relatively small such as b = 20 Nms 2 /rad, changing k s produces a similar effect to that observed for L2 (b = 0). Hence, the width of the zero-shimmy region reaches its maximum just before the δ-shimmy occurs. When b is increased to 30 Nms 2 /rad, the first qualitative change is observed -the δ-shimmy curve appears 270 as a narrow-strip below the boundary point BP1 (recall the points BP1-BP4 from Fig. 10(a)) and moves upwards for larger k s . Consequently, for b = 30 Nms 2 /rad, the zero-shimmy region is bounded by two points BP1 and BP3 and is maximised just before the upper boundary of the δ-shimmy region BP2 u exceeds BP1. Increasing b to 70 Nms 2 /rad produces a similar variation as that seen for b = 30 Nms 2 /rad case but now the shape of the δ-shimmy region when it appears is 275 a small circle. When b equals 100 Nms 2 /rad, the upper boundary of the δ-shimmy region appears higher than BP1. Similar variation trends of the δ-shimmy region are observed for b = 130 and 150 Nms 2 /rad. Note that for the b ∈ [70, 150] Nms 2 /rad cases, an additional (third) ψ-shimmy dominated region is observed, as illustrated in Figs. 12(c)-(f). However, as this region is much lower than the point BP1 vertically in the (V , c d )-plane, it will not affect the zero-shimmy region.

Effects of layout LI1
occurs. Here the widest zero-shimmy region is still determined by the two points BP1 and BP3. In contrast, for b = 150 Nms 2 /rad, the upper boundary of the third ψ-shimmy region grows higher than BP1 before the appearance of the δ-shimmy region. Therefore in this situation, the zero-285 shimmy region is determined by BP3 and the upper boundary of the third ψ-shimmy region, and reaches the maximum just before the δ-shimmy emerges.  The boundary points in the (V , c d )-plane are plotted in Figs. 13(a1) and (a2) as k s and b vary. The zero-shimmy regions are shaded in grey and the vertical lines indicate the maximum width in damping of the zero-shimmy region occurs. In contrast to the results for the layout L2, the lower 290 bounds of the grey areas are not just dependent on BP1 and the maximum width in damping of the zero-shimmy region is not always determined by the appearance of the δ-shimmy. The variation of this width with respect to b (vertical lines in Figs. 13(a1) and (a2)) is plotted in Fig. 13(b). It can be observed that when b = 140 Nms 2 /rad and k s = 6.38 × 10 5 Nm/rad, the width of the zero-shimmy region is maximised to 1.12 × 10 4 Nms/rad (13561.6 − 2319.1), with a maximal 9.8% 295 improvement over that with L2 (1.02 × 10 4 Nms/rad). One notable benefit of the LI1 layout over L2 is that the lower bound of this region, 2319.1 Nms/rad, is significantly smaller than that of L2, 3813.7 Nms/rad. This is potentially beneficial from a manufacturing perspective.

Effects of layout LI2
The same approach is now taken to investigate the influence of layout LI2. Fixing b of LI2 to 300 the specific values, 50 and 100 Nms 2 /rad, and then varying k s , the variations of the bifurcation diagrams in the (V , c d )-plane are obtained and shown in Fig. 14. The zero-shimmy parameter regions are plotted in the (k s , c d )-plane in Fig. 15(a). As for the b = 0 case (layout L2), when b = 50 Nms 2 /rad the maximum width in damping of the zero-shimmy region arises just before the occurrence of δ-shimmy. If b is increased further, such as to 100 Nms 2 /rad, a relatively small k s 305 will result in a δ-shimmy as shown in Fig. 14(b1). This changes rapidly with the variation of k s and disappears if k s is allowed to increase sufficiently. The zero-shimmy region is affected by this δ-shimmy. Hence, as illustrated by the zero-shimmy areas bounded by the magenta and black lines of Fig. 15(a), the left boundaries are not formed by smooth curves. A further increased k s leads to the appearance of another δ-shimmy region. It can be seen that the width in damping of the zero-310 shimmy region reaches its maximum value just before this δ-shimmy arises. As with in Fig. 13(b), Fig. 15(b) gives the width in damping of the widest zero-shimmy region of Fig. 15(a) for each b. It can be calculated that the maximum width in damping of the zero-shimmy region increases with b, however, a much larger k s is required to obtain this maximum value. When comparing with the layout LI1, this growth trend ends with b = 140 Nms 2 /rad and gives the maximum width in

Conclusions
This paper investigates the effects of several shimmy-suppression devices on a MLG system and 320 proposes a method of selecting the device parameter values to prevent shimmy for any forward speed within the operating region. The forward velocity V and the device damping c d are used as continuation parameters to obtain two-dimensional bifurcation diagrams. It is found that, to ensure the MLG is free of shimmy instability for any forward velocity over the considered range, the damping of the conventional shimmy damper (L1) needs to be designed within a narrow range of 325 values and hence a tight tolerance. This device damping range is defined as the width in damping of the zero-shimmy region and needs to be enlarged to provide a larger parametric operating range for the damping device. By adding an additional spring k s in series with the shimmy-damper layout, a spring-damper layout (L2) is proposed. The effective torsional stiffness of the system is altered, resulting in an enlarged zero-shimmy region. The maximum width in damping of this zero-shimmy 330 region is about 38 times that obtained by the shimmy damper. Consequently, two inerter-combined layouts LI1 and LI2 are proposed. It is found that for LI1 a maximal 9.8% improvement over L2 can be obtained. One notable benefit of the LI1 layout over L2 is that the lower bound of this maximum zero-shimmy region is significantly smaller than that of L2. For LI2, the maximum width in damping of the zero-shimmy region increases significantly with growing b and k s . where i = L/R and N δ 5i1 = Λ i sin ψ(sin φ cos δ sin θ + sin δ cos θ) + cos φ cos δ sin ψ, N δ 5i2 = Λ i (− sin φ sin δ sin θ + cos δ cos θ) − cos φ sin δ, N δ 5i3 = −Λ i cos ψ(sin φ cos δ sin θ + sin δ cos θ) − cos δ cos ψ cos φ.
For Eq. (34), N z 11 = m D , N z 12 = m D l D cos φsinδ, N z 21 = m D l D cos φ cos δ.