Stochastic Stability of Coupled Viscoelastic Systems Excited by Real Noise

The moment stochastic stability and almost-sure stochastic stability of two-degree-of-freedom coupled viscoelastic systems, under the parametric excitation of a real noise, are investigated through the moment Lyapunov exponents and the largest Lyapunov exponent, respectively. The real noise is also called the Ornstein-Uhlenbeck stochastic process. For small damping and weak random fluctuation, the moment Lyapunov exponents are determined approximately by using the method of stochastic averaging and a formulated eigenvalue problem. The largest Lyapunov exponent is calculated through its relation with moment Lyapunov exponents. The stability index, the stability boundaries, and the critical excitation are obtained analytically. The effects of various parameters on the stochastic stability of the system are then discussed in detail. Monte Carlo simulation is carried out to verify the approximate results of moment Lyapunov exponents. As an application example, the stochastic stability of a flexural-torsional viscoelastic beam is studied.


Introduction
In the study of dynamic stability of structures, Lyapunov exponents play an important role in characterizing sample or almost-sure stability of stochastic systems [1]. If the largest Lyapunov exponent is negative, the trivial solution of system is stable with probability 1; otherwise, it is unstable almost surely. Hence, the vanishing of the largest Lyapunov exponent indicates the almost-sure stability boundaries in parameter space.
It is known that the almost-sure stability cannot assure the moment stability due to the probabilistic theory of large deviation [1]. Because the moment Lyapunov exponents are not responsible to recognize only the moment stability and the almost-sure stability, it is important to take into consideration both the almost-sure stability and moment stability and to determine both the Lyapunov exponents and the moment Lyapunov exponents. If the th moment Lyapunov exponent is negative, then the th moment of the system is asymptotically stable. The vanishing of the th moment Lyapunov exponent indicates the th moment stability boundaries in parameter space.
For stability of one-degree-of-freedom (1DOF) systems under both white noise and real noise excitations, many interesting results have been obtained in the past [2,3]. For 2DOF coupled systems, however, relatively fewer studies can be found in the literature [4,5]. The reason is that the computational process is much more complicated even if only one more degree is considered, especially in the case when the viscoelasticity is involved. Due to the presence of parametric coupling, vibration initiated in one or more of the modes of motion would disturb other modes originally at rest under certain conditions.
For viscoelastic systems, Ariaratnam [6] was among the first who studied the stochastic almost-sure stability by evaluating the largest Lyapunov exponent and the rotation number. Ling et al. [7] discussed the response and almost-sure stability of a 1DOF viscoelastic system with strongly nonlinear stiffness under wideband noise excitations. Moment Lyapunov exponents of a viscoelastic system under the excitation of a wideband noise were investigated by using the averaging method of both first order and second order [8]. Abdelrahman [9] extended Ariaratnam's method from 1DOF to 2DOF systems, but still only the Lyapunov exponent was calculated.

Mathematical Problems in Engineering
Sufficient conditions for almost-sure stability were obtained for both elastic and viscoelastic columns under the excitation of a random wideband stationary process using Lyapunov's direct method [10]. Later, Potapov [11] described the behavior of stochastic viscoelastic systems evaluating the Lyapunov exponents numerically.
The objective of this paper is to study the moment stability and almost-sure stability of 2DOF coupled viscoelastic systems under the real noise excitations. This research is motivated by problems in the dynamic stability of viscoelastic systems subjected to stochastically fluctuating loads. The paper is different from [12], where a perturbation approach was used, and is unlike [4], where the Girsanov theorem and Feynman-Kac formula were applied and the viscoelasticity was not considered. The study in this paper is an extension of [13] from white noise process to real noise process, that is, Ornstein-Uhlenbeck process. Furthermore, this study carries out Monte Carlo simulation of moment Lyapunov exponents for coupled systems.

