Abstract
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.
1. 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 (DOF) systems under both white noise and real noise excitations, many interesting results have been obtained in the past [2, 3]. For DOF 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 DOF 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 DOF to DOF systems, but still only the Lyapunov exponent was calculated. 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 DOF 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.
2. 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 and the autocorrelation function , where denotes expectation and presents the Dirac delta function. The formal relationship iswhere is a standard Wiener process. The cosine power spectral density is a constant:with the sine power spectral density over the entire frequency range. The variance of a white noise process, denoting the average power, is . 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).

(a) White noise

(b) Real noise
Another important wideband noise is called real noise, which is often characterized by an Ornstein-Uhlenbeck process and is constructed asIf the initial condition is normally distributed , then is a stationary Gaussian process with mean, , and the autocorrelation function in exponential form . 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 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 , the real noise is reduced to Gaussian white noise with constant spectral density . 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.
3. Formulation
The systems considered are described by Stratonovich stochastic differential equations of the formwhere 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: . presents a viscoelastic operator given bywhere 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 yieldsSubstituting (8) into (6) yieldswhereThus, (6) can be written in amplitudes, , and phase, , by solving (9) and (10):
whereFor 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 formwhere and are -dimensional vectors of independent standard Wiener processes and and are drift coefficients and the diffusion matrices and , in which , , , and , are given bywhere 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 asWhen applying the averaging operator, the integration is performed over explicitly appearing only. Applying the transformation and changing the order of integration lead toSimilarly, it is found thatwhereare the sine and cosine transformations of the viscoelastic kernel function , 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 DOF systems (), one can obtain Stratonovich stochastic differential equations aswhere the cases for and are called symmetric coupling and skew-symmetric coupling, respectively. Equations (22a) and (22b) can be approximated by ItĂ´ stochastic differential equations:whereIt 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 () 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
4. 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, , is a periodic function in of period , and and denote the first and second derivatives of with respect to , respectively, andwhere and , are given in (25). Substituting the mathematical transformations [13] and into (29), one may find that the coefficients , , and 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 of the formHere only cosine functions are adopted because of appearance of the member . This Fourier expansion method is actually a method for solving partial differential equations [16].
Substituting this expansion and and into the eigenvalue problem (28), multiplying both sides by , , and performing integration with respect to from to yield a set of equations for the unknown coefficients :where
The eigenvalue can be obtained by solution of a polynomial equation as follows. Rearranging (31) leads to
The third-order submatrix, , is listed here. For convenience, is inserted in .where ,  , is the pseudodamping defined in (21), and the upper sign is taken when and the lower sign is taken when .
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:where are coefficients and are functions of ; denotes the approximate value of the 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 . One may approximate the moment Lyapunov exponent of the system by

4.1. Determination of the Largest Lyapunov Exponent
The th moment Lyapunov exponent is a convex analytic function in which passes through the origin and the slope at the origin is equal to the largest Lyapunov exponent ; that is,which is obtained directly from (35).
For comparison, the largest Lyapunov exponent for system (23) can be analytically derived from invariant probability density by solving a Fokker-Plank equation in () of [1] as follows:
(1) If , that is, ,where .
(2) If , that is, ,where .
(3) If , that is, ,The constants and are defined as
5. 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 , the eigenfunction in (30) is ; then, from , the moment Lyapunov exponent is determined asIf viscoelasticity is not considered, that is, and , and the noise is taken as a white noise, that is, , then (42) is reduced to the moment Lyapunov exponent for DOF linear systems subjected to white noise parametric excitation, reported in in [5], where a perturbation method was applied. The moment stability boundary is then obtained asIt is important to note one more time that if the viscoelasticity and coupling 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 in [17]).
The largest Lyapunov exponent is given by
The almost-sure stability region is found to be
From (43) and (45), the moment stability and the almost-sure stability boundaries are both a straight line in the coordinates with and .
When , the eigenfunction is ; the moment Lyapunov exponent can be determined fromThis is a quadratic equation, which is , where the analytical expression for moment Lyapunov exponent can be obtained. The general Lyapunov exponent is given by
It is noted that only those values of the excitation spectrum at frequencies , , and 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 and [18]. Note that and are for symmetric coupling and skew-symmetric coupling, respectively.
Case 1 (coupled systems under band-limited noise excitation close to ). The corresponding largest Lyapunov exponents are
Case 2 (coupled systems under band-limited noise excitation close to ). 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 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, can be taken as damping coefficients as a whole.
Lyapunov exponents listed in (50), (53), and (55) are special cases of (47). For example, substituting , , , and into (47) yields (55).
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 , the eigenfunction is ; the moment Lyapunov exponent can be solved fromThis is a cubic equation, which is . 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 , no explicit expressions can be presented, as a quartic equation is involved.
6. Stability Index
The stability index is the nontrivial zero of the moment Lyapunov exponent. Hence, it can be determined as a root-finding problem such that . Here we consider the case for white noise. When the order of Fourier expansion , from (42), the stability index is given by
When , 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 , the larger the stability index and then the more stable the system.

7. Simulation
Monte Carlo simulation is carried out to determine the moment Lyapunov exponents and to check the accuracy of the approximate results from the stochastic averaging. DenotingEquations (22a) and (22b) can be converted to the system of first-order Stratonovich stochastic differential equations:where