Stochastic Stability Analysis of Coupled Viscoelastic Systems with Nonviscously Damping Driven by White Noise

Nonviscously damped structural system has been raised in many engineering fields, in which the damping forces depend on the past time history of velocities via convolution integrals over some kernel functions. This paper investigates stochastic stability of coupled viscoelastic system with nonviscously damping driven by white noise through moment Lyapunov exponents and Lyapunov exponents. Using the coordinate transformation, the coupled Itô stochastic differential equations of the norm of the response and angles process are obtained. Then the problem of the moment Lyapunov exponent is transformed to the eigenvalue problem, and then the second-perturbation method is used to derive the moment Lyapunov exponent of coupled stochastic system. Lyapunov exponent also can be obtained according to the relationship between moment Lyapunov exponent and Lyapunov exponent. Finally, the effects of various physical quantities of stochastic coupled system on the stochastic stability are discussed in detail. These results are validated by the direct Monte Carlo simulation technique.


Introduction
Viscoelastic materials are widely used in aerospace, construction, textile, and other industries, because they have a series of excellent properties, including light weight, high strength, wide source, and good shock absorption. These materials' stress depends on the past time history of stain; that is, the stress will increase correspondingly with the increase of time and the bucking will be more likely to occur. Therefore, the research of dynamical behavior of viscoelastic system has received a lot of interests in recent years [1][2][3][4].
To better understand the dynamical behavior of the viscoelastic system, many methods had been put forward. McTavish et al. [5] presented GHM method to analyze the linear multiple degree of freedom viscoelastic damping systems. Li et al. [6] studied the decay rate of energy functional of nonlinear viscoelastic full Marguerre-von Kármán shallow shell system. Later, Li [7] also showed the existence and uniqueness of weak solution of this system by applying the Galerkin finite element method. Nieto and Ahmad [8] used the generalized quasi-linearization method to obtain the explicit approximate solutions of an initial and terminal value problem for the forced Duffing equation with nonviscous damping. Recently, Li and Du [9] studied general energy decay of a degenerate viscoelastic Petrovsky-type plate equation with boundary feedback through a priori estimate and analysis of Lyapunovlike functional.
The dynamic stability as an important field of nonlinear dynamics has attached increasing attention in recent years. Guo et al. [10] derived a new simple sufficient condition of the global asymptotic stability of the integrodifferential systems with delay through constructing suitable Lyapunov functionals combined with the analytical technique. And then Meng et al. [11,12] obtained a sufficient condition of uniform stability for nonlinear integrodifferential system. However, when we consider the dynamical behavior of the system in a real environment, the stochastic excitation cannot be ruled out, such as wind gusts, earthquakes, and ocean waves [13][14][15]. Stochastic stability of the system under the stochastic excitation has attracted more and more attention [16][17][18][19][20]. The Lyapunov exponent specially has been researched by many scholars as an important indication for judging the stability of stochastic system. For example, Potapov [21] derived the sufficient condition of the almost-sure asymptotic stability of elastic and viscoelastic systems by the Lyapunov's direct method. He [22] also studied the stability of elastic and viscoelastic systems under the non-Gaussian excitation. For stochastic stability of high dimensional coupled system, Pavlović et al. [23] proposed the direct Lyapunov method to investigate the almost-sure stochastic stability of a viscoelastic double-beam system under parametric excitation. Then they used the same method to analyze the instability of coupled nanobeam systems subject to compressive axial loading [24]. These works only indicate that the system is almost-surely stable.
According to the theory of large deviation, which was first proposed by Baxendale and Stroock [25], the almost-sure stability of stochastic system does not mean the moment stability. Therefore, it is necessary to study the moment Lyapunov exponent of a stochastic system. The moment Lyapunov exponent can judge not only almost-sure stability but also the moment stability [26]. Nevertheless, analytic solutions of the moment Lyapunov exponent are very difficult to derive. To overcome this, many scholars developed different approximation methods to study the moment Lyapunov exponents of the stochastic system. For example, Sri Namachchivaya and Van Roessel [27] studied the moment Lyapunov exponents of two-degrees-of-freedom coupled elastic oscillators under real noise excitation by combining an asymptotic approximation for the moment Lyapunov exponents. Kozić and his associates [28][29][30] developed the first-order perturbation approach to obtain weak noise expansion of moment Lyapunov exponents and Lyapunov exponents for a stochastically coupled doublebeam system and Timoshenko beam system, respectively. Subsequently, Stojanović and Petković [31] used a perturbation method to study the moment Lyapunov exponents and the Lyapunov exponents of the three elastically connected Euler beams. More recently, Deng [32] applied the averaging method to establish the moment Lyapunov exponent of coupled gyroscopic stochastic system under bounded noise excitation. But it is rare to see studies of the moment Lyapunov exponent of a viscoelastic coupled system with nonviscous damping. Deng [33] investigated stochastic stability of two-degrees-of-freedom coupled viscoelastic system under white noise through moment Lyapunov exponent, but its damping term is viscous.
The purpose of this paper is to study stochastic stability of a viscoelastic coupled system with nonviscous damping subject to Gaussian white noise excitation through moment Lyapunov exponents and Lyapunov exponents, in which the nonviscous damping term is expressed by Boltzmanns superposition integral with a hereditary relaxation kernel. The paper is organized as follows. Section 2 derives the governing equations of motion of a stochastic viscoelastic system with a nonviscous damping. And then solving the problem of the moment Lyapunov exponent is changed into the eigenvalue and eigenfunction problem through a suitable transformation. In Section 3, the zeroth-order and the first-order and second-order perturbation solutions of the moment Lyapunov exponents are derived, respectively, based on the second-order perturbation method. The effects of different physical quantities on the stochastic stability of the coupled system are discussed under the given parameters.
Correspondingly, the results are verified by means of the direct Monte Carlo simulation in Section 4. Finally, conclusions will be drawn in Section 5.

