# Stochastic Stability of Coupled Viscoelastic Systems Excited by Real Noise

Authors: Jian Deng
Journal: Mathematical Problems in Engineering
 (2018)

---

## 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.


---

## Body

## 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 pth moment Lyapunov exponent is negative, then the pth moment of the system is asymptotically stable. The vanishing of the pth moment Lyapunov exponent indicates the pth 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. 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 of2DOF 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ξ(t) on -∞<t<+∞ with the mean Eξt=0 and the autocorrelation function Eξtξt+τ=Rτ=S0δ(τ), where E denotes expectation and δ(τ) presents the Dirac delta function. The formal relationship is(1)ξt=σW˙t,σ=S0,where W(t) is a standard Wiener process. The cosine power spectral density S(ω) is a constant:(2)Sω=∫-∞+∞cos⁡ωτRτdτ=∫-∞+∞e-iωτS0δτdτ=S0,with the sine power spectral density Ψω=∫-∞+∞sin⁡ωτRτdτ=0 over the entire frequency range. The variance of a white noise process, denoting the average power, is Eξ2t=R0=∞. 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).Figure 1
Autocorrelation functionR(τ) and corresponding power spectral density function S(ω).
(a)
White noise (b)
Real noiseAnother important wideband noise is called real noise, which is often characterized by an Ornstein-Uhlenbeck process and is constructed as(3)dξt=-αξtdt+σdWt,ξt0=ξ0.If the initial condition ξ0 is normally distributed ξ0~N(0,σ2/(2α)), then ξ(t) is a stationary Gaussian process with mean, Eξt=0, and the autocorrelation function in exponential form Rτ=E[ξ(t)ξ(t+τ)]=e-α|τ|σ2/(2α). The cosine and sine power spectral density are denoted as S(ω) and Ψ(ω), respectively:(4)Sω=∫-∞+∞cos⁡ωτe-ατσ22αdτ=σ2α2+ω2=σ2α211+ω/α2,(5)Ψω=∫-∞+∞sin⁡ωτe-ατσ22αdτ=ωσ2αα2+ω2.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 σ=αS0→∞, the real noise is reduced to Gaussian white noise with constant spectral density Sω=S0. 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 form(6)q¨i+2ε∑j=1nβijq˙j+ωi2qi-εωi2∑j=1nτijεHqj+ε1/2ωiξt∑j=1nkijtqj=0,i=1,2,…,n,where qi are generalized coordinates, βij are damping constants, ωi are natural frequencies, τijε are constant coefficients, and kij are constants. ξ(t) 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(7)Hψt=∫0tHt-τψτdτ,0⩽∫0∞Hθdθ<1,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:(8)qi=aicos⁡Φi,q˙i=-ωiaisin⁡Φi,Φi=ωit+ϕi,i=1,2,…,n.Differentiating qi with respect to t and comparing with q˙i yields(9)a˙icos⁡Φi-aiϕ˙isin⁡Φi=0.Substituting (8) into (6) yields(10)a˙isin⁡Φi+aiϕ˙icos⁡Φi=1ωiFia,Φ,where(11)Fia,Φ=-2ε∑j=1nβijωjajsin⁡Φj-εωi2∑j=1nτijεHajcos⁡Φj+ε1/2ωiξt∑j=1nkijajcos⁡Φj.Thus, (6) can be written in amplitudes, ai, and phase, ϕi, by solving (9) and (10):(12)a˙i=εFa,i1+ε1/2Fa,i0,ϕ˙i=εFϕ,i1+ε1/2Fϕ,i0,where(13)Fa,i0=sin⁡Φiξt∑j=1nkijajcos⁡Φj,Fa,i1=-2ωi∑j=1nβijajωjsin⁡Φisin⁡Φj-sin⁡Φiωi∑j=1nτijεHajcos⁡Φj,Fϕ,i0=1aicos⁡Φiξt∑j=1nkijajcos⁡Φj,Fϕ,i1=-2aiωi∑j=1nβijajωjcos⁡Φisin⁡Φj-1aicos⁡Φiωi∑j=1nτijεHajcos⁡Φj.For small values of ε, the relaxation times of the ai(t) and ϕi(t) processes are much larger than the correlation time of the excitation processes ξ(t). 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(14)dai=εmiadt+ε1/2∑j=12σijadWja,(15)dϕi=εmiϕdt+ε1/2∑j=12σijϕdWjϕ,i=1,2,…,n,where Wa={W1a,W2a,…,Wna}T and Wϕ={W1ϕ,W2ϕ,…,Wnϕ}T are n-dimensional vectors of independent standard Wiener processes and εmia and εmiϕ are drift coefficients and the n×n diffusion matrices εba=εσa(σa)T and εbϕ=εσϕ(σϕ)T, in which σa=σija, σϕ=σijϕ, ba=bija, and bϕ=bijϕ, are given by(16)mia=MtFa,i1+∫-∞0E∑j=12∂Fa,i0∂ajFa,jτ0+∂Fa,i0∂ϕjFϕ,jτ0dτ=ai-βii-ωiτiiεMtIisc+316kii2S2ωi+18∑j=1j≠inkijkjiSij-+116∑j=1j≠inaj2aikij2Sij+,miϕ=MtFϕ,i1+∫-∞0E∑j=12∂Fϕ,i0∂ajFϕ,jτ0+∂Fϕ,i0∂ϕjFϕ,jτ0dτ=-ωiτiiεMtIicc+18∑j=1nkijkjiΨij-1+ajai-14∑j=1nkij2Ψij+aj2ai2,biia=Mt∫-∞∞EFa,i0Fa,iτ0dτ=18kii2ai2Sii2ωi+∑j=1j≠inkij2aj2Sij+,bija=Mt∫-∞∞EFa,i0Fa,jτ0dτ=18aiajkijkjiSij-,i≠j,biiϕ=Mt∫-∞∞EFϕ,i0Fϕ,iτ0dτ=18kii22S0+S2ωi+18∑j=1j≠inaj2ai2kij2Sij+,bijϕ=Mt∫-∞∞EFϕ,i0Fϕ,jτ0dτ=14kiikjjS0+18kijkjiSij+,i≠j,Fjτ0=Fj0a,ϕ,ξt+τ,t+τ,i,j=1,2,…,n,Sij±=Sωi+ωj±Sωi-ωj,Ψij±=Ψωi+ωj±Ψωi-ωj,where S(ω) and Ψ(ω) are cosine and sine power spectral density of process ξ(t) at ω, respectively. The coupled oscillators are assumed to have noncommensurable frequencies; that is, ωi≠ωj,i≠j. Both the drift coefficients and the diffusion coefficients are obtained by using the method of stochastic averaging [1]. The averaging operator is defined as(17)Mτ·=limT→∞1T∫ττ+T·dτ.When applying the averaging operator, the integration is performed over explicitly appearing τ only. Applying the transformation s=t-τ and changing the order of integration lead to(18)MtIisc=limT→∞1T∫t=0T∫s=0tHt-ssin⁡Φitcos⁡Φisdsdt=12∫0∞Hτsin⁡ωiτdτ=12Hsωi.Similarly, it is found that(19)MtIicc=Mtcos⁡Φit∫0tHt-scos⁡Φisds=12Hcωi,where(20)Hsω=∫0∞Hτsin⁡ωτdτ,Hcω=∫0∞Hτcos⁡ωτdτare the sine and cosine transformations of the viscoelastic kernel function H(t), respectively. The term Ei,(21)Ei=βi+ωiτiiεMtIisc,i=1,2,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 (n=2), one can obtain Stratonovich stochastic differential equations as(22a)q¨1+2εβ1q˙1+ω121-εHq1+ε1/2ω1k11q1+k12q2ξt=0,(22b)q¨2+2εβ2q˙2+ω221-εHq2+ε1/2ω2k21q1+k22q2ξt=0,where the cases for k12=k21=k and k12=-k21=k are called symmetric coupling and skew-symmetric coupling, respectively. Equations (22a) and (22b) can be approximated by Itô stochastic differential equations:(23)dai=εmiadt+ε1/2∑j=12σijadWja,(24)dϕi=εmiϕdt+ε1/2∑j=12σijϕdWjϕ,i=1,2,where(25)mia=ai-βi-ωiMtIisc+316kii2S2ωi+18kijkjiS-+116aj2aikij2S+,biia=18kii2ai2S2ωi+kij2aj2S+,bija=18aiajkijkjiS-,S±=Sω1+ω2±Sω1-ω2.It is important to note that both the averaged amplitude ai and phase angle equations ϕi do not involve the phase angles and the amplitude equations are advantageously decoupled from the phase angle equations. Hence, the averaged amplitude vector (a1,a2) is a two-dimensional diffusion process. In this study, the viscoelastic kernel function is supposed to follow ordinary Maxwell model [15]:(26)Ht=γe-ηt,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(27)Hsω=ωγη2+ω2,Hcω=γηη2+ω2.