Real Noise Process
In many engineering applications, the excitations of dynamical systems can be described as wideband stochastic processes with relatively flat power spectral density functions over a large frequency range. The equations of motion of many physical systems under the excitations of wideband stochastic processes can be approximated by Stratonovich stochastic differential equations. Gaussian white noise process and real noise process are two typical wideband stochastic processes.
Gaussian white noise process is mathematically defined as a stationary process ( ) on −∞ < < +∞ with the mean E[ ( )] = 0 and the autocorrelation function E[ ( ) ( + )] = ( ) = 0 ( ), where E denotes expectation and ( ) presents the Dirac delta function. The formal relationship is where ( ) is a standard Wiener process. The cosine power spectral density ( ) is a constant: with the sine power spectral density Ψ( ) = ∫ +∞ −∞ sin( ) ( )d = 0 over the entire frequency range. The variance of a white noise process, denoting the average power, is E[ 2 ( )] = (0) = ∞. White noise process presents a stochastic excitation that is characterized to have equal intensity at different frequencies, thus not physically realizable and only an idealization of practical excitations. The total average power of a white noise process is infinite, as shown in Figure 1(a).
Another important wideband noise is called real noise, which is often characterized by an Ornstein-Uhlenbeck process and is constructed as If the initial condition 0 is normally distributed 0 ∼ (0, 2 /(2 )), then ( ) is a stationary Gaussian process with mean, E[ ( )] = 0, and the autocorrelation function in exponential form ( ) = E[ ( ) ( + )] = e − | | 2 /(2 ). The cosine and sine power spectral density are denoted as ( ) and Ψ( ), respectively: It is found that the parameter characterizes the bandwidth of the process and ( / ) 2 indicates the spectral density of the process. For large values of , a real noise's power will spread over a wide frequency band. Specifically when = √ 0 → ∞, the real noise is reduced to Gaussian white noise with constant spectral density ( ) = 0 . Thus, by suitably selecting parameter , as shown in Figure 1(b), a real noise may be used as the mathematical model of a wideband noise. When the excitation is a wideband process or the correlation function of random excitation decays to zero fast enough, a non-Markov process of physical responses can converge to a Markov process, whose governing Itô stochastic differential equations are obtained by using the method of stochastic averaging.

Formulation
The systems considered are described by Stratonovich stochastic differential equations of the form where are generalized coordinates, are damping constants, are natural frequencies, are constant coefficients, and are constants. ( ) is the excitation load which is a real noise process, an ergodic stochastic process with zero mean and sufficiently small correlation time. is a small parameter introduced to show small damping, light viscosity, and weak random fluctuation. These assumptions have been used by Xie and other researchers [1], which would make the analysis more convenient: 0 < | | ≪ 1. H presents a viscoelastic operator given by where H( ) is the viscoelastic kernel function.  Equation (6) describes approximately the motion of certain continuous viscoelastic structures whose equations of motion have been discretized from partial differential equations by some suitable technique such as Rayleigh-Ritz, Galerkin, finite differences, or finite elements [14]. Equation (6) can be expressed in polar coordinates by using the polar transformation: Differentiating with respect to and comparing withẏ ieldṡc Substituting (8) into (6) yieldṡ where Thus, (6) can be written in amplitudes, , and phase, , by solving (9) and (10):̇= , + 1/2 (0) , , Mathematical Problems in Engineering For small values of , the relaxation times of the ( ) and ( ) processes are much larger than the correlation time of the excitation processes ( ). The solution of the set of equations (12) converges in the weak sense and up to first order in to a diffusive Markov process whose governing Itô equations are of the form , , where ( ) and Ψ( ) are cosine and sine power spectral density of process ( ) at , respectively. The coupled oscillators are assumed to have noncommensurable frequencies; that is, ̸ = , ̸ = . Both the drift coefficients and the diffusion coefficients are obtained by using the method of stochastic averaging [1]. The averaging operator is defined as When applying the averaging operator, the integration is performed over explicitly appearing only. Applying the Mathematical Problems in Engineering 5 transformation = − and changing the order of integration lead to Similarly, it is found that where are the sine and cosine transformations of the viscoelastic kernel function H( ), respectively. The term , may be called pseudodamping, because it plays the role of damping but includes viscoelasticity as well. Now consider a special case of (6). For 2DOF systems ( = 2), one can obtain Stratonovich stochastic differential equations as̈1 where the cases for 12 = 21 = and 12 = − 21 = are called symmetric coupling and skew-symmetric coupling, respectively. Equations (22a) and (22b) can be approximated by Itô stochastic differential equations: It is important to note that both the averaged amplitude and phase angle equations do not involve the phase angles and the amplitude equations are advantageously decoupled from the phase angle equations. Hence, the averaged amplitude vector ( 1 , 2 ) is a two-dimensional diffusion process. In this study, the viscoelastic kernel function is supposed to follow ordinary Maxwell model [15]: which can be used as an approximation to most linear viscoelastic behaviour as closely as possible if enough number of Maxwell units are arranged in parallel. Its sine and cosine transformations in (20) are given, respectively, by

