Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation

ARTICLE IN PRESS

JID: FI

[m1+;October 11, 2017;9:14]

Available online at www.sciencedirect.com

Journal of the Franklin Institute 000 (2017) 1–29 www.elsevier.com/locate/jfranklin

Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation Jian Deng Department of Civil Engineering, Lakehead University, Thunder Bay, Ontario P7B 5E1, Canada Received 24 November 2016; received in revised form 1 September 2017; accepted 20 September 2017 Available online xxx

Abstract The moment Lyapunov exponents and stochastic stability of a single-degree-of-freedom (SDOF) fractional viscoelastic system under bounded noise excitation are studied by using the method of higher-order stochastic averaging. A realistic example of such a system is the transverse vibration of a viscoelastic column under the excitation of stochastic axial compressive load. The excitation is modeled as a bounded noise, which is a realistic model of stochastic fluctuation in engineering applications. The viscoelastic material is assumed to follow a fractional Kelvin–Voigt constitutive relation. The method of higher-order stochastic averaging is used to approximate the fractional stochastic differential equation of motion, and then moment Lyapunov exponents are determined for the system with small damping and weak random fluctuation. The approximate results are confirmed by Monte-Carlo simulations. It is found that convergence of moment Lyapunov exponents depends on the width of power spectral density of the bounded noise process. For this viscoelastic structure, second-order averaging analysis is adequate for stability analysis. The effects of various parameters on the stochastic stability of the system are discussed and possible explanations are explored. © 2017 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.

1. Introduction Viscoelasticity is referred to as a mechanical behavior that combines liquidlike and solidlike characteristics. Viscoelasticity has been observed in many engineering materials, such E-mail address: [email protected]. https://doi.org/10.1016/j.jfranklin.2017.09.019 0016-0032/© 2017 The Franklin Institute. Published by Elsevier Ltd. All rights reserved.

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

2

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

as polymers, concrete, metals at elevated temperatures, and rocks, which have been used extensively in industrial world [1]. Most theories on viscoelasticity are devoted to traditional binary viscoelastic models allowing for only two possible states for its components: pure elasticity and pure viscosity, where Hoek’s spring and Newton’s dash-pot are used to represent the purely elastic and purely viscous behaviors, respectively. Several conventional mechanical models for viscoelasticity were proposed. For example, the Maxwell model is represented by a dash-pot and a spring connected in series, and the Kelvin–Voigt model consists of a dash-pot and a spring in parallel. However, many real-world materials are composed of multi-state components, which have various performance levels between elasticity and viscosity, with various effects on the entire system performance. By using the theory of fractional calculus (Appendix A), another component can be formulated with property between spring and dash-pot, which is introduced in Appendix B. Application of fractional calculus in viscoelasticity [2] and stability analysis of fractional-order systems can be found in theoretical and practical studies [3]. However, loadings from earthquakes, explosion, wind, and ocean waves can be described satisfactorily only in probabilistic domain, which results in the fact that the equation of motion of the viscoelastic system under such excitations is usually governed by stochastic differential equation and the response and stability of the system are difficult to be obtained exactly. The method of stochastic averaging has been widely used to approximately solve stochastic differential equations containing a small parameter. Under certain conditions, the stochastic averaging is able to reduce the dimension of problems and then greatly simplify the solution [4]. The popularity of stochastic averaging can be felt from the overwhelming number of papers in the literature [5,6]. The basic idea of stochastic averaging method is to approximate the original Stratonovich stochastic system by an averaged Itô stochastic system (diffusive Markov process). The Itô system is presumably easier to study, and we can infer properties of the dynamics of the original system by the understanding of the dynamics of the averaged system. It has been shown that the ordinary stochastic averaging is a first order approximation method [4]. Deng et al. [7] investigated stochastic stability of a fractional viscoelastic column under bounded noise excitation using a first-order stochastic averaging method. A systematic and unified approach of second-order stochastic averaging based on the Stratonovich–Khasminskii limit theorem was presented [8]. Huang and Xie found that the second-order averaging results agree better with those estimated by Monte Carlo simulation than the first-order averaging [9]. More recently, moment Lyapunov exponents of Hill equation were studied by the method of stochastic averaging, both the first-order and the second-order [10]. The aim of the present paper is to study the stability of fractional viscoelastic systems under bounded noise excitation by using the method of higher-order stochastic averaging. The paper is an extention of [7] from stochastic averaging of first-order to higher-order. It also extends [9–11] from white noise to bounded noise. The paper is structured as follows. Section 2 derives the stochastic equation of motion and uses the method of higher order stochastic averaging to obtain Itô stochastic differential equations. Results and discussions are given in Section 3. Conclusions are presented in Section 4.

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

ARTICLE IN PRESS

JID: FI

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

3

2. Stochastic averaging 2.1. Formulation Consider the stability problem of a column of uniform cross section under dynamic axial compressive load F(t). The study of such as vibrating column would yield interesting insights on the behavior of a much wider class of structures [12]. The equation of motion is derived as [4] ∂ 2 M(x, t ) ∂ 2 v(x, t ) ∂v(x, t ) ∂ 2 v(x, t ) = ρA + β0 + F (t ) , (2.1) 2 2 ∂x ∂t ∂t ∂x 2 where L is the length of the column, ρ is the mass density per unit volume of the column, A is the cross-sectional area, v(x, t ) is the transverse displacement of the central axis, β 0 is the damping constant. The moment M(x, t) at the cross-section x and the geometry relation are  ∂ 2 v(x, t ) M(x, t ) = σ (x , t )z dA, ε(x , t ) = − z. (2.2) ∂x 2 A The viscoelastic material is supposed to follow fractional constitutive model in Eq. (B.1), μ RL μ which can be recast as σ (t ) = [E + η RL 0 Dt ]ε(t ). Here 0 Dt is an operator for Riemann– Liouville(RL) fractional derivative of order μ, which is defined in Eq. (A.8). Substituting into Eq. (2.2) and then to Eq. (2.1) yields  4  ∂ v(x, t ) ∂ 2 v(x, t ) ∂v(x, t ) ∂ 4 v(x, t ) ∂ 2 v(x, t ) RLDμ + F (t ) ρA + β + E I + I η = 0. (2.3) 0 0 t 2 4 4 ∂t ∂t ∂x ∂x ∂x 2 If the column is simply supported, the solution of Eq. (2.3) can be written in the following form: πx v(x, t ) = q(t ) sin . (2.4) L Substitute Eq. (2.4) into Eq. (2.3), suppose the damping, viscoelastic effect, and the amplitude of load are all small, and introduce a small parameter 0 < ε  1. If the function F(t) is taken to be a stochastic process, then the equation of motion is given by  μ q¨ (t ) + 2εβ q˙ (t ) + ω2 1 − εξ (t ) + ετε RL (2.5) 0 Dt q(t ) = 0, where β=

 π 2 β0 E I  π 4 F (t ) η , ω2 = , P = EI , ξ (t ) = , τε = . 2ρA ρA L L P E

(2.6)

Mathematically, random excitations can be described as stochastic processes, either wideband noise processes or narrow-band noise processes. A typical wide-band noise process is White noise which is not a physically realizable process. A bounded noise η(t) is a realistic and versatile narrow-band stochastic fluctuation in engineering applications and is normally represented as   η(t ) = ζ cos νt + ε 1/2 σW (t ) + θ , (2.7) where ζ is the noise amplitude, σ is the noise intensity, W(t) is the standard Wiener process, and θ is a random variable uniformly distributed in the interval [0, 2π ]. It overcomes the Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

4

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

shortcoming of Ornstein–Uhlenbeck process, which is not bounded. The inclusion of the phase angle θ makes the bounded noise η(t) a stationary process. Eq. (2.7) may be written as η(t ) = ζ cos Z (t ), dZ (t ) = νt + σ dW (t ),

(2.8)

where the initial condition of Z(t) is Z (0) = θ . This process is of course bounded between −ζ and +ζ for all time t and hence is a bounded stochastic process. The auto-correlation function of η(t) is given by

σ2 1 R(τ ) = E[η(t )η(t + τ )] = ζ 2 cos ντ exp − |τ | , (2.9) 2 2 and the spectral density function of η(t) is

