© Wuhan University 2024

0 Introduction

As the operating speed of high-speed trains increases, the issue of hunting instability in train systems becomes increasingly prominent. Hunting motion is an inherent characteristic of train systems. Investigating the causes of hunting motion in trains through research on the nonlinear dynamics of the wheelset system is beneficial for optimizing the design and parameters of trains, thereby improving the overall operational stability, ride comfort, and safety of trains.

In practice, train systems exhibit numerous nonlinear factors^[1-3], including nonlinear constraints of wheel-rail geometry, wheel-rail creep theory, and suspension components and so on. From the perspective of hunting stability issues, the wheelset system can represent a typical subsystem of the overall train system's hunting stability characteristics to a considerable extent. Furthermore, the wheelset system has fewer degrees of freedom and parameters, making it amenable to analytical methods for in-depth analysis, aiding in a better understanding of vehicle dynamic response characteristics.

At present, both domestical and international scholars have undertaken significant research to explore the role of the wheelset system in the hunting motion of high-speed trains. These studies hold promising implications for addressing derailment issues caused by instability motion and improving the overall performance of the entire train system. Zhang et al^[4] addressed the first passage problem of the wheelset system, the influence of noise intensity on mean first passage time and the sensitivity of the yaw damper to the wheelset system at different operating speeds are analyzed. Guo et al^[5] discussed the influence of the nonlinear geometry of the wheel-rail contact on the bifurcation of Hopf and the bifurcation of the limit cycle. The result shows that a larger suspension stiffness would increase the running stability under wheel wear. Li et al^[6] proposed an improved railway wheelset system, the chaotic behaviors of the motion of the railway wheelset and the dynamical behaviors of the railway wheelset system under geometric constraint are studied. Xiong et al^[7] proposed a method to simultaneously diagnose multiple faults such as axle cracks, wheel flats, wheel out-of-roundness, and their coupling faults in the wheelset system. Miao et al^[8] considered the rigid impact model and soft impact model of the railway wheelset system, studying the bifurcation behavior of these two models separately.

However, the role of time delay was not considered in the studies of the above literature. In real situation, the response of the primary suspension of the wheelset system is in place, but it does not reach the corresponding elastic force^[9]. In order to be more in line with the actual situation, this article considers the effect of the time delay in the primary suspension on the system, a dynamic model for a time delay wheelset system under Gaussian white noise excitation conditions is established. The combination of theoretical derivations and numerical simulations to discuss the conditions for the occurrence of hunting instability and stochastic bifurcation in the system. The influence of equivalent conicity and time delay on the occurrence of hunting motion instability in wheelset systems is deeply investigated, and the dynamic behavior before and after the occurrence of hunting instability is analyzed, hoping to provide some theoretical value for train safety and design.

1 Model Description

Considering the influences of lateral wheel-rail forces, wheel-rail normal forces, wheel flange forces, and a series of suspension forces on the wheelset system, we establish a nonlinear wheel-rail contact relationship in a dynamic model of the wheelset, as illustrated in Fig.1. The model neglects the vertical motion of the wheelset and the impact of uneven track. The creep force and creep coefficient between the wheel and rail are based on Kaller linear creep theory^[10], while the nonlinear wheel-rail contact relationship is represented by a fifth-degree polynomial. The springs exhibit linear elasticity, with the wheelset's gravity stiffness and gravity angle stiffness representing the wheel-rail normal forces.

Fig. 1 Model of the wheelset system

Establishing the deterministic motion equations for a wheelset system based on Newton's Second Law yields the following expressions