Moment Lyapunov Exponents
By using mathematical transformations, the following eigenvalue problem for the moment Lyapunov exponent can be obtained from (23) [13]: where Λ is the Lyapunov exponent of the th moment of system (22a) and (22b), that is, Λ = Λ q( ) ( ), ( ) is a periodic function in of period , and and denote 6 Mathematical Problems in Engineering the first and second derivatives of ( ) with respect to , respectively, and where and , , = 1, 2, are given in (25). Substituting the mathematical transformations [13] 1 = cos and 2 = sin into (29), one may find that the coefficients 0 , 1 , and 2 are function of and only. The coefficients in (28) are periodic with period , so it is reasonable to consider a Fourier cosine series expansion of the eigenfunction ( ) in the order + 1 of the form Here only cosine functions are adopted because of appearance of the member 1/sin (2 ). This Fourier expansion method is actually a method for solving partial differential equations [16]. Substituting this expansion and 1 = cos and 2 = sin into the eigenvalue problem (28), multiplying both sides by cos(2 ), = 0, 1, . . . , , and performing integration with respect to from 0 to /2 yield a set of equations for the unknown coefficients , = 0, 1, . . . , : (32) The eigenvalue can be obtained by solution of a polynomial equation as follows. Rearranging (31) The third-order submatrix, = 2, is listed here. For convenience, Λ is inserted in .
To have a nontrivial solution of , it is required that the determination of the coefficient matrix of (33) equals zero, from which the eigenvalue Λ( ) can be obtained: moment Lyapunov exponent under the assumption that the expansion of eigenfunction ( ) is up to th order Fourier cosine series. The set of approximate eigenvalues obtained by this procedure converges to the corresponding true eigenvalues as → ∞. However, as shown in Figure 2, the approximate eigenvalues converge so quickly that the approximations almost coincide after the order ≥ 1. One may approximate the moment Lyapunov exponent of the system by Λ q( ) ( ) ≈ Λ ( ) ( ) . (36)

Stability Boundary
Moment Lyapunov exponents can be numerically determined from (35). To obtain analytical stability boundaries, some cases of lower order of are discussed in the following. When = 0, the eigenfunction in (30) is ( ) = 0 ; then, from 00 − Λ (0) = 0, the moment Lyapunov exponent is determined as If viscoelasticity is not considered, that is, = 0 and = , and the noise is taken as a white noise, that is, (2 1 ) = (2 2 ) = ( ± ) = ( ∓ ), then (42) is reduced to the moment Lyapunov exponent for 2DOF linear systems subjected to white noise parametric excitation, reported in (3.19) in [5], where a perturbation method was applied. The moment stability boundary is then obtained as It is important to note one more time that if the viscoelasticity and coupling ( = 0) are neglected and a white noise is assumed, then (43) can be reduced to already obtained one by other approximate methods such as asymptotic expansion of integrals and stochastic averaging (see (40) in [17]). The largest Lyapunov exponent is given by The almost-sure stability region is found to be From (43) and (45) It is noted that only those values of the excitation spectrum at frequencies 2 1 , 2 2 , and 1 ± 2 have an effect on the stability condition. To simplify the expressions for moment Lyapunov exponents, one can consider some cases of real noise (or band-limited noise) with special power spectral density ( ), where 0 − Δ 0 /2 < < 0 + Δ 0 /2 and Δ 0 ≪ 0 [18]. Note that Λ (1)+ ( ) and Λ (1)− ( ) are for symmetric coupling and skew-symmetric coupling, respectively.
The corresponding largest Lyapunov exponents are .
The corresponding largest Lyapunov exponents are Case 3 (coupled systems under band-limited noise excitation close to ( − )).
The corresponding largest Lyapunov exponents are The analytically obtained moment Lyapunov exponents and Lyapunov exponents in (49) to (55) provide invaluable tools for stability characterization with respect to the variation of the spectral density. For example, qualitatively (51) for symmetric coupling can give a stability index and provide stability boundaries because it is a polynomial of . However, (52) for skew symmetric coupling has no stability index, which is a consequence of a linear function of . Quantitatively, the second-order expansion = 1 is almost as accurate as higher-order expansion, as shown in Figure 2. It is not surprisingly found that (49) to (54) can be reduced to corresponding cases in [4], where viscoelasticity is not considered. However, it is of importance to note that viscoelasticity plays the same role in stability assessment as the damping does. Therefore, ( = 1, 2) can be taken as damping coefficients as a whole.
It is noted that the analytical expressions from (38) to (40) seem not to agree with the expression given in (47). This is a direct consequence of the discrepancy between exact solutions and the sequence of approximations. When = 2, the eigenfunction is ( ) = 0 + 1 cos 2 + This is a cubic equation, which is (Λ (2) ) 3 + (2) 2 (Λ (2) ) 2 + (2) 1 Λ (2) + (2) 0 = 0. The coefficients are too long to list here. The analytical expression for moment Lyapunov exponents can then be obtained by solving this cubic equation. However, for ⩾ 3, no explicit expressions can be presented, as a quartic equation is involved.