1  +∞ ζ 2 σ 2 ω2 + ν 2 + σ 4 4  . S(ω) = R(τ )e−iωτ dτ =  (2.10) 1 1 4 −∞ 2 4 2 2 (ω + ν) + σ (ω − ν) + σ 4 4 When the noise intensity σ is small, the bounded noise can be used to model a process about frequency ν. Hence, the bounded noise is a very good realistic model of narrow-band noise and is widely used in many engineering applications. In the limit as σ approaches zero, the bounded noise reduces to a deterministic sinusoidal function. On the other hand, in the limit as σ approaches infinite, the bounded noise becomes a white noise of constant spectral density. However, since the mean-square value is fixed at 21 , this constant spectral density level reduces to zero in the limit. The bounded noise process has been widely applied in engineering [4], such as wind turbulence and earthquake load. The moment stability of a stochastic dynamical system with state vector x(t) is determined by moment Lyapunov exponents, which are defined by Xie [4]   1 (p) = lim log E x(t ) p , (2.11) t →∞ t where E[·] denotes the expected value and  ·  represents a suitable vector norm. The pth moments are asymptotically stable if (p) < 0. Moreover, (p) is a convex function of p and  (p) is equal to the largest Lyapunov exponent λ, which is defined by 1 λ = lim log x(t ), (2.12) t →∞ t and describes the almost-sure or sample stability of the system. As shown in Fig. 1, the almost-sure stability or instability cannot always guarantee the moment stability, which can be explained by the theory of large deviation [4]. Thus, it is important to obtain the moment Lyapunov exponents of stochastic systems so that the complete properties of dynamic stability can be described. However, it is difficult to obtain the exact moment Lyapunov exponents of the stochastic system concerned in Eq. (2.5). In this paper, the method of higher-order stochastic averaging will be used to obtain the differential equations governing the pth moment. The moment stability of viscoelastic system (2.5) can then be determined by solving the averaged equations. 2.2. First-order stochastic averaging Because higher-order stochastic averaging is based on lower-order, for completeness, this subsection introduces first-order stochastic averaging, which was used in [7]. To apply the Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

ARTICLE IN PRESS

JID: FI

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

5

Fig. 1. Almost sure stability and moment stability.

averaging method, the transformation is made to the amplitude and phase variable a and ϕ by means of the relations ν q(t ) = a(t ) cos (t ), q˙ (t ) = −ω a(t ) sin (t ), (t ) = t + ϕ(t ). (2.13) 2 Substituting Eq. (2.13) into Eq. (2.5) yields a˙ cos (t ) − aϕ˙ sin (t ) = −εa sin (t ), a˙ sin (t ) + aϕ˙ cos (t ) = εa cos (t ) − 2εaβ sin (t ) μ +εωξ (t )a cos (t ) + εωτε RL 0 Dt q(s),

(2.14)

where ε = ω − 21 ν. Solving Eq. (2.14) and assuming P = a p , the pth norm of system (2.5), one can obtain   1 P˙ (t ) = ε pP −2β sin2 (t ) + ωξ (t ) sin 2(t ) + ωτε U ss , 2   ϕ˙ (t ) = ε  − β sin 2(t ) + ω cos2 (t )ξ (t ) + ωτε U cs , (2.15) where

  ω sin (t )  t P (s) 1/p sin (s) ds, (1 − μ) 0 P (t ) (t − s)μ   ω cos (t )  t P (s) 1/p sin (s) =− ds, (1 − μ) 0 P (t ) (t − s)μ

U ss = − U cs

and for ease of presentation, the fractional derivative in Eq. (A.8) is recast as  t   g(t ) f  (τ ) RLDμ f (τ )g(t ) = dτ. 0 t (1 − μ) 0 (t − τ )μ

(2.16)

The bounded noise is recast as, by assuming that the magnitude is small and then introducing a small parameter ε1/2 , ξ (t ) = ζ cos [νt + ψ (t )],

ψ (t ) = ε 1/2 ϑW (t ) + θ .

Substituting ξ (t) in Eq. (2.15) yields Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

6

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29



ζω P˙ (t ) = ε pP −β[1 − cos 2] + cos[νt + ψ (t )] sin 2 + ωτε U ss , 2   ϕ˙ (t ) = ε  − β sin 2 + ωζ cos[νt + ψ (t )] cos2  + ωτε U cs , ψ˙ (t ) = ε 1/2 ϑ W˙ (t ).

(2.17)

Eqs. (2.17) are exactly equivalent to Eq. (2.5) and still cannot be accurately solved. However, it is observed the presence of a small parameter ε means both P and ϕ change slowly. Therefore, one can expect to obtain reasonably accurate results by averaging the response over one period. This may be done by applying the averaging operator given by  1 τ +T M(·) = lim (·)dτ, τ T →∞ T τ where ( · ) represents the equation of motion. When applying the averaging operator the integration is performed over explicitly appearing τ only. The averaging method can be applied to obtain the averaged equations as follows, without distinction between the averaged and the original non-averaged variables P and ϕ. When applying the averaging operation, P(t) and ϕ(t) are treated as unchanged, i.e., they are replaced by the corresponding averaged variables directly.    ss  1 ˙ , P (t ) = ε pP −β + ζ ω sin (2ϕ − ψ ) + ωτε M U t 4     1 ϕ˙ (t ) = ε  + ζ ω cos(2ϕ − ψ ) + ωτε M U cs , (2.18) t 4 where some averaged identities have been used, M(cos 2) = M(sin 2) = 0, t t   1 M cos[νt + ψ (t )] sin 2 = sin (2ϕ − ψ ), t 2   1 M cos[νt + ψ (t )] cos2  = cos(2ϕ − ψ ). (2.19) t 4 Applying the transformation τ = t − s and changing the order of integration lead to     ω ν  ω ν  , M U cs = H s , M U ss = − H c (2.20) t t 2 2 2 2 where  ∞    ν μ−1 1 ντ μπ c ν = H τ −μ cos dτ = sin , 2 (1 − μ) 0 2 2 2  ν   ν μ−1 ∞ 1 ντ μπ = Hs τ −μ sin dτ = cos . (2.21) 2 (1 − μ) 0 2 2 2 Substituting Eqs. (2.20) into Eq. (2.18) yields   1 P˙ (t ) = ε pP −βˆ + ζ ω sin (2ϕ − ψ ) , 4   ˆ + 1 ζ ω cos(2ϕ − ψ ) , ϕ˙ (t ) = ε  4

(2.22)

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

7

where

ν    1 ˆ =  + 1 ωτε H s ν . ,  βˆ = β + ωτε H c 2 2 2 2

(2.23)

Introducing the transformation  = ϕ − 21 ψ and Itô Lemma results in three Itô stochastic differential equations (SDE) dP (t ) = mP P dt,

ϑ d W (t ), 2 1/2 dψ (t ) = ε ϑ d W (t ),

d(t ) = m dt − ε 1/2

(2.24)

where

    1 ˆ + 1 ζ ω cos(2) . mP = ε p −βˆ + ζ ω sin (2) , m = ε  4 4

(2.25)

The eigenvalue problem governing the pth moment Lyapunov exponent can be formulated from these SDEs. Applying a linear stochastic transformation [13] S = T ()P, P = T −1 ()S,

0 ≤  ≤ π,

the Itô equation for the transformed pth norm process S is given by,

εϑ 2 ε 1/2 ϑ dS = mP T + m T + T P dt − T P d W (t ) 8 2

1 εϑ 2 ε 1/2 ϑ T mP T + m T + = T S dt − S d W (t ), T 8 2 T

(2.26)

(2.27)

where T and T are the first and second derivatives of T with respect to , respectively. For bounded and non-singular transformation T(), both processes P and S are expected to have the same stability behavior, as Eq. (2.26) is a linear transformation. Therefore, T() is chosen so that the drift term of the Itô Eq. (2.27) is independent of , i.e.,

1 εϑ 2 mP T + m T + T = , T 8 and ε 1/2 ϑ T S d W (t ). (2.28) 2 T Comparing the drift terms in Eqs. (2.27) and (2.28), such a transformation T() is given by the following differential equation dS = S dt −

ϑ2 ˆ , T + m1 T + m p1 T = T 8 where   m ˆ + 1 ζ ω cos(2), m p1 = m p = p −βˆ + 1 ζ ω sin (2) . m1 = = ε 4 ε 4

(2.29)

(2.30)

ˆ T() is a periodic function in  of period π , and q(t ) (p) = (p) = ε (p) . Eq. (2.29) deˆ fines an eigenvalue problem for a second-order differential operator with (p) being the eigenvalue and T() the associated eigenfunction. Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

ARTICLE IN PRESS

JID: FI

8

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

Taking the expected value of Eq. (2.28) leads to dE[S] = E[S] dt, from which it can be seen that  is the Lyapunov exponent of the transformed pth moment E[S] of Eqs. (2.24) or Eq. (2.5). Since both processes P and S have the same stability behavior,  is the Lyapunov exponent of pth moment E[P]. The remaining task for determining the moment Lyapunov exponents is to solve the eigenvalue problem. 2.3. Determination of moment Lyapunov exponents Since the coefficients in Eq. (2.29) are periodic with period π , it is reasonable to consider a Fourier series expansion of the eigenfunction T() in the form T () = C0 +

K 

(Ck cos 2k + Sk sin 2k).

(2.31)

k=1

To solve Eq. (2.29), substituting the expansion into eigenvalue problem (2.29) and equating the coefficients of like trigonometric terms sin 2k and cos 2k, k = 0, 1, . . . , K yield a set of homogeneous linear algebraic equations with infinitely many equations for the unknown coefficients C0 , Ck and Sk , Constant : cos 2 :

ˆ C0 + b0 S1 = 0, − (c0 + ) ˆ C1 + d1 S1 + b1 S2 = 0, − (c1 + )