Problem Formulation
Considering a coupled viscoelastic system with nonviscous damped structure driven by white noise excitation, the governing equations can be expressed as where 1 and 2 are generalized displacements. 1 and 2 are natural frequencies. is a small parameter. 1 and 2 are nonviscous damping coefficients. and are constants. ( ) is a Gaussian white noise with zero mean and noise intensity 2 . ( ) is a relaxation function of nonviscous damping, which can be expressed by in which is the relaxation parameter. Note that the limit case → 0, the nonviscous damping force reduces to classical viscous damping.
Let 1 = 1 ,̇1 = 1 2 , 2 = 3 ,̇2 = 2 4 ; system (1) can be represented aṡ  where represents the norm of the response. 0 ≤ 1 ≤ 2 and 0 ≤ 2 ≤ 2 are the angles process of the coupled stochastic system. 0 ≤ ≤ /2 denotes the coupling or exchange of energy between the coupled system. Equation (3) can be rewritten aṡ Therefore, the solution forms of can be expressed as = + ( ) from (6), where ( ) is quickly varying process relative to the process due to ≪ 1. Therefore, the following approximations will be obtained: Substituting (7) into the nonviscous damping, one obtains Similarly, the following results also can be derived: Substituting (8) and (9) into (6) and letting = = ( 2 1 + 2 2 + 2 3 + 2 4 ) /2 (−∞ < < +∞), (6) can be written as the following Itô stochastic differential equations: where ( ) is the standard Wiener process; the drift and diffusion coefficients can be rewritten as ⋅ (1 − cos 4 2 ) cos 4 }⟩ , 11 ( 1 , 2 , ) = −√ 1 cos 2 1 , To obtain the eigenvalue problem for the th moment Lyapunov exponents of a four-dimensional linear Itô stochastic system, a linear stochastic transformation [34] is introduced where is a new norm process depending on the transformation function ( 1 , 2 , ). The transformation function ( 1 , 2 , ) is defined on the independent stationary phase processes 1 , 2 , and in the ranges If the transformation function ( 1 , 2 , ) is bounded and nonsingular, both processes and have the same stability behavior. Therefore, transformation function ( 1 , 2 , ) is chosen so that the drift term of the Itô differential equation (14) does not depend on the phase processes 1 , 2 , and , so one can obtain = Λ ( ) Comparing (14) and (15), such transformation function ( 1 , 2 , ) can be written by the following equation: Mathematical Problems in Engineering 5 where 0 and 1 are the following first-order and secondorder differential operators: in which 6

Mathematical Problems in Engineering
In the following investigation, the perturbation theory will be used to obtain the solution Λ( ) of (16).