## 4. Moment Lyapunov Exponents

By using mathematical transformations, the following eigenvalue problem for the moment Lyapunov exponent can be obtained from (23) [13]:(28)LpT=λ2Tφφ′′+λ1Tφ′+λ0T=ΛT,0⩽φ<π,where Λ is the Lyapunov exponent of the pth moment of system (22a) and (22b), that is, Λ=Λq(t)(p), T(φ) is a periodic function in φ of period π, and Tφ′ and Tφφ′′ denote the first and second derivatives of T(φ) with respect to φ, respectively, and(29)λ2=12b11aa12+a22sin2⁡φ-b12a+b21aa12+a22sin⁡φcos⁡φ+b22aa12+a22cos2⁡φ,λ1=m2acos⁡φa12+a22-m1asin⁡φa12+a22+b11asin⁡2φ2a12+a22-b12a+b21a2cos2⁡φ-12a12+a22-b22asin⁡2φ2a12+a22+p-b11aa12+a22cos⁡φsin⁡φ+b12aa12+a22cos2⁡φ-sin2⁡φ+b22aa12+a22cos⁡φsin⁡φ,λ0=pm1aa12+a22cos⁡φ+m2aa12+a22sin⁡φ+12b11aa12+a22p-2cos2⁡φ+1+b12a+b21a2a12+a22p-2sin⁡2φ+b22aa12+a22p-2sin2⁡φ+1,where mia and bija,i,j=1,2, are given in (25). Substituting the mathematical transformations [13] a1=rcos⁡φ and a2=rsin⁡φ into (29), one may find that the coefficients λ0, λ1, and λ2 are function of φ and p only.The coefficients in (28) are periodic with period π, so it is reasonable to consider a Fourier cosine series expansion of the eigenfunction T(φ) in the order K+1 of the form(30)Tφ=∑i=0KCicos⁡2iφ.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 anda1=rcos⁡φ and a2=rsin⁡φ into the eigenvalue problem (28), multiplying both sides by cos⁡(2jφ), j=0,1,…,K, and performing integration with respect to φ from 0 to π/2 yield a set of equations for the unknown coefficients Ci,i=0,1,…,K:(31)∑i=0KaijCi=ΛCj,j=0,1,2,…,where(32)aij=4π∫0π/2Lpcos⁡2iφcos⁡2jφdφ,j=0,1,…,K,∫0π/2cos⁡2iφcos⁡2jφdφ=π4,i=j,0,i≠j.The eigenvalue can be obtained by solution of a polynomial equation as follows. Rearranging (31) leads to(33)a00-Λa01a02a03⋯a10a11-Λa12a13⋯a20a21a22-Λa23⋯a30a31a32a33-Λ⋯⋯⋯⋯⋯⋱C0C1C2C3⋯=0.The third-order submatrix,K=2, is listed here. For convenience, Λ is inserted in aij.(34)a00=11283k112S2ω1+3k222S2ω2+4k2Sω±p2+1645k112S2ω1+5k222S2ω2+12k2Sω±-32E1+E2p-Λ,a01=164k112S2ω1-k222S2ω2p2+132k112S2ω1-k222S2ω2-8E1-E2p,a02=1256k112S2ω1+k222S2ω2-4k2Sω±p2-1128k112S2ω1+k222S2ω2-8k2Sω±p,a10=164k112S2ω1-k222S2ω2p2+116k112S2ω1-k222S2ω2-4E1-E2p⁡+116k112S2ω1-k222S2ω2-8E1-E2,a11=15127k112S2ω1+7k222S2ω2+4k2Sω±p2+125611k112S2ω1+11k222S2ω2+20k2Sω±-64E1+E2p-164k112S2ω1+k222S2ω2+12k2Sω±+16k2Sω∓-Λ2,a12=1128k112S2ω1-k222S2ω2p2-18E1-E2p-132k112S2ω1-k222S2ω2-8E1-E2,a20=1256k112S2ω1+k222S2ω2-4k2Sω±p2+3128k112S2ω1+k222S2ω2-4k2Sω±p+132k112S2ω1+k222S2ω2-20k2Sω±-16k2Sω∓,a21=1128k112S2ω1-k222S2ω2p2+1643k112S2ω1-3k222S2ω2-8E1-E2p+116k112S2ω1-k222S2ω2-8E1-E2,a22=12563k112S2ω1+3k222S2ω2+4k2Sω±p2+11285k112S2ω1+5k222S2ω2+12k2Sω±-32E1+E2p-116k112S2ω1+k222S2ω2+8k2Sω±+12k2Sω∓-Λ2,where ω∓=ω1∓ω2,  ω±=ω1±ω2, Ei is the pseudodamping defined in (21), and the upper sign is taken when k12=k21=k and the lower sign is taken when k12=-k21=k.To have a nontrivial solution ofCk, it is required that the determination of the coefficient matrix of (33) equals zero, from which the eigenvalue Λ(p) can be obtained:(35)eK+1KΛKK+1+eKKΛKK+⋯+e1KΛK+e0K=0,where ekK(k=0,0,…,K+1) are coefficients and are functions of aij; Λ(K) denotes the approximate value of the moment Lyapunov exponent under the assumption that the expansion of eigenfunction T(φ) is up to Kth order Fourier cosine series. The set of approximate eigenvalues obtained by this procedure converges to the corresponding true eigenvalues as K→∞. However, as shown in Figure 2, the approximate eigenvalues converge so quickly that the approximations almost coincide after the order K≥1. One may approximate the moment Lyapunov exponent of the system by(36)Λqtp≈ΛKp.Figure 2
Moment Lyapunov exponents for differentKth order Fourier expansion.