ˆ S1 = 0, sin 2 : 2a1C0 + d1C1 + b1C2 + (c1 + ) ˆ Ck = 0, cos 2k : ak Sk−1 + dk Sk + bk Sk+1 − (ck + ) ˆ Sk = 0, sin 2k : ak Ck−1 + dk Ck + bk Ck+1 + (ck + ) ˆ CK = 0, cos 2K  : aK SK−1 + dK SK − (cK + ) ˆ SK = 0, sin 2K  : aK CK−1 + dK CK + (cK + )

(2.32)

where k = 2, 3, . . . , K − 1, 1 1 an = ζ ω(2n − 2 − p), bn = ζ ω(2n + 2 + p), 8 8 1 2 2 ˆ dn = 2n, ˆ cn = σ n + pβ, n = 0, 1, 2, . . . , K. 2 To have a non-trivial solution of the C0 , Ck and Sk , it is required that the determinant of the coefficient matrix of Eqs. (2.32) equal zero to yield a polynomial equation of degree ˆ (K ) (p), 2K + 1 for   2K+1  2K (K ) (K ) ˆ (K ) ˆ (K ) ˆ (K ) + e(K ) = 0, e2K+1  + e2K  + · · · + e1(K )  (2.33) 0 ˆ (K ) denotes the approximate moment Lyapunov exponent. Then, one may approximate where  the moment Lyapunov exponent of the system by ˆ (K ) (p). q(t ) (p) ≈ ε 

(2.34)

For K = 1, Eq. (2.33) reduces to a cubic equation with Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

ARTICLE IN PRESS

JID: FI

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

9

ˆ e3(1) = 1, e2(1) = σ 2 + 3 pβ, 4 ˆ 2 + σ − 1 ζ 2 ω2 p(p + 2) + 3βˆ 2 p2 + 2σ 2 βˆ p, e1(1) = 4 4 32   4 ˆ 2 + σ − 1 ζ 2 ω2 p(p + 2) − 1 ζ 2 ω2 σ 2 p(p + 2), e0(1) = βˆ 3 p3 + σ 2 βˆ 2 p2 + pβˆ 4 4 32 64

and an approximate expression of the moment Lyapunov exponent is given by 2 ˆ (1) (p) = A2 − 2(3e1 − e2 ) − e2 ,  6 3A2 3  1/3 A2 = 12A1 − 108e0 + 36e1 e2 − 8e32 ,  2 1/2 A1 = 81e0 + 12e31 − 54e0 e1 e2 + 12e0 e32 − 3e21 e22 ,

(2.35)

in which the superscript “(1)” in the coefficients e s is dropped for clarity of presentation. When K > 1, there are no analytical solutions for the polynomial Eq. (2.33) and numerical approaches must be applied to solve it. 2.4. Second-order stochastic averaging Rewrite the transformed equation of motion in Eq. (2.17) P˙ (t ) = ε f1 (P, ϕ, t ), ϕ˙ (t ) = ε f2 (P, ϕ, t ), ψ˙ (t ) = ε 1/2 σ W˙ (t ), where

(2.36)



ζω sc f1 (P, ϕ, t ) = pP −β[1 − cos 2(t )] + cos [νt + ψ (t )] sin 2 + ωτε U , 2

f2 (P, ϕ, t ) =  − β sin 2 + ωζ cos [νt + ψ (t )] cos2  + ωτε U cc .

(2.37)

To formulate a second-order stochastic averaging, the near-identity transformation is introduced [8]     P (t ) = P¯ (t ) + ε P1 P¯ , ϕ¯ , t , ϕ(t ) = ϕ¯ (t ) + ε ϕ1 P¯ , ϕ¯ , t , (2.38) where P¯ (t ) and ϕ¯ (t ) stand for nonoscillatory amplitude and phase angle, respectively. Differentiating Eq. (2.38) and equating each result with the corresponding functions from Eq. (2.36) give ⎧ ⎫ ⎧  ⎫   ∂P1 ˙¯ ∂P1 ˙ ∂P1 ⎪ ⎪ ⎪ ⎪ ¯ P + ϕ ¯ + ⎨ ⎬ ⎨ ⎬ ˙ P f + ε P , ϕ ¯ + ε ϕ , t 1 1 1 ∂ ϕ¯ ∂t ∂ P¯ P¯ +ε =ε (2.39)   . ⎪ ⎪ ϕ˙¯ ⎩ ∂ϕ1 P˙¯ + ∂ϕ1 ϕ˙¯ + ∂ϕ1 ⎪ ⎭ ⎩ f2 P¯ + ε P1 , ϕ¯ + ε ϕ1 , t ⎪ ⎭ ∂ ϕ¯ ∂t ∂ P¯ Substituting the Taylor expansion   ∂ fj ∂ fj f j P¯ + ε P1 , ϕ¯ + ε ϕ1 , t = f j (P¯ , ϕ¯ , t ) + ε P1 + ε ϕ1 + o(ε), j = 1, 2, ¯ ∂ ϕ¯ ∂P and the inverse matrix ⎡ ⎤ ⎡ ⎤ 1 1 1 1 1 + ε ∂P ε ∂P 1 − ε ∂P −ε ∂P ∂ ϕ¯ ∂ ϕ¯ ∂ P¯ ∂ P¯ ⎦, A−1 = ⎣ ⎦ + o(ε) A=⎣ ∂ϕ1 ∂ϕ1 ∂ϕ1 1 ε ∂ P¯ 1 + ε ∂ ϕ¯ −ε ∂ϕ 1 − ε ∂ ϕ¯ ∂ P¯

(2.40)

(2.41)

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

10

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

into Eq. (2.39) yields     P¯˙ = ε u1 + ε 2 u2 + o ε 2 , ϕ˙¯ = ε v1 + ε 2 v2 + o ε 2 ,

(2.42)

where ∂ f1 ∂P1 , u2 = P1 + ∂t ∂ P¯ ∂ f2 ∂ϕ1 v1 = f 2 − , v2 = P1 + ∂t ∂ P¯

u1 = f 1 −

∂ f1 ϕ1 − ∂ ϕ¯ ∂ f2 ϕ1 − ∂ ϕ¯

∂P1 ∂P1 ∂ϕ1 ∂P1 ∂P1 ∂P1 f1 − f2 + + , ¯ ¯ ∂ ϕ ¯ ∂t ∂t ∂ ϕ¯ ∂P ∂P ∂P1 ∂ϕ1 ∂ϕ1 ∂ϕ1 ∂ϕ1 ∂ϕ1 f1 − f2 + + . ∂ ϕ¯ ∂t ∂ P¯ ∂t ∂ ϕ¯ ∂ P¯

(2.43)

Now efforts are being made to obtain u1 , u2 , v1 and v2 , and their averaged values. Firstly, to obtain u1 and v1 , one can assume that ∂P1 = f 1 − M{ f 1 } , t ∂t

∂ϕ1 = f 2 − M{ f 2 } , t ∂t

(2.44)

then the terms f1 − (∂ P1 /∂ t ) and f2 − (∂ ϕ1 /∂ t ) contain only nonoscillatory terms. It also suggests that P1 and ϕ 1 are oscillatory functions. Substituting Eqs. (2.20) into the averaged Eq. (2.37) yields   ζω u1 = M{ f1 } = pP −βˆ + sin (2ϕ − ψ ) , t 4 ζ ω ˆ + v 1 = M{ f 2 } =  cos(2ϕ − ψ ), (2.45) t 4 where P and ϕ are treated as a constant under the averaging operation. Substituting Eqs. (2.37) and (2.45) into Eq. (2.44) yields ν  ∂P1 1  sin (νt + 2ϕ¯ ) = pP¯ 4β cos(νt + 2ϕ¯ ) + ωζ sin (ψ − 2ϕ¯ ) + 2ωτε H s ∂t 4 2 ν   cos(νt + 2ϕ¯ ) + 2ωζ cos(νt + ψ ) sin (νt + 2ϕ¯ ) , + 2ωτε H c 2 ν  ν  ∂ϕ1 1 1 cos (νt + 2ϕ¯ ) − ωτε H c sin (νt + 2ϕ¯ ) = −β sin (νt + 2ϕ¯ ) + ωτε H s ∂t 2 2 2 2 1 1 1 − ωζ cos(ψ − 2ϕ¯ ) + ζ ω cos(νt + ψ ) + ωζ cos(νt + ψ ) cos(νt + 2ϕ¯ ). (2.46) 4 2 2 During the derivation of Eq. (2.46), the following equations must be used,      t sin ν t + ϕ¯ sin ν s + ϕ¯ ν 2 2 U sc = − ds 2(1 − μ) 0 (t − s)μ

  νt ντ   ∞ sin νt2 + ϕ¯ sin − + ϕ¯ ν 2 2 =− dτ 2(1 − μ) 0 τμ



νt νt ντ   ∞ sin + ϕ¯ sin − + ϕ¯ 2 2 2 − dτ τμ t

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

ARTICLE IN PRESS

JID: FI

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

  ν νt cos − τ + 2 ϕ ¯ ν 1 ν  2 = dτ − H c μ 4(1 − μ) 0 τ 2 2 ν   ν  1 c ν  1  + cos(νt + 2ϕ¯ )H c + sin (νt + 2ϕ¯ )H s . =− H 2 2 2 2 2 