$\{\begin{array}{l}m\ddot{y}+\frac{\mathrm{2}{f}_{\mathrm{22}}}{V}\dot{y}+K_{\mathrm{1}y}y+\frac{W{\lambda }_{\mathrm{e}}}{b}y-\mathrm{2}{f}_{\mathrm{22}}\psi +{\delta }_{\mathrm{1}}{y}^{\mathrm{3}}+{\delta }_{\mathrm{2}}{y}^{\mathrm{5}}=\mathrm{0},\\ J\ddot{\psi }+\frac{\mathrm{2}{f}_{\mathrm{11}}{b}^{\mathrm{2}}}{V}\dot{\psi }+\frac{\mathrm{2}{f}_{\mathrm{11}}b{\lambda }_{\mathrm{e}}}{r_{\mathrm{0}}}y+(K_{\mathrm{1}x}{b_{\mathrm{0}}}^{\mathrm{2}}-Wb{\lambda }_{\mathrm{e}})\psi =\mathrm{0}\end{array}$(1)

Considering the increase in mileage, corresponding wear occurs on both the wheelset tread and rail surface, leading to variations in the equivalent conicity at contact points. The stochastic perturbation associated with the equivalent conicity is considered as Gaussian white noise. Simultaneously, considering that the system response has reached the primary suspension but the elastic force has not reached the required level, we account for the presence of lateral stiffness time delay in the primary suspension. At this point, the system equations transform into

$\{\begin{array}{l}m\ddot{y}+\frac{\mathrm{2}{f}_{\mathrm{22}}}{V}\dot{y}+K_{\mathrm{1}y}y(t-\tau )+{\epsilon }^{\frac{\mathrm{1}}{\mathrm{2}}}\xi (t)y+\frac{W{\lambda }_{\mathrm{e}}}{b}y-\mathrm{2}{f}_{\mathrm{22}}\psi \\ +{\delta }_{\mathrm{1}}{y}^{\mathrm{3}}+{\delta }_{\mathrm{2}}{y}^{\mathrm{5}}=\mathrm{0},\\ J\ddot{\psi }+\frac{\mathrm{2}{f}_{\mathrm{11}}{b}^{\mathrm{2}}}{V}\dot{\psi }+\frac{\mathrm{2}{f}_{\mathrm{11}}b{\lambda }_{\mathrm{e}}}{r_{\mathrm{0}}}y+(K_{\mathrm{1}x}{b_{\mathrm{0}}}^{\mathrm{2}}-Wb{\lambda }_{\mathrm{e}})\psi =\mathrm{0}\end{array}$(2)

In the given equations: $\xi (t)$ represents Gaussian white noise with a mean of 0 and a noise intensity of $\mathrm{2}D$; $m$ denotes the mass of the wheelset system; $y$ is the lateral displacement of the wheelset; $\psi$ signifies the yaw angle of the wheelset; $V$ stands for the forward speed of the wheelset; $K_{\mathrm{1}y}$ represents lateral stiffness of the primary suspension; $W$ is the axle load; ${\lambda }_{\mathrm{e}}$ denotes the equivalent conicity; $b$ is half of the lateral displacement of the rolling circle; ${\delta }_{\mathrm{1}},{\delta }_{\mathrm{2}}$ correspond to the nonlinear fitting parameters of the wheel flange force; $J$ is the yaw rotational inertia; $f_{\mathrm{11}},f_{\mathrm{22}}$ represents the longitudinal and lateral creep coefficients of the wheel-rail contact; $r_{\mathrm{0}}$ is the rolling radius of the wheel; $K_{\mathrm{1}x}$ signifies longitudinal stiffness of the primary suspension; and $b_{\mathrm{0}}$ is half of the primary yaw spring arm.

2 Derivation Process

2.1 Model Processing

Assuming the equilibrium point of the wheelset system (2) is denoted as $(y_{\mathrm{1}}^{\mathrm{*}},y_{\mathrm{2}}^{\mathrm{*}},y_{\mathrm{3}}^{\mathrm{*}},y_{\mathrm{4}}^{\mathrm{*}})$, we discuss the dynamic behavior of the wheelset system (2) in the vicinity of the equilibrium point, suppose $y=y_{\mathrm{1}}^{\mathrm{*}}+\epsilon x_{\mathrm{1}},\dot{y}=y_{\mathrm{2}}^{\mathrm{*}}+\epsilon x_{\mathrm{2}},\psi =$$y_{\mathrm{3}}^{\mathrm{*}}+\epsilon x_{\mathrm{3}},\dot{\psi }=y_{\mathrm{4}}^{\mathrm{*}}+\epsilon x_{\mathrm{4}}$. The wheelset system can be expressed as follows

$\dot{\mathit{X}}=\mathit{A}\mathit{X}+\mathit{B}\mathit{X}(t-\tau )+\mathit{F}(\mathit{X})$(3)

where $\mathit{X}=\{x_{\mathrm{1}},x_{\mathrm{2}},x_{\mathrm{3}},x_{\mathrm{4}}{\}}^{\mathrm{T}}$

$\mathit{A}=\left(\begin{array}{cccc}\mathrm{0}& \mathrm{1}& \mathrm{0}& \mathrm{0}\\ -a& -\frac{b}{V}& -d& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{1}\\ -e& \mathrm{0}& -f& -\frac{g}{V}\end{array}\right),\mathit{B}=\left(\begin{array}{cccc}\mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ -h& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right),\mathit{F}(\mathit{X})=\left(\begin{array}{c}\mathrm{0}\\ -{\epsilon }^{\frac{\mathrm{1}}{\mathrm{2}}}\frac{x_{\mathrm{1}}\xi (t)}{m}+\epsilon (m_{\mathrm{1}}{x}_{\mathrm{1}}^{\mathrm{3}}+m_{\mathrm{2}}{x}_{\mathrm{1}}^{\mathrm{5}})\\ \mathrm{0}\\ \mathrm{0}\end{array}\right)$

$a=\frac{W{\lambda }_{\mathrm{e}}}{mb},c=\frac{\mathrm{2}{f}_{\mathrm{22}}}{m},d=-\frac{\mathrm{2}{f}_{\mathrm{22}}}{m},e=\frac{\mathrm{2}{f}_{\mathrm{11}}b{\lambda }_{\mathrm{e}}}{J{r}_{\mathrm{0}}},f=\frac{K_{\mathrm{1}x}{b_{\mathrm{0}}}^{\mathrm{2}}-Wb{\lambda }_{\mathrm{e}}}{J},g=\frac{\mathrm{2}{f}_{\mathrm{11}}{b}^{\mathrm{2}}}{J},h=\frac{K_{\mathrm{1}y}}{m},m_{\mathrm{1}}=-\frac{{\delta }_{\mathrm{1}}\epsilon }{m},m_{\mathrm{2}}=-\frac{{\delta }_{\mathrm{1}}{\epsilon }^{\mathrm{3}}}{m}.$

The characteristic equation of the linear part of the wheelset system (3) is expressed as follows:

$\mathrm{\Delta }(\lambda ,V)={\lambda }^{\mathrm{4}}+\frac{g+c}{V}{\lambda }^{\mathrm{3}}+(\frac{cg}{V^{\mathrm{2}}}+f+a+h{\mathrm{e}}^{-\lambda \tau }){\lambda }^{\mathrm{2}}+\frac{cf+ag+hg{\mathrm{e}}^{-\lambda \tau }}{V}\lambda +(af-de+hf{\mathrm{e}}^{-\lambda \tau })=\mathrm{0}$(4)

When the speed $V$ reaches the critical speed $V_{\mathrm{c}}$, the characteristic Eq. (4) has a pair of complex conjugate eigenvalues $\pm \mathrm{i}\omega$ with zero real parts at the equilibrium point. Additionally, all other eigenvalues of the characteristic Eq. (4) have negative real parts. Substituting $\lambda =\mathrm{i}\omega$ into the characteristic Eq. (4) yields

${(\mathrm{i}\omega )}^{\mathrm{4}}+\frac{g+c}{V}{(\mathrm{i}\omega )}^{\mathrm{3}}+(\frac{cg}{V^{\mathrm{2}}}+f+a+h{\mathrm{e}}^{-\mathrm{i}\omega \tau }){(\mathrm{i}\omega )}^{\mathrm{2}}+\frac{cf+ag+hg{\mathrm{e}}^{-\mathrm{i}\omega \tau }}{V}(\mathrm{i}\omega )+(af-de+hf{\mathrm{e}}^{-\mathrm{i}\omega \tau })=\mathrm{0}$(5)

After separating the real and imaginary components, the following results can be obtained

$(\frac{hg}{V}\omega )\mathrm{s}\mathrm{i}\mathrm{n}\omega \tau +(hf-h{\omega }^{\mathrm{2}})\mathrm{c}\mathrm{o}\mathrm{s}\omega \tau =-{\omega }^{\mathrm{4}}+(\frac{cg}{V^{\mathrm{2}}}+f+a){\omega }^{\mathrm{2}},(\frac{hg}{V}\omega )\mathrm{c}\mathrm{o}\mathrm{s}\omega \tau -(hf-h{\omega }^{\mathrm{2}})\mathrm{s}\mathrm{i}\mathrm{n}\omega \tau =\frac{g+c}{V}{\omega }^{\mathrm{3}}-\frac{cf+ag}{V}\omega$(6)

Squaring the two equations of Eq. (6) and summing them, we get

${(\frac{hg}{V}\omega )}^{\mathrm{2}}+(hf-h{\omega }^{\mathrm{2}}{)}^{\mathrm{2}}=(-{\omega }^{\mathrm{4}}+(\frac{cg}{V^{\mathrm{2}}}+f+a){\omega }^{\mathrm{2}}{)}^{\mathrm{2}}+(\frac{g+c}{V}{\omega }^{\mathrm{3}}{-\frac{cf+ag}{V}\omega )}^{\mathrm{2}}$(7)

According to the implicit function theorem, differentiating both sides of Eq. (4) with respect to $V$, we get the transversality condition

${\left(\frac{\mathrm{d}\lambda }{\mathrm{d}V}\right)|}_{\lambda =\mathrm{i}\omega }=\frac{\frac{g+c}{V^{\mathrm{2}}}{\lambda }^{\mathrm{3}}+\frac{\mathrm{2}cg}{V^{\mathrm{3}}}{\lambda }^{\mathrm{2}}+\frac{cf+ag+hg{\mathrm{e}}^{-\lambda \tau }}{V^{\mathrm{2}}}\lambda }{\left(\mathrm{4}{\lambda }^{\mathrm{3}}+\left[\mathrm{3}\left(\frac{g+c}{V}\right)-h\tau {\mathrm{e}}^{-\lambda \tau }\right]{\lambda }^{\mathrm{2}}+\left[\mathrm{2}\left(\frac{cg}{V^{\mathrm{2}}}+f+a+h{\mathrm{e}}^{-\lambda \tau }\right)-\frac{hg\tau {\mathrm{e}}^{-\lambda \tau }}{V}\right]\lambda +\frac{cf+ag+hg{\mathrm{e}}^{-\lambda \tau }}{V}-hf\tau {\mathrm{e}}^{-\lambda \tau }\right)}$(8)

2.2 Dimensionality Reduction Processing

Next, we will analyze the stochastic dynamics of the wheelset system (3) near the critical speed $V=V_{\mathrm{c}}$ and investigate the hunting stability of the system. Since the solution space of the wheelset system (3) is infinite-dimensional, by mapping the interval $[-\tau ,\mathrm{0}]$ onto a continuous space function ${\mathbb{R}}^{\mathrm{2}}$, the wheelset system (3) can be expressed in the following form

$\dot{\mathit{X}}=\mathit{L}(\phi (\theta ),V)+\epsilon \mathrm{\Delta }\mathit{f}(\phi (\theta ),V,\epsilon ),\mathit{L}(\phi (\theta ),V)=\mathit{A}\mathit{X}+\mathit{B}\mathit{X}(t-\tau )$(9)

where $\mathit{L}(\phi (\theta ),V)$ and $\mathit{f}(\phi (\theta ),V,\epsilon )$ represent the linear and nonlinear parts of the wheelset system (3), respectively.

Let ${\mathit{\Psi }}_j(s)\in \overline{C}([\mathrm{0},\tau ],{\mathbb{R}}^{\mathrm{2}}),{\mathit{\Phi }}_k\in ([-\tau ,\mathrm{0}],{\mathbb{R}}^{\mathrm{2}})$, and define the bilinear operator

$<{\mathit{\Psi }}_j(s),{\mathit{\Phi }}_k(\theta )>={\mathit{\Psi }}_j(s)\cdot {\mathit{\Phi }}_k(\theta )-h{\int }_{-\tau }^{\mathrm{0}}{\mathit{\Psi }}_{j\mathrm{2}}(\chi +\tau )\cdot {\mathit{\Phi }}_{k\mathrm{1}}(\chi )\mathrm{d}\chi ,j,k=\mathrm{1,2}$(10)

where ${\mathit{\Psi }}_{j\mathrm{2}}$ and ${\mathit{\Phi }}_{k\mathrm{1}}$ represent the second element in ${\mathit{\Psi }}_j$ and the first element in ${\mathit{\Phi }}_k$, respectively.

When the characteristic value $\lambda =\pm \mathrm{i}\omega$, the corresponding basis function $\mathit{\Phi }(\theta )\in C$ of the characteristic Eq. (4) and the associated basis function $\mathit{\Psi }(s)\in \overline{C}$ of the adjoint equation are as follows:

$\begin{array}{l}\mathit{\Phi }(\theta )=\left(\begin{array}{cc}\mathrm{c}\mathrm{o}\mathrm{s}\omega \theta & \mathrm{s}\mathrm{i}\mathrm{n}\omega \theta \\ -\omega \mathrm{s}\mathrm{i}\mathrm{n}\omega \theta & \omega \mathrm{c}\mathrm{o}\mathrm{s}\omega \theta \\ b_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\omega \theta -c_{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\omega \theta & b_{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\omega \theta +c_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\omega \theta \\ b_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\omega \theta -c_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\omega \theta & b_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\omega \theta +c_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\omega \theta \end{array}\right)={\left(\begin{array}{c}{\mathit{\Phi }}_{\mathrm{1}}(\theta )\\ {\mathit{\Phi }}_{\mathrm{2}}(\theta )\end{array}\right)}^{\mathrm{T}},-\tau \le \theta \le \mathrm{0},\\ \mathit{\Psi }(s)={\left(\begin{array}{cc}\mathrm{c}\mathrm{o}\mathrm{s}\omega s& \mathrm{s}\mathrm{i}\mathrm{n}\omega s\\ -\omega \mathrm{s}\mathrm{i}\mathrm{n}\omega s& \omega \mathrm{c}\mathrm{o}\mathrm{s}\omega s\\ b_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s-c_{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s& b_{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s+c_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s\\ b_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s-c_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s& b_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s+c_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s\end{array}\right)}^{\mathrm{T}}=\left(\begin{array}{c}{\mathit{\Psi }}_{\mathrm{1}}(s)\\ {\mathit{\Psi }}_{\mathrm{2}}(s)\end{array}\right),\mathrm{0}\le s\le \tau ,\end{array}$(11)

where

$b_{\mathrm{2}}=\frac{a+h\mathrm{c}\mathrm{o}\mathrm{s}\omega \tau -{\omega }^{\mathrm{2}}}{d},c_{\mathrm{2}}=\frac{\frac{c}{V}\omega -h\mathrm{s}\mathrm{i}\mathrm{n}\omega \tau }{d},b_{\mathrm{3}}=-\omega c_{\mathrm{2}},c_{\mathrm{3}}=\omega b_{\mathrm{2}}$

Substituting the base function (11) into the bilinear expression (10) yields the following

$\left(\begin{array}{cc}<{\mathit{\Psi }}_{\mathrm{1}}(s),{\mathit{\Phi }}_{\mathrm{1}}(\theta )>& <{\mathit{\Psi }}_{\mathrm{1}}(s),{\mathit{\Phi }}_{\mathrm{2}}(\theta )>\\ <{\mathit{\Psi }}_{\mathrm{2}}(s),{\mathit{\Phi }}_{\mathrm{1}}(\theta )>& <{\mathit{\Psi }}_{\mathrm{2}}(s),{\mathit{\Phi }}_{\mathrm{2}}(\theta )>\end{array}\right)={(\mathit{\Psi },\mathit{\Phi })}_{nsg}=\left(\begin{array}{cc}{a}_{\mathrm{11}}& a_{\mathrm{12}}\\ a_{\mathrm{21}}& a_{\mathrm{22}}\end{array}\right)$(12)

where

$\begin{array}{l}{a}_{\mathrm{11}}=\mathrm{1}+b_{\mathrm{2}}^{\mathrm{2}}+b_{\mathrm{3}}^{\mathrm{2}}+\frac{h\tau \omega \mathrm{s}\mathrm{i}\mathrm{n}\omega \tau }{\mathrm{2}},a_{\mathrm{12}}=b_{\mathrm{2}}{c}_{\mathrm{2}}+b_{\mathrm{3}}{c}_{\mathrm{3}}-h(\frac{\mathrm{s}\mathrm{i}\mathrm{n}\omega \tau }{\mathrm{2}}-\frac{\tau \omega \mathrm{c}\mathrm{o}\mathrm{s}\omega \tau }{\mathrm{2}}),\\ a_{\mathrm{21}}=b_{\mathrm{2}}{c}_{\mathrm{2}}+b_{\mathrm{3}}{c}_{\mathrm{3}}-h(\frac{\mathrm{s}\mathrm{i}\mathrm{n}\omega \tau }{\mathrm{2}}+\frac{\tau \omega \mathrm{c}\mathrm{o}\mathrm{s}\omega \tau }{\mathrm{2}}),a_{\mathrm{22}}={\omega }^{\mathrm{2}}+c_{\mathrm{2}}^{\mathrm{2}}+c_{\mathrm{3}}^{\mathrm{2}}+\frac{h\tau \omega \mathrm{s}\mathrm{i}\mathrm{n}\omega \tau }{\mathrm{2}}.\end{array}$(13)

Normalize the basis function $\mathit{\Psi }(s)$ to obtain a new basis function $\overline{\mathit{\Psi }}(s)$,

$\overline{\mathit{\Psi }}(s)={(\overline{\mathit{\Psi }},\mathit{\Phi })}_{nsg}^{-\mathrm{1}}=\left(\begin{array}{cccc}{A}_{\mathrm{11}}& A_{\mathrm{12}}& A_{\mathrm{13}}& A_{\mathrm{14}}\\ A_{\mathrm{21}}& A_{\mathrm{22}}& A_{\mathrm{23}}& A_{\mathrm{24}}\end{array}\right),\mathrm{0}\le s\le \tau$(14)

where

$\begin{array}{l}{A}_{\mathrm{11}}(s)=(a_{\mathrm{11}}{a}_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}}{)}^{-\mathrm{1}}(a_{\mathrm{22}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s-a_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s),\\ A_{\mathrm{12}}(s)=-(a_{\mathrm{11}}{a}_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}}{)}^{-\mathrm{1}}\omega (a_{\mathrm{22}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s+a_{\mathrm{12}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s),\\ A_{\mathrm{13}}(s)=(a_{\mathrm{11}}{a}_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}}{)}^{-\mathrm{1}}[a_{\mathrm{22}}(b_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s-c_{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s)-a_{\mathrm{12}}(b_{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s+c_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s)],\\ A_{\mathrm{14}}(s)=(a_{\mathrm{11}}{a}_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}}{)}^{-\mathrm{1}}[a_{\mathrm{22}}(b_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s-c_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s)-a_{\mathrm{12}}(b_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s+c_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s)],\\ A_{\mathrm{21}}(s)=(a_{\mathrm{11}}{a}_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}}{)}^{-\mathrm{1}}(-a_{\mathrm{21}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s+a_{\mathrm{11}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s),\\ A_{\mathrm{22}}(s)=(a_{\mathrm{11}}{a}_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}}{)}^{-\mathrm{1}}\omega (a_{\mathrm{21}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s+a_{\mathrm{11}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s),\\ A_{\mathrm{23}}(s)=(a_{\mathrm{11}}{a}_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}}{)}^{-\mathrm{1}}[-a_{\mathrm{21}}(b_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s-c_{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s)+a_{\mathrm{11}}(b_{\mathrm{2}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s+c_{\mathrm{2}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s)],\\ A_{\mathrm{24}}(s)=(a_{\mathrm{11}}{a}_{\mathrm{22}}-a_{\mathrm{12}}{a}_{\mathrm{21}}{)}^{-\mathrm{1}}[-a_{\mathrm{21}}(b_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s-c_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s)-a_{\mathrm{11}}(b_{\mathrm{3}}\mathrm{s}\mathrm{i}\mathrm{n}\omega s+c_{\mathrm{3}}\mathrm{c}\mathrm{o}\mathrm{s}\omega s)].\end{array}$(15)

Assuming the unique solution of the system (3) can be represented as ${\mathit{x}}_t(\varphi (\theta ),V,\epsilon )={\mathit{x}}_t^{\mathit{P}}(\varphi (\theta ),V,\epsilon )+{\mathit{x}}_t^{\mathit{Q}}(\varphi (\theta ),V,\epsilon )$, where $\varphi (\theta )={\varphi }^{\mathit{P}}(\theta )+{\varphi }^{\mathit{Q}}(\theta )$ lies in the semigroup $J(t,\beta )$ and its infinitesimal generator $A(\theta ,V)$, and ${\mathit{x}}_t^{\mathit{P}}(\varphi (\theta ),V,\epsilon )$$\in \mathit{P},{\mathit{x}}_t^{\mathit{Q}}(\varphi (\theta ),V,\epsilon )\in \mathit{Q}$ are the projections onto ${\mathit{x}}_{\mathit{t}}(\varphi (\theta ),V,\epsilon )$ in two subspaces $\mathit{P},\mathit{Q}$, respectively. According to the relation $\mathit{A}(\theta ,p)\mathit{\Phi }(\theta )=\mathit{\Phi }(\theta )\mathit{B}$, we computed $\mathit{B}=[{\left(\mathrm{0},\omega \right)}^{\mathrm{T}},{\left(-\omega ,\mathrm{0}\right)}^{\mathrm{T}}]$. The approximate expression for transformation ${\varphi }^P(\theta )=\mathit{\Phi }(\theta )\mathit{B}\in \mathit{C}$ on the central manifold ${\mathit{M}}_{\mathit{V}}\in \mathit{C}([-\tau ,\mathrm{0}],{\mathbb{R}}^{\mathrm{2}})$ and variable changes ${\mathit{x}}_t^{\mathit{P}}(\varphi (\theta ),V,\epsilon )=\varphi (\theta )\mathit{u}(t)+{\mathit{x}}_{\mathit{t}}^{\mathit{Q}}($$\varphi (\theta ),V,\epsilon )$ are written as follows:

$\begin{array}{l}{x}_{\mathrm{1}}(t)=u_{\mathrm{1}}(t),x_{\mathrm{2}}(t)=\omega u_{\mathrm{2}}(t),\\ x_{\mathrm{3}}(t)=b_{\mathrm{2}}{u}_{\mathrm{1}}(t)+c_{\mathrm{2}}{u}_{\mathrm{2}}(t),x_{\mathrm{4}}(t)=b_{\mathrm{3}}{u}_{\mathrm{1}}(t)+c_{\mathrm{3}}{u}_{\mathrm{2}}(t).\end{array}$(16)

The expression form of the system (3) on the central manifold^[11] is as follows:

${\dot{u}}_{\mathrm{1}}(t)=\omega u_{\mathrm{2}}(t)+\epsilon [(m_{\mathrm{1}}{u}_{\mathrm{1}}^{\mathrm{3}}+m_{\mathrm{2}}{u}_{\mathrm{1}}^{\mathrm{5}})A_{\mathrm{12}}]-{\epsilon }^{\frac{\mathrm{1}}{\mathrm{2}}}[A_{\mathrm{12}}\frac{u_{\mathrm{1}}(t)}{m}]\xi (t),{\dot{u}}_{\mathrm{2}}(t)=-\omega u_{\mathrm{1}}(t)+\epsilon [(m_{\mathrm{1}}{u}_{\mathrm{1}}^{\mathrm{3}}+m_{\mathrm{2}}{u}_{\mathrm{1}}^{\mathrm{5}})A_{\mathrm{22}}]-{\epsilon }^{\frac{\mathrm{1}}{\mathrm{2}}}[A_{\mathrm{22}}\frac{u_{\mathrm{1}}(t)}{m}]\xi (t)$(17)

Introducing a polar coordinate transformation^[12]: $u_{\mathrm{1}}(t)=A(t)\mathrm{c}\mathrm{o}\mathrm{s}\theta ,u_{\mathrm{2}}(t)=-A(t)\mathrm{s}\mathrm{i}\mathrm{n}\theta ,\theta =\omega t+\varphi (t)$. The wheelset system (3) is transformed into an equation concerning the amplitude $A(t)$ and phase $\varphi (t)$, as follows:

$\{\begin{array}{l}\dot{A}(t)=\epsilon [A_{\mathrm{12}}(m_{\mathrm{1}}{A}^{\mathrm{3}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{4}}\theta +m_{\mathrm{2}}{A}^{\mathrm{5}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{6}}\theta )-A_{\mathrm{22}}(m_{\mathrm{1}}{A}^{\mathrm{3}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{3}}\theta \mathrm{s}\mathrm{i}\mathrm{n}\theta +m_{\mathrm{2}}{A}^{\mathrm{5}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{5}}\theta \mathrm{s}\mathrm{i}\mathrm{n}\theta )]+{\epsilon }^{\frac{\mathrm{1}}{\mathrm{2}}}A\mathrm{c}\mathrm{o}\mathrm{s}\theta \frac{-A_{\mathrm{12}}\mathrm{c}\mathrm{o}\mathrm{s}\theta +A_{\mathrm{22}}\mathrm{s}\mathrm{i}\mathrm{n}\theta }{m}\xi (t),\\ \dot{\phi }(t)=-\epsilon [A_{\mathrm{12}}(m_{\mathrm{1}}{A}^{\mathrm{2}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{3}}\theta \mathrm{s}\mathrm{i}\mathrm{n}\theta +m_{\mathrm{2}}{A}^{\mathrm{4}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{5}}\theta \mathrm{s}\mathrm{i}\mathrm{n}\theta )-A_{\mathrm{22}}(m_{\mathrm{1}}{A}^{\mathrm{2}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{4}}\theta +m_{\mathrm{2}}{A}^{\mathrm{4}}\mathrm{c}\mathrm{o}{\mathrm{s}}^{\mathrm{6}}\theta )]+{\epsilon }^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{c}\mathrm{o}\mathrm{s}\theta \frac{A_{\mathrm{12}}\mathrm{s}\mathrm{i}\mathrm{n}\theta +A_{\mathrm{22}}\mathrm{c}\mathrm{o}\mathrm{s}\theta }{m}\xi (t).\end{array}$(18)

According to the stochastic averaging method^[13,14], Itô stochastic differential equation for the wheelset system (3) can be derived

$\mathrm{d}A=m(A)\mathrm{d}t+\sigma (A)\mathrm{d}B(t)$(19)

where $B(t)$ is a Wiener process with unit intensity, while $m(A)$ and $\sigma (A)$ represent the drift and diffusion coefficients of the system (18), respectively.

$\begin{array}{c}m(A)={\mu }_{\mathrm{1}}A+{\mu }_{\mathrm{2}}{A}^{\mathrm{3}}+{\mu }_{\mathrm{3}}{A}^{\mathrm{5}},\sigma {\sigma }^{\mathrm{T}}={\mu }_{\mathrm{4}}{A}^{\mathrm{2}}\\ {\mu }_{\mathrm{1}}=(\frac{\mathrm{3}{A}_{\mathrm{12}}^{\mathrm{2}}+A_{\mathrm{22}}^{\mathrm{2}}}{\mathrm{8}{m}^{\mathrm{2}}}+\frac{A_{\mathrm{12}}^{\mathrm{2}}+A_{\mathrm{22}}^{\mathrm{2}}}{\mathrm{4}{m}^{\mathrm{2}}})\mathrm{\pi }D,{\mu }_{\mathrm{2}}=\frac{\mathrm{3}}{\mathrm{8}}{A}_{\mathrm{12}}{m}_{\mathrm{1}},{\mu }_{\mathrm{3}}=\frac{\mathrm{5}}{\mathrm{16}}{A}_{\mathrm{12}}{m}_{\mathrm{2}},{\mu }_{\mathrm{4}}=\frac{\mathrm{3}{A}_{\mathrm{12}}^{\mathrm{2}}+A_{\mathrm{22}}^{\mathrm{2}}}{\mathrm{4}{m}^{\mathrm{2}}}\mathrm{\pi }D.\end{array}$(20)

3 Stochastic Stability Analysis

3.1 The Local Stochastic Stability

According to the definition of the Lyapunov exponent, the maximum Lyapunov exponent for the Itô stochastic differential equation corresponding to the wheelset system (3) is calculated as follows

$\lambda =\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\mathrm{1}}{t}\Vert A(t)\Vert =\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}(A{(t))}^{\frac{\mathrm{1}}{\mathrm{2}}}=\frac{\mathrm{1}}{\mathrm{2}}\left(m^{\text{'}}(\mathrm{0})-\frac{\mathrm{1}}{\mathrm{2}}({\sigma }^{\text{'}}{(\mathrm{0}))}^{\mathrm{2}}\right)=\frac{\mathrm{1}}{\mathrm{2}}({\mu }_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}}{\mu }_{\mathrm{4}}).$(21)