Stability Index
The stability index is the nontrivial zero of the moment Lyapunov exponent. Hence, it can be determined as a rootfinding problem such that Λ q( ) ( q( ) ) = 0. Here we consider the case for white noise. When the order of Fourier expansion = 0, from (42), the stability index is given by When = 1, the stability index is the nontrivial solution of (46), which is hard to express analytically. Typical results of the stability index are shown in Figure 3. The stability index decreases from positive to negative values with the increase of the amplitude of power spectrum, which suggests that the noise destabilizes the system. The larger the pseudodamping coefficient 1 , the larger the stability index and then the more stable the system.

Simulation
Monte Carlo simulation is carried out to determine the moment Lyapunov exponents and to check the accuracy Equations (22a) and (22b) can be converted to the system of first-order Stratonovich stochastic differential equations:  ] . (60) The associated Itô stochastic differential equations are the same as the Stratonovich stochastic differential equations. Thus, the discretized equations using explicit Euler scheme are with Δ being the time step and denoting the th iteration.

Application
As an application, the flexural-torsional stability of a simply supported, uniform, narrow, rectangular, viscoelastic beam of length subjected to a stochastically varying concentrated load ( ) acting at the center of the beam cross section is considered. Both nonfollower and follower loading cases are studied.
For a nonfollower force case, the only difference is that = 0 in (64). The equations of motion are of the same form of (66), but the parameters are different: It is seen that (66) has the same form of (22a) and (22b), except that 11 = 22 = 0. By introducing the polar transformation and using the method of stochastic averaging, Substituting (69) and 1 = cos and 2 = sin into (28) yields an eigenvalue problem, from which moment Lyapunov exponents can be determined by solving (35).
The almost-sure and th moment stability boundaries for systems under real noise excitations are shown in Figure 5. The almost-sure stability region is not always larger than that of moment stability region, which is different from systems under white noise excitations. This may be due to the fact that the spectral density function is constant over all frequencies for the white noise but concentrated in a small area around ( 1 + 2 ) for real noise, changeable with the property of boundaries. Figure 6 illustrates that the moment stability is conservative than the almost-sure stability, and the higherorder moment stability is more conservative than the lowerorder moment stability. Figures 7 and 8 illustrate the variation of moment Lyapunov exponents with respect to the parameters of real noise processes. The larger noise parameter leads to the following: the narrower the stability region for > 0, the more unstable the system. On the contrary, parameter plays a stabilizing role; that is, with the increase of , the stability region of the th moment (for > 0) increases and the system is more stable, as shown in Figure 8. This can be explained from (4), where it is suggested that smaller or larger reduces the power cosine spectral density. Then the power of noise would spread over a wider frequency band, which in turn reduces the power of the noise in the neighborhood of resonance and stabilizes the system. The good agreement between analytical approximate results and numerical simulation results shows that the tedious derivation in this paper is correct.

Conclusions
The stochastic stability of coupled viscoelastic systems described by stochastic integrodifferential equations of 2DOF was investigated. The system was parametrically excited by a real noise of small intensity and small damping. The equations of motion were firstly decoupled into two-dimensional Itô stochastic differential equations with the method of stochastic averaging. The moment Lyapunov exponents were determined approximately from a formulated eigenvalue problem and were verified by results from the Monte Carlo simulation. Analytical expressions for the stability boundaries and the stability index were then obtained. It is found that they are in good agreement with other analytical results. As an application, the flexural-torsional stability of a simply supported rectangular viscoelastic beam subjected to a stochastically varying concentrated load acting at the center of the beam cross section was considered. Under the real noise excitation, the parameters of damping and the noise intensity have stabilizing effects on both the moment stability and almost-sure stability. However, the real noise parameter plays a destabilizing role. These results are useful in engineering applications.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The author declares that there are no conflicts of interest regarding the publication of this paper.