11

∞

(2.47)

The third equation exists when t → ∞. Similarly, one may obtain 1 s ν  − H 2 2 1 ν  + U ss = H s 2 2   ν 1 + U cs = H c 2 2

U cc =

ν   ν  1 − cos(νt + 2ϕ¯ )H s sin (νt + 2ϕ¯ )H c , 2 2 2 ν   ν  1 − cos(νt + 2ϕ¯ )H s sin (νt + 2ϕ¯ )H c , 2 2 2      ν ν 1 + sin (νt + 2ϕ¯ )H s cos(νt + 2ϕ¯ )H c . 2 2 2

(2.48)

Solving these ordinary first-order differential equations in Eqs. (2.46) results in  ν  ν  pP¯  4 2β + ωτε H s cos(νt + 2ϕ¯ ) P1 = sin (νt + 2ϕ¯ ) − 4ωτε H s 8ν 2 2

−ωζ cos(2νt + ψ + 2ϕ¯ ) ,   ν  1 4 2β + ωτε H c ϕ1 = cos(νt + 2ϕ¯ ) + 4ωζ sin (νt + ψ ) 8ν 2

ν  sin (νt + 2ϕ¯ ) + ωζ sin (2νt + ψ + 2ϕ¯ ) . +4ωτε H s 2

(2.49)

Secondly, to obtain u2 and v2 , the following derivatives of f1 and f2 in Eqs. (2.37) are required ∂ f1 ∂ P¯ ∂ f1 ∂ ϕ¯ ∂ f2 ∂ P¯ ∂ f2 ∂ ϕ¯

  p ωζ cos(νt + ψ ) sin 2 + ωτε pU sc + P¯ V sc , 2   = pP¯ −2β sin 2 + ωζ cos(νt + ψ )cos2 + ωτε (U cc − U ss ) , = −2β p sin2  +

=

ωτε V cc , p

= −2β cos 2 − ωζ cos(νt + ψ ) sin 2 − ωτε (U cs + U sc ),

(2.50)

where 

P¯ (s) 1/p cos(ωt + ϕ¯ ) cos(ωs + ϕ¯ ) , P¯ (t ) 

P¯ (s) 1/p μ ss RL U = 0 Dt sin (ωt + ϕ¯ ) sin (ωs + ϕ¯ ) , P¯ (t )  1/p

P¯ (s) μ cs RL U = 0 Dt cos(ωt + ϕ¯ ) sin (ωs + ϕ¯ ) , P¯ (t )

U cc =

RLDμ 0 t

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

ARTICLE IN PRESS

JID: FI

12

U

sc

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

=

RLDμ 0 t



1/p

sin (ωt + ϕ¯ ) cos(ωs + ϕ¯ ) ,

1 1  V = − cos(ωt + ϕ¯ ) cos(ωs + ϕ¯ ) , P¯ (s) P¯ (t )  1/p 

P¯ (s) 1 1  μ sc RL V = 0 Dt − sin (ωt + ϕ¯ ) cos(ωs + ϕ¯ ) , P¯ (t ) P¯ (s) P¯ (t )  1/p 

P¯ (s) 1 1  μ ss RL V = 0 Dt − sin (ωt + ϕ¯ ) sin (ωs + ϕ¯ ) , P¯ (t ) P¯ (s) P¯ (t ) cc

RLDμ 0 t

μ V cc = RL 0 Dt



P¯ (s) P¯ (t )



and ∂U sc V sc = , p ∂ P¯

P¯ (s) P¯ (t )

P¯ (s) P¯ (t )

1/p 

1/p 

1 1  − cos(ωt + ϕ¯ ) cos(ωs + ϕ¯ ) , P¯ (s) P¯ (t )

∂U cc V cc = , p ∂ P¯

∂U cc = −U sc − U cs , ∂ ϕ¯

∂U sc = U cc − U ss . ∂ ϕ¯

(2.51)

(2.52)

Substituting P1 and ϕ 1 into Eq. (2.43) leads to the expressions of u2 and v2 , but they are too tedious to list here. Our main interest is the averaged values. Performing averaging on them gives:

∂ f1 ∂ f1 ∂P1 ∂P1 M{ u 2 } = M P1 + ϕ1 − M{ f 1 } − M{ f 2 } t t ∂ ϕ¯ ∂ ϕ¯ t ∂ P¯ ∂ P¯ t

     ¯ pPωζ c ν s ν sin (ψ − 2ϕ¯ ) =− 2β + ωτε H cos(ψ − 2ϕ¯ ) − ωτε H 4ν 2 2  pP¯ ζ ω  ˆ =− β cos(ψ − 2ϕ¯ ) + αˆ sin (ψ − 2ϕ¯ ) , 2ν ∂ f  ∂ f2 ∂ϕ1 ∂ϕ1 2 M{ v 2 } = M P1 + ϕ1 − M{ f 1 } − M{ f 2 } t t ∂ ϕ¯ ∂ ϕ¯ t ∂ P¯ ∂ P¯ t    ν   ν 1 8ζ ω2 τε H s cos(ψ − 2ϕ¯ ) + 8ωζ 2β + ωτε H c =− sin (ψ − 2ϕ¯ ) 32ν 2 2

  2   2   2 2 s ν c ν c ν 2 + 32β + ω ζ + 8 ωτε H + 8 ωτε H + 32ωβτε H 2 2 2  1 ˆ sin (ψ − 2ϕ¯ ) + δˆ , = ζ ωαˆ cos(ψ − 2ϕ¯ ) − βωζ (2.53) 2ν where ν  1 , αˆ = − ωτε H s 2 2

  2 1   2   c ν c ν 2 ˆδ = − 1 ω2 ζ 2 + 1 ωτε H s ν + 2β . + ωτε H + 2ωβτε H 16 2 2 2 2 2 Applying the transformation  = ϕ − 21 ψ and Itô Lemma results in three Itô stochastic differential equations Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

  dP (t ) = ε m p1 + ε 2 m p2 P dt,   σ d(t ) = ε m1 + ε 2 m2 dt − ε 1/2 d W (t ), 2 dψ (t ) = ε 1/2 σ d W (t ),

13

(2.54)

where

  1   ζω M u1 = p −βˆ + sin (2) , P t 4   1   pζ ω ˆ m p2 = M u2 = − β cos(2) − αˆ sin (2) , P t 2ν   ˆ + ζ ω cos(2), m1 = M v1 =  t 4    1 ˆ sin (2) + δˆ . m2 = M v2 = ζ ωαˆ cos(2) + βωζ (2.55) t 2ν It is observed that Eqs. (2.54) are also coupled. Introducing a linear stochastic transformation m p1 =

S = T ()P, P = T −1 ()S,

0 ≤  ≤ π,

the Itô equation for the transformed pth norm process S is given by, again using Itô Lemma in the vector case       εσ 2 ε 1/2 σ 2 2 dS = ε m p1 + ε m p2 T + ε m1 + ε m2 T + T P dt − T P d W (t ) 8 2      1  εσ 2 ε 1/2 σ T ε m p1 + ε 2 m p2 T + ε m1 + ε 2 m2 T + = T S dt − S d W (t ). T 8 2 T (2.56) The eigenvalue problem governing the pth moment Lyapunov exponent is   σ2 ˆ , T + (m1 + εm2 )T + m p1 + εm p2 T = T 8

(2.57)

ˆ where T() is a periodic function in  of period π , and q(t ) (p) = (p) = ε (p) . Substituting the expansion in Eq. (2.31) into eigenvalue problem (2.57) and equating the coefficients of like trigonometric terms sin 2k and cos 2k, k = 0, 1, . . . yield a set of homogeneous linear algebraic equations with infinitely many equations for the unknown coefficients C0 , Ck and Sk , k = 1, 2, . . . Constant terms:     ˆ 1 − αS ˆ C0 − b0 S1 + ε 2b0 βC c0 +  ˆ 1 = 0, υ Coefficients for cos 2:     ˆ 0 + 2βb ˆ 1C2 = 0, ˆ C1 − d1 S1 − b1 S2 + ε 1 −e1 S1 − 2αb c1 +  ˆ 1 S2 − 4a1 βC υ Coefficients for sin 2:     ˆ 0 + 2βb ˆ 1 S2 = 0, ˆ S1 + ε 1 e1C1 + 2αb 2a1C0 + d1C1 + b1C2 + c1 +  ˆ 1C2 + 4a1 βC υ Coefficients for cos 2k: Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

ARTICLE IN PRESS

JID: FI