### 4.1. Determination of the Largest Lyapunov Exponent

Thepth moment Lyapunov exponent Λq(t)(p) is a convex analytic function in p which passes through the origin and the slope at the origin is equal to the largest Lyapunov exponent λq(t); that is,(37)λqtp=limp→0∂ΛKp∂p=-limp→0e0Kpe1K,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 (8.6.14) of [1] as follows:(1) Ifk112S2ω1+k222S2ω2>4|k12k21|S(ω1±ω2), that is, Δ0>0,(38)λ=12λ1+λ2+λ1-λ2coth⁡λ1-λ2Δ0α+18k12k21S-,where α=cosh-1⁡Q/(2k12k21S+).(2) Ifk112S2ω1+k222S2ω2<4k12k21S(ω1±ω2), that is, Δ0<0,(39)λ=12λ1+λ2+λ1-λ2cothλ1-λ2-Δ0α+18k12k21S-,where α=cos-1Q/(2k12k21S+).(3) Ifk112S2ω1+k222S2ω2=4k12k21S(ω1±ω2), that is, Δ0=0,(40)λ=12λ1+λ2+λ1-λ2coth4λ1-λ2k12k21S++18k12k21S-.The constants Q and Δ0 are defined as(41)Q=k112S2ω1+k222S2ω2-2k12k21S+,Δ0=164Q2-4k122k212S+2,λ1=-β1-ω1MtI1sc+18k112S2ω1,λ2=-β2-ω2MtI2sc+18k222S2ω2.