When ${\mu }_{\mathrm{1}}-{\mu }_{\mathrm{4}}/\mathrm{2}<\mathrm{0}$, the equilibrium solution of the wheelset system (3) is locally asymptotically stable in a probabilistic sense. Conversely, when ${\mu }_{\mathrm{1}}-{\mu }_{\mathrm{4}}/\mathrm{2}>\mathrm{0}$, the equilibrium solution of the wheelset (3) is probabilistically locally asymptotically unstable.

3.2 The Global Stochastic Stability

According to the theory of singular boundaries, the right boundary at $A=+\mathrm{\infty }$ is a second-type singular boundary. Analysis of Itô stochastic differential equations yields

${\alpha }_R=\mathrm{2},{\beta }_R=\mathrm{5},c_R=-\underset{A\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\mathrm{2}m(A)|A{|}^{{\alpha }_R-{\beta }_R}}{{\sigma }^{\mathrm{2}}(A)}=-\frac{\mathrm{2}{\mu }_{\mathrm{3}}}{{\mu }_{\mathrm{4}}}.$(22)

The left boundary at $A=\mathrm{0}$ is a first-type singular boundary. Analysis of Itô stochastic differential equations yields

${\alpha }_L=\mathrm{2},{\beta }_L=\mathrm{1},c_L=\underset{A\to \mathrm{0}}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\mathrm{2}m(A)|A{|}^{{\beta }_L-{\alpha }_L}}{{\sigma }^{\mathrm{2}}(A)}=\frac{\mathrm{2}{\mu }_{\mathrm{1}}}{{\mu }_{\mathrm{4}}}.$(23)

