Response analysis of the vibro-impact system under fractional-order joint random excitation

Abstract

As a kind of good damping material, viscoelastic material is widely used in machinery, civil engineering, and other fields. In this paper, the viscoelasticity of the system is described by fractional differentiation. The dynamic response of a unilateral vibro-impact system with a viscoelastic oscillator under joint random excitation is studied, in which joint random excitation is composed of additive and multiplicative white noise. The fractional-order derivative was calculated based on Caputo’s definition, and the fractional derivative was equivalent to the corresponding linear damping force and linear restoring force. As a result, a new random system without fractional-order terms was obtained. A non-smooth transformation was introduced, which was equivalent to the original system to a new system without a velocity jump. The steady-state probability density functions of fractional-order vibro-impact systems under joint random excitation are solved by using the random average method and non-smooth transformation. In addition, the effects of parameters on the steady-state response of the system are analyzed.

Introduction

As a typical non-smooth dynamic model, a vibro-impact system is involved in several fields, such as civil engineering and machinery, and has been closely related to daily life. Random factors will always appear in practical engineering, and adding random factors on this basis will make its dynamic characteristics rich and complex. For instance, Feng et al.^1,2 used the quasi-conservative averaging method to study the steady-state response of a vibro-impact Duffing oscillator system under joint random excitation and single additive stochastic excitation, respectively. In the same study, the researchers also discussed the problem of stochastic bifurcation. Meanwhile, Li et al.^3,4 used the modified quasi-conservative averaging method to calculate the steady-state response function of the vibro-impact system under the joint-colored noise excitation. They compared it with the numerical solution, and the stochastic bifurcations were also explored. Xu et al.^5–7 adopted the modified Hertzian contact model. The nonlinear factors were equated to linear factors, and the stochastic response of an inelastic vibro-impact system was analyzed by the energy dissipation balance method and the stochastic averaging method. Then, the joint probability density function (PDF) of displacement velocity was obtained. Recently, Xu et al. ⁸ analyzed the steady-state response of an inelastic impact system under additive and multiplicative joint random excitation. The results showed that the joint random excitation has less damage to the stability of the system than the single random excitation. Huang et al. ⁹ expressed the multi-degree-of-freedom vibro-impact system as a dissipative Hamiltonian system, and the steady-state response of the system was studied by using the quasi-Hamiltonian stochastic averaging method. Furthermore, Zhuravlev and Iourtchenko et al.^10–16 conducted numerical and analytical calculations on the dynamic behavior of a classical vibro-impact system under stochastic excitation. The steady-state response of the system was studied by using the energy envelope stochastic averaging method. In addition, the energy dissipation balance method was used to study the mean value of the system, and the energy dissipation of the system was analyzed. Meanwhile, the path integral method was used to study the high energy dissipation of the system. Yang et al. ¹⁷ studied the stochastic response of the Rayleigh-Van der Pol vibro-impact system under joint parametric excitation. They analytically obtained the steady-state PDF of the system by using the non-smooth transformation and stochastic averaging method, based on which the stochastic bifurcation of the system was analyzed.

As a kind of good damping material, viscoelastic materials can achieve good vibration and noise reduction in a wide frequency band and are widely used in machinery, civil engineering, and other fields, such as magnetorheological dampers, metal rubber, and air springs. When the fractional-order is set to 0∼1, the fractional differential behaves as viscoelastic. Therefore, fractional differentiation can be used to describe the viscoelasticity of the system, to further study the dynamic response of vibro-impact systems with viscoelastic materials under joint random excitation. Some scholars extend fractional calculus to the field of random non-smoothing and study its rich dynamic characteristics. Zhao et al. ¹⁸ studied the stochastic response of a vibro-impact system under fractional additive and multiplicative stochastic excitation and found that fractional coefficient and fractional-order would induce the stochastic bifurcation of the system. Chen et al. ¹⁹ studied the stochastic response of the Duffing oscillator under the combined excitation of fractional harmonic and white noise. They used the reduced Fokker-Plank-Kolmogorov (FPK) equation to predict the steady-state response of the original system. According to the FPK equation, a stochastic jump would appear in the system, and a bifurcation phenomenon would occur when the fractional parameters were changed. Xiao et al. ²⁰ applied the generalized stochastic averaging method to study the strongly nonlinear vibro-impact system under fractional random excitation. Their study showed that the fractional derivative term, impact conditions, and stochastic excitation would affect the response of the fractional vibro-impact dynamic system. Meanwhile, Yang et al. ²¹ applied non-smooth transformation and the stochastic averaging method to study the stochastic bifurcation of the vibro-impact system under fractional random excitation. The results showed that the change of fractional-order would change the number of peaks in the image of the steady-state PDF of the system. Furthermore, Li et al. ²² studied the system with a tri-stable van der Pol oscillator under fractional derivative excitation. Song et al. ²³ investigated the bifurcation problem of a system with a Van der Pol oscillator under stochastic excitation through experiments. The study showed that with the change of random excitation intensity, the joint probability density of the system changed the topology shape. Yang et al. ²⁴ studied the stochastic response of bistable systems under fractional stochastic excitation. They found that the fractional-order value and the noise intensity would cause stochastic P-bifurcation. The fractional-order might cause the steady-state PDF to turn from a single-peak mode to a double-peak mode, but the noise intensity might change the steady-state PDF from a double-peak mode to a single-peak mode. Moreover, Sun et al. ²⁵ studied the stochastic P-bifurcation problem of a fractional damping system under joint excitation. Li et al. ²⁶ studied the stochastic bifurcation problem of the fractional derivative Duffing-van der Pol system under color noise excitation. They obtained the critical conditions for the stochastic bifurcation of the system by using singularity theory. They also found that the change in noise intensity and fractional derivative would induce the stochastic bifurcation of the system. Tao et al. ²⁷ solved traditional differential equations and fractional differential equations by Aboodh transform-based homotopy perturbation method. They found that this approach has been shown to have the potential to solve both linear and nonlinear problems. Tao et al. ²⁸ solved traditional differential equations and fractional differential equations by Aboodh transformation-based homotopy perturbation method. They found that this approach has been shown to have the potential to solve both linear and nonlinear problems. Ain et al. ²⁹ analyzed the effect of time MERS-CoV disease transmission explored using a nonlinear fractional order in the sense of the Caputo operator. Anjum et al. ³⁰ couple fractional complex transformation with HPM to obtain a fairly accurate analytic solution of the nonlinear time fractional CH equation. They found the high accuracy and efficiency of this proposed coupling for the solution. Zhao et al.^31,32 by comparing the performance of the single-layer linear and that of the double-layer one. They found that the double-layer liner can increase the damping effect at a higher frequency range. In the follow-up study, they found that porous concrete samples are found be more effective in absorbing noise by comparing the noise damping performances of different materials. Zhao et al. ³³ established a time-domain numerical model of cylindrical lined duct. They verified that the model can be used to simulate the propagation and dissipation of acoustic plane waves in a lined duct in real-time.