## 4.1. Determination of the Largest Lyapunov Exponent

Thepth moment Lyapunov exponent Λq(t)(p) is a convex analytic function in p which passes through the origin and the slope at the origin is equal to the largest Lyapunov exponent λq(t); that is,(37)λqtp=limp→0∂ΛKp∂p=-limp→0e0Kpe1K,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 (8.6.14) of [1] as follows:(1) Ifk112S2ω1+k222S2ω2>4|k12k21|S(ω1±ω2), that is, Δ0>0,(38)λ=12λ1+λ2+λ1-λ2coth⁡λ1-λ2Δ0α+18k12k21S-,where α=cosh-1⁡Q/(2k12k21S+).(2) Ifk112S2ω1+k222S2ω2<4k12k21S(ω1±ω2), that is, Δ0<0,(39)λ=12λ1+λ2+λ1-λ2cothλ1-λ2-Δ0α+18k12k21S-,where α=cos-1Q/(2k12k21S+).(3) Ifk112S2ω1+k222S2ω2=4k12k21S(ω1±ω2), that is, Δ0=0,(40)λ=12λ1+λ2+λ1-λ2coth4λ1-λ2k12k21S++18k12k21S-.The constants Q and Δ0 are defined as(41)Q=k112S2ω1+k222S2ω2-2k12k21S+,Δ0=164Q2-4k122k212S+2,λ1=-β1-ω1MtI1sc+18k112S2ω1,λ2=-β2-ω2MtI2sc+18k222S2ω2.