According to the three-exponential method, when $m(+\mathrm{\infty })<\mathrm{0}$ or $m(-\mathrm{\infty })>\mathrm{0}$, the right boundary $A=+\mathrm{\infty }$ is an entry into the natural boundary. In this case, we only need to consider the left boundary of the Itô stochastic differential equation.

Through analysis, it can be obtained that when $c_L<\mathrm{1}$, the left boundary $A=\mathrm{0}$ is an attractive natural boundary, and the equilibrium solution of the wheelset system is globally asymptotically stable in a probabilistic sense. Conversely, when $c_L>\mathrm{1}$, the left boundary $A=\mathrm{0}$ is a repulsive natural boundary, and the equilibrium solution of the wheel system is globally unstable in a probabilistic sense. In the case of $c_L=\mathrm{1}$, the left boundary $A=\mathrm{0}$ becomes a strict natural boundary, situated at a critical state between repulsion and attraction, where stochastic bifurcations may occur.

4 Stochastic Bifurcation Analysis

4.1 Stochastic D-Bifurcation

The occurrence of stochastic D-bifurcations in the wheelset system can be inferred based on the maximum Lyapunov exponent. When ${\mu }_{\mathrm{4}}=\mathrm{0}$, Itô stochastic differential equation to the wheelset system is a deterministic system that does not bifurcate. Subsequently, we will explore the stochastic D-bifurcation of the wheelset system under the condition ${\mu }_{\mathrm{4}}\ne \mathrm{0}$.