Moment Lyapunov Exponents
The method of deriving the eigenvalue problem for the moment Lyapunov exponents of a two-dimensional linear Itô stochastic system was first applied by Wedig [34]. The eigenvalue problem for a differential operator of three independent variables 1 , 2 , and will be identified from (16), in which Λ( ) is the eigenvalue and ( 1 , 2 , ) is the associated eigenfunction. Meanwhile, the eigenvalue Λ( ) is seen to be the Lyapunov exponent of the th moment of system (1) from (15). Applying the perturbation theory, both the moment Lyapunov exponent Λ( ) and the eigenfunction ( 1 , 2 , ) are expanded in power series of ; that is, Substituting (19) into (16) and equating terms of the equal powers of lead to the following equations: where ( 1 , 2 , ) must be positive and periodic in the range 0 ≤ 1 ≤ 2 , 0 ≤ 2 ≤ 2 . (17) into (20), the zeroth-order perturbation equation is

The First-Order Perturbation.
From (20), the first-order perturbation equation is Based on the analysis results in Section 3.1, one obtains Mathematical Problems in Engineering 7 Then the solvability condition of (30) is Applying the sense of the arbitrary of function ( ), (31) will become where ⋅ sin 2 − ( 2 1 + 2 2 ) 2 sin 2 cos 2 } , then (32) is simplified as the following equation: Then the boundary conditions of (34) are determined by considering the adjoint equation for the case of = 0, * ̃( ) = 0, where * is the Fokker-Planck operator and = 0 and = /2 are the entrance boundaries [27]. The eigenfunction 0 ( ) satisfies zero Neumann boundary condition, and Λ 1 ( ) is the largest eigenvalue of (34) with zeros Neumann boundary. Therefore, the solution of (34) can be calculated from an orthogonal expansion [34]. Under zeros Neumann boundary conditions, 0 ( ) can be expressed by a Fourier cosine series [27]; that is, Substituting (36) into (34), multiplying by cos 2 in both sides and integrating for , we can obtain In order to guarantee the existence of the nontrivial solution for each , the coefficient matrix = ( ) must equal zero. Therefore, the problem of calculating Λ 1 ( ) is translated into calculating the leading eigenvalue of . The sequence of approximations will be constructed by the eigenvalues of a sequence of the following submatrices: The comparison results of approximate analytical moment Lyapunov exponents obtained by truncating and the direct Monte Carlo simulation results are shown in Figure 1. It is easily found that the approximate results agree well with the simulation results. Therefore, the first-order perturbation will be reduced to Without loss of generality, there is a relationship between two frequencies of the form 1 1 = 2 2 , in which 1 and 2 are integers. The second frequency can be rewritten as 2 = 1 . Then the function ( 1 , 2 , ) in (41) can be rewritten as where function ( 1 , 2 ) is the periodic function on 1 , 2 and given as Note that the general solution 1 ( 1 , 2 , ) of (41) can not be obtained except in some special cases. Considering the nature Mathematical Problems in Engineering 9 of the coefficients of (41), the series expansion of function 1 ( 1 , 2 , ) can be presented in the following form: Substituting (44) into (41) and equating terms of the equal trigonometry function to give a set of partial differential equations, one obtains where the function 1 ( 1 , 2 ) can be written as in which ( 2 − 1 ) are arbitrary functions of two variables, and we make assumptions as follows: Here, the unknown constants 1 , 1 , 1 ( = 0, 1, 2, 3) will be determined by using the following conditions: That is,

The Second-Order Perturbation.
To further investigate the effects of parameter 2 on stochastic stability of system, the second-order perturbation equation of (20) will be used; that is, Based on the solvability condition, one obtains The solution has the following form: where the values Λ 21 ( ), Λ 22 ( ), Λ 23 ( ) are given in Appendix A. Substituting (40) and (52) into (19), the second-order approximate solution of the moment Lyapunov exponent is obtained as (3 + 10) − 4 ( 1 1 + 2 2 )] The corresponding Lyapunov exponents can be expressed as where the values 21 , 22 , and 23 are given in Appendix B.