## 5. Stability Boundary

Moment Lyapunov exponents can be numerically determined from (35). To obtain analytical stability boundaries, some cases of lower order of K are discussed in the following.WhenK=0, the eigenfunction in (30) is Tφ=C0; then, from a00-Λ0=0, the moment Lyapunov exponent is determined as(42)Λqtp≈Λ0p=11283k112S2ω1+3k222S2ω2+4k2Sω±p2+1645k112S2ω1+5k222S2ω2+12k2Sω±-32E1+E2p.If viscoelasticity is not considered, that is, τε=0 and Ei=βi, and the noise is taken as a white noise, that is, S2ω1=S2ω2=Sω±=S(ω∓), 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(43)3k112S2ω1+3k222S2ω2+4k2Sω±p+25k112S2ω1+5k222S2ω2+12k2Sω±-32E1+E2=0.It is important to note one more time that if the viscoelasticity and coupling (k=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(44)λ0p=limp→0Λ0pp=1645k112S2ω1+5k222S2ω2+12k2Sω±-32E1+E2.The almost-sure stability region is found to be(45)5k112S2ω1+5k222S2ω2+12k2Sω±-32E1+E2<0.From (43) and (45), the moment stability and the almost-sure stability boundaries are both a straight line in the coordinates with E1 and E2.WhenK=1, the eigenfunction is Tφ=C0+C1cos⁡2φ; the moment Lyapunov exponent can be determined from(46)a00-Λ1a01a10a11-Λ1=0.This is a quadratic equation, which is (Λ1)2+e11Λ1+e01=0, where the analytical expression for moment Lyapunov exponent can be obtained. The general Lyapunov exponent is given by(47)λ1=-limp→0e0Kpe1K=YZ,(48)Y=-5k124+k214+74k122k212S+2+-36k12k21k122+k212S--46k122+30k212k112S2ω1-30k122+46k212k222S2ω2+288E1+160E2k122+160E1+288E2k212S++12k122k212S-2+4k112k12k21S2ω1+4k222k12k21S2ω2-64E1+E2k12k21S--13k114S2ω12-13k224S2ω22+6k112k222S2ω2+160E1-96E2k112S2ω1-96E1-160E2k222S2ω2-512E1-E22,Z=642k12k21S--k112S2ω1-k222S2ω2-7k122+k212S+.It is noted that only those values of the excitation spectrum at frequencies2ω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 S(ω), where ω0-Δω0/2<ω<ω0+Δω0/2 and Δω0≪ω0 [18]. Note that Λ1+p and Λ1-p are for symmetric coupling and skew-symmetric coupling, respectively.Case 1 (coupled systems under band-limited noise excitation close toS(2ωi),i=1,2).
(49) Λ 1 + p = Λ 1 - p = 13 512 p 2 + 21 256 p - 1 64 K i i 2 S 2 ω i - p 2 E 1 + E 2 + 1 512 64 + 992 p + 1524 p 2 + 772 p 3 + 129 p 4 K i i 4 S 2 ω i 2 - 4096 p p + 2 2 E 1 - E 2 2 K i i 2 S 2 ω i + 32768 p p + 2 E 1 - E 2 2 1 / 2 .
The corresponding largest Lyapunov exponents are(50)λ1+p=λ1-p=16413k114S2ωi2-k112S2ωi160E1-96E2+512E1-E22k112S2ωi.Case 2 (coupled systems under band-limited noise excitation close toS(ω+)).
(51) Λ 1 + p = 3 128 p 2 + 11 64 p - 3 16 k 2 S ω + - p 2 E 1 + E 2 + 1 128 576 + 96 p + 52 p 2 + 4 p 3 + p 4 k 4 S ω + 2 + 2048 p p + 2 E 1 - E 2 2 1 / 2 , (52) Λ 1 - p = - 1 4 k 2 S ω + - p 2 E 1 + E 2 + 1 4 k 4 S ω + 2 + 2 p p + 2 E 1 - E 2 2 1 / 2 .
The corresponding largest Lyapunov exponents are(53)λ1+p=-1489k4Sω+2-24k2Sω+E1+E2+32E1-E22k2Sω+,λ1-p=-12k2Sω+E1+E2-E1+E22k2Sω+.Case 3 (coupled systems under band-limited noise excitation close toS(ω-)).
(54) Λ 1 + p = - 1 4 k 2 S ω - - p 2 E 1 + E 2 + 1 4 k 4 S ω - 2 + 2 p p + 2 E 1 - E 2 2 1 / 2 , Λ 1 - p = 3 128 p 2 + 11 64 p - 3 16 k 2 S ω - - p 2 E 1 + E 2 + 1 128 576 + 96 p + 52 p 2 + 4 p 3 + p 4 k 4 S ω - 2 + 2048 p p + 2 2 E 1 - E 2 2 1 / 2 .
The corresponding largest Lyapunov exponents are(55)λ1+p=-12k2Sω-E1+E2-E1-E22k2Sω-,λ1-p=-1489k4Sω-2-24k2Sω-E1+E2+32E1-E22k2Sω-.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 p. However, (52) for skew symmetric coupling has no stability index, which is a consequence of a linear function of p. Quantitatively, the second-order expansion K=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, Ei(i=1,2) 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 k12=k21=k, S±=S(ω+)±S(ω-), S2ω1=S2ω2=0, and Sω+=0 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.WhenK=2, the eigenfunction is Tφ=C0+C1cos⁡2φ+C2cos⁡4φ; the moment Lyapunov exponent can be solved from(56)a00-Λ2a01a02a10a11-Λ2a12a20a21a22-Λ2=0.This is a cubic equation, which is (Λ2)3+e22(Λ2)2+e12Λ2+e02=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 K⩾3, 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Λqtδqt=0. Here we consider the case for white noise. When the order of Fourier expansion K=0, from (42), the stability index is given by(57)δqt=25k112S2ω1+5k222S2ω2+12k2Sω±-32E1+E23k112S2ω1+3k222S2ω2+4k2Sω±.WhenK=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 E1, the larger the stability index and then the more stable the system.Figure 3
Stability index for system under real noise.


## 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. Denoting(58)x1t=q1t,x2t=q˙1t,x3t=q2t,x4t=q˙2t,x5t=∫0tγe-ηt-sq1sds,x6t=∫0tγe-ηt-sq2sds,x7t=ξt,dξt=-αξtdt+σdWt,Equations (22a) and (22b) can be converted to the system of first-order Stratonovich stochastic differential equations:(59)dx=Axdt+BdW,where x={x1,x2,x3,x4,x5,x6,x7}T, and(60)A=0100000-ω12-2εβ100-εω120-ε1/2ω1k11x1+k12x30001000-ω22-2εβ2000-εω22-ε1/2ω2k21x1+k22x3γ000-η0000γ00-η0000000-α,B=000000σ.The associated Itô stochastic differential equations are the same as the Stratonovich stochastic differential equations. Thus, the discretized equations using explicit Euler scheme are(61)x1k+1=x1k+x2k·Δt,x2k+1=x2k+-ω12x1k-2εβ1x2k+εω12x5k-ε1/2ω1k11x1k+k12x3kx7kΔt,x3k+1=x3k+x4k·Δt,x4k+1=x4k+-ω22x3k-2εβ2x4k+εω22x6k-ε1/2ω2k21x1k+k22x3kx7kΔt,x5k+1=x5k+γx1k-ηx5kΔt,x6k+1=x6k+γx3k-ηx6kΔt,x7k+1=x7k+-αx7k·Δt+σ·ΔWk,with Δt being the time step and k denoting the kth iteration.


## 8. Application

As an application, the flexural-torsional stability of a simply supported, uniform, narrow, rectangular, viscoelastic beam of lengthL subjected to a stochastically varying concentrated load P(t) acting at the center of the beam cross section is considered. Both nonfollower and follower loading cases are studied.For the elastic beam under dynamic loading as shown in Figure4, the flexural and torsional equations of motion are given by [1](62)EIy∂4u∂z4+∂2Mxψ∂z2-∂2My∂z2+m∂2u∂t2+Du∂u∂t=0,-GJ∂2ψ∂z2+∂Mx∂z∂u∂t+Mx∂2u∂z2+∂Mz∂z+mr2∂2ψ∂t2+Dψ∂ψ∂t=0,with boundary conditions(63)u0,t=uL,t=∂2u0,t∂z2=∂2uL,t∂z2=0,ψ0,t=ψL,t=0,where Mx, My, and Mz are the bending moments in the x, y, and z directions, respectively. u(t) is the lateral deflection and ψ(t) is the angle of twist. E is Young’s modulus of elasticity and G is the shear modulus. Iy is the moment of inertia in y direction. J is Saint-Venant’s torsional constant. EIy is the flexural rigidity, GJ is the torsional rigidity, m is the mass per unit length of the beam, r is the polar radius of gyration of the cross section, Dψ and Du are the viscous damping coefficients in the ψ and u directions, respectively, t is the time, and z is the axial coordinate.Figure 4
Flexural-torsional vibration of a rectangular beam.The governing equation for a viscoelastic deep beam can be obtained from (62) as(64)EIy1-H∂4u∂z4+∂2Mxψ∂z2-∂2My∂z2+m∂2u∂t2+Du∂u∂t=0,-GJ1-H∂2ψ∂z2+∂Mx∂z∂u∂z+Mx∂2u∂z2+∂Mz∂z+mr2∂2ψ∂t2+Dψ∂ψ∂t=0,where H is the relaxation operator defined in (7). These equations are stochastic partial differential equations and are difficult to analyze. An approximate solution can be obtained by using the Galerkin method [19]. Seek solutions of the form(65)uz,t=q1sin⁡πzL,ψz,t=q2sin⁡πzL,in which q1=um=u(1/2L,t) and q2=ψm=ψ(1/2L,t). Substituting into (64), multiplying by sin(πz/L), and integrating with respect to z from 0 to L yield [1](66)Q¨1+2εβ1Q˙1+ω121-εHQ1+ε1/2ω1k12ξtQ2=0,Q¨2+2εβ2Q˙2+ω221-εHQ2+ε1/2ω2k21ξtQ1=0,where(67)K=ω212-π2ω14+π2,β1=Du2m,β2=Dψ2mr2,ω12=πL4EIym,ω22=πL2GJmr2,ξt=PtPcr,Pcr=4mrLω12-ω2212-π24+π21/2,k12=-k21=ω12-ω222ω1ω2=kF,where β1 and β2 are reduced viscous damping coefficients; Pcr is the critical force for the simply supported narrow rectangular beam.For a nonfollower force case, the only difference is thatMy=0 in (64). The equations of motion are of the same form of (66), but the parameters are different:(68)K=-ω2ω1,2β1=Dum,2β2=Dψmr2,ω12=πL4EIym,ω22=πL2GJmr2,ξt=PtPcr,Pcr=8mrLω1ω24+π2,k12=k21=ω1ω2=kN.It is seen that (66) has the same form of (22a) and (22b), except that k11=k22=0. By introducing the polar transformation and using the method of stochastic averaging, (66) can be approximated in amplitude by Itô stochastic differential equations in (23), where the drift and diffusion terms are given by(69)mia=ai-βi-ωiMtIisc+18kijkjiS-+116aj2aikij2S+,biia=18kij2aj2S+,bija=18aiajkijkjiS-,S±=Sω1+ω2±Sω1-ω2.Substituting (69) and a1=rcos⁡φ and a2=rsin⁡φ into (28) yields an eigenvalue problem, from which moment Lyapunov exponents can be determined by solving (35).The almost-sure andpth 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 S(ω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 higher-order moment stability is more conservative than the lower-order moment stability.Figure 5
Almost-sure andpth moment stability boundaries under real noise.Figure 6
Critical excitation and pseudodamping under real noise.Figures7 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 p>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 pth moment (for p>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.Figure 7
Effect ofσ on the moment Lyapunov exponents.Figure 8
Effect ofα on the moment Lyapunov exponents.


## 9. Conclusions

The stochastic stability of coupled viscoelastic systems described by stochastic integrodifferential equations of2DOF 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.



---

*Source: 4725148-2018-06-28.xml*