In scenario one, when ${\mu }_{\mathrm{2}}={\mu }_{\mathrm{3}}=\mathrm{0}$, the corresponding linear Itô stochastic differential equation to the wheelset system is

$\mathrm{d}A={\mu }_{\mathrm{1}}A\mathrm{d}t+({\mu }_{\mathrm{4}}{A}^{\mathrm{2}}{)}^{\frac{\mathrm{1}}{\mathrm{2}}}\mathrm{d}B(t).$(24)

According to the relationship between the Stratonovich stochastic differential equation and the Itô stochastic differential equation, let $m(A)=({\mu }_{\mathrm{1}}-{\mu }_{\mathrm{4}}/\mathrm{2})A,\sigma (A)=\sqrt[]{{\mu }_{\mathrm{4}}}A$, the continuous dynamic system generated by equation (24) is given by

$\eta (t)x=x{\int }_{\mathrm{0}}^{t}m(\eta (s)x)\mathrm{d}s+{\int }_{\mathrm{0}}^t\sigma (\eta (s)x)\circ \mathrm{d}B(s).$(25)

Eq. (25) is the unique strong solution of Eq. (24) with x as the initial value, where $\circ \mathrm{d}B(s)$ is the stochastic term of the Stratonovich stochastic differential equation^[13,14]. When $m(\mathrm{0})=\mathrm{0},\sigma (\mathrm{0})=\mathrm{0}$, 0 is as a fixed point for $\eta (t)$. If $m(A)$ is a bounded function, then when $\sigma (A)\ne \mathrm{0}$, the stationary probability density function for all $A\ne \mathrm{0}$ is unique. The FPK (Fokker-Planck-Kolmogorov) equation corresponding to amplitude $A(t)$ is given as follows:

$\frac{\partial P}{\partial A}=-\frac{\partial }{\partial A}\left(\left(({\mu }_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}}{\mu }_{\mathrm{4}})A\right)P\right)+\frac{\mathrm{1}}{\mathrm{2}}\frac{{\partial }^{\mathrm{2}}}{\partial A^{\mathrm{2}}}\left(({\mu }_{\mathrm{4}}{A}^{\mathrm{2}}{)}^{\frac{\mathrm{1}}{\mathrm{2}}}P\right)$(26)

Let $\partial P/\partial A=\mathrm{0}$, the stationary probability density function of Eq. (26) can be computed as follows:

$P(A)=\frac{C}{{\sigma }^{\mathrm{2}}(A)}\mathrm{e}\mathrm{x}\mathrm{p}\left({\int }_{\mathrm{0}}^A\frac{\mathrm{2}m(x)}{{\sigma }^{\mathrm{2}}(x)}\mathrm{d}x\right)$(27)

where $C$ is the normalization constant, and Eq. (25) has two stationary states: fixed point stationary state and nontrivial stationary state, where $w$ represents the invariant measure of the nontrivial stationary state, and its probability density function is Eq. (27), the invariant measure of the fixed point stationary state is represented by ${\delta }_{\mathrm{0}}$, and the stationary probability density function is represented by ${\delta }_x$. Solving the linear Itô stochastic differential equation (24), we get