14

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

  ˆ Ck + dk Sk + bk Sk+1 ak Sk−1 − ck +   1 ˆ k+1 = 0, + ε 2ak αS ˆ k−1 + ek Sk + 2bk αS ˆ k+1 + 2ak αC ˆ k−1 − 2bk βC υ k = 2, 3, . . . , K − 1, Coefficients for sin 2k:

  ˆ Sk ak Ck−1 + dk Ck + bk Ck+1 + ck +   1 ˆ k+1 = 0, + ε 2ak αC ˆ k−1 + ek Ck + 2bk αC ˆ k+1 − 2ak αS ˆ k−1 + 2bk βS υ k = 2, 3, . . . , K − 1, Coefficients for cos 2K:     ˆ K−1 = 0, ˆ CK + dK SK + ε 1 2ak αS aK SK−1 − cK +  ˆ K−1 + ek SK + 2ak βC υ Coefficients for sin 2K:     ˆ K−1 = 0, ˆ SK + ε 1 2ak αC aK CK−1 + dK CK + cK +  ˆ K−1 + ek CK + 2ak βS υ where 1 1 ak = ζ ω(2k − 2 − p), bk = ζ ω(2k + 2 + p), 8 8 1 2 2 ˆ dk = 2k , ˆ k = 0, 1, 2, . . . . ˆ ck = σ k + pβ, ek = k δ, 2

(2.58)

(2.59)

As in the first-order stochastic averaging, Eqs. (2.58) is reduced to a set of 2K + 1 homogeneous linear equations for C0 , Ck and Sk , k = 1, 2, . . . , K . The determinant is set to zero to ˆ K (p). Solving the polynomial equation yield a polynomial equation of degree 2K + 1 for    for the largest root yields the approximate moment Lyapunov exponents. The terms H s 21 ν   and H c 21 ν depend upon the material’s constitutive model. 2.5. Higher-Order stochastic averaging Following the procedure of Section 2.4, higher-order stochastic averaging analysis can be carried out. Let’s take the three-order stochastic averaging as an example of higher-order method. Introducing a near-identity transformation as in Eq. (2.38), but up to order ε2 P (t ) = P¯ (t ) + ε P1 + ε 2 P2 , ϕ(t ) = ϕ¯ (t ) + ε ϕ1 + ε 2 ϕ2 ,

(2.60)

where P1 = P1 (P¯ , ϕ¯ , t ), P2 = P2 (P¯ , ϕ¯ , t ), ϕ1 = ϕ1 (P¯ , ϕ¯ , t ), and ϕ2 = ϕ2 (P¯ , ϕ¯ , t ) are unknown functions. When P2 = 0 and ϕ2 = 0, it is reduced to second-order averaging. Differentiating Eq. (2.60) and equating each result with the corresponding functions from Eq. (2.36) gives ⎧ ⎫ ⎧ ⎫   ∂P1 ˙¯ ∂P1 ˙ ∂P1 ⎪ ∂P2 ˙¯ ∂P2 ˙ ∂P2 ⎪ ⎪ ⎪ P P + ϕ ¯ + + ϕ ¯ + ⎨ ⎬ ⎨ ˙ ¯ ¯ ∂ ϕ¯ ∂t ∂ ϕ¯ ∂t ⎬ ∂P ∂P P¯ +ε + ε2 ⎪ ⎪ ϕ˙¯ ⎩ ∂ϕ1 P˙¯ + ∂ϕ1 ϕ˙¯ + ∂ϕ1 ⎪ ⎭ ⎩ ∂ϕ2 P˙¯ + ∂ϕ2 ϕ˙¯ + ∂ϕ2 ⎪ ⎭ ∂ ϕ¯ ∂t ∂ ϕ¯ ∂t ∂ P¯ ∂ P¯    f1 P¯ + ε P1 + ε 2 P2 , ϕ¯ + ε ϕ1 + ε 2 ϕ2 , t =ε (2.61)   , f2 P¯ + ε P1 + ε 2 P2 , ϕ¯ + ε ϕ1 + ε 2 ϕ2 , t Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

ARTICLE IN PRESS

JID: FI

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

or in matrix form      f1 P¯ + ε P1 + ε 2 P2 , ϕ¯ + ε ϕ1 + ε 2 ϕ2 , t − P˙¯ A =ε   ϕ˙¯ f2 P¯ + ε P1 + ε 2 P2 , ϕ¯ + ε ϕ1 + ε 2 ϕ2 , t − where ⎡ A=⎣

1 2 1 + ε ∂P + ε 2 ∂P ∂ P¯ ∂ P¯

1 2 ε ∂P + ε 2 ∂P ∂ ϕ¯ ∂ ϕ¯

1 2 ε ∂ϕ + ε 2 ∂ϕ ∂ P¯ ∂ P¯

1 2 1 + ε ∂ϕ + ε 2 ∂ϕ ∂ ϕ¯ ∂ ϕ¯

∂P1 ∂t ∂ϕ1 ∂t



15

 ∂P2  −ε

2

∂t

∂ϕ2 ∂t

,

(2.62)

⎤ ⎦.

(2.63)

Using the expansion of 1/|A| to the order of ε2 , the inverse matrix of A is obtained as ⎡ ⎤ 1 2 1 2 1 + ε ∂ϕ + ε 2 ∂ϕ −ε ∂P − ε 2 ∂P   ∂ ϕ¯ ∂ ϕ¯ ∂ ϕ¯ ∂ ϕ¯ 1 ⎣ ⎦ = A11 A12 , A−1 = (2.64) A21 A22 |A| ∂P1 2 ∂ϕ2 2 ∂P2 1 −ε ∂ϕ − ε 1 + ε + ε ¯ ¯ ¯ ¯ ∂P ∂P ∂P ∂P where

$

% ∂P1 2 ∂P1 ∂P2 2 ∂P1 ∂ϕ1 A11 =1 − ε +ε − + , ∂ ϕ¯ ∂ P¯ ∂ P¯ ∂ P¯ ∂ P¯   ∂P1 ∂ϕ1 ∂P1 ∂P1 ∂P1 ∂P2 , A12 = − ε + ε2 − + ∂ ϕ¯ ∂ ϕ¯ ∂ ϕ¯ ∂ ϕ¯ ∂ P¯ ∂ ϕ¯   ∂ϕ1 ∂ϕ1 ∂ϕ1 ∂P1 ∂ϕ1 ∂ϕ2 , A21 = − ε + ε2 − + ∂ P¯ ∂ P¯ ∂ ϕ¯ ∂ P¯ ∂ P¯ ∂ P¯ $

% ∂ϕ1 2 ∂ϕ1 ∂ϕ2 2 ∂P1 ∂ϕ1 A22 =1 − ε +ε − + . ∂ ϕ¯ ∂ ϕ¯ ∂ P¯ ∂ ϕ¯ ∂ ϕ¯ Considering the Taylor expansion of a function with multi-variables         ∂ f j P¯ , ϕ¯ , t ∂ f j P¯ , ϕ¯ , t 2 2 ¯ ¯ f j P + ε P1 + ε P2 , ϕ¯ + ε ϕ1 + ε ϕ2 , t = f j P, ϕ¯ , t + ε P1 + ϕ1 ∂ ϕ¯ ∂ P¯       2 ∂ f j P¯ , ϕ¯ , t ∂ f j P¯ , ϕ¯ , t 1 ∂ f j P¯ , ϕ¯ , t + ε 2 P2 + ϕ2 + P12 ∂ ϕ¯ 2 ∂ P¯ ∂ P¯ 2    

  ∂ 2 f j P¯ , ϕ¯ , t 1 ∂ 2 f j P¯ , ϕ¯ , t + o ε 2 , j = 1, 2. + ϕ12 + P ϕ (2.65) 1 1 2 ¯ 2 ∂ ϕ¯ ∂ P∂ ϕ¯ Substituting Eqs. (2.64) and (2.65) into Eq. (2.62) yields     P˙¯ = ε μ1 + ε 2 μ2 + ε 3 μ3 + o ε 3 , ϕ˙¯ = ε ν 1 + ε 2 ν 2 + ε 3 ν 3 + o ε 3 ,

(2.66)

where





∂P1 ∂P2 ∂ f1 ∂ f1 ∂P1 ∂P1 ∂P1 ∂ϕ1 f1 − − f2 − − μ1 = f 1 − , μ2 = P1 + ϕ1 − , ¯ ¯ ∂t ∂ ϕ ¯ ∂t ∂ ϕ ¯ ∂t ∂t ∂P ∂P ∂ f1 ∂ f1 P2 ∂ 2 f1 ϕ12 ∂ 2 f1 ∂ 2 f1 μ3 = P2 + ϕ2 + 1 + + P ϕ 1 1 ∂ ϕ¯ 2 ∂ P¯ 2 2 ∂ ϕ¯ 2 ∂ P¯ ∂ P¯ ∂ ϕ¯

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

16

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29





∂P1 ∂ f1 ∂ f1 ∂ f2 ∂ f2 ∂P1 ∂P2 ∂ϕ2 P1 − P1 + ϕ1 − + ϕ1 − ∂ ϕ¯ ∂t ∂ ϕ¯ ∂ ϕ¯ ∂t ∂ P¯ ∂ P¯ ∂ P¯

$

2 %



∂P1 ∂P1 ∂ϕ1 ∂P2 ∂P1 ∂P1 ∂P1 ∂ϕ1 ∂P2 ∂P1 ∂ϕ1 , + f1 − − + + f2 − − + ∂t ∂ ϕ¯ ∂ P¯ ∂ P¯ ∂t ∂ ϕ¯ ∂ ϕ¯ ∂ ϕ¯ ∂ P¯ ∂ ϕ¯ ∂ P¯ (2.67)