Akbar et al. ³⁴ used the sine-Gordon expansion approach and generalized Kudryashov scheme for the Boussinesq equation, and the dynamical behavior of the waves is analyzed. Ahmad et al.^35,36 propose a new method for finding the numerical solutions of nonlinear noninteger order partial differential equations, which is called the fractional iteration algorithm-I, the research shows that this method can ensure a rapid convergence. In subsequent studies,^37–40 they proposed Variational iteration algorithm-I and Variational iteration algorithm-II, the research shows that this method can reduce the amount of computation, and has the advantages of fast convergence, high computational efficiency, accurate results, and better robustness. In addition, there are also generalized auxiliary equation technique, ⁴¹ Riccatti transformation method, ⁴² meshless techniques,^43,44 modified (G’/G)-expansion method, ⁴⁵ Lattice Boltzmann method,^46,47 and other methods that can be applied to the nonlinear systems analyzed.

Although many achievements have been made in existing research of fractional-order vibro-impact systems,^48–50 there is a lack of steady-state response research under joint excitation. Most of these results have studied the fractional-order vibro-impact system under a single excitation, while the system is often affected by more than one kind of excitation in real works. It is often not enough to consider only one kind of excitation when studying the dynamics behavior of the system. Therefore, considering the complexity of the excitation of the actual system, the steady-state response of the fractional-order vibro-impact system under joint random excitation is studied, in which the joint excitation is composed of additive and multiplicative white noise.

The following is the format of the rest of the article: In this paper, the dynamic response of a unilateral vibro-impact system with a viscoelastic oscillator under joint random excitation is studied, in which joint random excitation is composed of additive and multiplicative white noise. The steady-state probability density functions of fractional-order vibro-impact systems under joint random excitation are solved by using the random average method and non-smooth transformation, and the effects of different parameters on the steady-state response of the system are analyzed.

Vibro-impact system

In practical scenarios, systems can be affected by multiple types of excitations, and it is inadequate to consider only one kind of excitation when studying the system. Therefore, considering the complexity of the excitation of the actual system, the steady-state response of the fractional-order collision system under the joint random excitation is studied, in which the joint random excitation is composed of additive and multiplicative random excitation. The vibro-impact system model is illustrated in Figure 1. Where $F_1$ and $F_2$ are the additive and multiplicative white noise, respectively. $m$ is the mass of the colliding oscillator; $k$ is the spring stiffness coefficient; $c$ is the linear damping coefficient. In addition, $p$ and $a$ are the fractional damping coefficient and fractional-order coefficient of the impact oscillator, respectively. A single rigid barrier is installed at the static equilibrium position of the system ( $h=0$ ).

Figure 1.

Vibro-impact System model.

The equation of the vibro-impact system under fractional-order joint random excitation can be expressed as:

$\begin{array}{l}m\ddot{x}+kx+c\dot{x}+p{D}^{\alpha }\left[x\right]=f_1{\zeta }_1\left(t\right)+f_{2}x{\zeta }_2\left(t\right),x>0\\ {\dot{x}}^{+}=-e{\dot{x}}^{-},x=0,\end{array}$

(1)

where ${\dot{x}}^{-}$ and ${\dot{x}}^{+}$ are impact and rebound velocities of the mass, respectively. $e$ is the collision recovery coefficient, and $0<e\le 1$ . The value of $e$ reflects the degree of energy loss of the oscillator. The smaller the value is, the more energy is lost. $f_1$ and $f_2$ are the amplitudes of two white Gaussian noises, respectively. ${\zeta }_1\left(t\right)$ and ${\zeta }_2\left(t\right)$ are two independent Gaussian white noises, whose expressions are

$\begin{array}{ll}E\left[{\zeta }_1\left(t\right)\right]=0& E\left[{\zeta }_1\left(t+\tau \right){\zeta }_1\left(t\right)\right]=2D_1\mathrm{\delta }\left(\tau \right)\\ E\left[{\zeta }_2\left(t\right)\right]=0& E\left[{\zeta }_2\left(t+\tau \right){\zeta }_2\left(t\right)\right]=2D_2\mathrm{\delta }\left(\tau \right)\end{array}.$

System (1) may be rescaled as follows:

${\omega }_0=\sqrt{\frac{k}{m}},2{\mu }_1=\frac{c}{m},K=\frac{p}{m},F_1=\frac{f_1}{m},F_2=\frac{f_2}{m}.$

The following equation is derived:

$\begin{array}{l}\ddot{x}+{\omega }_0^{2}x+2{\mu }_1\dot{x}+K{D}^{\alpha }\left[x\right]=F_1{\zeta }_1\left(t\right)+F_{2}x{\zeta }_2\left(t\right),x>0\\ {\dot{x}}^{+}=-e{\dot{x}}^{-},x=0,\end{array}$

(2)

There are several definitions for fractional-order derivatives, such as Grunwald-Letnikov’s definition, Riemann-Liouville’s definition, and Caputo’s definition, and they are equivalent under some conditions for a wide class of functions. Here we adopt Caputo’s definition. $D^{\alpha }\left[x\right]$ represents the fractional derivative in the Caputo sense, which can be expressed in the following form:

$D^{\alpha }\left[x\right]=\frac{1}{\mathrm{\Gamma }\left(1-\alpha \right)}{\int }_0^t\frac{\dot{x}\left(u\right)}{{\left(t-u\right)}^{\alpha }}du,\left(0<\alpha \le 1\right).$

(3)

where $\mathrm{\Gamma }\left(z\right)$ is the Gamma function and satisfies $\mathrm{\Gamma }\left(z+1\right)=z\mathrm{\Gamma }\left(z\right)$ . According to the conclusion obtained in Refs. [46,47], the fractional derivative defined by Caputo can be replaced by classical damping force and linear restoring force:

$D^{\alpha }\left[x\right]={\omega }_0^{\alpha -1}\hspace{0.17em}\sin \frac{\alpha \pi }{2}\dot{x}+{\omega }_0^{\alpha }\hspace{0.17em}\cos \frac{\alpha \pi }{2}x.$

(4)

Substituting equation (4) into equation (2), the original differential equation of the system can be expressed as

$\begin{array}{l}\ddot{x}+{\omega }_0^{2}x+2{\mu }_1\dot{x}+\\ K\left({\omega }_0^{\alpha -1}\hspace{0.17em}\sin \frac{\alpha \pi }{2}\dot{x}+{\omega }_0^{\alpha }\hspace{0.17em}\cos \frac{\alpha \pi }{2}x\right)=F_1{\zeta }_1\left(t\right)+F_{2}x{\zeta }_2\left(t\right),x>0\\ {\dot{x}}^{+}=-e{\dot{x}}^{-},x=0\end{array}.$

(5)

Non-smooth transformation

According to Refs. [51,52], the system is difficult to examine directly due to the non-smooth characteristic. Therefore, the non-smooth transformation is performed on the system first.

$x=\left|y\right|,\dot{x}=\dot{y}\text{sgn}\left(y\right),\ddot{x}=\ddot{y}\text{sgn}\left(y\right)$