$A(t)=A(\mathrm{0})\mathrm{e}\mathrm{x}\mathrm{p}\left({\int }_{\mathrm{0}}^t\left(m^{\text{'}}(A)+\frac{\sigma (A){\sigma }^{\text{'}}(A)}{\mathrm{2}}\right)\mathrm{d}x+{\int }_{\mathrm{0}}^t{\sigma }^{\text{'}}(A)\mathrm{d}B\right).$(28)

Below, we calculate the maximum Lyapunov exponents concerning the invariant measures ${\delta }_{\mathrm{0}}$ and $w$. The maximum Lyapunov exponent for system $\eta$ with respect to measure $\varsigma$ is defined as follows:

${\lambda }_{\eta }(\varsigma )=\underset{x\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\mathrm{1}}{t}\mathrm{l}\mathrm{n}\Vert A(t)\Vert .$(29)

Due to $\sigma (\mathrm{0})=\mathrm{0},{\sigma }^{\text{'}}(\mathrm{0})=\mathrm{0}$, the maximum Lyapunov exponent concerning the invariant measure ${\delta }_{\mathrm{0}}$ is as follows:

${\lambda }_{\eta }({\delta }_{\mathrm{0}})=\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\mathrm{1}}{t}(\mathrm{l}\mathrm{n}\Vert m(\mathrm{0})\Vert +m^{\text{'}}(\mathrm{0}){\int }_{\mathrm{0}}^t\mathrm{d}x+{\sigma }^{\text{'}}(\mathrm{0}){\int }_{\mathrm{0}}^t\mathrm{d}{B}_A(x)={\mu }_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}}{\mu }_{\mathrm{4}}.$(30)

Substituting Eq. (28) into Eq. (29), we obtain the maximum Lyapunov exponent concerning the invariant measure $w$

$\begin{array}{l}{\lambda }_{\eta }(w)=\underset{t\to +\mathrm{\infty }}{\mathrm{l}\mathrm{i}\mathrm{m}}\frac{\mathrm{1}}{t}{\int }_{\mathrm{0}}^t\left(m^{\text{'}}(A)+\sigma (A){\sigma }^{\text{'}}(A)\right)\mathrm{d}x={\int }_R\left(m^{\text{'}}(A)+\frac{\sigma (A){\sigma }^{\text{'}}(A)}{\mathrm{2}}\right)P(A)\mathrm{d}A\\ \begin{array}{cc}& \end{array}=-\mathrm{2}{m}^{\text{'}}(\mathrm{0}){{\int }_R\left(\frac{m(A)}{\sigma (A)}\right)}^{\mathrm{2}}P(A)\mathrm{d}A=-{\mu }_{\mathrm{4}}^{\frac{\mathrm{3}}{\mathrm{2}}}\left({\mu }_{\mathrm{1}}-\frac{{\mu }_{\mathrm{4}}}{\mathrm{2}}\right)\mathrm{e}\mathrm{x}\mathrm{p}\left(\frac{\mathrm{2}}{{\mu }_{\mathrm{4}}}\left({\mu }_{\mathrm{1}}-\frac{{\mu }_{\mathrm{4}}}{\mathrm{2}}\right)\right).\end{array}$(31)

Let $\alpha ={\mu }_{\mathrm{1}}-{\mu }_{\mathrm{4}}/\mathrm{2}$, by analyzing Eq. (30) and Eq. (31), we obtain that when $\alpha <\mathrm{0}$, the invariant measure ${\delta }_{\mathrm{0}}$ is stable, while the invariant measure $w$ is unstable, when $\alpha >\mathrm{0}$, the invariant measure ${\delta }_{\mathrm{0}}$ becomes unstable, while the invariant measure $w$ is stable. Consequently, a stochastic D-bifurcation occurs in the wheelset system when $\alpha =\mathrm{0}$.

In scenario two, when ${\mu }_{\mathrm{2}}\ne \mathrm{0},{\mu }_{\mathrm{3}}=\mathrm{0}$, let $\varphi =\sqrt[]{-{\mu }_{\mathrm{2}}}A,{\mu }_{\mathrm{2}}<\mathrm{0}$, Eq. (19) transforms into

$\mathrm{d}\varphi =\left(\left({\mu }_{\mathrm{1}}-\frac{\mathrm{1}}{\mathrm{2}}{\mu }_{\mathrm{4}}\right)\varphi -{\varphi }^{\mathrm{3}}\right)\mathrm{d}t+{{\mu }_{\mathrm{4}}}^{\frac{\mathrm{1}}{\mathrm{2}}}\varphi \circ \mathrm{d}B(s).$(32)

When ${\mu }_{\mathrm{4}}=\mathrm{0}$, the system (32) has a deterministic fork bifurcation. When ${\mu }_{\mathrm{4}}\ne \mathrm{0},\alpha <\mathrm{0}$, Eq. (32) only has an invariant measure ${\delta }_{\mathrm{0}}$ of the fixed point stationary state, with a probability density function ${\delta }_x$. When ${\mu }_{\mathrm{4}}\ne \mathrm{0},\alpha >\mathrm{0}$, Eq. (32) has three invariant measures of stationary states: the invariant measure ${\delta }_{\varphi }$ of the stationary state of the fixed point, and the invariant measure $w_{\varphi }^{\pm }$ of the two nontrivial stationary states. Their probability density functions are as follows

$P_w^{+}(\varphi )=\{\begin{array}{c}{N}_w{\varphi }^{\frac{\mathrm{2}\alpha }{{\mu }_{\mathrm{4}}}-\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{{\varphi }^{\mathrm{2}}}{{\mu }_{\mathrm{4}}}\right),\begin{array}{ccc}& & \varphi >\mathrm{0}\end{array},\\ \begin{array}{ccc}\begin{array}{ccc}\begin{array}{cc}\begin{array}{ccc}\mathrm{0},& & \end{array}& \end{array}& & \begin{array}{cc}& \end{array}\end{array}& & \varphi <\mathrm{0}\end{array},\end{array}$(33)

$P_w^{-}(\delta )=P_w^{+}(-\delta )$(34)

where $N_w^{-\mathrm{1}}=\mathrm{\Gamma }(\alpha /{\mu }_{\mathrm{4}}){\sqrt[]{{\mu }_{\mathrm{4}}}}^{-\mathrm{2}\alpha /{\mu }_{\mathrm{4}}}$, the Lyapunov function index for the invariant measure ${\delta }_{\varphi }$ can be obtained from Eq. (29) and Eq. (30) as follows:

${\lambda }_{\eta }({\delta }_{\varphi })={\mu }_{\mathrm{1}}-{\mu }_{\mathrm{4}}/\mathrm{2}$(35)

According to the ergodic theorem^[15,16], the Lyapunov exponent concerning the invariant measure $w_{\varphi }^{\pm }$ can be derived from Eq. (31), Eq. (33), and Eq. (35),

${\lambda }_{\eta }(w_{\varphi }^{\pm })=-\mathrm{2}\alpha =-\mathrm{2}\left({\mu }_{\mathrm{1}}-{\mu }_{\mathrm{4}}\right)/\mathrm{2}$(36)

After analyzing Eq. (35) and Eq. (36), it can be inferred that a qualitative change occurs in the dynamic behavior of the wheelset system when $\alpha =\mathrm{0}$. Consequently, the wheelset system occurs a stochastic D-bifurcation when $\alpha =\mathrm{0}$.

4.2 P-Bifurcation

Next, the combination of the stationary probability density method and numerical simulation is used to discuss the stochastic P-bifurcation of the wheelset system. FPK equation for the one-dimensional diffusion process of the wheelset system (3)

$\frac{\partial P}{\partial t}=-\frac{\partial }{\partial A}\left[({\mu }_{\mathrm{1}}A+{\mu }_{\mathrm{2}}{A}^{\mathrm{3}}+{\mu }_{\mathrm{3}}{A}^{\mathrm{5}})P(A)\right]+\frac{\mathrm{1}}{\mathrm{2}}\frac{{\partial }^{\mathrm{2}}}{\partial A^{\mathrm{2}}}[{\mu }_{\mathrm{4}}{A}^{\mathrm{2}}P(A)]$(37)

The initial value of the FPK equation satisfies the condition: $P(A,t|A_{\mathrm{0}},t_{\mathrm{0}})\to \delta (A-A_{\mathrm{0}}),t\to t_{\mathrm{0}}$, where $P(A,t|A_{\mathrm{0}},t_{\mathrm{0}})$ is the transition probability density in the diffusion process and the invariant measure of $A(t)$ is the stationary probability density $P(A)$. Let $\partial P/\partial t=\mathrm{0}$, we get

$-\frac{\partial }{\partial A}\left[({\mu }_{\mathrm{1}}A+{\mu }_{\mathrm{2}}{A}^{\mathrm{3}}+{\mu }_{\mathrm{3}}{A}^{\mathrm{5}})P(A)\right]+\frac{\mathrm{1}}{\mathrm{2}}\frac{{\partial }^{\mathrm{2}}}{\partial A^{\mathrm{2}}}[{\mu }_{\mathrm{4}}{A}^{\mathrm{2}}P(A)]=\mathrm{0}$(38)

Upon solving Eq. (38), the stationary probability density function is obtained:

$P(A)=\frac{C}{{\mu }_{\mathrm{4}}}{A}^{\frac{\mathrm{2}{\mu }_{\mathrm{1}}}{{\mu }_{\mathrm{4}}}-\mathrm{2}}\mathrm{e}\mathrm{x}\mathrm{p}(\frac{{\mu }_{\mathrm{2}}}{{\mu }_{\mathrm{4}}}{A}^{\mathrm{2}}+\frac{{\mu }_{\mathrm{3}}}{\mathrm{2}{\mu }_{\mathrm{4}}}{A}^{\mathrm{4}}),$(39)

where C is the normalization constant, such that ${\int }_{\mathrm{0}}^{+\infty }P(A)\mathrm{d}A=\mathrm{1}$. Let $A=\sqrt[]{u_{\mathrm{1}}^{\mathrm{2}}+u_{\mathrm{2}}^{\mathrm{2}}}$, then the joint probability density function of the system is obtained as follows

$P(u_{\mathrm{1}},u_{\mathrm{2}})=\frac{C}{u_{\mathrm{4}}}(u_{\mathrm{1}}^{\mathrm{2}}+u_{\mathrm{2}}^{\mathrm{2}}{)}^{\frac{{\mu }_{\mathrm{1}}}{{\mu }_{\mathrm{4}}}-\mathrm{1}}\mathrm{e}\mathrm{x}\mathrm{p}[\frac{{\mu }_{\mathrm{2}}}{{\mu }_{\mathrm{4}}}(u_{\mathrm{1}}^{\mathrm{2}}+u_{\mathrm{2}}^{\mathrm{2}})+\frac{{\mu }_{\mathrm{3}}}{\mathrm{2}{\mu }_{\mathrm{4}}}(u_{\mathrm{1}}^{\mathrm{2}}+u_{\mathrm{2}}^{\mathrm{2}}{)}^{\mathrm{2}}]$(40)