−





∂ϕ1 ∂ϕ2 ∂ f2 ∂ f2 ∂ϕ1 ∂ϕ1 ∂P1 ∂ϕ1 f1 − − f2 − − ν 1 = f2 − , ν2 = P1 + ϕ1 − , ¯ ¯ ∂t ∂ ϕ¯ ∂t ∂ ϕ¯ ∂t ∂t ∂P ∂P ∂ f2 ∂ f2 P2 ∂ 2 f2 ϕ12 ∂ 2 f2 ∂ 2 f2 ν3 = P2 + ϕ2 + 1 + + P ϕ 1 1 ∂ ϕ¯ 2 ∂ P¯ 2 2 ∂ ϕ¯ 2 ∂ P¯ ∂ P¯ ∂ ϕ¯



∂ϕ1 ∂ f1 ∂ f1 ∂ f2 ∂ f2 ∂ϕ1 ∂P2 ∂ϕ2 P1 − P1 − + ϕ1 − + ϕ1 − ∂ ϕ¯ ∂t ∂ ϕ¯ ∂ ϕ¯ ∂t ∂ P¯ ∂ P¯ ∂ P¯

$

2 % ∂ϕ1 ∂P1 ∂ϕ1 ∂ϕ1 ∂ϕ2 + f2 − − + ¯ ∂t ∂ ϕ¯ ∂ P ∂ ϕ¯ ∂ ϕ¯



∂P1 ∂ϕ1 ∂ϕ1 ∂ϕ1 ∂P1 ∂ϕ2 . + f1 − − + (2.68) ∂t ∂ P¯ ∂ P¯ ∂ P¯ ∂ P¯ ∂ ϕ¯ Since third-order averaging can be reduced to the second order averaging, it is natural to assume that P1 and ϕ 1 have the same forms as in Eqs. (2.49). Following the same idea of deriving Eqs. (2.49), one can integrate μ2 with respect to t, solve the result for P2 , and remove any terms from P2 that are linear in ϕ¯ (in order for the near-identity transformation (2.60) to be uniformly valid in the large t limit). It can be assumed that   ∂ϕ2   ∂P2 = u2 − M u2 , = v2 − M v2 . (2.69) t t ∂t ∂t Substituting u2 and v2 in Eq. (2.43) and their averaging values in Eq. (2.53) into Eq. (2.69) and solving the differential equations yields the functions P2 and ϕ 2 , whose expressions, however, are too long to list here. To obtain μ¯ 3 and ν¯3 in Eq. (2.66), the derivatives of f1 and f2 in Eq. (2.50) are needed and some second order derivatives must be used,  sc  ∂ 2 f1 ∂  p 2 ¯ V sc P pU = −2β p sin  + ωζ cos (νt + ψ ) sin 2 + ωτ + ε 2 ∂ P¯ 2 ∂ P¯ ∂V sc sc = ωτε V + ωτε P¯ , ∂ P¯   ∂ 2 f1 ∂  p = −2β p sin2  + ωζ cos(νt + ψ ) sin 2 + ωτε pU sc + P¯ V sc 2 ∂ P¯ ∂ ϕ¯ ∂ ϕ¯ = − 2β p sin 2 + pωζ cos(νt + ψ ) cos 2 + ωτε p(U cc − U ss ) + ωτε P¯ (V cc − V ss ),   ∂ 2 f1 ∂ = pP¯ −2β sin 2 + ωζ cos(νt + ψ )cos2 + ωτε (U cc − U ss ) 2 ∂ ϕ¯ ∂ ϕ¯   = − 2pP¯ 2β cos 2 + ωζ cos(νt + ψ )sin2 + ωτε (U sc + U cs ) ,   ∂ 2 f2 ∂ ωτε V cc = p ∂ P¯ 2 ∂ P¯ Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

ARTICLE IN PRESS

JID: FI

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

17

%



1 1 1 1 1 2 −1 − + +1 2 2 p p P¯ (t )P¯ (s) p P¯ (s) P¯ (t )

cos(ωt + ϕ¯ ) cos(ωs + ϕ¯ ) ,

ωτε RL μ = D p 0 t



P¯ (s) P¯ (t )

1/p $

  ∂ 2 f2 ∂ ωτε V cc ωτε sc =− = (V + V cs ), ∂ ϕ¯ p p ∂ P¯ ∂ ϕ¯  ∂ 2 f2 ∂  −2β cos 2 − ωζ cos(νt + ψ ) sin 2 − ωτε (U cs + U sc ) = ∂ ϕ¯ 2 ∂ ϕ¯ = 4β sin 2 − 2ωζ cos(νt + ψ ) cos 2 − 2ωτε (U cc − U ss ),

(2.70)

where Eq. (2.52) and ∂U ss = U cs + U sc , ∂ ϕ¯

∂U cs = −U ss + U cc , ∂ ϕ¯

∂V sc = V cc − V ss , ∂ ϕ¯

∂V cc = −V sc − V cs , (2.71) ∂ ϕ¯

are used, and U sc , U cc , U cs , U ss , V ss , V cc , V sc , and V cs are defined in Eq. (2.51). Substituting these derivatives and P1 , ϕ 1 , P2 , and ϕ 2 into Eqs. (2.67) and (2.68) yields the values of μ3 and ν 3 . Their averaged values are

  pP¯ ζ ω ˆ cos(ψ − 2ϕ¯ ) , ˆ M μ3 = γ sin (ψ − 2 ϕ ¯ ) − 128 α ˆ β t 256ν 2

  −ζ ω ˆ ˆ M v3 = γˆ cos(ψ − 2ϕ¯ ) + 128αˆ β sin (ψ − 2ϕ¯ ) − λ , (2.72) t 256ν 2 where

ν  1 , γˆ = 7ω2 ζ 2 − 128αˆ 2 , αˆ = − ωτε H s 2    2     ν 2 2 c ν ˆλ = 128ωτε  ˆ ωτε H c ν + ωτε H s + 4βH 2 2 2     2 2 2 2 2 2 − 8αˆ 64β − 3ω ζ + 8 64β + ω ζ . Applying the transformation  = ϕ − 21 ψ and Itô Lemma results in three Itô stochastic differential equations   dP (t ) = ε m p1 + ε 2 m p2 + ε 3 m p3 P dt,   σ d(t ) = ε m1 + ε 2 m2 + ε 3 m3 dt − ε 1/2 d W (t ), 2 dψ (t ) = ε 1/2 σ d W (t ), (2.73) where

  1   ζω M u1 = p − βˆ + sin (2) , P t 4  1   pζ ω  ˆ = M u2 = − β cos(2) − αˆ sin (2) , P t 2ν

1   pζ ω ˆ = M μ3 = γˆ sin (ψ − 2ϕ¯ ) − 128αˆ β cos(ψ − 2ϕ¯ ) , P t 256ν 2

m p1 = m p2 m p3

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

18

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

  ˆ + ζ ω cos(2), m1 = M v1 =  t 4    1 ˆ sin (2) + δˆ , m2 = M v2 = ζ ωαˆ cos(2) + βωζ t 2ν

  −ζ ω ˆ ˆ m3 = M v3 = γˆ cos(ψ − 2ϕ¯ ) + 128αˆ β sin (ψ − 2ϕ¯ ) − λ . t 256ν 2

(2.74)

It is observed that Eqs. (2.73) are also coupled. Introducing a linear stochastic transformation S = T ()P, P = T −1 ()S,

0 ≤  ≤ π,

the Itô equation for the transformed pth norm process S is given by, again using Itô Lemma in the vector case      εσ 2 dS = ε m p1 + ε 2 m p2 + ε 3 m p3 T + ε m1 + ε 2 m2 + ε 3 m3 T + T P dt 8 ε 1/2 σ − T P d W (t ) 2     1  εσ 2 ε m p1 + ε 2 m p2 + ε 3 m p3 T + ε m1 + ε 2 m2 + ε 3 m3 T + = T S dt T 8 ε 1/2 σ T − S d W (t ). (2.75) 2 T The eigenvalue problem governing the pth moment Lyapunov exponent is     σ2 ˆ , T + m1 + εm2 + ε 2 m3 T + m p1 + εm p2 + ε 2 m p3 T = T 8

(2.76)