$\text{sgn}\left(y\right)=\{\begin{array}{ll}1& y>0\\ 0& y=0\\ -1& y<0.\end{array}$

(6)

Substituting equation (6) in equation (5) and using $\left|y\right|=y\text{sgn}\left(y\right)$ , the following equation is derived:

$\begin{array}{l}\ddot{y}+{\omega }_0^{2}y+2{\mu }_1\dot{y}+K\left({\omega }_0^{\alpha -1}\hspace{0.17em}\sin \frac{\alpha \pi }{2}\dot{y}+{\omega }_0^{\alpha }\hspace{0.17em}\cos \frac{\alpha \pi }{2}y\right)=\\ F_1{\zeta }_1\left(t\right)\text{sgn}\left(y\right)+F_{2}y{\zeta }_2\left(t\right),t\ne t_k\\ {\dot{y}}^{+}=e{\dot{y}}^{-},t=t_k.\end{array}$

(7)

where $t_k$ is the moment of collision in equation (7). $x=0$ is equivalent to $t=t_k$ . Before and after the collision, the transition quantity of velocity in the original state variable is $\mathrm{\Delta }\dot{x}={\dot{x}}^{+}-{\dot{x}}^{-}=\left(-e-1\right){\dot{x}}^{-}$ , while the transition quantity of velocity in the new state variable is $\mathrm{\Delta }\dot{y}={\dot{y}}^{+}-{\dot{y}}^{-}=\left(e-1\right){\dot{y}}^{-}$ . With reference to the Dirac delta function, the time of collision can be regarded as an impulse term applied to the system at time $t_k$ . By Dirac delta function property:

$\mathrm{\delta }\left(y\right)\left|\dot{y}\right|=\mathrm{\delta }\left(t-t_k\right).$

Using the continuity at the time of collision and combining with the Dirac delta function, the impulse term is shown by equation (8):

$\left({\dot{y}}^{-}-{\dot{y}}^{+}\right)\mathrm{\delta }\left(t-t_k\right)=\left(1-e\right)\dot{y}\left|\dot{y}\right|\mathrm{\delta }\left(y\right).$

(8)

Substituting equation (8) in equation (7), the following equation is derived:

$\ddot{y}+{\omega }_1^{2}y+b\dot{y}+\left(1-e\right)\dot{y}\left|\dot{y}\right|\mathrm{\delta }\left(y\right)=F_1{\zeta }_1\left(t\right)\text{sgn}\left(y\right)+F_{2}y{\zeta }_2\left(t\right),$

(9)

Where

${\omega }_1^2={\omega }_0^2+K{\omega }_0^{\alpha }\hspace{0.17em}\cos \left(\frac{\alpha \pi }{2}\right),b=2{\mu }_1+K{\omega }_0^{\alpha -1}\hspace{0.17em}\sin \left(\frac{\alpha \pi }{2}\right).$

Stochastic averaging method

According to Ref. [53], the solution of equation (9) can be expressed as equation (10).

$\begin{array}{l}y\left(t\right)=A\hspace{0.17em}\cos \hspace{0.17em}\varphi \\ \dot{y}\left(t\right)=-A{\omega }_1\hspace{0.17em}\sin \hspace{0.17em}\varphi \\ \varphi ={\omega }_{1}t+\theta \\ U\left(A\right)={\int }_0^A{\omega }_1^{2}udu=\frac{1}{2}{\omega }_1^2{A}^2\end{array}.$

(10)

Take the derivative of time $t$ in equation (10) can be obtained:

$\{\begin{array}{l}\dot{y}\left(t\right)=\dot{A}\cos \hspace{0.17em}\varphi -A\left({\omega }_1+\dot{\theta }\right)\sin \hspace{0.17em}\varphi \\ \ddot{y}\left(t\right)=-\dot{A}{\omega }_1\hspace{0.17em}\sin \hspace{0.17em}\varphi -A{\omega }_1\left({\omega }_1+\dot{\theta }\right)\cos \hspace{0.17em}\varphi .\end{array}$

(11)

Substituting equation (11) in equation (9), the following equation is derived:

$\{\begin{array}{l}\dot{A}\cos \hspace{0.17em}\varphi -A\dot{\theta }\sin \hspace{0.17em}\varphi =0\\ -\dot{A}{\omega }_1\hspace{0.17em}\sin \hspace{0.17em}\varphi -A{\omega }_1\left({\omega }_1+\dot{\theta }\right)\cos \hspace{0.17em}\varphi +{\omega }_1^{2}A\hspace{0.17em}\cos \hspace{0.17em}\varphi +b\left(-A{\omega }_1\hspace{0.17em}\sin \hspace{0.17em}\varphi \right)\\ +\left(1-e\right)\left(-A{\omega }_1\hspace{0.17em}\sin \hspace{0.17em}\varphi \right)\left|-A{\omega }_1\hspace{0.17em}\sin \hspace{0.17em}\varphi \right|\delta \left(A\hspace{0.17em}\cos \hspace{0.17em}\varphi \right)=\\ F_1{\zeta }_1\left(t\right)\text{sgn}\left(A\hspace{0.17em}\cos \hspace{0.17em}\varphi \right)+F_2{\zeta }_2\left(t\right)A\hspace{0.17em}\cos \hspace{0.17em}\varphi ,\end{array}$

(12)

$\dot{A}$ and $\dot{\theta }$ in equation (12) are taken as unknowns to solve the bivariate primary equation:

$\{\begin{array}{l}\frac{\mathrm{d}A}{\mathrm{d}t}=M_1+G_{11}{\zeta }_1\left(t\right)+G_{12}{\zeta }_2\left(t\right)\\ \frac{\mathrm{d}\theta }{\mathrm{d}t}=M_2+G_{21}{\zeta }_1\left(t\right)+G_{22}{\zeta }_2\left(t\right),\end{array}$

(13)

Where

$\begin{array}{l}{M}_1=-\frac{\sin \hspace{0.17em}\varphi }{{\omega }_1}\left\{\begin{array}{l}b\left(A{\omega }_1\hspace{0.17em}\sin \hspace{0.17em}\varphi \right)\\ -\left(1-e\right)\left|-A{\omega }_1\hspace{0.17em}\sin \hspace{0.17em}\varphi \right|\left(-A{\omega }_1\hspace{0.17em}\sin \hspace{0.17em}\varphi \right)\mathrm{\delta }\left(A\hspace{0.17em}\cos \hspace{0.17em}\varphi \right)\end{array}\right\}\\ G_{11}=-\frac{\sin \hspace{0.17em}\varphi }{{\omega }_1}{F}_1\text{sgn}\left(A\hspace{0.17em}\cos \hspace{0.17em}\varphi \right)\\ G_{12}=-\frac{\sin \hspace{0.17em}\varphi }{{\omega }_1}{F}_{2}A\hspace{0.17em}\cos \hspace{0.17em}\varphi ,\end{array}$

(14)

$\begin{array}{l}{M}_2=-\frac{\cos \hspace{0.17em}\varphi }{A{\omega }_1}\left\{\begin{array}{l}b\left(A{\omega }_1\hspace{0.17em}\sin \hspace{0.17em}\varphi \right)\\ -\left(1-e\right)\left|-A{\omega }_1\hspace{0.17em}\sin \hspace{0.17em}\varphi \right|\left(-A{\omega }_1\hspace{0.17em}\sin \hspace{0.17em}\varphi \right)\mathrm{\delta }\left(A\hspace{0.17em}\cos \hspace{0.17em}\varphi \right)\end{array}\right\}\\ G_{21}=-\frac{\cos \hspace{0.17em}\varphi }{A{\omega }_1}{F}_1\text{sgn}\left(A\hspace{0.17em}\cos \hspace{0.17em}\varphi \right)\\ G_{22}=-\frac{\cos \hspace{0.17em}\varphi }{A{\omega }_1}{F}_{2}A\hspace{0.17em}\cos \hspace{0.17em}\varphi .\end{array}$

(15)

According to the stochastic averaging method and the Itô law of differentiation, the amplitude is a slow variable, and the phase is a fast variable. Averaging the phases over time gives a deterministic stochastic mean value of the amplitude and phase difference. The Itô differential equation can be obtained as follows:

$\begin{array}{l}\mathrm{d}A=a_1\left(A,\varphi ,{\omega }_{0}t\right)\mathrm{d}t+{\sigma }_{1r}\left(A,\varphi \right)\mathrm{d}{B}_r\left(t\right)\\ \mathrm{d}\theta =a_2\left(A,\varphi ,{\omega }_{0}t\right)\mathrm{d}t+{\sigma }_{2r}\left(A,\varphi \right)\mathrm{d}{B}_r\left(t\right),r=1,2,\end{array}$

(16)

where ${\alpha }_1$ and ${\alpha }_2$ are the drift coefficient. ${\sigma }_{11}$ , ${\sigma }_{12}$ , ${\sigma }_{21}$ , and ${\sigma }_{22}$ are the diffusion coefficient. $B\left(t\right)$ is a unit Wiener process.

$\begin{array}{l}{a}_i=M_i+D_l\left(G_{1l}\frac{\partial G_{il}}{\partial A}+G_{2l}\frac{\partial G_{il}}{\partial \varphi }\right)\\ b_{ij}={\displaystyle \sum }_{r,s=1}^2{\sigma }_{ir}{\sigma }_{js}=2D_l{G}_{il}{G}_{jl},i,j,l=1,2.\end{array}$

(17)

By averaging $\varphi$ in equation (17) within a period, the mean drift coefficient and diffusion coefficient can be expressed as follows:

$\{\begin{array}{l}{\overline{a}}_1\left(A\right)=\frac{1}{T}{\int }_0^T{M}_1+D_1\left(G_{11}\frac{\partial G_{11}}{\partial A}+G_{21}\frac{\partial G_{11}}{\partial \varphi }\right)+D_2\left(G_{12}\frac{\partial G_{12}}{\partial A}+G_{22}\frac{\partial G_{12}}{\partial \varphi }\right)\mathrm{d}\varphi \\ {\overline{b}}_{11}=\frac{1}{T}{\int }_0^{T}2D_1{G}_{11}{G}_{11}+2D_2{G}_{12}{G}_{12}\mathrm{d}\varphi \end{array}.$

(18)

The termination time of integration is set as $T=2\pi$ , which means the integral interval is $\left(0,2\pi \right)$ . By integrating equation (18), we can obtain

$\{\begin{array}{l}{\overline{a}}_1\left(A\right)=-\frac{bA}{2}-\frac{\left(1-e\right){\omega }_{1}A}{\pi }+\frac{D_1{F_1}^2}{2A{\omega }_1^2}+\frac{3D_2{F_2}^{2}A}{8{\omega }_1^2}\\ {\overline{b}}_{11}=\frac{{\mathrm{D}}_1{F_1}^2}{{\omega }_1^2}+\frac{D_2{F_2}^2{A}^2}{4{\omega }_1^2}\end{array}.$

(19)

Therefore,

$\begin{array}{l}\mathrm{d}A=\left(-\frac{bA}{2}-\frac{\left(1-e\right){\omega }_{1}A}{\pi }+\frac{D_1{F_1}^2}{2A{\omega }_1^2}+\frac{3D_2{F_2}^{2}A}{8{\omega }_1^2}\right)\mathrm{d}t\\ +\sqrt{\frac{D_1{F_1}^2}{{\omega }_1^2}+\frac{D_2{F_2}^2{A}^2}{4{\omega }_1^2}}\mathrm{d}B\left(t\right)\end{array}.$

(20)

where $B\left(\mathrm{t}\right)$ is an independent unit Wiener process. Equation (20) shows that amplitude A does not depend on the change of $\varphi$ . The FPK equation about amplitude A can be obtained from equation (20) as follows:

$\begin{array}{l}\frac{\partial p\left(A,t\right)}{\partial t}=-\frac{\partial }{\partial A}\left[\left(-\frac{bA}{2}-\frac{\left(1-e\right){\omega }_{1}A}{\pi }+\frac{D_1{F_1}^2}{2A{\omega }_1^2}+\frac{3D_2{F_2}^{2}A}{8{\omega }_1^2}\right)p\right]\\ +\frac{1}{2}\frac{{\partial }^2}{\partial A^2}\left[\left(\frac{D_1{F_1}^2}{{\omega }_1^2}+\frac{D_2{F_2}^2{A}^2}{4{\omega }_1^2}\right)p\right].\end{array}$

(21)

When $A=0,p=c,c\in \left(-\infty ,+\infty \right)$ , when $A\to \infty ,p\to 0,\frac{\partial p}{\partial A}\to 0$ , the expression of steady-state probability density in equation (21) can be obtained:

$p\left(A\right)=\frac{C}{{\overline{b}}_{11}\left(A\right)}\exp \left[\int \frac{2{\overline{a}}_1\left(A\right)}{{\overline{b}}_{11}\left(A\right)}\mathrm{d}A\right],$

(22)

where $C$ is the normalization constant. The expression is

$C={\left[{\int }_0^{\infty }\left(\frac{1}{{\overline{b}}_{11}\left(A\right)}\exp \left[\int \frac{2{\overline{a}}_1\left(A\right)}{{\overline{b}}_{11}\left(A\right)}\mathrm{d}A\right]\right)\mathrm{d}A\right]}^{-}.$

Substituting equation (19) in equation (22), a complete expression of steady-state probability density can be obtained:

$p\left(A\right)=\frac{CA\mathrm{M}{\left(A\right)}^{\mathrm{N}\left(A\right)}}{{\overline{b}}_{11}\left(A\right)},$

(23)

where

$\{\begin{array}{l}\mathrm{M}\left(A\right)=4D_1{F}_1^2+D_2{F}_2^2{A}^2\\ \mathrm{N}\left(A\right)=1-\frac{2{\omega }_1^2\left(b\pi -2{\omega }_1\left(e-1\right)\right)}{D_2{F}_2^2\pi }\end{array}.$

The steady-state PDF of the system Hamiltonian $H=U\left(A\right)$ can be obtained from equation (22).

$p_{Y,\dot{Y}}\left(\mathrm{y},\dot{\mathrm{y}}\right)={\frac{p\left(A\right){\omega }_1}{2\pi {\omega }_0^{2}A}|}_{A=U^{-1}\left(H\right)}.$

(24)

In equation (24), $U^{-1}$ represent the inverse function of $U$ . According to the fourth equation of equation (10), we can obtain that $U^{-1}\left(H\right)=\sqrt{2H/{\omega }_1^2}$ . Substitute $H=y^2/2+U\left(y\right)$ and $U\left(y\right)={\omega }_1^2{y}^2/2$ into $U^{-1}\left(H\right)$ the following equation is derived

$A=U^{-1}\left(H\right)=\sqrt{\frac{y^2}{{\omega }_1^2}+y^2}.$

According to the non-smooth transformation, the PDF of displacement and velocity of the system are as follows:

$p\left(x\right)={\int }_{-\infty }^{\infty }{p}_{X,\dot{X}}\left(x,\dot{x}\right)\mathrm{d}\dot{x}=2{\int }_{-\infty }^{\infty }{p}_{\mathrm{Y},\dot{\mathrm{Y}}}\left(x,\dot{x}\right)\mathrm{d}\dot{x},$

(25)

$p\left(\dot{x}\right)={\int }_0^{\infty }{p}_{X,\dot{X}}\left(x,\dot{x}\right)\mathrm{d}x=2{\int }_0^{\infty }{p}_{\mathrm{Y},\dot{\mathrm{Y}}}\left(x,\dot{x}\right)\mathrm{d}x.$

(26)

Numerical modeling

The numerical scheme proposed above is used to simulate the fractional derivative. Gaussian white noise is generated by random numbers. Finally, the Runge-Kutta method is used to simulate the original system (1) to obtain numerical results, which can prove the correctness of the analytical results. The simulation results are shown in Figure 2, where Figure 2(a)–(c) presents the steady-state PDF of amplitude, displacement, and velocity, respectively. The figure also shows that the numerical and analytical results agree well with each other within the allowed range of errors. Thus, the effectiveness of the stochastic averaging method is verified.

Figure 2.

The fitting diagram of the analytical and numerical solution of the system ( $\alpha =0.5,{\omega }_0=1,e=0.9,K=1,F_1=1,F_2=1,{\mu }_1=0.3,D_1=0.5$ and ${\mathrm{D}}_2=0.5$ ). (a) PDF of amplitude A; (b) PDF of displacement $x_1$ ; (c) PDF of velocity $x_2$ .

Parametric analysis

In this section, the effects of different parameters on the steady-state response of the system are analyzed. In the following analysis, the meanings of system parameters are shown in Table 1.

Table 1.

The parameters of the system.

Item	Notation	Value
Natural frequency	${\omega }_0$	1
Fractional-order	$a$	0.5
Fractional-order coefficient	$K$	1
Restoring force	$e$	0.9
Linear damping force	$u_1$	0.3
Additive noise amplitude	$F_1$	1
Multiplicative noise amplitude	$F_2$	1
Additive noise intensity	$D_1$	0.5
Multiplicative noise intensity	$D_2$	0.5

Influence of fractional-order terms on vibration characteristics of the system

The effect of fractional-order on the system

When other parameters remain unchanged and fractional-order $\alpha$ is set at 0.1, 0.3, and 0.5, the steady-state probability density curves of the amplitude, displacement, and velocity of the system are shown in Figure 3(a)–(c), respectively.

Figure 3.

Probability density function of the system ( ${\omega }_0=1,e=0.9,K=1,{\mu }_1=0.3,F_1=1,F_2=1,{\mathrm{D}}_1=0.5$ and ${\mathrm{D}}_2=0.5$ ). (a) PDF of amplitude A; (b) PDF of displacement $x_1$ ; (c) PDF of velocity $x_2$ .

Figure 3 demonstrates that with the gradual decrease of fractional-order, the peak value of the steady-state PDF curve of the system’s amplitude, displacement, and velocity also decreases, which indicates that the system will gradually deviate from the stable state as the fractional-order decreases. In addition, when the fractional-order decreases from 0.5 to 0.3, the peak value of the PDF curves of the amplitude, displacement, and velocity of the system decreases slightly. When the fractional-order decreases from 0.3 to 0.1, the peak value of the PDF curves of the amplitude, displacement, and velocity decreases greatly. This correlation indicates that the closer the fractional-order is to 0, the more sensitive the system is to its change. In addition, Figure 3(a) shows that with the increase of fractional-order $\alpha$ , the peaks of these probability density function curves gradually shift to the direction with a smaller amplitude. The increase of fractional-order $\alpha$ leads to the increase of fractional-order equivalent linear damping and the decrease of equivalent linear stiffness. This phenomenon indicates that a larger fractional-order will lead to a weaker response of the system.

The time displacement diagram and time velocity of the system are shown in Figure 4 and Figure 5, respectively. Comparing Figures 4(a) and Figure5(a) we can know, when fractional-order $\alpha =0.5$ , the displacement response of the system is mostly below 2. Meanwhile, when fractional-order $\alpha =0.1$ , the displacement is more than 2 in many cases. Similarly, comparing Figures 4(b) and Figure 5(b) shows that the velocity amplitude of fractional-order $\alpha =0.1$ is larger than that of fractional-order $\alpha =0.5$ . In addition, Figure 3 demonstrates that with the increase of fractional-order, not only does the peak value increase but the probability of larger amplitude, displacement, and velocity value decreases as well.

Figure 4.

Time displacement diagram and time velocity diagram ( $a=0.1,{\omega }_0=1,e=0.9,K=1,{\mu }_1=0.3,F_1=1,F_2=1,{\mathrm{D}}_1=0.5$ and ${\mathrm{D}}_2=0.5$ ). (a) Time displacement diagram; (b) Time velocity diagram.

Figure 5.

Time displacement diagram and time velocity diagram ( $a=0.5,{\omega }_0=1,e=0.9,K=1,{\mu }_1=0.3,F_1=1,F_2=1,{\mathrm{D}}_1=0.5$ and ${\mathrm{D}}_2=0.5$ ). (a) Time displacement diagram; (b) Time velocity diagram.

The effect of the fractional coefficient on the system

With other parameters unchanged, the fractional-order is fixed $\alpha =0.5$ . Figure 6 shows three probability density graphs of the system with different fractional-order coefficients K. Furthermore, equation (5) demonstrates that the equivalent linear stiffness and equivalent linear damping of the system increase with the increase in fractional-order coefficient K. Comparing Figure 6 with Figure 3 shows that the fractional-order coefficient has a greater effect on the system variation amplitude than the numeral-order coefficient. However, when the fractional-order coefficient $K$ is from 1 to 2 or 2 to 3, the variation range of the three PDF graphs is almost the same. When the fractional-order coefficient $K$ decreases, the peak value of the PDF curve of the amplitude, displacement, and velocity of the system also decreases, thus making the response of the system relatively strong. Furthermore, the probability of a strong response from the system increases. Comparing Figure 5 with Figure 7 shows that when $K=3$ , the maximum value of both displacement and velocity is much smaller than the value when $K=1$ . Moreover, the maximum value of amplitude and displacement can be regarded as less than 1. Those larger than 1 can be ignored, and the maximum value of velocity can be regarded as between −2 and 2. Meanwhile, the rest can be ignored. Thus, the steady-state response interval of the system can be estimated relatively accurately.

Figure 6.

Probability density function of the system ( $a=0.5,{\omega }_0=1,e=0.9,{\mu }_1=0.3,F_1=1,F_2=1,{\mathrm{D}}_1=0.5$ and ${\mathrm{D}}_2=0.5$ ). (a) PDF of amplitude A; (b) PDF of displacement $x_1$ ; (c) PDF of velocity $x_2.$

Figure 7.

Time displacement diagram and time velocity diagram ( $K=3,a=0.5,{\omega }_0=1,e=0.9,{\mu }_1=0.3,F_1=1,F_2=1,{\mathrm{D}}_1=0.5$ and ${\mathrm{D}}_2=0.5$ ). (a) Time displacement diagram; (b) Time velocity diagram.

Effect of restoring force e on vibration characteristics of system

Selecting an appropriate collision restitution coefficient in practical engineering can be difficult. Hence, the influence of the collision restitution coefficient on the system is analyzed in this section. With the other parameters unchanged, restoring force e is set at 0.3, 0.6, and 0.9, respectively. The steady-state probability density diagram of the amplitude, displacement, and velocity of the system is shown in Figure 8. When the restoring force $e$ gradually increases, the peak value of the PDF curve of the system amplitude, displacement, and velocity decreases. However, the probability of the larger values of the amplitude, displacement, and velocity increases. This phenomenon indicates that within a certain range, the greater the restoring force e is, the stronger the vibration response of the system will be. The time displacement diagram in Figure 9(a) demonstrates that when $e=0.3$ , most of the system displacements are less than 1.5. Meanwhile, the time displacement diagram in Figure 5(a) shows that when $e=0.9$ , the displacement has a value greater than 6, and relatively many cases have a displacement greater than 4. Similarly, the change in the recovery coefficient has the same effect on the system velocity. In theory, when $e=1$ , the system does not consume energy when it collides. When $e=0$ , it completely consumes energy. Therefore, the conclusion observed from the figure follows the collision recovery theory, which also verifies the correctness of the analytical method used in this study. The choice of the size of the restitution coefficient should not be limited to it being larger or smaller. A specific system should have a specific restitution coefficient to satisfy the use requirement according to the requirements in the practical application of the reasonable selection of the crash restitution coefficient.

Figure 8.

Probability density function of the system ( $a=0.5,{\omega }_0=1,K=1,{\mu }_1=0.3,F_1=1,F_2=1,{\mathrm{D}}_1=0.5$ and ${\mathrm{D}}_2=0.5$ ). (a) PDF of amplitude A; (b) PDF of displacement $x_1$ ; (c) PDF of velocity $x_2$ .

Figure 9.

Time displacement diagram and time velocity diagram ( $a=0.5,{\omega }_0=1,K=1,e=0.3,{\mu }_1=0.3,F_1=1,F_2=1,{\mathrm{D}}_1=0.5$ and ${\mathrm{D}}_2=0.5$ ). (a) Time displacement diagram; (b) Time velocity diagram.

Effect of linear damping force on the system

As in the previous sections, other parameters remain unchanged. The linear damping coefficient ${\mu }_1$ is set as 0.3, 0.5, and 0.7, respectively, and the steady-state probability density diagram of the amplitude, displacement, and velocity of the system is shown in Figure 10. The increase in the damping coefficient ${\mu }_1$ can increase the peaks of the three curves. Figures 5(a) and Figure 11(a) show that when ${\mu }_1=0.3$ and ${\mu }_1=0.7$ , the change range of the displacement of the system decreases when the damping of the system increases, thus indicating that the increase of damping inhibits the motion of the system. Meanwhile, Figures 5(b) and Figure 11(b) demonstrate that when ${\mu }_1=0.3$ and ${\mu }_1=0.7$ , with the increase of damping, the velocity of the system not only decreases in amplitude but also increases in frequency. This outcome indicates that the larger the damping is, the smaller the amplitude and the faster the vibration velocity will be. The above conclusion is in line with the characteristics of general vibration and proves the correctness of the analytical solution.

Figure 10.

Probability density function of the system ( $a=0.5,{\omega }_0=1,K=1,e=0.9,{\mu }_1=0.3,F_1=1,F_2=1,{\mathrm{D}}_1=0.5$ and ${\mathrm{D}}_2=0.5$ ). (a) PDF of amplitude A; (b) PDF of displacement $x_1;$ (c) PDF of velocity $x_2$ .

Figure 11.

Time displacement diagram and time velocity diagram ( $a=0.5,{\omega }_0=1,K=1,e=0.3,{\mu }_1=0.7,F_1=1,F_2=1,{\mathrm{D}}_1=0.5$ and ${\mathrm{D}}_2=0.5$ ). (a) Time displacement diagram; (b) Time velocity diagram.

Effect of excitation on the system

Based on the impact of noise intensity on the system, Figure 12 presents that the peaks of all three curves tend to decrease with the increase in noise intensity. Combined with Figure 5 and 12, and Figure 13, the peaks of the amplitude, displacement, and velocity probability density curves of the system all show a decreasing trend regardless of whether the single noise is enhanced or both noises are enhanced. The decrease in the peak value indicates that the probability of the system response increasing will increase, which means that the system response is enhanced. In addition, when $D_1$ increases and $D_2$ remains unchanged, the change of the three curves is gentler than that when $D_2$ increases and $D_1$ is unchanged. This phenomenon indicates that additive noise has a greater influence than multiplicative noise in this system. Similarly, it can be seen from Figures 14 and 15 that the trend of excitation amplitude contrast curve is consistent with that of noise intensity contrast curve.

Figure 12.

Probability density function of the system ( $a=0.5,{\omega }_0=1,K=1,e=0.9,{\mu }_1=0.3,F_1=1$ and $F_2=1$ ). (a) PDF of amplitude A; (b) PDF of displacement $x_1;$ (c) PDF of velocity $x_2$ .

Figure 13.

Time displacement diagram and time velocity diagram ( $a=0.5,{\omega }_0=1,K=1,e=0.9,{\mu }_1=0.3,F_1=1,F_2=1,{\mathrm{D}}_1=0.5$ and ${\mathrm{D}}_2=0.5$ ). (a) Time displacement diagram; (b) Time velocity diagram.

Figure 14.

Probability density function of the system ( $a=0.5,{\omega }_0=1,K=1,e=0.3,{\mu }_1=0.3,{\mathrm{D}}_1=0.5$ and ${\mathrm{D}}_2=0.5$ ). (a) PDF of amplitude A; (b) PDF of displacement $x_1;$ (c) PDF of velocity $x_2.$

Figure 15.

Time displacement diagram and time velocity diagram ( $a=0.5,{\omega }_0=1,K=1,e=0.3,{\mu }_1=0.3,F_1=1.5,F_2=1.5,{\mathrm{D}}_1=0.5$ and ${\mathrm{D}}_2=0.5$ ). (a) Time displacement diagram; (b) Time velocity diagram.

Discussion and conclusions

In practical scenarios, systems can be affected by multiple types of excitations. Therefore, considering the complexity of the excitation of the actual system, the steady-state response of the fractional-order collision system under the joint random excitation is studied, in which the joint random excitation is composed of additive and multiplicative random excitation. The original system is transformed into a smooth system by non-smooth transformation, and the steady-state probability density of the system is derived by the stochastic averaging method. Based on the probability density curve, time displacement curve, and time velocity curve of the system, the influence of parameters on the response of the system is analyzed. The research reveals the following:

(1) When the fractional-order decreases, the system will gradually deviate from the stable state. The closer the fractional-order is to 0, the more sensitive the system changes will be. In contrast, when the fractional-order increases, the peak value will increase, and the probability of large values of amplitude, displacement, and velocity will also decrease. However, a larger fractional-order will lead to a weaker response of the system. When the fractional-order coefficient decreases, the peak value of the PDF curve of the amplitude, displacement, and velocity of the system also declines, which leads to an increase in the probability of a strong response of the system.

(2) When the restoring force e gradually increases, the peak value of the PDF curve of the system amplitude, displacement, and velocity decreases. However, the probability of larger values of the amplitude, displacement, and velocity rises. This phenomenon indicates that within a certain range, the greater the restoring force e is, the stronger the vibration response of the system will be.

(3) When the system damping increases, the change in amplitude of the system displacement decreases. Not only does the change in amplitude decrease but the frequency of change increases as well. Increased damping inhibits the motion of the system. The larger the damping, the smaller the amplitude of the system will be, and the faster the vibration velocity becomes.

(4) Based on the influence of noise intensity on the system, the vibration response of the system will be enhanced regardless of whether a single noise is enhanced or both noises are enhanced. However, with the increase in noise intensity, the peak value of the PDF curve of the amplitude, displacement, and velocity of the system will decrease. The results show that the additive noise has a greater influence on the vibration response of the system than the multiplicative noise. Therefore, we should pay more attention to the effect of additive noise on the system in engineering practice.