Research on the Stability and Bifurcation Characteristics of a Landing Gear Shimming Dynamics System

Abstract

A dynamic model is established to investigate the shimmy instability of a landing gear system, considering the influence of nonlinear damping. The stability criterion is utilized to determine the critical speed at which the landing gear system becomes unstable. The central manifold theorem and canonical method are employed to simplify the dynamic model of the landing gear. The first Lyapunov coefficient of the system is theoretically derived and verified using numerical simulation. Further investigation on the Hopf bifurcation characteristics and stability of the shimmy in the landing gear system is conducted. The results indicate that above a certain threshold speed, with a tire stability distance greater than half the tire length in contact with the ground plus the slack length, the aircraft remains stable during taxiing. At critical speeds, a shimmy system with higher-order nonlinear damping will undergo supercritical Hopf bifurcation. Quantitative analysis suggests an increase in the linear damping coefficient within a range that ensures a stability margin to mitigate undesired oscillation, while the nonlinear damping coefficient should be designed within a reasonable range to decrease the amplitude of the limit cycle.

1. Introduction

The phenomenon of shimmy is a self-excited lateral yaw oscillation that can occur during the takeoff, landing, or taxiing of an aircraft’s landing gear. It is caused and sustained by the transfer of kinetic energy from the aircraft to the landing gear. In order to effectively address significant issues in the design of anti-swing mechanisms for aircraft landing gears, it is crucial to have a comprehensive understanding of the mechanism behind shimmy and investigate how various factors impact its stability [1]. The issue of shimmy has always posed a formidable challenge in the design and operation of aircraft landing gears. Mastering the accurate and appropriate analysis of shimmy dynamics is pivotal in addressing the oscillation problems associated with shimmy [2].

Liu [3] found that shimmy is primarily caused by the torsional mode of the wheel’s swinging part around its directional axis and the lateral bending mode of the landing gear around its longitudinal axis. The longitudinal mode, on the other hand, has minimal impact. At the same time, the tire and ground contact areas experience alternating deformation, and the oscillation of the wheel can cause the landing gear strut and the fuselage to shake, potentially resulting in the entire fuselage trembling, on the basis of the mechanics of the landing gear and linear models of the tire elasticity. Somieski [4] applied well-known linear mathematical methods to analyzing the shimmy of a basic nose landing gear model and investigated the impact of various structural parameters on the shimmy dynamics. Tartaruga [5] performed sensitivity and uncertainty analysis using singular value decomposition to investigate the impact of various structural parameters on shimmy onset. The results indicate that the proposed method reduces the computational time compared to full-scale full element simulation. The study conducted by Thota et al. [6] examined a seven-dimensional landing gear model encompassing torsional, lateral, and longitudinal degrees of freedom. Their analysis focused on the overall shimmy behavior and frequency, revealing that variations in the forward velocity and vertical force influence the interaction between the torsional mode and lateral bending mode, resulting in diverse types of shimmy phenomena. Conversely, minimal involvement was observed from the longitudinal mode regarding any potential occurrence of shimmy.

The landing gear is a complex structural system characterized by the high nonlinearity arising from the nonlinear deformation of tires, elastic tire deformation, and side slip [7]. The nonlinearities in the landing gear include shock absorbers, dampers, joint free clearance, and coulomb friction components [8,9], which contribute significantly to nonlinear vibrations. From a dynamics perspective, bifurcation phenomena primarily encompass saddle node bifurcation, transcritical bifurcation, pitchfork bifurcation, period-doubling bifurcation, and Hopf bifurcation [10]. Among these types of bifurcation, Hopf bifurcation is closely associated with the limit cycles and self-excited vibration within the system; thus, it represents a dynamic form of bifurcation. Shimmy motion corresponds to periodic solutions in nonlinear dynamics that typically emerge due to Hopf bifurcations. Consequently, it holds substantial research value.

Thota et al. [11] distinguished the stable region from the shimmy oscillation region by means of the Hopf bifurcation curve. They classified different types of shimmy as torsional, lateral, and quasi-periodic shimmy. Ghadami [12] proposed an effective method for predicting Hopf bifurcation and analyzing the complete post-bifurcation behavior of nonlinear oscillatory systems. By measuring the system’s response to a disturbance in its pre-bifurcation state, a three-dimensional bifurcation diagram was generated for prediction purposes. The advantages of this prediction method were demonstrated by utilizing the predicted bifurcation diagram of a fluid structure system, which accurately predicted both supercritical and subcritical bifurcations. Thota et al. [13] investigated a nose gear model with a two-wheel configuration. They selected forward velocity and vertical load as the bifurcation parameters and analyzed how structural parameters influenced two types of shimmy. Additionally, they provided transformations between different two-parameter bifurcation diagrams. The research revealed that besides stable torsional shimmy and stable lateral vibration, it could also trigger quasi-periodic shimmy. Li [14] studied the Hopf bifurcation of an aircraft nose landing gear model based on two pairs of continuous parameters. The results showed that higher-order terms caused deviations between the original system’s bifurcation curve and truncated amplitude system’s bifurcation line. This deviation led to a contraction in the bistable region compared to in previous findings.

In recent years, extensive research has been conducted on Hopf bifurcation and the central manifold of nonlinear dynamical systems. Zhang [15] proposed a model for a piecewise smooth suspended single wheelset, investigating the influence of the vehicle parameters on the characteristics of Hopf bifurcation. By employing the central manifold theorem, they successfully reduced the dimensionality of the wheelset equation and derived an expression for the first-order fine focus. Dong [16] utilized the normal form method to demonstrate that high-speed EMUs exhibit subcritical and supercritical Hopf bifurcations under simplified wheel–rail contact relationships. Additionally, Dong identified the factors contributing to the distinct types of bifurcations observed in two types of bogies. Sandor Beregi [17] employed a brush tire model to analyze the nonlinear dynamics of traction wheels by reducing its center manifold into a normal form comprising linear and piecewise smooth second-order terms. This normal form was used to establish stability at non-hyperbolic equilibrium points within the system and provided an estimate for the limit cycles appearing at linear stability boundaries. It has been proven that subcritical Hopf bifurcations in non-smooth delay systems generate parameter ranges with bistability that is undetectable using standard tire models. Wang [18], using the vehicle speed as a bifurcation parameter, applied central manifold theory to obtain a two-dimensional central manifold at critical vehicle speeds in order to comprehensively analyze the characteristics of Hopf bifurcation in shimmy systems while obtaining an approximate periodic solution.

The current research on the stability of aircraft nose landing gear has advanced from linear stability theory to nonlinear bifurcation theory [19,20]; however, there remains a lack of comprehensive explanations and solutions for certain fundamental shimmy problems encountered in practical applications. Furthermore, the analysis of shimmy problems has not yet considered the central manifold method. This article employs a combination of theoretical derivation and numerical simulation to analyze a landing gear system that incorporates higher-order nonlinear factors. The center manifold simplification method is utilized to streamline the nose landing gear system, while the position of the equilibrium point is examined using classical linear theory. The Hopf bifurcation of the equilibrium point is investigated using classical bifurcation theory, with its specific type determined based on the first Lyapunov coefficient. Both numerical simulation and theoretical derivation mutually validate each other’s findings. The research results presented in this article can serve as valuable theoretical references for optimizing the structural design of aircraft landing gear shimmy dynamics systems.

2. Dynamic Modeling of Nose Landing Gear Shimmy

2.1. Coordinate System and Degrees of Freedom