ˆ where T() is a periodic function in  of period π , and q(t ) (p) = (p) = ε (p) . Substituting the expansion in Eq. (2.31) into eigenvalue problem (2.76) and equating the coefficients of like trigonometric terms sin 2k and cos 2k, k = 0, 1, . . . yield a set of homogeneous linear algebraic equations with infinitely many equations for the unknown coefficients C0 , Ck and Sk , k = 1, 2, . . . Constant terms:    a1  ˆ ˆ 1 − αS ˆ C0 − b0 S1 + ε 2b0 βC (c0 + ) ˆ 1 + ε 2 2 λC 1 − γˆ S1 = 0, υ 8υ Coefficients for cos 2:   ˆ 0 + 2 βb ˆ 1C2 ˆ C1 − d1 S1 − b1 S2 + ε 1 − e1 S1 − 2 αb (c1 + ) ˆ 1 S2 − 4a1 βC υ  1  ˆ 0 + 4a3 λC ˆ 2 = 0, − 4a3 γˆ S2 − p λC + ε2 2 32υ Coefficients for sin 2:   ˆ 0 + 2 βb ˆ 1 S2 ˆ S1 + ε 1 e1C1 + 2 αb 2a1C0 + d1C1 + b1C2 + (c1 + ) ˆ 1C2 + 4a1 βC υ   2 1 ˆ 2 + p γˆ C0 + 4a3 γˆ C2 = 0, +ε 4a3 λS 32υ 2 Coefficients for cos 2k: Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

ARTICLE IN PRESS

JID: FI

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

19

ˆ Ck + dk Sk + bk Sk+1 ak Sk−1 − (ck + )  1 ˆ k+1 + ε 2ak αS ˆ k−1 + ek Sk + 2bk αS ˆ k+1 + 2ak αC ˆ k−1 − 2bk βC υ  2 1 − bk−2 γˆ Sk−1 − ak+2 γˆ Sk+1 +ε 8 υ2  ˆ k−1 + ak+2 λC ˆ k+1 = 0, − bk−2 λC k = 2, 3, . . . , K − 1, Coefficients for sin 2k: ˆ Sk ak Ck−1 + dk Ck + bk Ck+1 + (ck + )  1 ˆ k+1 + ε 2ak αC ˆ k−1 + ek Ck + 2bk αC ˆ k+1 − 2ak αS ˆ k−1 + 2bk βS υ 2 1 ˆ k−1 + ak+2 λS ˆ k+1 − bk−2 λS +ε 8υ 2

+ bk−2 γˆ Ck−1 + ak+2 γˆ Ck+1 = 0, k = 2, 3, . . . , K − 1, Coefficients for cos 2K:

  ˆ K−1 ˆ CK + dK SK + ε 1 2ak αS aK SK−1 − (cK + ) ˆ K−1 + ek SK + 2ak βC υ  1  ˆ K−1 = 0, + ε 2 bk−2 γˆ SK−1 + bk−2 λC 8υ Coefficients for sin 2K:   ˆ K−1 ˆ SK + ε 1 2ak αC aK CK−1 + dK CK + (cK + ) ˆ K−1 + ek CK + 2ak βS υ  1  ˆ K−1 = 0, + ε 2 bk−2 γˆ CK−1 − bk−2 λS 8υ where

(2.77)

1 1 ζ ω(2k − 2 − p), bk = ζ ω(2k + 2 + p), 8 8 1 2 2 ˆ dk = 2k , ˆ k = 0, 1, 2, . . . . ˆ ck = σ k + pβ, ek = k δ, 2

ak =

As in the first-order and second-order stochastic averaging, Eq. (2.77) is reduced to a set of 2K + 1 homogeneous linear equations for C0 , Ck and Sk , k = 1, 2, . . . , K . The determinant ˆ K (p). Solving the is set to zero to yield a polynomial equation of degree 2K + 1 for  polynomial equation   for thelargest  root yields the approximate moment Lyapunov exponents. The terms H s 21 ν and H c 21 ν depend upon the material’s constitutive model. 3. Results and discussions The moment Lyapunov exponent is the nontrivial solution of Eq. (2.33), where the Fourier series expansion for T() is truncated to the Kth order. Theoretically, the set of approximate moment Lyapunov exponents obtained by this procedure converges to the corresponding true exponent values as K → ∞. Practically, Fig. 2 illustrates that moment Lyapunov moments converge so quickly with the increase of Fourier expansion order K that results of K = 3 are Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

20

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

Fig. 2. Effect of Fourier series expansion on moment Lyapunov exponents.

Fig. 3. Fourier series expansion on moment Lyapunov exponents.

accurate enough. However, this convergence depends on the width of power spectral density of the bounded noise process. Fig. 3 shows when ζ = 5, even K = 12 is not adequate. This result is different from [14], where real noise was considered. From Fig. 4, it is seen that for small σ or large ζ the bounded noise is a narrow-band process; whereas when σ is increased or ζ is decreased, the power spectral density function becomes flat. Hence, if the power spectral density function of the bounded noise is relatively sharp, which means it is a narrow band process, a large number of sinusoidal terms in the Fourier series Eq. (2.33) must be taken to obtain a good approximation of the moment Lyapunov exponent (p). On the other hand, if the power spectral density function of the bounded noise is relatively flat, which means it is a wide band process, fewer sinusoidal terms in the Fourier series Eq. (2.33) can yield a satisfactory approximation of (p). Effect of damping β on moment Lyapunov exponents is shown in Fig. 5, in which the parameter β plays a stabilizing role in the system stability. The approximate moment Lyapunov Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

21

Fig. 4. Power spectral density of a bounded noise process.

Fig. 5. Effect of damping β on moment Lyapunov exponents.

exponents are compared with the results from Monte-Carlo simulation [7], which shows that the two results agree well, although discrepancy is found when p is large. Effect of fractional orders μ on moment Lyapunov exponents is shown in Fig. 6, which indicates that the parameter μ plays an important role in the system stability. The increase of the fractional order μ would enhance the system stability. When μ changes from 0 to 1, the system gradually becomes stabilized. As mentioned in Appendix B, when μ shifts from 0 to 1, the fractional Newton element in a fractional Kelvin–Voigt model moves from a Hooke element to an integer-value Newton element, i.e., from pure elasticity to pure viscosity. Therefore, the fractional order μ enables a subtle and smooth transition from elasticity to viscosity in one element, which abandons the conventional binary theory that only allow an element either to function totally elastic (a spring) or perform perfect viscosity (a dashpot). Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

22

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

Fig. 6. Effect of fractional orders μ on moment Lyapunov exponents.

Fig. 7. Effect of σ of bounded noise on moment Lyapunov exponents.

This multi-state fractional component has various performance levels and creates various viscoelastic modes. Fig. 6 clearly shows fractional components have vital effects on the entire system’s performance. Fig. 7 shows that with the increase of the noise intensity parameter σ , the stability of the system increases. Three-dimensionally, Fig. 8 show that the smaller the noise intensity σ , the more significant the parametric resonance. Effect of noise amplitude (ζ ) on moment Lyapunov exponents is shown in Fig. 9, which indicates the de-stabilizing effect of noise amplitude on the system stability. Great effect of ζ on parametric resonance is observed in Fig. 10. Larger amplitude would result in prominent parametric resonance. One possible explanation is that from the power spectrum density function of bounded noise in Eq. (2.10), Fig. 4 shows the larger of σ or the smaller of ζ , the wider the frequency band of power spectrum, which Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

23

Fig. 8. Parametric resonance of SDOF fractional viscoelastic system.

Fig. 9. Effect of noise amplitude ζ on moment Lyapunov exponents.

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

24

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

Fig. 10. Parametric resonance of a SDOF fractional viscoelastic system.

Fig. 11. Effect of retardation τ ε on moment Lyapunov exponents.

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

ARTICLE IN PRESS

JID: FI

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

25

suggests that the power of the noise not be concentrated in the neighborhood of the central frequency ν, which reduces the effect of the primary parametric resonance. Fig. 11 illustrates that the system retardation factor (τ ε ) plays a stabilizing role in moment stability, for the stability index increases with the retardation factor. The stability index is the non-trivial zero δ of the moment stability (p), i.e., (δ) = 0, which characterizes the probability with which the response of an almost-surely system exceeds a threshold. The above figures illustrate that there is some difference between first-order and higherorder averaging. However, the third-order averaging is almost identical to the second-order averaging, which suggests that for viscoelastic structures, second-order averaging analysis is adequate. 4. Conclusions The stochastic stability of a fractional viscoelastic column under the excitation of a bounded noise is studied by using the method of higher order stochastic averaging. The higher order stochastic averaging is based on lower order stochastic averaging and is formulated by introducing a near identity transformation. We investigate the stochastic stability of the system through moment Lyapunov moment, which is an ideal avenue to characterize moment stability. The moment Lyapunov exponent is determined numerically by solving a set of homogeneous linear equations, which is obtained from an eigenvalue problem by using linear stochastic transformation and Itô Lemma. It is found that there is some difference between the results from the first-order and higher-order averaging. However, Results from the third-order averaging are almost identical to those of the second-order averaging, which suggests that for this viscoelastic structure, second-order averaging analysis is adequate. The approximate results are compared well to Monte Carlo simulations, which shows that the tedious derivation for higher-order stochastic averaging is correct in this paper. The convergence issue and the effects of various parameters on the stochastic stability of the system are discussed in detail and possible explanations are explored. These results pave a road for utilizing or controlling parametric resonances in engineering applications. Appendix A. Fractional calculus Consider the first three repeated integration and reduce to a single integral by using the usual change of variables formula  t (−1) f (t ) = f (τ )dτ, 0

f

(−2)

 (t ) =

t

 dτ2

0

f

(−3)

 (t ) = 

0



