Stochastic stability of a fractional viscoelastic plate driven by non-Gaussian colored noise

In this paper, the moment Lyapunov exponent and stochastic stability of a fractional viscoelastic plate driven by non-Gaussian colored noise is investigated. Firstly, the stochastic dynamic equations with two degrees of freedom are established by piston theory and Galerkin approximate method. The fractional Kelvin–Voigt constitutive relation is used to describe the material properties of the viscoelastic plate, which leads to that the fractional derivation term is introduced into the stochastic dynamic equations. And the noise is simplified into an Ornstein–Uhlenbeck process by utilizing the path-integral method. Then, via the singular perturbation method, the approximate expansions of the moment Lyapunov exponent are obtained, which agree well with the results obtained by the Monte Carlo simulations. Finally, the effects of the noise, viscoelastic parameters and system parameters on the stochastic dynamics of the viscoelastic plate are discussed.


Introduction
In the past several decades, the flutter problem of plates has attracted more and more attention with the development of supersonic aircraft. Flutter instability will lead to the decrease in structure life and even cause structural failure. Viscoelastic materials are widely used in various parts of aircraft and high-speed trains because of their excellent vibration damping properties, which can suppress the fluttering of plates. On the other hand, random factors exist in nature, such as high winds, ocean waves, earthquakes, atmospheric turbulence, and so on. Therefore, determining the dynamic behaviors of the viscoelastic systems under random excitation is of great help for engineering applications, which also has been investigated by many researchers [1][2][3][4][5].
Although plenty of achievements have been made in the study of stochastic dynamics of elastic [6][7][8] and viscoelastic systems [9][10][11], the research on stochastic flutter of the viscoelastic plates is still limited. Potapov [12] investigated the stochastic flutter of a viscoelastic plate under white noise excitation through the stability conditions in the first and second statistical moments. Based on the results of largest Lyapunov exponent determined by the stochastic averaging method, Ariaratnam and Abdelrahman [13] examined the almost sure stability of a viscoelastic plate in the supersonic gas flow. Later, by applying numerical evaluation of the largest Lyapunov exponent, Potapov [14] studied the stability of elastic and viscoelastic plates under a stochastic excitation. The largest Lyapunov exponent also has been used to determine the almost sure stability of other viscoelastic systems, such as viscoelastic beam [15], viscoelastic column [16], non-gyroscopic viscoelastic system [17], and so on.
According to the large deviation theory [18], even for the stochastic dynamical system with almost sure stable, the pth moment of response can be unstable and grow exponentially. Hence, it is necessary to investigate the pth moment Lyapunov exponent of a stochastic dynamical system. The pth moment Lyapunov exponent is defined as follows: where x t; x 0 ð Þ is a solution process to a stochastic dynamical system. If K p; x 0 ð Þ\0, then E ½ xðt; x 0 Þ k k p ! 0 as t ! þ1, which implies that the pth moment is stable; otherwise, it is unstable. In Ref. [19], it has been disclosed that the limit of Eq. (1) exists and is independent of x 0 under the specified conditions. Thus K p; x 0 ð Þ can be expressed as K p ð Þ, which is a convex analytic function in p 2 R. And the largest Lyapunov exponent k is just the first-order term coefficient of moment Lyapunov exponent with respect to p [20], i.e., Based on the pth moment Lyapunov exponent, the stochastic stability, D-bifurcation and P-bifurcation [18,21] of a random system can be determined. Therefore, the pth moment Lyapunov exponent can describe the stochastic dynamic behavior of a system more profoundly and completely.
Although the pth moment Lyapunov exponent is very significant in the research of stochastic dynamics of a system, it is especially difficult to calculate the moment Lyapunov exponent. In recent years, researchers have studied the moment Lyapunov exponent for different stochastic dynamical systems, and achieved gratifying results [22][23][24][25][26][27]. However, the research on the viscoelastic plates is still limited. For a wide-band noise or non-Gaussian colored noise excited viscoelastic plate, Huang [28,29] studied the pth moment Lyapunov exponent and obtained an asymptotic expansion of the moment Lyapunov exponent using the stochastic averaging method. By determining the moment Lyapunov exponent, Deng [30] examined the stochastic stability of a viscoelastic plate under bounded noise excitation. Recently, Wu [31] discussed the moment stability of a viscoelastic plate under non-Gaussian colored noise excitation by calculating the pth moment Lyapunov exponent.
With the development of fractional calculus theory [32][33][34][35], it is found that the integer-order viscoelastic model cannot meet the actual needs, while the fractional viscoelastic model is more accurate and more in line with practical engineering. The constitutive relation of viscoelastic materials can be described by the fractional Kelvin-Voigt model, which also has been adopted in a lot of researches on viscoelastic systems [36][37][38]. For a fractional viscoelastic plate excited by non-Gaussian colored noise, the stochastic stability and pth moment Lyapunov exponent are investigated in the present paper. The fractional Kelvin-Voigt model is introduced to describe the material properties of the viscoelastic plates. And via the path-integral method [39][40][41], the non-Gaussian colored noise is reduced to an Ornstein-Uhlenbeck (O-U) process. By applying the singular perturbation method and Fourier series expansion, the approximate analytical results of the pth moment Lyapunov exponent and largest Lyapunov exponent are given. Furthermore, the numerical results of the original system are also obtained by using Monte Carle simulation. Finally, based on the pth moment Lyapunov exponent and largest Lyapunov exponent, the impacts of the noise, viscoelastic parameters and system parameters on the stochastic dynamics of the fractional viscoelastic plate are studied and discussed in detail.