The stability of shimmy is influenced by various factors, including linear parameters such as stability distance, tire characteristics, and aircraft taxiing speed and non-linear factors such as damping, friction, and clearance. In order to investigate the impact of key parameters on the shimmy behavior, this study employs a simplified model of the nose landing gear in a two-wheeled aircraft equipped with an oil damper under a specific vertical load, as shown in Figure 1:

The landing gear follows the fuselage at a constant speed V, with the reverse direction of the taxiing defined as the X-axis of the landing gear coordinate system. The Y-axis represents horizontal movement to the right and the Z-axis represents vertical movement upward. The landing gear strut torsion is in a neutral position, so any lateral motion of the fuselage is not taken into consideration. During smooth taxiing, the aircraft does not experience any vertical displacement, and the stroke of the landing gear buffer remains constant. Therefore, the vertical displacement of the landing gear is not taken into account.

(1)

The landing gear strut can rotate around the directional axis with an angle of $\psi$.

(2)

In the X-axis direction, it is considered that the stiffness of the landing gear strut is sufficient, and the deformation has a small impact on the shimmy.

(3)

The lateral displacement $\lambda$ of the left and right tires is used to describe the deformation of the tires.

(4)

The degree of freedom $\phi$ of the axle torsion deformation indicates that the two wheels do not rotate together.

Due to the fact that the strut does not have a forward inclination angle, the tire angle $\theta$ is the same as the landing gear strut angle. Under a small angle rotation, assuming that there is no change in the vertical load on the two wheels, it can be considered that the lateral and torsional deformations of the two tires are the same. For the convenience of subsequent calculations, the stiffness of the two tires is superimposed and set as the stiffness damping parameter of a single tire.

2.2. Tire Dynamics Model

In Smiley’s tire mechanics theory, the forces and moments exerted by the ground on the tire can be expressed as:

$$F_{\lambda }=K_{\lambda }\lambda +C_{\lambda }\dot{\lambda }$$

(1)

$$M_{\phi }=K_{\phi }\phi +C_{\phi }\dot{\phi }$$

(2)

where $K_{\lambda }$ and $C_{\lambda }$ are the lateral stiffness and damping coefficients of the tire, respectively; $K_{\phi }$ and $C_{\phi }$ are the torsional stiffness and damping coefficient of the tire, respectively; and $\lambda$ and $\phi$ are the lateral elastic deformation and torsional elastic deformation of the tire, respectively.

As shown in Figure 2, $\lambda$ and $\phi$ can be expressed in the following form:

$$\lambda =y-e\theta$$

(3)

$$\phi =\dot{y}/V-\theta$$

(4)

$y$ is the lateral displacement of the tire’s touchdown center point while $e$ is the landing gear stability distance. Smiley’s theory uses Taylor series, ignoring the higher-order terms in the series, and can form different tire approximate constraint equations. This article adopts Smiley’s second-order approximation formula:

$$L_2\ddot{y}+V{L}_1\dot{y}+V^2\lambda +V^2{L}_1\theta =0$$

(5)

In the equation, $L_2=h(h+2\sigma )/2$,$L_1=h+\sigma$. $h$ is the half length of the tire touching the ground; $\sigma$ is the tire slack length.

We calculate the tire formula using equation coordination:

$$\dot{\lambda }+e\dot{\theta }+\mathrm{V}(\theta +\phi )=0$$

(6)

$$\dot{\phi }+\dot{\theta }-V\lambda /L_2+L_{1}V\phi /L_2=0$$

(7)

2.3. Nonlinear Analysis of Pendulum Reducers

The torsional damping moment provided by anti-swing dampers typically exhibits a nonlinear relationship with the swing speed of the column’s torsional degrees of freedom. During the front wheel shimmy process, the oil in the stabilizer is confined to enclosed circuits such as damping pipelines and throttling devices. The damping pipeline of a pendulum reducer generally adopts a circular cross-section. Ideally, the damping torque of the pendulum reducer demonstrates a quadratic power function characteristic. In reality, however, factors such as temperature and fluid boundary layer effects may result in high-power damping as well. The torsional anti-swing damping torque provided by the pendulum reducer can be roughly described as follows:

$$M_C=C\dot{\psi }+C_2{\dot{\psi }}^{2}sgin(\dot{\psi })+C_3{\dot{\psi }}^3+\mathrm{h}.\mathrm{o}.\mathrm{t}.$$

(8)

$C$, $C_2$, and $C_3$ are the linear term torsional damping coefficient, square term torsional damping coefficient, and cubic term torsional damping coefficient, respectively.

When considering the square term in nonlinear damping, its inclusion of a sign function poses challenges for analysis. Within the range of small angle vibrations, cubic damping also exhibits an increasing trend with angular velocity. In terms of the nonlinear influence properties, studying cubic damping can serve as a substitute for investigating square term damping.

2.4. Dynamic Equations of the System

The dynamic equilibrium equation for the torsion direction of the pillar is:

$$I\ddot{\psi }+C\dot{\psi }+C_3{\dot{\psi }}^3=e\cdot F_{\lambda }+M_{\phi }$$

(9)

Make $\mathit{x}={[x_1,x_2,x_3,x_4]}^T={[y,\dot{y},\lambda ,\phi ]}^T$. The equation can be represented as a system of first-order differential equations, as follows:

$$\dot{\mathit{x}}=\mathit{A}\mathit{x}=\mathit{F}(\mathit{x})$$

(10)

where:

$$\mathit{A}=\left[\begin{array}{cccc}0& 1& 0& 0\\ 0& -\frac{C}{I}& \frac{e{K}_{\lambda }}{I}& \frac{K_{\phi }}{I}\\ -V& -e& 0& -V\\ 0& -1& \frac{V}{L_2}& -\frac{V{L}_1}{L_2}\end{array}\right]$$

(11)

$$\mathit{F}(\mathit{x})={[0,f_2(x_2),0,0]}^T$$

(12)

The parameters of the nose landing gear of the aircraft used in this study are shown in

Table 1. Structural model parameter table.

Symbolic	Description	Value	Unit
$e$	Stability distance	0.08	[m]
$I$	Moment of inertia around the axis of the pillar	0.25	[kg∙m2]
$\sigma$	Tire slack length	0.1376	[m]
$h$	Half length of tire touching the ground	0.0517	[m]
$K_{\lambda }$	Tire lateral stiffness	2 × 2 × 105	[N/m]
$K_{\phi }$	Tire torsional stiffness	2 × 1.2 × 103	[Nm/rad]
$C$	Linear damping coefficient	10	[Nms]
$C_3$	Nonlinear damping coefficient	1	[Nms3]

:

3. Hopf Bifurcation Theory of the Oscillation Model

3.1. Linear Stability Theory

Before conducting stability analysis on the nose landing gear system, it is imperative to initially acquire the eigenvalues of the linearized system at the equilibrium point. The characteristic equation of the linear component within the dynamics equation pertaining to the landing gear can be expressed as follows:

$$|\mathit{A}-\lambda {\mathit{I}}_4|=0$$

(13)

In the equation, ${\mathit{I}}_4$ is the fourth-order identity matrix, and λ represents the characteristic root of the landing gear system. Expanding the equation yields:

$${\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4=0$$

(14)

In the equation, $a_{1~4}$ is the coefficient of each order of the characteristic equation, expressed as follows:

$$\{\begin{array}{l}{a}_1=(C\cdot L_2+I\cdot V\cdot L_1)/(I\cdot L_2)\hfill \\ a_2=(V\cdot C\cdot L_1+K_{\phi }\cdot L_2+I\cdot V^2+K_{\lambda }\cdot e^2\cdot L_2)/(I\cdot L_2)\hfill \\ a_3=(V\cdot e\cdot K_{\phi }+V^2\cdot C+V\cdot K_{\lambda }\cdot e^2\cdot L_1)/(I\cdot L_2)\hfill \\ a_4=(V^2\cdot K_{\lambda }\cdot e\cdot L_1+V^2\cdot K_{\phi })/(I\cdot L_2)\hfill \end{array}$$

(15)

When calculating the linear critical speed of the landing gear system, there are typically two methods available: the root trajectory method and the Hurwitz stability criterion method. In comparison to the root locus method, the Hurwitz stability criterion method eliminates the need to solve the characteristic roots of the system equation. Instead, it solely relies on determining the system motion stability based on the coefficients within the characteristic equation [21]. The Liénard–Chipart stability criterion adopted in this paper further simplifies the calculation process based on Hurwitz [22].

The Hurwitz matrix is as follows:

$${\Delta }_4=\left|\begin{array}{llll}{a}_1& 1& 0& 0\\ a_3& a_2& a_1& 1\\ 0& a_4& a_3& a_2\\ 0& 0& 0& a_4\end{array}\right|$$

(16)

According to the Liénard–Chipart stability criterion, the necessary and sufficient conditions for the characteristic equation to have a pair of pure imaginary roots and the other two roots to have negative real parts are:

$$\{\begin{array}{l}{a}_{1~4}>0,\hfill \\ {\Delta }_1>0,{\Delta }_3\ge 0\hfill \end{array}$$

(17)

When ${\Delta }_3=0$, the system is in a critical stable state, and when ${\Delta }_3>0$, the system is in a convergence state.

3.2. Hopf Bifurcation Theory

According to the Hopf bifurcation criterion, the necessary conditions for Hopf bifurcation are as follows:

$${\Delta }_3=\left|\begin{array}{ccc}{a}_1& 1& 0\\ a_3& a_2& a_1\\ 0& a_4& a_3\end{array}\right|=a_1{a}_2{a}_3-a_3^2-a_1^2{a}_4=0$$

(18)

According to the practical significance of the structural parameters of the landing gear, it is easy to verify $a_{1~4}>0$; from this, it can be seen that when ${\Delta }_3=0$ is satisfied, it is the critical state of landing gear system equilibrium point instability, and the univariate quantic equation of speed can be obtained:

$$b_0{v}^5+b_1{v}^4+b_2{v}^3+b_3{v}^2+b_4{v}^1+b_5=0$$

(19)

where $b_{0~5}$ is the coefficient of each order, and the expression is as follows:

$$\{\begin{array}{l}{b}_0=I^2\cdot C\cdot L_1/(I^3\cdot L_2^3)\hfill \\ b_1=(I\cdot C^2\cdot L_1^2+I^2\cdot L1\cdot e\cdot K\phi +I^2\cdot K\lambda \cdot e^2\cdot L_1^2-I^2\cdot e\cdot K\lambda \cdot L_1^3-I^2\cdot K\phi \cdot L_1^2)/(I^3\cdot L_2^3)\hfill \\ \begin{array}{l}{b}_2=(C^3\cdot L_1\cdot L_2+I\cdot C\cdot K\phi \cdot L_1\cdot L_2-I\cdot e\cdot C\cdot K\phi \cdot L_2+I\cdot C\cdot e\cdot K\phi \cdot L_1^2+I\cdot C\cdot e^2\cdot K\lambda \cdot L_1^3\\ -I\cdot 2\cdot e\cdot C\cdot K\lambda \cdot L_1^2\cdot L_2-I\cdot 2\cdot C\cdot K\phi \cdot L_1\cdot L_2)/(I^3\cdot L_2^3)\end{array}\hfill \\ \begin{array}{l}{b}_3=(C^2\cdot K\phi \cdot L_2^2+C^2\cdot e^2\cdot K\lambda \cdot L_2^2+(e\cdot K\phi +e^2\cdot K\lambda \cdot L_1)\cdot (I\cdot K\phi \cdot L_1\cdot L_2+C^2\cdot L_1\cdot L_2-I\cdot e\cdot K\phi \cdot L_2)\\ -C^2\cdot L_2^2\cdot (K\lambda \cdot e\cdot L_1+K\phi ))/(I^3\cdot L_2^3)\\ b_4=(e\cdot K\phi +e^2\cdot K\lambda \cdot L_1)\cdot (C\cdot K\phi \cdot L_2^2+C\cdot K\lambda \cdot e^2\cdot L_2^2)/(I^3\cdot L_2^3)\\ b_5=0\end{array}\hfill \end{array}$$

(20)

The smallest positive real root of all solutions of the equation is the linear critical speed of the landing gear system. At the Hopf bifurcation point of the nose landing gear system, the characteristic equation has a pair of pure imaginary roots ${\lambda }_{1,2}=\pm i{\omega }_0$, where ${\omega }_0>0$, because ${\lambda }_{1,2}=\pm i{\omega }_0$ is the root of the characteristic equation. We bring ${\lambda }_{1,2}=\pm i{\omega }_0$ into the characteristic equation to obtain:

$$\left({\omega }_0^4-a_2{\omega }_0^2+a_4\right)-i\left(a_1{\omega }_0^3-a_3{\omega }_0\right)=0$$

(21)

Since the real part is zero, we have:

$${\omega }_0^2=\frac{a_3}{a_1}=\frac{V^2\cdot C+V\cdot K_{\phi }\cdot e+V\cdot K_{\lambda }\cdot e^2\cdot L_1}{V\cdot I\cdot L_1+C\cdot L_2}$$

(22)

So, the characteristic equation is decomposed into:

$$\left({\lambda }^2+{\omega }_0^2\right)\left({\lambda }^2+p\lambda +q\right)=0$$

(23)

According to the equal coefficient of the same term:

$$p=a_1,q=\frac{a_1{a}_4}{a_3}$$

(24)

When the discriminant satisfies $\Delta =p^2-4q<0$, a pair of eigenvalues with negative real parts of the characteristic equation can be obtained, ${\lambda }_{3,4}=\alpha \pm i{\omega }_1$, where ${\omega }_1>0$. The specific expression is as follows:

$$\alpha =-p/2,{\omega }_1=\frac{\sqrt{4q-p^2}}{2}$$

(25)

When the discriminant satisfies $\Delta =p^2-4q\ge 0$, a pair of real root eigenvalues of the characteristic equation can be obtained: ${\lambda }_{3,4}=\frac{-p\pm \sqrt{\Delta }}{2}$.

3.3. Dimension Reduction in the Central Manifold Theorem

If the bifurcation characteristics of the high-dimensional system are directly analyzed, the process is cumbersome, and the dimensions of the high-dimensional system can be reduced to the plane system using the central manifold theorem [23]. Then, the analysis method of the plane system is used to complete the research on the high-dimensional system. When the current landing gear equation is at a critical speed, let ${\mathit{\beta }}_1\in {ℂ}^4$ be the eigenvector of matrix $\mathit{A}$ corresponding to eigenvalue $i{\omega }_0$. Take two eigenvalues containing imaginary parts as an example: let ${\mathit{\beta }}_2\in {ℂ}^4$ be the eigenvector of matrix $\mathit{A}$ corresponding to eigenvalue $\alpha \pm i{\omega }_1$. Then, ${\mathit{\beta }}_1,{\mathit{\beta }}_2$ satisfies the following equation:

$$\{\begin{array}{l}\mathit{A}{\mathit{\beta }}_1=i{\omega }_0{\mathit{\beta }}_1\hfill \\ \mathit{A}{\mathit{\beta }}_2=(\alpha +i{\omega }_1){\mathit{\beta }}_2\hfill \end{array}$$

(26)

where the expression of $B_{ij}$ is as follows:

$$\{\begin{array}{l}{B}_{31}=(-K_{\phi }^2\cdot e\cdot {\omega }_0^2+K_{\lambda }\cdot K_{\phi }\cdot V^2\cdot e-C\cdot K_{\phi }\cdot V\cdot {\omega }_0^2-K_{\lambda }\cdot I\cdot V^2\cdot e\cdot {\omega }_0^2)/(K_{\lambda }^2\cdot V^2\cdot e^2+K_{\phi }^2\cdot {\omega }_0^2)\\ B_{41}=(-K_{\lambda }^2\cdot V^2\cdot e^2+C\cdot K_{\lambda }\cdot V\cdot e\cdot {\omega }_0^2+K_{\phi }\cdot K_{\lambda }\cdot e^2\cdot {\omega }_0^2-K_{\phi }\cdot I\cdot {\omega }_0^4)/(K_{\lambda }^2\cdot V^2\cdot e^2+K_{\lambda }^2\cdot {\omega }_0^2)\\ B_{32}=(K_{\phi }^2\cdot V\cdot {\omega }_0+K_{\lambda }\cdot K_{\phi }\cdot V\cdot e^2\cdot {\omega }_0-I\cdot K_{\phi }\cdot V\cdot {\omega }_0^3+C\cdot K_{\lambda }\cdot V^2\cdot e\cdot {\omega }_0)/(K_{\lambda }^2\cdot V^2\cdot e^2+K_{\phi }^2\cdot {\omega }_0^2)\\ B_{42}=(-K_{\lambda }^2\cdot V\cdot e^3\cdot {\omega }_0+I\cdot K_{\lambda }\cdot V\cdot e\cdot {\omega }_0^3-K_{\phi }\cdot K_{\lambda }\cdot V\cdot e\cdot {\omega }_0+C\cdot K_{\phi }\cdot {\omega }_0^3)/(K_{\lambda }^2\cdot V^2\cdot e^2+K_{\phi }^2\cdot {\omega }_0^2)\\ B_{33}=(-K_{\phi }^2\cdot V\cdot a-K_{\phi }^2\cdot a^2\cdot e-K_{\phi }^2\cdot e\cdot {\omega }_1^2+K_{\lambda }\cdot K_{\phi }\cdot V^2\cdot e-I\cdot K_{\phi }\cdot V\cdot a^3-C\cdot K_{\phi }\cdot V\cdot a^2+K_{\lambda }\cdot K_{\phi }\cdot V\cdot a\cdot e^2\\ -I\cdot K_{\phi }\cdot V\cdot a\cdot {\omega }_1^2-C\cdot K_{\phi }\cdot V\cdot {\omega }_1^2+K_{\lambda }\cdot I\cdot V^2\cdot a^2\cdot e+C\cdot K_{\lambda }\cdot V^2\cdot a\cdot e-K_{\lambda }\cdot I\cdot V^2\cdot e\cdot {\omega }_1^2)/(K_{\lambda }^2\cdot V^2\cdot e^2-\\ 2\cdot K_{\lambda }\cdot K_{\phi }\cdot V\cdot a\cdot e+K_{\phi }^2\cdot a^2+K_{\phi }^2\cdot {\omega }_1^2)\\ B_{43}=(K_{\phi }^2\cdot V\cdot {\omega }_1-I\cdot K_{\phi }\cdot V\cdot a^2\cdot {\omega }_1+K_{\lambda }\cdot K_{\phi }\cdot V\cdot e^2\cdot {\omega }_1-I\cdot K_{\phi }\cdot V\cdot {\omega }_1^3+\\ K_{\lambda }\cdot I\cdot V^2\cdot a\cdot e\cdot {\omega }_1\cdot 2+C\cdot K_{\lambda }\cdot V^2\cdot e\cdot {\omega }_1)/(K_{\lambda }^2\cdot V^2\cdot e^2-2\cdot K_c\cdot K_{\phi }\cdot V\cdot a\cdot e+K_{\phi }^2\cdot a^2+K_{\phi }^2\cdot {\omega }_1^2)\\ B_{34}=(-K_{\lambda }^2\cdot V^2\cdot e^2-K_{\lambda }^2\cdot V\cdot a\cdot e^3-I\cdot K_{\lambda }\cdot V\cdot a^3\cdot e-C\cdot K_{\lambda }\cdot V\cdot a^2\cdot e+3\cdot I\cdot K_{\lambda }\cdot V\cdot a\cdot e\cdot {\omega }_1^2+\\ K_{\phi }\cdot K_{\lambda }\cdot V\cdot a\cdot e+C\cdot K_{\lambda }\cdot V\cdot e\cdot {\omega }_1^2+K_{\phi }\cdot K_c\cdot a^2\cdot e^2+K_{\phi }\cdot K_{\lambda }\cdot e^2\cdot {\omega }_1^2+K_{\phi }\cdot I\cdot a^4+C\cdot K_{\phi }\cdot a^3+\\ C\cdot K_{\phi }\cdot a\cdot {\omega }_1^2-K_{\phi }\cdot I\cdot {\omega }_1^4)/(K_c^2\cdot V^2\cdot e^2-2\cdot K_{\lambda }\cdot K_{\phi }\cdot V\cdot a\cdot e+K_{\phi }^2\cdot a^2+K_{\phi }^2\cdot {\omega }_1^2)\\ B_{44}=(-K_{\lambda }^2\cdot V\cdot e^3\cdot {\omega }_1-I\cdot K_{\lambda }\cdot V\cdot a^2\cdot e\cdot {\omega }_1\cdot 3-C\cdot K_{\lambda }\cdot V\cdot a\cdot e\cdot {\omega }_1\cdot 2+I\cdot K_{\lambda }\cdot V\cdot e\cdot {\omega }_1^3-K_{\phi }\cdot K_{\lambda }\cdot V\cdot e\cdot {\omega }_1+\\ K_{\phi }\cdot I\cdot a^3\cdot {\omega }_1\cdot 2+C\cdot K_{\phi }\cdot a^2\cdot {\omega }_1+K_{\phi }\cdot I\cdot a\cdot {\omega }_1^3\cdot 2+C\cdot K_{\phi }\cdot {\omega }_1^3)/(K_{\lambda }^2\cdot V^2\cdot e^2-2\cdot K_{\lambda }\cdot K_{\phi }\cdot V\cdot a\cdot e+K_{\phi }^2\cdot a^2+\\ K_{\phi }^2\cdot {\omega }_1^2)\end{array}$$

(27)

Before applying the central manifold theorem, first, the real eigenbasis matrix $\mathit{B}$ is constructed by using the eigenvector ${\mathit{\beta }}_1,{\mathit{\beta }}_2$, and then the landing gear equation is transformed into the Jordan canonical form using matrix $\mathit{B}$. Matrix $\mathit{B}$ is as follows:

$$B=[\mathrm{Re}{\mathit{\beta }}_1,-\mathrm{Im}{\mathit{\beta }}_1,\mathrm{Re}{\mathit{\beta }}_2,-\mathrm{Im}{\mathit{\beta }}_2]$$

(28)

For the convenience of symbol derivation, matrix $\mathit{B}$ is recorded as:

$$B=\left[\begin{array}{cccc}1& 0& 1& 0\\ 0& -{\omega }_0& \alpha & -{\omega }_1\\ B_{31}& B_{32}& B_{33}& B_{34}\\ B_{41}& B_{42}& B_{43}& B_{44}\end{array}\right]$$

(29)

Let $\mathit{y}={\left[y_1,y_2,y_3,y_4\right]}^{\mathrm{T}}\in {ℂ}^4$ where $\mathit{x}=\mathit{B}\mathit{y}$ be substituted into the landing gear equation and finally transformed into the Jordan standard:

$$\left[\begin{array}{c}{\dot{y}}_1\\ {\dot{y}}_2\\ {\dot{y}}_3\\ {\dot{y}}_4\end{array}\right]=\left[\begin{array}{cccc}0& -{\omega }_0& 0& 0\\ {\omega }_0& 0& 0& 0\\ 0& 0& \alpha & -{\omega }_1\\ 0& 0& {\omega }_1& {\alpha }_1\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ y_4\end{array}\right]+\left[\begin{array}{c}{b}_{12}f(\mathit{B}\mathit{y})\\ b_{22}f(\mathit{B}\mathit{y})\\ b_{32}f(\mathit{B}\mathit{y})\\ b_{42}f(\mathit{B}\mathit{y})\end{array}\right]$$

(30)

Including:

$$\begin{array}{cc}\left|B\right|& =a·B_3{}_4·B_4{}_2-a·B_3{}_2·B_4{}_4+{\omega }_0·B_3{}_4·B_4{}_1-{\omega }_0·B_3{}_4·B_{43}\hfill \\ & -{\omega }_0·B_3{}_1·B_{44}+{\omega }_0·B_3{}_3·B_{44}-{\omega }_1·B_3{}_2·B_{41}\hfill \\ & +{\omega }_4·B_3{}_4·B_{45}-{\omega }_1·B_3{}_3·B_{42}+{\omega }_1·B_3{}_2·B_{43}\hfill \end{array}$$

(31)

$$\{\begin{array}{l}{b}_1{}_2=(-B_3{}_4\cdot B_4{}_2+B_3{}_2\cdot B_4{}_4)/|B|\hfill \\ b_2{}_2=(-B_3{}_4\cdot B_4{}_1+B_3{}_4\cdot B_4{}_3+B_3{}_1\cdot B_4{}_4-B_3{}_3\cdot B_4{}_4)/|B|\hfill \\ b_{32}=(B_{34}\cdot B_{42}-B_{32}\cdot B_{44})/|B|\hfill \\ b_4{}_2=(B_3{}_2\cdot B_4{}_1-B_3{}_1\cdot B_4{}_2+B_3{}_3\cdot B_4{}_2-B_3{}_2\cdot B_4{}_3)/|B|\hfill \end{array}$$

(32)

For the critical point of Hopf bifurcation, there are two eigenvalues of the Jordan type matrix with zero real parts and one with negative real parts.

And because $f(\mathit{B}\mathit{y})\in {ℂ}^2$, then there is a central manifold $y_3=h_1(y_1,y_2)$, $y_4=h_2(y_1,y_2)$, $h(0,0)=\mathit{D}h(0,0)=0$, and we define the operator:

$$\begin{array}{c}N(h(y_1,y_2))=h^{\prime }(y_1,y_2)[A(y_1,y_2)+f(y_1,y_2+h(y_1,y_2))]\\ -Bh(y_1,y_2)-g(y_1,y_2,h(y_1,y_2))]\end{array}$$

(33)

It can be seen from $f(\mathit{B}\mathit{y})=-C_3{(-{\omega }_0{y}_2+\alpha y_3-{\omega }_1{y}_4)}^3$ that $h(y_1,y_2)$ is a nonlinear function of no less than order 3, and the Taylor expansion of the central manifold is as follows:

$$\begin{array}{l}{h}_1\left(y_1,y_2\right)=k_{30}{y}_1^3+k_{21}{y}_1^2{y}_2+k_{12}{y}_1{y}_2^2+k_{03}{y}_2^3+O(\Vert \mathit{y}{\Vert }^4)\hfill \\ h_2\left(y_1,y_2\right)=l_{30}{y}_1^3+l_{21}{y}_1^2{y}_2+l_{12}{y}_1{y}_2^2+l_{03}{y}_2^3+O\left(\Vert \mathit{y}{\Vert }^4\right)\hfill \end{array}$$

(34)

We bring the equation into the operator and merge the terms with the same degree. Since the coefficient is 0, eight equations and eight unknowns are obtained, which can be solved into the $h(y_1,y_2)$ coefficient. Considering that only quadratic and cubic terms of the state variables are needed to calculate the first Lyapunov coefficient, the higher-order term system of the central manifold is ignored, and the final dimension reduction is as follows:

$$\left[\begin{array}{c}{\dot{y}}_1\\ {\dot{y}}_2\end{array}\right]=\left[\begin{array}{cc}0& -{\omega }_0\\ {\omega }_0& 0\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right]+\left[\begin{array}{c}{\omega }_0^3{b}_{12}{C}_3{y}_2{}^3\\ {\omega }_0^3{b}_{22}{C}_3{y}_2{}^3\end{array}\right]$$

(35)

In order to judge the Hopf bifurcation type of the shimmy system and investigate the limit cycle stability of the system, the stability parameter of the central manifold is introduced:

$$\begin{array}{ll}\partial =& \frac{1}{16}\left(\frac{{\partial }^3{f}_1}{\partial y_1^3}+\frac{{\partial }^3{f}_1}{\partial y_1\partial y_2^2}+\frac{{\partial }^3{f}_2}{\partial y_1^2\partial y_2}+\frac{{\partial }^3{f}_2}{\partial y_2^3}\right)+\frac{1}{16{\omega }_0}[\frac{{\partial }^2{f}_1}{\partial y_1\partial y_2}(\frac{{\partial }^2{f}_1}{\partial y_1^2}+\frac{{\partial }^2{f}_1}{\partial y_2^2})\\ & -\frac{{\partial }^2{f}_2}{\partial y_1\partial y_2}(\frac{{\partial }^2{f}_2}{\partial y_1^2}+\frac{{\partial }^2{f}_2}{\partial y_2^2})-\frac{{\partial }^2{f}_1}{\partial y_1^2}\frac{{\partial }^2{f}_2}{\partial y_1^2}+\frac{{\partial }^2{f}_1}{\partial y_2^2}\frac{{\partial }^2{f}_2}{\partial y_2^2}]\end{array}$$

(36)

The central popular stability parameter is obtained as $\partial =\frac{3}{8}{\omega }_0^3{C}_3{b}_{22}$.

3.4. Analysis of the Time-Delay Effect under Linear Damping

Hopf bifurcation refers to the phenomenon of periodic motion that arises from the equilibrium state when the bifurcation parameters cross the bifurcation point. In order to distinguish the type of Hopf bifurcation in the landing gear system, this section utilizes canonical theory to derive the discriminant for identifying the type of Hopf bifurcation in the landing gear system [24].

When Hopf bifurcation occurs in the landing gear system, there is a pair of pure imaginary roots in the equation $\pm i{\omega }_0$. Let $\mathit{q},\mathit{p}\in {ℂ}^2$ be the right and left eigenvectors of matrix ${\mathit{J}}_c$: then, $\mathit{q},\mathit{p}$ satisfies the following equation:

$$\{\begin{array}{l}{\mathit{J}}_c\mathit{q}=i{\omega }_0\mathit{q}\hfill \\ {\mathit{J}}_c^{\mathrm{T}}\mathit{p}=-i{\omega }_0\mathit{p}\hfill \end{array}$$

(37)

Normalize $\mathit{p}$ to $\mathit{q}$; then, at this time:

$$\{\begin{array}{l}\langle \mathit{p},\mathit{q}\rangle =1\hfill \\ \langle \mathit{p},\overline{\mathit{q}}\rangle =0\hfill \end{array}$$

(38)

Among them,

$$<\mathit{p},\mathit{q}>=\overline{\mathit{p}}\cdot \mathit{q}={\overline{p}}_1{q}_1+{\overline{p}}_2{q}_2$$

(39)

By solving the above equations together, it can be obtained that:

$$\mathit{q}={\left[\frac{1}{2},-\frac{i}{2}\right]}^{\mathrm{T}},\mathit{p}={\left[1,-i\right]}^{\mathrm{T}}$$

(40)

Introducing complex variables $z$, any vector in the equation $y_c\in {ℂ}^2$ can be uniquely represented as $y_c=zq+\overline{z} \overline{\mathit{q}}$: we bring it into the available:

$$z=i\omega z+g(z,\overline{z},v_c)$$

(41)

Among them,

$$g(z,\overline{z},v_c)=\langle \mathit{p},{\mathit{\Gamma }}_c(z\mathit{q}+\overline{z} \overline{\mathit{q}},v_c)\rangle =-\frac{C_3(b_{22}+i{b}_{12}){(z+\overline{z})}^3}{2}$$

(42)

$g(z,\overline{z},v_c)$ is expressed in the Taylor series form of conjugate complex variable $(z,\overline{z})$ as follows:

$$\{\begin{array}{l}g(z,\overline{z},v_c)=g_2(z,\overline{z},v_c)+g_3(z,\overline{z},v_c)+O({|z|}^4)\hfill \\ g_r(z,\overline{z},v_c)={\displaystyle \sum _{k+l=i}\frac{1}{k!l!}}{g}_{kl}{z}^k{\overline{z}}^l,(k,l=0,1,\cdot \cdot \cdot \cdot )\hfill \end{array}$$

(43)

From the equality of the power coefficients of the same degree, it can be concluded that:

$$g_{11}=0,g_{20}=0,g_{21}=\frac{3}{4}{\omega }_0^3{C}_3(b_{22}+i{b}_{12})$$

(44)

The first-order fine focus $c_1$ of the system is as follows:

$$c_1=\frac{1}{2}{\tilde{g}}_{21}=\frac{i}{2{\omega }_0}(g_{20}{g}_{11}-2|g_{11}{|}^2-\frac{1}{3}|g_{02}{|}^2)+\frac{1}{2}{g}_{21}$$

(45)

The first Lyapunov coefficient can be obtained as:

$$l_1(0)=\frac{\mathrm{Re}{c}_1}{{\omega }_0}=\frac{1}{2{\omega }_0}\mathrm{Re}(\frac{3}{4}{\omega }_0^3{C}_3(b_{22}+i{b}_{12}))=\frac{3}{8}{\omega }_0^2{C}_3{b}_{22}$$

(46)

When $l_1(0)<0$, supercritical Hopf bifurcation occurs in the landing gear system. When $l_1(0)>0$, the undercarriage system undergoes subcritical Hopf bifurcation. When $l_1(0)=0$, the first Lyapunov coefficient degenerates, and the landing gear system degenerates the Hopf bifurcation. Since it is greater than zero, the central popular stability parameter and Lyapunov coefficient have the same sign, so it is accurate to judge the type of Hopf bifurcation.

4. Stability Analysis of Oscillation System

4.1. Numerical Verification of the Theoretical Results

To verify the accuracy of all the derived formulas, given the parameter values, in addition to the parameters in Table 1, we set the linear damping coefficient C = 10 Nms, substitute it into the formula, and calculate it to obtain $b_0$ = 1.2375 × 10⁷, $b_1$ = −1.0603 × 10⁹, $b_2$ = −2.6470 × 10⁶, $b_3$ = 4.8429 × 10¹⁰, and $b_4$ = 2.5335 × 10¹¹. By simplifying the equation, it can be concluded that:

$$v(b_0{v}^4+b_1{v}^3+b_2{v}^2+b_3{v}^1+b_4)=0$$

(47)

Taking the roots of the equation yields the following five roots: 85.1091 + 0.0000i, 8.9839 + 0.0000i, −4.2065 + 3.0134i, −4.2065 − 3.0134i, 0. The smallest positive real roots 8.9839 and 85.1091 in the solution are the linear critical speeds of the landing gear system. We select speed points near the speed for analysis as follows:

The stability analysis of time-delay systems, similar to their non-time-delay counterparts, can be primarily categorized into two methods: time-domain and frequency-domain approaches. The key advantage of the time-domain analysis method lies in its ability to effectively handle system nonlinearity and uncertainty. By setting the time delay as 0 ms, 3 ms, and 4 ms, the simulation results are depicted in the Figure 3 and Figure 4 below.

At the critical speed point, the swing angle of the landing gear exhibits constant amplitude vibration. When the speed value is slightly changed, the system’s state either becomes divergent or convergent, thus verifying the accuracy of the derived formula.

Shimmy primarily involves the stability and frequency of a system, with stability being just one factor in determining the criticality of the system. Consistency in frequency can ensure the accuracy of the theoretical results. According to the analysis in the previous section, the frequency expression of oscillation is as follows:

$$f=\frac{{\omega }_0}{2\pi }=\frac{1}{2\pi }\sqrt{\frac{V^2\cdot C+V\cdot K_{\phi }\cdot e+V\cdot K_{\lambda }\cdot e^2\cdot L_1}{V\cdot I\cdot L_1+C\cdot L_2}}$$

(48)

Therefore, the relationship between frequency and speed is shown in the following Figure 5:

According to the calculations, the oscillation frequencies corresponding to the two critical speeds are 18.483 Hz and 28.266 Hz, respectively.

The FFT algorithm is a fast algorithm for the discrete Fourier transform, which transforms the time-domain signals into frequency-domain signals. By analyzing the frequencies of time-domain, it can be observed that the frequencies of the two are 18.5 Hz and 28.3 Hz, respectively. It can be seen that the values of the two are similar. In summary, the numerical analysis results are consistent with the aforementioned analytical analysis results, confirming the accuracy of the model and the effectiveness of the analytical method.

4.2. The Influence of the Stability Distance on the Stability of Shimmy

The landing gear, being a complex mechanical dynamic system, exhibits significant coupling effects among its various components and subsystems. Consequently, there exist numerous factors that influence the aircraft shimmy, along with diverse engineering approaches to mitigating landing gear shimmy. However, parameters such as the tire radius, lateral stiffness, and equivalent mass of the entire nose landing gear pose challenges in terms of adjustment and optimization within engineering practice. The stability distance holds paramount importance in shimmy analysis as it plays a crucial role in suppressing shimmy by generating a restoring torque on the landing gear strut through the lateral force exerted on the tire due to its stable distance from the ground. In light of this significance, this article primarily focuses on investigating the geometric stability distance of the landing gear, which can be easily adjusted. The objective is to provide technical guidance for enhancing anti-shimmy design in engineering practice.

Due to the large overall mass of the tire, changes in the stability distance can lead to changes in the inertia around the axis. Assuming a wheel and axle mass of 10 kg, the total rotational inertia of the piston rod, torsion arm, and other rotating components around the pillar $I_0$ = 0.186 kgm², and we have:

$$I=I_0+Mw\cdot e^2$$

(49)

When the damping of the anti-swing damper and the structural damping are not considered, then $b_0=b_2=b_4=b_5=0$. Therefore, it can be obtained that:

$$b_1{v}^4+b_3{v}^2=0$$

(50)

In addition to having two zero roots in the speed, there are also two other ones:

$$\begin{array}{ll}v& =\pm \sqrt{-\frac{b_3}{b_1}}=\pm \sqrt{\frac{(e\cdot K\phi +e^2\cdot K\lambda \cdot L_1)\cdot (K\phi \cdot L_1\cdot L_2-e\cdot K\phi \cdot L_2)}{I\cdot L_1(e\cdot K\lambda \cdot L_1^2+K\phi \cdot L_1-e\cdot K\phi -K\lambda \cdot e^2\cdot L_1)}}\\ & =\pm \sqrt{\frac{e\cdot K\phi \cdot L_2}{I\cdot L_1}\frac{(e\cdot K\lambda \cdot L_1+K\phi )}{(e\cdot K\lambda \cdot L_1+K\phi )}\frac{(L_1-e)}{(L_1-e)}}=\pm \sqrt{\frac{e\cdot K\phi \cdot L_2}{I\cdot L_1}}\end{array}$$

(51)

Due to all the parameters being greater than zero, there is a constant critical velocity in the system without considering damping dissipation, which stabilizes the system. We name the sum of the half length of the tire touching the ground and the slack length as the critical stable distance or climbing distance. When the stable distance is less than the crawling distance, we have $b_1=(L_1-e)(K_{\phi }+e\cdot K_{\lambda }\cdot L_1)/(I\cdot L_1^2)>0$, and the stability interval of the system is $0<v_{\mathrm{cr}}\le \sqrt{e\cdot K_{\phi }\cdot L_2/(I\cdot L_1)}$; when the stable distance is equal to the creep distance, we have $b_1=(L_1-e)(K_{\phi }+e\cdot K_{\lambda }\cdot L_1)/(I\cdot L_1^2)=0$, and the stability interval of the system is $v\in R$; and when the stable distance is less than the crawling distance, we have $b_1=(L_1-e)(K_{\phi }+e\cdot K_{\lambda }\cdot L_1)/(I\cdot L_1^2)<0$, and the stability interval of the system is $v_{cr}\ge \sqrt{e\cdot K_{\phi }\cdot L_2/(I\cdot L_1)}$. Therefore, it is possible to draw a stability zone map corresponding to the landing gear stability distance and speed.

From Figure 6, it can be observed that the aircraft remains stable during taxiing at a certain speed when the distance between the tires is greater than the combined length of half of the tire in contact with the ground and the slack length. Therefore, when the problem of aircraft shimmy is extremely prominent, lengthening the landing gear stability distance or switching to smaller-sized tires can be considered.

4.3. Analysis of Shimmy Stability

After the completion of a general landing gear design, modifying its structure poses significant challenges. This modification process entails addressing various aspects, including the dynamic stability performance, static strength concerns, retraction and extension issues, as well as the spatial and fuselage layout of the landing gear cabin. To effectively resolve nose landing gear shimmy while minimizing any impact on the other components, the installation of anti-shimmy dampers is deemed the most efficient approach.

In order to determine the maximum critical damping coefficient, Shen Tianheng was inspired by the formula for finding the root of a one-dimensional quartic equation. He simplified the formula and provided a more convenient criterion for determining the number of real number solutions and the case of multiple roots of the equation. Moreover, Tianheng’s formula does not involve imaginary number square roots, which simplifies the operation and makes it more convenient. This also allows for simplification of the code in computer root-finding tools. For equations $b_0{v}^4+b_1{v}^3+b_2{v}^2+b_3{v}^1+b_4=0$, the double root discriminant is:

$$\{\begin{array}{l}D=3b_1^2-8b_0{b}_2\\ E=-b_1^3+4b_0{b}_1{b}_2-8b_0^2{b}_3\\ F=3b_1^4+16b_0^2{b}_2^2-16b_0{b}_1^2{b}_2+16b_0^2{b}_1{b}_3-64b_0^3{b}_4\\ A=D^2-3F\\ B=DF-9E^2\\ C=F^2-3D{E}^2\\ \Delta =B^2-4AC\end{array}$$

(52)

When $ABC\ne 0$, $\Delta =0$, the equation has a pair of double real roots; if $AB<0$, the other two are conjugate imaginary roots, $x_1=x_2=\frac{-b_1-\frac{2AE}{B}}{4b_0}$ and $x_{3,4}=\frac{-b_1+\frac{2AE}{B}\pm i\sqrt{-\frac{2B}{A}}}{4b_0}$.

$\Delta$ is a complex expression containing the damping coefficient $\mathrm{C}$, which contains the high-order coefficient. In order to obtain the critical value concisely and quickly, the dichotomy method is used to change the damping coefficient value. By judging the sign of $\Delta$, the critical damping coefficient is obtained. Finally, after iteration, $\mathrm{C}$ is 19.8389 Nms. Substituting this value into the characteristic equation, its root can be calculated as y = 21.4322 + 0.0000i, 21.3154 + 0.0000i, −5.4031 + 3.9523i, −5.4031 − 3.9523i. Therefore, in the full speed domain, the maximum required damping corresponds to a speed of approximately 21.4 m/s. As long as the damping coefficient of the stabilizer is greater than 19.8389 Nms, the landing gear oscillation is stable.

To validate the aforementioned conclusion, it is imperative to initially perform a dual parameter analysis of the critical stable state of the landing gear system in order to conduct linear stability analysis of the shimmy. This analytical process entails modifying the values of both the speed and damping coefficient within the system. By means of time-domain analysis and numerical computation, the following can be deduced.

As shown in Figure 7, the analysis reveals that the required anti-swing damping force of the system initially increases and subsequently decreases as the speed escalates. Notably, at a velocity of 21.4 m/s, the maximum damping is achieved with an anti-swing damping value of 19.9 Nms, which aligns consistently with the previous findings. Consequently, it can be inferred that this velocity damping coefficient curve represents the critical damping coefficient essential for ensuring the landing gear’s oscillation stability. Remarkably, the peak critical damping value corresponds to the minimum requirement for overall system stability across all speeds.

4.4. Hopf Bifurcation Analysis

The forward taxiing speed of an aircraft is also an important factor affecting shimmy. Usually, the landing gear system is prone to inducing shimmy within a certain speed range during forward sliding, known as shimmy speed. According to the critical speeds of 8.9839 m/s and 85.1091 m/s calculated in the previous section, the two parameters are substituted into the nonlinear analysis: for the critical speed of 8.9839 m/s, their corresponding eigenvalues are ${\lambda }_{1,2}=\pm 116.13i$ and ${\lambda }_{3,4}=-120.41\pm 96.86i$. The calculated $b_{22}$ is −0.0077. We obtain the first Lyapunov coefficient of the landing gear system according to the formula $l_1(0)=-39.0325<0$; from this, it can be determined that when the speed exceeds 8.9839 m/s, the landing gear system undergoes supercritical Hopf bifurcation. Similarly, for speeds of 85.1091 m/s, their corresponding characteristic values are ${\lambda }_{1,2}=\pm 177.6i$, ${\lambda }_{3,4}=-80.79,-113.46$, and the calculated $b_{22}$ is −0.0061. We obtain the first Lyapunov coefficient of the landing gear system according to the formula $l_1(0)=-71.6482<0$, and from this, it can be determined that a supercritical Hopf bifurcation occurs in the landing gear system when the speed is below 85.1091 m/s. In order to verify the correctness of the theoretical calculation results, the differential equation of the landing gear system was numerically solved to obtain the curve of the maximum real part of the eigenvalues of the landing gear system under the parameters in Table 1 as a function of the sliding speed, as shown in the Figure 8:

From the Figure 8, it can be observed that the maximum real part of the eigenvalues initially decreases as the sliding speed increases. It reaches the minimum value at a speed of 4.65 m/s. Then, it starts to increase again. At a speed of 8.9839 m/s, the maximum real part of the characteristic value of the landing gear system crosses the zero point, changing the stability of the equilibrium point of the landing gear system. This indicates that the minimum critical speed of the landing gear system is 8.9839 m/s. As the speed further increases to 23.95 m/s, the maximum value is reached. Finally, the real part of the system’s eigenvalues gradually decreases. At a speed of 85.1091 m/s, the maximum real part of the eigenvalues of the landing gear system crosses the zero point again, this time from positive to negative, indicating another change. In terms of the stability of the landing gear system, the maximum critical speed is determined to be 8.9839 m/s. The simulation results align with the theoretical calculations, thus confirming the accuracy of the theoretical analysis method presented in this paper.

During the forward taxiing process of an aircraft, shimmy occurs when the parameter values are located at the Hopf bifurcation point of the shimmy dynamics equation. Strictly speaking, it must be a supercritical Hopf bifurcation. In this case, any minor external excitation or disturbance can cause the system to oscillate.

We set the sliding speed to 30 m/s, the linear damping coefficient $C$ = 10 Nms, the nonlinear damping coefficient $C_3$ = 1 Nms³, and the initial swing angle to 0.001 rad. The following curve can be obtained using numerical simulation:

As can be seen in Figure 9, at the current parameter value, the system will gradually approach a state of constant amplitude vibration due to the non-zero initial conditions, leading to the loss of stability in the original neutral equilibrium state. The axes in the phase space represent different state parameters of the system. A point in the phase space signifies the system’s state, while its trajectory forms a phase orbit depicting the evolution of the system’s state over time. Drawing a phase diagram involves representing spatial coordinates for each degree of freedom at a specific time point, enabling analysis of the dynamic motion trends across different spatial states.

Based on the phase diagram of Figure 10, it is evident that the incorporation of non-linear damping into the system results in the convergence of the landing gear swing angle’s amplitude to a constant value during motion progression, irrespective of the initial conditions. Furthermore, all the trajectories depicted in the system phase diagram gradually approach a limit cycle. The presence of nonlinear damping in the landing gear system can induce limit cycles, consequently leading to sustained oscillation.

4.5. The Influence of Parameters on Amplitude

When the system undergoes Hopf bifurcation and generates limit cycles, it is necessary to restrict the range of the oscillation amplitude in order to mitigate aircraft sliding accidents. Theoretically, small constant amplitude vibrations can be effectively suppressed by Coulomb friction. Furthermore, considering the significant influence of the parameters on the limit cycle, it becomes imperative to investigate how parameter variations affect the amplitude of the limit cycle. Firstly, let us examine the impact of the aircraft taxiing speed on the limit cycle.

Observation Figure 11, in systems that are unstable or potentially divergent, the occurrence of a stable limit cycle during motion is expected, as it helps maintain the system stability through constant oscillations. The bifurcation diagram reveals two supercritical Hopf bifurcations, but in opposite directions. Specifically, as the speed increases, a limit cycle begins to emerge, and its amplitude gradually increases until reaching its maximum value. However, with a further increase in speed, the amplitude of the limit cycle starts to decrease until it eventually disappears. Notably, for the first bifurcation value, the stable limit cycle appears on the right side, while for the second bifurcation value, it appears on the left side. These findings align with previous theoretical analysis.

From the bifurcation diagram of Figure 12, it can be seen that as the linear damping increases, the amplitude of the limit cycle begins to decrease. At a linear damping value of 19 Nms, the limit cycle disappears, and the system tends to stabilize. By comparing the results of the linear analysis, it can be observed that when the speed is 30 m/s and the damping coefficient is 19 Nms, all the real parts of the eigenvalues are negative, indicating that the system is stable and convergent.

The primary factor leading to limit cycles in landing gear systems is the magnitude of the nonlinear damping coefficient. In the absence of nonlinear damping, linearly unstable systems would result in an unbounded increase in the landing gear swing angle, which contradicts the current scenario. However, it is observed that many landing gears experience structural damage or even aircraft instability due to excessive torsional deformation. Hence, a comprehensive investigation into the nonlinear damping coefficient across the entire amplitude range of oscillation limit cycle becomes imperative.

As illustrated in the Figure 13, an increase in the nonlinear damping coefficient leads to a gradual reduction in the amplitude of the limit cycle within the system, thereby maintaining the stable state of the limit cycle. When the nonlinear damping coefficient is small, the amplitude of the limit cycle becomes significantly large. However, as the nonlinear damping increases, there is no direct proportionality between its increment and a decrease in amplitude of the limit cycle. Hence, it can be inferred that the presence of nonlinear damping alone ensures sustained low-amplitude vibration without necessitating deliberate adjustment of its coefficient.

We distribute the nonlinear damping coefficient and linear damping coefficient within a specified range and analyze their impact on the amplitude of the limit cycle, as illustrated in the three-dimensional diagram above Figure 14. The system can only stabilize and converge by increasing the linear damping coefficient. Moreover, during Hopf bifurcation, the amplitude of the limit cycle remains relatively small. In terms of nonlinear coefficients, a limit cycle is always generated when there is linear instability; however, this helps mitigate the aircraft sliding risks caused by continuous divergence of the system. Therefore, during the stabilizer design stage, it is crucial to consider both the existence and significance of the damping coefficients.

5. Conclusions

The focus of this article lies in investigating the nose landing gear system of an aircraft. It examines the influence of nonlinear damping, establishes the dynamic equation for the model of the aircraft’s nose landing gear, simplifies the landing gear model using both the central manifold theorem and canonical method, and conducts theoretical analysis on the shimmy stability and Hopf bifurcation characteristics. The following conclusions have been derived.

• Using theoretical analysis and numerical simulation, the accuracy of the derived formula has been validated for linear systems. When the stability distance exceeds the combined length of half of the tire in contact with the ground and the relaxation length, the stable shimmy performance is maintained by the aircraft while sliding at a specific speed. However, when the stability distance is small, only a limited portion of the speed range exhibits stability. Therefore, adjusting either the stability distance or the tire size can alter the stability performance of aircraft shimmy.
• The critical speed of the landing gear system’s linear oscillation was determined based on the stability criterion. Furthermore, a theoretical derivation was conducted to obtain the expression for the first Lyapunov coefficient of the landing gear system. By analyzing its sign, the first Lyapunov coefficient can determine the type of Hopf bifurcation that will occur in the landing gear system without requiring a solution for its bifurcation diagram. It has been demonstrated that higher-order nonlinear damping leads to supercritical Hopf bifurcation in the landing gear system.
• Using quantitative analysis of the parameters’ influence on the limit cycle, the research findings indicate that an increase in velocity induces shimmy motion in the landing gear system through Hopf bifurcation from its equilibrium point. Initially, the amplitude of the limit cycle increases and subsequently decreases, ultimately converging to the equilibrium point via Hopf bifurcation, resulting in cessation of shimmy motion. To mitigate undesired oscillations, it is recommended to enhance the linear damping coefficient within a stability margin range. Similarly, for reducing the limit cycle amplitude, designing the nonlinear damping coefficient within a reasonable range is advisable.