t

 dτ3

 f (τ1 )dτ1

τ3

f (τ1 )dτ1

 dτ2 

τ2



0 t

τ1

t

f (τ1 )dτ1 =

(t3 − τ1 )dτ3 =

 dτ2 =



t

(t − τ1 ) f (τ1 )dτ1 ,

0



dτ3

0

1 2

t

τ1

0

0 t

t

f (τ1 )dτ1 =

0

0

=

τ2

t3

(t3 − τ1 ) f (τ1 )dτ1

0 t

(t − τ1 )2 f (τ1 )dτ1 .

(A.1)

0

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

ARTICLE IN PRESS

JID: FI

26

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

One can continue in this way to find that  t d−n f (t ) 1 f (−n) (t ) = = (t − τ )n−1 f (τ )dτ, n ∈ N + , dt (n − 1)! 0

(A.2)

for almost all x with −∞ ≤ 0 < t ≤ ∞, and N + is the set of positive integers. By replacing the factorial expression for its gamma function equivalent, one can generalize the formula for repeated integration in Eq. (A.2) for all n ∈ R+ ,  t  1 −μ  (−μ) RL f (t ) = 0 Dt f (t ) = (t − τ )μ−1 f (τ )dτ, μ ∈ R+ , (A.3) (μ) 0   −μ · is an operator where μ is the fractional coefficient, R+ is positive real numbers, RL 0 Dt for Riemann–Liouville fractional integral of order μ, (μ) is the gamma function,  ∞ (μ) = e−t t μ−1 d t. (A.4) 0

If μ is a positive integer, then (μ) = (μ − 1)! , so this definition of f (−μ) (t ) agrees with (A.2) when μ is a positive integer. Eq. (A.3) is the definition of Riemann–Liouville (RL) fractional integral [15]. In a similar way, fractional derivative can be defined from repeated derivative of a function. However, it is more convenient to formulate a definition of fractional derivative directly using the fractional integral. For a fixed n ≥ 1 and integer m ≥ n the (m − n)th derivative of the function f(t) can be written as  t  (−n)  1 (m−n) (m) m f f (t ) = f (t ) = D (t − τ )n−1 f (τ )dτ, m ∈ N + , (A.5) (n) 0 where the symbol Dm (m ≥ 0) denotes mth iterated differentiations. Substituting Eq. (A.3) into Eq. (A.5) leads to      1 dm t RLDm−α f (t ) = RLDm RLD−α f (t ) = (t − τ )α−1 f (τ )dτ, 0 < α ≤ 1, m ∈ N + . 0 t 0 t 0 t (α) dt m 0 (A.6) Denoting μ = m − α, one can write fractional derivative as    1 dm t f (τ ) RLDμ f (t ) = dτ, m − 1 ≤ μ < m, (A.7) 0 t (m − μ) d t m 0 (t − τ )−(m−μ)+1 μ  where RL for Riemann–Liouville fractional derivative of order μ. It a Dt · is an operator    0 m (m) f = can be seen that RL D (t ) f (t ) and RL (t ). Specifically, for 0 < μ ≤ 1, 0 t 0 Dt f (t ) = f assuming f(t) is absolutely continuous for t > 0 , one has (see Eq. (2.1.28) of [15])   f (0)  t f  (τ )    1 d t f (τ ) 1 RLDμ f (t ) = d τ = + d τ . (A.8) 0 t μ (1 − μ) d t 0 (t − τ )μ (1 − μ) t μ 0 (t − τ ) Further assuming f (0) = 0 gives the fractional derivative in Caputo (C) form  t  t    1 f  (τ ) 1 f (t − τ ) CDμ f (t ) = d τ = dτ. (A.9) 0 t (1 − μ) 0 (t − τ )μ (1 − μ) 0 τμ Fig. A.12 gives a diagram for a simple function f (t ) = t and its fractional calculus as μ  1−μ a showcase RL (2)/(2 − μ). The fractional calculus of arbitrary real order μ 0 Dt t = t can be considered as an interpolation of the sequence of integer-order calculus. Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

ARTICLE IN PRESS

JID: FI

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

27

Fig. A.12. Fractional calculus of a power function.

Fig. B.13. Schematic of a fractional element between a spring (μ = 0) and a dashpot (μ = 1).

Appendix B. Fractional viscoelasticity Fractional calculus provides a perfect tool to construct a fractional viscoelastic component shown in Fig. B.13. With the various values of μ from 0 to 1, property of the fractional element changes from elastic to viscous. Fractional Viscoelasticity may be formulated by combining Hook springs, fractional pots, and dash pots in series or in parallel. The fractional Kelvin–Voigt viscoelastic mechanical model consists of a fractional Newton dash-pot and a Hooke spring in parallel, which is mathematically defined as  dμ ε(t ) μ = E ε(t ) + η RL (B.1) 0 Dt ε(t ) , d tμ where E is the Young’s elastic modulus and η is the viscosity modulus or coefficient of viscosity. σ (t ) = E ε(t ) + η

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

28

ARTICLE IN PRESS

[m1+;October 11, 2017;9:14]

J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

Fig. B.14. Creep and recovery for a fractional viscoelastic material.

When μ = 0, the fractional Newton part becomes a constant η, which is the case for the Hooke model. It is further observed that, only when μ = 0, this model exhibits transient elasticity at t = 0. When μ = 1, the relaxation modulus of the fractional Newton part becomes a Dirac delta function, i.e., σ (t ) = η δ(t ), which is the case of the integer Newton model. Hence, the fractional Newton element is a continuum between the spring model and the Newton model (Fig. B.13). A typical relationship of creep and recovery for a fractional viscoelastic material is displayed in Fig. B.14, which clearly shows that the extra degreeof-freedom from the fractional order can improve the performance of traditional viscoelastic elements. Specifically, when μ = 1, this fractional model reduces to the ordinary integer Kelvin–Voigt model. References [1] R. Mainardi, Fractional Calculus and Waves in Linear Viscoelasticity: An Introduction to Mathematical Models, Imperial College Press, London, 2010. [2] H. Xu, X. Jiang, Creep constitutive models for viscoelastic materials based on fractional derivatives, Comput. Math. Appl. 73 (6) (2017) 1377–1384. [3] Q. Song, X. Yang, C. Li, T. Huang, X. Chen, Stability analysis of nonlinear fractional-order systems with variable-time impulses, J. Frankl. Inst. 354 (2017) 2959–2978. [4] W.-C. Xie, Dynamic Stability of Structures, Cambridge University Press, Cambridge, 2006. [5] Q. Ling, X.L. Jin, Z.L. Huang, Response and stability of SDOF viscoelastic system under wideband noise excitations, J. Frankl. Inst. 348 (2011) 2026–2043. [6] Y. Xu, B. Pei, J.L. Wu, Stochastic averaging principle for differential equations with non-Lipschitz coefficients driven by fractional Brownian motion, Stoch. Dyn. 17 (2017) 1750013. [7] J. Deng, W.-C. Xie, M.D. Pandey, Stochastic stability of a fractional viscoelastic column under bounded noise excitation, J Sound Vib. 333 (6) (2014) 1629–1643. [8] M. Hijawi, N. Moshchuk, R.A. Ibrahim, Unified second-order stochastic averaging approach, ASME J. Appl. Mech. 64 (2) (1997) 281–291. [9] Q. Huang, W.-C. Xie, Stability of SDOF linear viscoelastic system under the excitation of wideband noise, ASME J. Appl. Mech. 75 (2) (2008) 021012. [10] G. Janevski, P. Kozic, I. Pavlovic, Moment Lyapunov exponents of the parametrical Hill’s equation under the excitation of two correlated wideband noises, Struct. Eng. Mech. 52 (3) (2014) 525–540.

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019

JID: FI

ARTICLE IN PRESS J. Deng / Journal of the Franklin Institute 000 (2017) 1–29

[m1+;October 11, 2017;9:14]

29

[11] J. Deng, W.-C. Xie, M.D. Pandey, Higher-order stochastic averaging to study stability of a fractional viscoelastic column, J Sound Vib. 333 (23) (2014) 6121–6139. [12] C. Floris, Stochastic stability of a viscoelastic column axially loaded by a white noise force, Mech. Res. Comm. 38 (2011) 57–61. [13] W. Wedig, Analysis and simulation of nonlinear stochastic systems, in: W. Schiehlen (Ed.), Nonlinear Dynamics in Engineering Systems, Springer-Verlag, Berlin, 1989, pp. 337–344. [14] N.S. Namachchivaya, H.J.V. Roessel, Moment Lyapunov exponent and stochastic stability of two coupled oscillators driven by real noise, ASME J. Appl. Mech. 68 (6) (2001) 903–914. [15] A.A. Kilbas, H.M. Srivastava, J.J. Trujillo, Theory and Applications of Fractional Differential Equations, Elsevier, New York, 2006.

Please cite this article as: J. Deng, Higher-order stochastic averaging for a SDOF fractional viscoelastic system under bounded noise excitation, Journal of the Franklin Institute (2017), https://doi.org/10.1016/j.jfranklin.2017.09.019