Stochastic Stability Analysis
Through the above-mentioned analysis, the analytic expressions of moment Lyapunov exponents Λ( ) and Lyapunov exponents have been derived. Next, we use two indices to go on stochastic stability analysis. The th moment Lyapunov exponent is defined as where Ε[⋅] denotes the expected value and ‖ ‖ is a solution process of a random system. If Λ( ) < 0, Ε[‖ ‖ ] → 0 as → ∞, and it means that the th moment is stable. Lyapunov exponent is the slope of Λ( ) at the origin; that is, If < 0, it indicates that the system is almost-surely stable. The th moment and almost-sure stability boundary can be determined from Λ( ) = 0 and = 0, respectively. By comparing the second-order perturbation with the first-order perturbation, we find that the moment Lyapunov exponents Λ( ) and the Lyapunov exponents under the second-order perturbation contain all the parameter information of the coupled stochastic system. Therefore, we only consider the second-order perturbation in the subsequent research.
The boundaries of the th moment and the almost-sure stability with respect to nonviscous damping parameters 1 and 2 are shown in Figure 2. We find that the moment stability boundaries are more conservative than the almost-sure stability boundaries. This means that the almost-sure stability cannot assure the moment stability. And the boundaries of stability for th moment would increase with the increase of the order .
The effects of relaxation parameter on moment Lyapunov exponents are shown in Figure 3. It describes that the stability index and the root of Λ( ) = 0 will decrease and the system becomes more unstable. Correspondingly, the analytical results and the direct Monte Carlo results are shown in Figure 3(b) for the given . And Lyapunov exponent is negative; it means the system is almost-surely stable.
The moment Lyapunov exponents Λ( ) for different nonviscous damping parameter 2 are discussed in Figure 4. It is clear that the bigger values of 2 can enhance the system's stability for > 0. Correspondingly, the analytical results and the direct Monte Carlo results are shown in Figure 4(b) for the given different 2 . It illustrates that the system also is almost-surely stable. However, when becomes sufficiently large, Λ( ) > 0, the system becomes more unstable.
Next, the effects of on the moment Lyapunov exponents are studied in Figure 5. One easily finds the system is moment stability when 0 < < 3, and the stability of the system increases with the increase of from 0.01 to 0.1. Afterward, the moment Lyapunov exponents are changed from negative to positive; the system is varied from moment stability to moment instability with the increase of . Numerical results are verified by the direct Monte Carlo in Figure 5(b).
Then, in order to discuss the influences of the noise intensity on stochastic stability of coupled stochastic system, the moment Lyapunov exponents Λ( ) are shown in Figure 6 with change of values of . They describe that the coupled stochastic system could be almost-surely stable. Meanwhile, the moment stability of the coupled system would reduce with the increase of the noise intensity .
The variations of the moment Lyapunov exponents for different 2 are also analyzed. From Figure 7, we find that the coupled stochastic system may be almost-surely stable due to the value of slope Λ( ) less than 0 at the origin. However, the moment stability of system is not assured for < 0 or is sufficiently large due to the moment Lyapunov exponents greater than 0. Besides, one can easily see that the moment stability of the system enhances with the increase of 2 .
Finally, the effects of 2 on the moment Lyapunov exponents are investigated and shown in Figure 8. From those figures, the system may be almost-surely stable in terms of the Lyapunov exponents. It is also indicated that the system has the characteristic of the moment stability. And the moment stability of the coupled system will increase with the increase of 2 under given parameters. The approximate analytical results agree well with numerical simulation results by the direct Monte Carlo.

Conclusions
In this paper, the moment stability and almost-sure stability of coupled viscoelastic system with nonviscous damping subject to Gaussian white noise excitation are investigated. The nonviscously damped structure is assumed to follow the exponential integral relation, which is possible to model the rate of energy dissipation. Motion equations are transformed into the coupled Itô stochastic differential equations by means of the method of stochastic averaging. Then the linear transformation is introduced to derive the eigenvalue problem. In order to obtain the analytical expressions of eigenvalue, the second-order perturbation method is used. Finally, the effects of Gaussian noise, the relaxation parameter, the nonviscous damping coefficient, and two physical quantities, that is 2 , =0.05

Data Availability
All the numerical calculated data used to support the findings of this study can be obtained by calculating the equations in the paper, and the white noise is generated randomly. The codes used in this paper are available from the corresponding author upon request.