Formulation
Considered an infinitely long viscoelastic plate with simply supported conditions is shown in Fig. 1. Based on the classical bending theory of the small deflection of thin plates, the governing equation for the viscoelastic plate with small oscillation is given as follows where w x; t ð Þ is the transversal displacement of plate; M x is the internal moment per unit length; r x ,e x are the stress and strain of the plate. b, l, h are the length, width and thickness of plate, respectively, and b ! 1. q is the material density; c is the viscous damping coefficient per unit length; N t ð Þ is the external random force inplane load per unit length. P x; t ð Þis aerodynamic force per unit area, which can be obtained by applying the piston theory.
where u,v, C 1 , and P 1 represent the flow velocity, gas flow polytropy index, undisturbed gas speed, and undisturbed gas pressure, respectively. The constitutive relation of viscoelastic materials can be given as where E is the elastic modulus;c 0 is the viscoelastic coefficient; D a t represents the Riemann-Liouville fractional derivative with a 2 0; 1 ð Þ, which is defined as With the nondimensional variables t 0 and x 0 , the corresponding nondimensional form of Eq. (7) is given as Bolotin [42] pointed out that quantitative results are reliably predicted with the use of the first two modes. Therefore, the transversal displacement w x 0 ; t 0 ð Þ can be expressed as which satisfies the boundary conditions: According to Eqs. (8) and (9), through using the Galerkin method and via a series of transformations and simplifications, the system (9) becomes, by introducing small parameter e (0\e ( 1), where assuming the damping term and viscoelastic term are both the second-order small quantity of parameter e, and the random term is a small quantity of parameter e.
, c Ã is a constant.

Approximation to the Markov process
Suppose that the process n t 0 ð Þ in system (10) is a non-Gaussian colored noise given as follows where s, D are the noise correlation time and noise intensity, respectively. r is the departure coefficient, denoting the departure from the Gaussian noise. g t 0 ð Þ is a Gaussian white noise and its statistical properties are For r ! 1, the process nðt 0 Þ is an O-U process with correlation time s, and its correlation function is expressed as nðt 0 ÞnðsÞ h i¼ ðD=sÞe À t 0 Às j j=s . And nðt 0 Þ is further simplified to a Gaussian white noise as s ! 0. As in Ref. [43], the stationary probability density P s n ð Þ of Eq. (11) can be normalized if and only if r 2 À1; 3 ð Þ, which is given as where Z is a normalization constant. According to Eq. (14), one can easily obtain the statistical properties of the process nðt 0 Þ: As r À 1 j j( 1, utilizing a path-integral method [43,44], one may get the following expression, with the effective noise correlation time defined as and the associated noise intensity defined as Hence, Eq. (10) is reduced to an O-U process with the associated noise intensity D 1 and correlation time s 1 , where Wðt 0 Þ is a wiener process with unit intensity, and the symbol '''' denotes the Stratonovich stochastic integral. Because of diffusion coefficient of the O-U process (19) is constant, the Wong-Zakai correction term of Eq. (19) is equal to zero, i.e., the Itô stochastic differential equation of Eq. (19) can be written as.
with the power spectral density.

Reduction in fractional differentiation
Applying the following transformation: Equation (10) is represented as where Based on the transformation Equation (24) can be written as According to Eq. (26), one may let where D a t 0 e q sin h cos / 2 ð Þ¼D a t 0 e q sin h cos u 2 cos x 2 t 0 ð Þ À D a t 0 e q sin h sin u 2 sin x 2 t 0 ð Þ : From Eqs. (26)(27), one can easily find that q, h, u 1 , and u 2 are slow variables. Thus in Eq. (28), they can be considered as constants [45,46]. Then Eq. (28) is written as Applying the fractional derivative formula [47][48][49], the following is obtained Based on the results of Rossikhin [47,48], as t 0 ! 1, Eq. (30) becomes Substituting Eqs. (29,31) into Eq. (26), we can get the follows

Moment Lyapunov exponent
Based on the fact of Eqs. (26,32), one is easy to find that the processes h; / 1 ; / 2 ; n ð Þare independent of the variable q. Therefore, the processes h; / 1 ; / 2 ; n ð Þ alone form a diffusive Markov process with the following generator: where The moment Lyapunov exponent K p ð Þ of system (26) is the principal simple eigenvalue of the operator L p ð Þ [19,50], i.e.,

Zeroth-order perturbation
From the definition of KðpÞ, we know that K 0 ðpÞ ¼ 0. Thus, Eq. (36) can be reduced to Using the method of separation of variables and letting w 0 ðh; / 1 ; / 2 ; nÞ ¼ Hð#ÞU 1 ð/ 1 ÞU 2 ð/ 2 ÞZðnÞ; One can easily obtain that U 1 ð/ 1 Þ ¼ C 1 e c 1 / 1 and U 2 ð/ 2 Þ ¼ C 2 e c 2 / 2 . Applying the periodic boundary conditions we have the constants The solution to Eq. (41) is ffiffiffiffi ffi a 0 p n=r 0 Þ, where erf Á ð Þ is the error function and j denotes the imaginary unit. Since Z n ð Þ is a bounded function as n ! AE1, it is required that C 4 ¼ 0, i.e., ZðnÞ ¼ C 3 . Therefore, we obtain that where F ðhÞ is an arbitrary function and P s n ð Þ is the stationary probability density of the process n t 0 ð Þ.

Second-order perturbation
Applying the above results in Sects. 5.1 and 5.2, Eq. (38) is reduced as follows The solvability condition to Eq. (48) is Since Eq. (49) always holds for an arbitrary function F ðhÞ, it leads to Z Making use of the correlation function of the process n and the power spectral density given, respectively, by where (53) is also presented as  ; and I is an identity matrix.
In order for Eq. (54) to have a nontrivial solution, the determinant of its corresponding coefficient matrix should be zero. Therefore, the calculation of K 2 p ð Þ is transformed into solving the eigenvalue of matrix C. We can obtain a series of approximations by calculating the eigenvalue of a series of submatrices. And the series of approximations will converge to the corresponding true eigenvalues as n ! 1. However, with the increase in n, the amount of calculation will increase drastically. Here, we get the approximations by the truncation of n, such as, K 2 p ð Þ ¼ c 00 as n ¼ 0.

Numerical results and discussions
In order to verify the accuracy of the results obtained by the perturbation method, Monte Carlo simulation is used to calculate the moment Lyapunov exponent of system (10). Here, the numerical algorithms presented in [52,53] are applied to evaluate the time series and pth moment Lyapunov exponent of system (10). For different cases, the approximate analytical and numerical results of the moment Lyapunov exponent K p ð Þ are plotted in Fig. 2, and the parameters are given as Fig. 2, it can be seen that along with the increase in n, the approximate analytical results converge (the moment Lyapunov exponents of the system are almost identical for the cases n ¼ 2 and n ¼ 3, which agree well with the numerical results. It shows that the result of the moment Lyapunov exponent K p ð Þ obtained by perturbation method is correct, and it is enough as n ¼ 3. Then, based on the approximate analytical results of moment Lyapunov exponent and largest Lyapunov exponent, the stochastic stability of the fractional viscoelastic plate is discussed. The effects of fractional order a on the stochastic stability are shown in Figs. 3, 4, 5, and 6. From Figs. 3 and 4, it is clear that the moment stability is enhanced and the viscoelastic plate gradually becomes stabilized with the increase in the parameter a. But the effect of fractional order a on the stochastic stability is not monotonic, which is related to the natural frequencies x 1 and x 2 of the system, as shown in Fig. 4. When the difference between two natural frequencies of the system is large enough or both of them are less than 1, the change curve has a minimum point on the interval (0, 1). The effect of the natural frequencies and fractional order on the stochastic stability is plotted in Fig. 5. In the view of Fig. 5, one can observe that with the introduction of fractional order factor, the stability of the system becomes very sensitive to the natural frequencies. According to Eq. (10), one easily finds that the natural frequencies of the system depend on the parameter a. Then, the effect of parameter a on the stochastic stability is discussed in Fig. 6. Along with increase in parameter a, the almost sure stability of the viscoelastic plate is reduced, and the D-bifurcation occurs when the parameter a increases to a critical value. From Fig. 6, one can also find that the fractional order a plays an important role in the stochastic stability and D-bifurcation behavior of the system, and the effect of parameter a on the stochastic stability is not monotonic. Therefore, the stochastic stability and D-bifurcation behavior of the system can be controlled by adjusting the fractional order a.
The moment Lyapunov exponent and largest Lyapunov exponent of the system are depicted in Figs. 7  Figs. 7 and 8, one easily finds that the noise intensity D has a significant effect on the stochastic stability of the viscoelastic plate. Along with the increase in noise intensity D, the system will be unstable. And when the noise intensity D is relatively high, the parameter a has a great influence on the stability of the viscoelastic plate. The departure  coefficient r can also reduce the stability of the system. However, the effect of the parameter r on the stability of the system is very small as r ! 1, which can be neglected. The noise correlation time s will enhance the stability of the system. And in view of Fig. 7b, it is easily found that the moment Lyapunov exponents of the system are almost identical for the cases s = 0.0001, s = 0.001, and s = 0.01, which implies the parameter s has little effect on the stochastic stability of the system if s is very small. Therefore, we should pay more attention to the noise and avoid strong noise as far as possible in practical engineering.
The effects of parameter b on the stochastic stability are shown in Figs. 9 and 10. It can be found that the stochastic stability of the viscoelastic plate is significantly enhanced with the increase in parameter b. According to Eqs. (9-10), the value of parameter b depends mainly on the damping coefficient c. Hence, the damping factor can effectively enhance the stochastic stability of the viscoelastic plate. Figures 11 and 12 plot the moment Lyapunov exponent and largest Lyapunov exponent of the system for different values of viscoelastic coefficient c. From Fig. 11, one easily finds that the stochastic stability of the viscoelastic plate can be enhanced with the increase in parameter c. But along with the increase in parameter a, the effect of viscoelastic coefficient c on the stochastic stability will not be monotonic, as shown in Fig. 12. Therefore, the effect of viscoelastic coefficient c on the stochastic stability of the viscoelastic plate is directly related to the natural frequency of the system.

Conclusion
Based on the moment Lyapunov exponents, the stochastic dynamics of fractional-order viscoelastic plate under non-Gaussian colored noise excitation is investigated in the present paper. The fractional Kelvin-Voigt constitutive relation is used to describe the material properties of the viscoelastic plate. And applying a path integral approach, the noise is simplified to an Ornstein-Uhlenbeck process. For weak noise excitations, the singular perturbation method is applied to evaluate the moment Lyapunov exponents. And via expanding the eigenfunction as a Fourier series, the approximate analytic solution of the moment Lyapunov exponent is given. Then, the largest Lyapunov exponent is obtained by employing the relationship with the moment Lyapunov exponent. Furthermore, in order to validate the accuracy of the approximate analytical results, the Monte Carlo simulation results are calculated, which agree well with the analytical results. Finally, the effects of the noise  10 Largest Lyapunov exponent for D ¼ 0:5; a ¼ 0:1; c ¼ 0:5 and the system parameters on the stochastic stability are discussed, and some results are presented. With the introduction of the fractional Kelvin-Voigt constitutive relation, the natural frequency will affect the stochastic stability of the viscoelastic plate.