Abstract

Complex fractional moment (CFM), which is defined as the Mellin transform of a probability density function (PDF), has been successfully employed to find the response PDF of a wide variety of integer-order nonlinear oscillators. In this paper, a CFM-based analysis is performed to determine the transient response PDF of nonlinear oscillators with fractional derivative elements under Gaussian white noise. First, an equivalent linear system is introduced for the purpose of deriving the Fokker–Planck (FP) equation for response amplitude. The equivalent natural frequency and equivalent damping coefficient of the system need to be determined, taking into account both the nonlinear and fractional derivative elements of the original oscillator. Moreover, to convert the FP equation into the governing equation of CFMs, these equivalent coefficients must be given in polynomial form of amplitude. This paper proposes formulas for appropriately determining the equivalent coefficients, based on an equivalent linearization technique. Then, applying stochastic averaging, the FP equation is derived from the equivalent linear system. Next, the Mellin transform converts the FP equation into coupled linear ordinary differential equations for amplitude CFMs, which are solved with a constraint corresponding to the normalization condition for a PDF. Finally, the inverse Mellin transform of the CFMs yields the amplitude PDF. The joint PDF of displacement and velocity is also obtained from the amplitude PDF. Three linear and nonlinear fractional oscillators are considered in numerical examples. For all cases, the analytical results are in good agreement with the pertinent Monte Carlo simulation results.

Graphical Abstract Figure
Graphical Abstract Figure
Close modal

1 Introduction

Mechanical and structural systems in various engineering applications are exposed to random excitations, such as seismic ground motion, fluctuating pressure in turbulent flow, and hydraulic force of ocean waves. In order to grasp the response characteristics of dynamic systems subjected to random excitation and to assess their reliability, it is important to accurately determine the probability density function (PDF) of the system response. For a nonlinear system under Gaussian white noise, the system response can usually be regarded as a Markov process, and the time evolution of the response PDF is governed by a partial differential equation called Fokker–Planck (FP) equation. Unfortunately, exact solutions of the FP equation are not available except for the stationary response of a very limited class of nonlinear systems [1]. Solutions for transient/non-stationary responses are even more difficult to obtain, owing to their time dependence. Therefore, the development of approximate analytical or numerical techniques for finding system response PDFs is a significant issue in the field of stochastic dynamics.

In the last few decades, fractional calculus has attracted a great deal of attention due to its usefulness as a mathematical modeling tool for non-local or long-memory phenomena. Various topics on fractional calculus in science and engineering can be found in the books [27] and the papers [810]. One successful application of fractional calculus to engineering is the modeling of viscoelastic materials. Indeed, many experimental results have proven that fractional derivative models with a small number of parameters can describe viscoelastic behavior with a high degree of accuracy [1118]. Therefore, in mechanical and structural dynamics, the fractional models have become widely used as models of viscoelastic dampers for seismic isolation and vibration control [16,1923].

Under such circumstances, there has been a growing interest in the response determination of dynamic systems with fractional derivative elements under stochastic excitation. In this regard, Monte Carlo simulation (MCS) comes to mind as a general-purpose method. However, computing the fractional derivative operator requires considerable computational time and resources, since this operator is defined by a convolution integral, i.e., depends on the function’s entire past time history. Hence, direct numerical integration of system equations with fractional-order derivatives is much more computationally demanding than that of systems with only integer-order derivatives. For this reason, in recent years, a lot of efforts have been put into developing analytical techniques for determining the stochastic response of fractional-order systems. Meanwhile, machines and structures including viscoelastic materials have often been modeled as single-degree-of-freedom (SDOF) systems with fractional derivative elements [2,4,8,9,1214,16]. Therefore, the development of analytical methods for SDOF systems with fractional derivatives is interesting and important from a practical point of view.

For moment analysis, an approach based on convolutional representation of the response was presented to obtain the response variance of SDOF linear systems with fractional derivatives [24,25]. More recently, approximate closed-form solutions of response moments were obtained for fractional linear oscillators to evolutionary stochastic excitation by expressing their impulse response with Prony series [26]. A moment equation approach has also been employed [2729], where fractional-order systems are approximately converted into coupled integer-order differential equations. In addition, frequency-domain approaches have been developed for stationary moment analysis [3032].

For the analysis of PDFs, stochastic averaging has been widely utilized. The displacement and velocity responses of fractional systems are no longer a Markov vector due to the history dependence of the fractional terms, even if the excitation is a white noise process. Nevertheless, by applying stochastic averaging, the response amplitude can be approximately treated as a Markov process, resulting in the FP equation governing the amplitude PDF [33]. As an additional benefit, system dimension reduction is also achieved. In this manner, the amplitude PDF of the stationary response of various fractional nonlinear oscillators has been accurately obtained [3341].

Compared to moment and stationary PDF analyses, the analysis of transient response PDF is a more difficult task. Although the number of studies is still limited at present, efforts are being made to seek an approximate solution of the transient amplitude PDF of fractional nonlinear systems using a path integral technique [42] or a Rayleigh distribution assumption [43]. However, these existing methods have some drawbacks in terms of computational cost and applicability. The path integration generally requires very short time-steps to obtain a highly accurate solution. The Rayleigh distribution assumption for amplitude PDF may not be appropriate for some types of nonlinear systems, such as oscillators possessing limit cycles. Therefore, it remains an important challenge to develop efficient and versatile analytical methods for determining the transient response PDF of nonlinear oscillators with fractional derivative elements.

Meanwhile, a novel statistical moment named complex fractional moment (CFM) was proposed for the characterization of random variables [44,45]. The CFM is a complex-order moment, which is calculated as the Mellin transform of a PDF. Similar to the Fourier transform relation between PDFs and characteristic functions, the inverse Mellin transform of the CFM reconstructs the original PDF. This striking feature of the CFM allows us to estimate PDFs via CFMs. Di Paola presented a CFM-based approach for solving the FP equation corresponding to first-order nonlinear systems excited by Gaussian white noise [46], in which by applying the Mellin transform to the FP equation, it is converted into a set of linear ordinary differential equations with respect to response CFMs. In other words, a partial differential equation is converted into linear ordinary differential equations much easier to solve. Then, solving the CFM equations and taking the inverse Mellin transform on the resulting CFMs, the response PDF is constructed. Through the above method, the PDFs of both transient and stationary responses were accurately determined in the whole domain, including the tail region. This method has also been applied successfully for white noise excitations of other types [4749]. Jin et al. developed the CFM-based approach for SDOF oscillators with nonlinear damping under Gaussian white noise by combining stochastic averaging [50]. Stochastic averaging leads to the FP equation for response amplitude, which is then solved via CFMs. The CFM method combined with stochastic averaging yielded transient and stationary response PDFs with high accuracy, and the computational time was significantly reduced compared to the corresponding MCS. In the subsequent studies, this method has been extended to a wide class of integer-order nonlinear oscillators, including systems with nonlinearities in both stiffness and damping [51], nonlinear vibro-impact systems [52] and nonlinear systems with time delay [53]. Such great success of the CFM-based method suggests that it can also be a powerful tool in the transient response analysis of fractional-order nonlinear oscillators.

In this paper, a CFM-based analysis is performed to determine the transient response PDF of nonlinear oscillators with fractional derivative elements, whose nonlinearities are included in both stiffness and damping. This analysis is mainly based on an analytical method for integer-order nonlinear oscillators recently developed by some of the authors [51]. First, we aim to derive the FP equation for response amplitude using stochastic averaging. For this purpose, the equivalent natural frequency and equivalent damping coefficient of the oscillator need to be determined. As an important point in the analysis of fractional oscillators, these equivalent coefficients have to include not only the effects of nonlinearities, but also the contributions of the fractional derivative terms. Furthermore, the coefficients must be given in polynomial form of response amplitude to utilize the CFM-based approach. This paper proposes formulas for determining the equivalent coefficients that satisfy the above requirements, based on an equivalent linearization technique. Then, resorting to stochastic averaging, the FP equation is derived. The solution of the FP equation, i.e., the amplitude PDF, is sought through the CFM approach. The joint PDF of displacement and velocity is also obtained from the amplitude PDF. The accuracy and versatility of the present analytical technique are demonstrated through numerical examples of three linear and nonlinear oscillators with fractional derivatives.

2 Preliminaries

2.1 Mellin Transform.

The Mellin transform Mf(γ1) of a real function f(x) with domain 0x< is defined as [45]
Mf(γ1)=0xγ1f(x)dx,γ=ρ+iη
(1)
where i=1 and ρ,ηR. The inverse Mellin transform is given in the form
f(x)=12πMf(γ1)xγdη,x>0
(2)
The integration in Eq. (2) is calculated with respect to the imaginary part η of γ while the real part ρ is fixed. The condition for the existence of the Mellin transform is
{p<ρ<qlimx0f(x)=O(xp),limxf(x)=O(xq)
(3)
where O() denotes the Landau symbol. The range of ρ in Eq. (3) is referred to as Fundamental Strip (FS) of the Mellin transform [45].

2.2 Complex Fractional Moment.

Let Y be a non-negative random variable with a PDF pY(y). The CFM of Y is defined as [45]
E[Yγ1]=0yγ1pY(y)dy=MpY(γ1)
(4)
where γ is a complex number γ=ρ+iη, E[] represents mathematical expectation. As can be seen from Eq. (4), the CFM is given as the Mellin transform of the PDF pY(y). Therefore, CFM is henceforth denoted as MpY(γ1). Using the inverse Mellin transform, the PDF pY(y) can be reconstructed.
pY(y)=12πMpY(γ1)yγdη,y>0
(5)
pY(y) is determined from Eq. (5) uniquely, irrespective of the choice of the real part ρ of γ, provided that ρ belongs to the FS [45]. Equation (5) is discretized in preparation for the analysis in Sec. 5.
pY(y)=12bs=mmMpY(γs1)yγs,y>0
(6)
where b=π/Δη and γs=ρ+isΔη. Δη is the discretization step of η. The value of m defines a cutoff value η¯=mΔη. CFM generally converges to zero as |η| [44], which means that the effect of truncating η in Eq. (6) is negligible if a sufficiently large η¯ is chosen.

3 Equation of Motion

Consider an SDOF nonlinear oscillator with a fractional derivative element. The system equation of motion is given by
X¨+ελD0,tαX+εh(X,X˙)+ω02X+εω02g(X)=ε12W(t)
(7)
where X denotes the system displacement. The dot symbol represents the (integer-order) derivative with respect to time t. h(X,X˙) and g(X) are nonlinear functions that account for the nonlinearities of damping and stiffness, respectively. W(t) is a zero-mean Gaussian white noise of power spectral density S0 and is suddenly applied to the system at t=0. λ is a coefficient, and D0,tαX denotes the Caputo fractional derivative of order 0<α<1 of X(t) defined as [2]
D0,tαX(t)=1Γ(1α)0tX˙(s)(ts)αds
(8)
where Γ() is the gamma function. The fractional derivative term with 0<α<1, which is commonly referred to as “springpot,” represents viscoelastic damping force that is intermediate between the elastic force of a spring (α=0) and the viscous damping force of a dashpot (α=1). This term behaves more like a spring as α approaches 0, whereas dashpot-like behavior prevails when α is close to 1. The parameter ε>0 is employed to quantify the magnitude of the fractional derivative, nonlinear damping, nonlinear stiffness, and excitation terms. ω0 is the natural frequency of the corresponding linear system (ε=0). The system (7) is quiescent for t0.
Assume that h(X,X˙) is a polynomial of the form
h(X,X˙)=i=1ciXniX˙2mi+1,ni,miZ0
(9)
where ci are coefficients, and Z0 denotes the set of non-negative integers. From Eq. (9), the order of X˙ in h(X,X˙) is odd. Equation (9) includes various nonlinear damping models for dynamic systems of engineering interest, e.g., damping of ship roll motion [54] and energy-dependent damping [55]. Furthermore, Eq. (9) also includes nonlinear damping for oscillators exhibiting limit-cycle behavior such as the van der Pol oscillator. Thus, Eq. (9) can represent the damping nonlinearities of diverse oscillatory systems. It is also assumed that g(X) is given by a polynomial of X.

The purpose of this paper is to obtain the transient response PDF of the oscillator expressed by Eq. (7). Some of the authors recently developed a CFM-based analytical approach for determining the transient response PDF of an integer-order nonlinear oscillator obtained by removing the fractional derivative term from Eq. (7) [51]. The present analysis for the fractional nonlinear oscillator (7) is also performed based on this analytical approach. An important point that must be addressed to accomplish this analysis is the appropriate determination of the equivalent natural frequency and equivalent damping coefficient for the fractional nonlinear oscillator.

4 Derivation of Fokker–Planck Equation for Response Amplitude

4.1 Pseudo-Harmonic Representation of Response.

Assuming that ε1, the response of the oscillator of Eq. (7) exhibits pseudo-harmonic behavior. Thus, X(t) and X˙(t) are described by the following equations [33,42,43]:
X(t)=A(t)cos(ωe(A)t+Φ(t))X˙(t)=A(t)ωe(A)sin(ωe(A)t+Φ(t))
(10)
where A(t) and Φ(t) are amplitude and phase, respectively, which are slowly varying functions of time t. ωe(A) is an equivalent natural frequency of the oscillator, which depends on A(t). In Ref. [51], ωe(A) was determined considering the effect of nonlinear stiffness. For the fractional oscillator under study, its fractional derivative (viscoelastic) element contributes to both elasticity and viscosity. Therefore, ωe(A) needs to be determined by taking account of both the stiffness nonlinearity and the elastic contribution of the fractional element. To this aim, we propose a formula for determining ωe(A) based on an equivalent linearization technique.

4.2 Equivalent Linearization.

An equivalent linear system for the original oscillator (7) is introduced, whose equation of motion is given as
X¨+εβe(A)X˙+ωe2(A)X=ε12W(t)
(11)
where ωe(A) is an equivalent natural frequency. βe(A) is an equivalent damping coefficient, depending on the amplitude A(t). Similarly to ωe(A), we determine βe(A) by taking into account both the nonlinear damping and the contribution of the fractional (viscoelastic) term to viscous damping. The equivalent system behaves in a pseudo-harmonic manner. In other words, the response of the system (11) is given by Eq. (10).
In the equivalent linearization technique, the equivalent coefficients ωe(A) and βe(A) are selected according to some statistical criterion. In this paper, the criterion of least mean square error is used. Thus, ωe(A) and βe(A) are determined in such a way that the mean square error between the original oscillator (7) and the equivalent linear system (11) is minimized. Performing a mean square minimization procedure [5658] on the error between Eqs. (7) and (11) yields ωe(A) and βe(A) in the form
ωe2(A)=ω02{1+επA02πg(Acosψ)cosψdψ}+ελπ02π(D0,tαcosψ)cosψdψ
(12)
and
βe(A)=1πωe(A)A02πh(Acosψ,Aωe(A)sinψ)sinψdψλπωe(A)02π(D0,tαcosψ)sinψdψ
(13)
where ψ=ωe(A)t+Φ. In deriving Eqs. (12) and (13), A(t) and Φ(t) were considered constant during one period of oscillation, and the fact that h(X,X˙) is given by Eq. (9) was employed. The integrals involving fractional derivatives in Eqs. (12) and (13) can be calculated as follows [43]:
02π(D0,tαcosψ)cosψdψ=πωeα(A)cos(απ2)
(14)
02π(D0,tαcosψ)sinψdψ=πωeα(A)sin(απ2)
(15)
Then, Eqs. (12) and (13) become
ωe2(A)=ω02{1+ε1πA02πg(Acosψ)cosψdψ}+ελωeα(A)cos(απ2)
(16)
and
βe(A)=1πωe(A)A02πh(Acosψ,Aωe(A)sinψ)sinψdψ+λωeα1(A)sin(απ2)
(17)

To solve the FP equation derived in the following section via the CFM-based approach, ωe2(A) and βe(A) must be given as polynomials in amplitude A [51]. Unfortunately, this requirement is not met, due to the presence of ωeα(A) and ωeα1(A) in Eqs. (16) and (17) (0<α<1). In order to overcome this problem, ωeα(A) in Eq. (16) and ωeα1(A) in Eq. (17) are approximated by ω0α and ω0α1, respectively. This analysis assumes ε1. In Eq. (16), the term to which the approximation is made includes ε. On the other hand, Eq. (17) does not include ε, but the damping term of the equivalent linear system in Eq. (11) is multiplied by ε. Therefore, the effect of the approximation is considered sufficiently small.

Finally, the formulas for determining ωe(A) and βe(A) in the present analysis are as follows:
ωe2(A)=ω02{1+ε1πA02πg(Acosψ)cosψdψ}+ελω0αcos(απ2)
(18)
and
βe(A)=1πωe(A)A02πh(Acosψ,Aωe(A)sinψ)sinψdψ+λω0α1sin(απ2)
(19)
The last constant terms in Eqs. (18) and (19) represent the elastic and viscous contributions of the fractional element, respectively. Substituting these equivalent coefficients into Eq. (11) yields an equivalent linear system for the original fractional nonlinear oscillator. Subsequent analysis proceeds in the same manner as in Ref. [51].

4.3 Stochastic Averaging.

Using Eqs. (10) and (11), the first-order differential equation for the response amplitude A(t) is derived as [51]
A˙=εβe(A)Asin2(ωe(A)t+Φ)ε12sin(ωe(A)t+Φ)ωe(A)W(t)
(20)
Herein, Eq. (20) is decoupled from the phase Φ(t) by utilizing stochastic averaging [59,60]. Then, A(t) is described by the following one-dimensional stochastic differential equation:
A˙=εβe(A)A2+επS02Aωe2(A)+επS0ωe2(A)ξ(t)
(21)
where ξ(t) is a zero-mean delta-correlated process of unit intensity (i.e., E[ξ(t)ξ(t+τ)]=δ(τ), in which δ() is the Dirac delta function).
Equation (21) is a stochastic differential equation in the Stratonovich sense. Considering the so-called Wong–Zakai correction term, the FP equation corresponding to Eq. (21) is derived as follows [51]:
pA(a,t)t=a[{εβe(a)a2+επS02aωe2(a)+επS04a(1ωe2(a))}pA(a,t)]+επS022a2(pA(a,t)ωe2(a))
(22)
where pA(a,t) is the PDF of the amplitude A(t). Since the initial conditions of the system are X(0)=0 and X˙(0)=0, the initial amplitude PDF pA(a,0) is given as
pA(a,0)=δ^(a)
(23)
where δ^(a) is the one-sided Dirac delta function, which satisfies
0δ^(a)da=1
(24)

Finally, the magnitude of the parameter ε of the oscillator in performing the above analysis is mentioned. In Refs. [61] and [62], the stationary response PDFs and transient response statistics of various integer-order nonlinear systems were analyzed by combining the equivalent linearization and the stochastic averaging. It was shown that these PDFs and statistics can be obtained correctly when ε is less than 0.1. Thus, the range of ε for which the combination of the equivalent linearization and the stochastic averaging works well is expected to be ε<0.1.

5 Complex Fractional Moment Approach

5.1 Conversion of Fokker–Planck Equation to Complex Fractional Moment Equations.

From Eq. (18), the equivalent natural frequency ωe(A) is represented in the form
ωe2(A)=ω02{1+εs(A)}
(25)
where
s(A)=1πA02πg(Acosψ)cosψdψ+λω0α2cos(απ2)
(26)
Since ε1, 1/ωe2(A) is approximated according to the Maclaurin expansion as follows [51]:
1ωe2(A)=1ω02{1+εs(A)}1ω02{1εs(A)}
(27)
Substituting Eq. (27) into Eq. (22) leads to
pA(a,t)t=a[{εβe(a)a2+επS0{1εs(a)}2aω02ε2πS04ω02s(a)a}pA(a,t)]+επS02ω022a2[{1εs(a)}pA(a,t)]
(28)
Because the system nonlinear damping h(X,X˙) and nonlinear stiffness g(X) are given in polynomial form, s(A) and βe(A), calculated by Eqs. (19) and (26), respectively, are represented as polynomials in A
s(A)=j=0nssjAj
(29)
βe(A)=j=0nββjAj
(30)
where ns and nβ are the highest degrees of the polynomials, and sj and βj are coefficients. As described in Sec. 4.2, s(A) and βe(A) have polynomial forms thanks to the approximation of Eqs. (16) and (17) to Eqs. (18) and (19).
The Mellin transform of both sides of Eq. (28) while considering Eqs. (29) and (30) yields
tMpA(γ1,t)=ε2{j=0nββj[aγ+jpA(a,t)]0(γ1)j=0nββj0aγ1+jpA(a,t)da}+επS02ω02{[aγ1apA(a,t)]0γ[aγ2pA(a,t)]0+(γ1)20aγ3pA(a,t)da}+ε2πS04ω02{j=0ns(2γj)sj[aγ2+jpA(a,t)]02j=0nssj[aγ1+japA(a,t)]0(γ1)j=0ns(2γ2+j)sj0aγ3+jpA(a,t)da}
(31)
where MpA(γ1,t) is the amplitude CFM. Herein, it is assumed that for any time t, the amplitude PDF pA(a,t) does not diverge as a0, and pA(a,t) and pA(a,t)/a converge to zero faster than power functions of a as a. These assumptions have been used successfully in CFM-based analyses for various integer-order nonlinear systems [5053]. Compared to the integer-order systems, the oscillator in Eq. (7) has an additional fractional derivative term, but this term is multiplied by the sufficiently small parameter ε. Thus, it is unlikely that the response of this fractional oscillator is far from the response of the integer-order oscillators. Furthermore, it was recently shown that the stationary amplitude PDF for a fractional linear oscillator under Gaussian white noise (i.e., the special case of Eq. (7) in which h(X,X˙)=0 and g(X)=0) is a Rayleigh distribution [33]. Therefore, the above assumptions about pA(a,t) seem to be valid even for the fractional nonlinear oscillator given by Eq. (7). Under the assumptions, choosing the real part ρ of γ from the range 2<ρ< determined from Eq. (3), the non-integral terms in Eq. (31) disappear, and then, Eq. (31) becomes the governing equation of the time evolution of MpA(γ1,t)
tMpA(γ1,t)=ε2(γ1)j=0nββjMpA(γ1+j,t)+επS02ω02(γ1)2MpA(γ3,t)ε2πS04ω02(γ1)j=0ns(2γ2+j)sj×MpA(γ3+j,t)
(32)

5.2 Solving Complex Fractional Moment Equations.

Unfortunately, Eq. (32) is not closed due to the existence of multiple CFMs with different real parts of order. As a closure technique, an approximate relation between CFMs with different real parts of order is utilized [46]. Let γk(l)=ρl+ikπ/b (l=1,2;k=m,m+1,,m1,m;b and m are defined in Eq. (6)) and the real-part values ρ1=ρ and ρ2=ρ+Δρ. Then, the following approximate relation between CFMs with different real parts of order is established [46]:
MpA(γs(1)1,t)=k=mmMpA(γk(2)1,t)cks(Δρ)cks(Δρ)=sin[π(ks)ibΔρ]π(ks)ibΔρ
(33)
where s=m,,m. Thus, MpA(γs(1)1,t) can be approximately represented by MpA(γk(2)1,t). By utilizing Eq. (33), a closed set of CFM linear equations is obtainable. First, MpA(γs1+j,t), MpA(γs3,t) and MpA(γs3+j,t) are rewritten as the linear combination of MpA(γk1,t)(k=m,,m) using Eq. (33). Next, substituting them into Eq. (32) discretized as γ=γs(=ρ+isπ/b,s=m,,m) leads to (2m+1) linear ordinary differential equations for MpA(γs1,t)(s=m,,m).
tMpA(γs1,t)=ε2(γs1)j=0nββjk=mmMpA(γk1,t)cks(j)+επS02ω02(γs1)2k=mmMpA(γk1,t)cks(2)ε2πS04ω02(γs1)j=0ns(2γs2+j)sj×k=mmMpA(γk1,t)cks(2j)
(34)
In solving Eq. (34), it is necessary to impose a constraint that corresponds to the normalization condition for a PDF, which is given by [46]
MpA(γ01,t)=1ρeb(1ρ)eb(1ρ)[2bk=mk0meb(1γk)eb(1γk)1γkMpA(γk1,t)]
(35)
Combining Eqs. (34) and (35), we anew obtain 2m coupled linear differential equations for MpA(γs1,t)(s=m,,1,1,,m), which are solved with the initial condition given by
MpA(γs1,0)=0pA(a,0)aγs1da=0
(36)
where pA(a,0)=δ^(a).

In this way, the solution of the FP equation (22) is attained in terms of CFMs by solving Eqs. (34)(36). The computation time is only a few seconds. When performing a stationary response analysis, linear algebraic equations derived by setting the time derivative term in Eq. (34) to zero are solved in conjunction with Eq. (35).

5.3 Calculation of Response Probability Density Function.

The amplitude PDF pA(a,t) is acquired by substituting the CFMs MpA(γs1,t)(s=m,,m) determined from Eqs. (34)(36) into Eq. (6). The joint PDF pXX˙(x,x˙,t) of displacement X(t) and velocity X˙(t) is also obtained as
pXX˙(x,x˙,t)=12πpA(a,t)|J|
(37)
where J is the Jacobian of variable transformation pertaining to the pseudo-harmonic representation (10), and |J| is given by [51]
|J|=1ωe(a)adda(1ωe(a))x˙2
(38)
Equations (37) and (38) include a. For this a, substitute the solution for a of the following equation:
a=x2+x˙2ωe2(a)
(39)
By way of example, let us consider the case of g(X)=X3. From Eq. (18), ωe2(A) is derived as
ωe2(A)=ω02(1+34εA2)+ελω0αcos(απ2)
(40)
Then, |J| and a are
|J|=4ωe3(a){4ωe4(a)+3εω02x˙2}aa=23ε{(C+34εx2)2+3εx˙2ω02+34εx2C}
(41)
where
C=1+ελω0α2cos(απ2)

The displacement PDF pX(x,t) and the velocity PDF pX˙(x˙,t) are determined as the marginal PDFs of pXX˙(x,x˙,t).

5.4 Summary of Analytical Procedure.

In Secs. 4 and 5, the derivation process for the equations used in the present analysis was described in detail. However, when applying the present approach, it is not necessary to perform all of the analyses in Secs. 4 and 5 in sequence. The actual analysis consists of fewer equations and procedures. The practical analytical procedure is summarized as follows:

  1. Calculate the equivalent natural frequency ωe(A) and the equivalent damping coefficient βe(A) using Eqs. (18) and (19). ωe2(A) is given in the form of Eqs. (25) and (29), while βe(A) is given in the form of Eq. (30).

  2. Substitute ωe(A) and βe(A) into Eq. (34) to get the linear ordinary differential equations for amplitude CFMs MpA(γs1,t)(s=m,,1,1,,m).

  3. Solve Eq. (34) in combination with Eq. (35) and the initial condition given by Eq. (36) to determine MpA(γs1,t)(s=m,,m).

  4. Apply the inverse Mellin transform to MpA(γs1,t), i.e., substitute MpA(γs1,t) into Eq. (6) to obtain the amplitude PDF pA(a,t).

  5. Calculate a and |J| in the joint PDF pXX˙(x,x˙,t) of X(t) and X˙(t) in Eq. (37) via Eqs. (38) and (39).

  6. Substitute pA(a,t), a, and |J| into Eq. (37) to obtain pXX˙(x,x˙,t).

6 Numerical Examples

The accuracy and versatility of the present analytical approach are assessed through numerical examples of three linear and nonlinear oscillators with fractional derivative elements. To the best of the author’s knowledge, no exact PDF solution has yet been found for oscillators with fractional derivatives, even for stationary responses. For this reason, the analytical results of response PDFs are compared with the results of the pertinent MCS.

6.1 Linear Oscillator With Fractional Derivative Element.

First, to verify the effectiveness of the present approach for dynamic systems involving fractional-order derivatives, a linear oscillator with a fractional derivative element is considered, once ignoring nonlinearities. The equation of motion is given by
X¨+ελD0,tαX+ω02X=ε12W(t)
(42)
that corresponds to the case of h(X,X˙)=0 and g(X)=0 in Eq. (7).
Using Eqs. (18) and (19), ωe(A) and βe(A) are calculated as
ωe2(A)=ωe2=ω02+ελω0αcos(απ2)βe(A)=βe=λω0α1sin(απ2)
(43)
Since there is no nonlinear element in the system, these equivalent coefficients do not depend on the amplitude A(t). The linear differential equations for amplitude CFMs to be solved are derived from Eqs. (34) and (43) as
tMpA(γs1,t)=ελω0α12sin(απ2)(γs1)MpA(γs1,t)+επS02ω02{1ελω0α2cos(απ2)}×(γs1)2k=mmMpA(γk1,t)cks(2)
(44)
where s=m,,1,1,m. Equation (44) is solved in combination with Eq. (35) and the initial condition given by Eq. (36). Then, the discretized inverse Mellin transform (Eq. (6)) of the resulting CFMs MpA(γs1,t)(s=m,,m) yields the amplitude PDF pA(a,t). Further, the joint PDF pXX˙(x,x˙,t) of X(t) and X˙(t) is obtained by Eq. (37). From Eqs. (38) and (39), a and |J| in Eq. (37) are
a=x2+x˙2ωe2;|J|=1ωea
(45)
where ωe is given by Eq. (43).

The reliability of the analytical method is validated with three different values of the fractional derivative order α: α=0.3,0.5, and 0.7. It is known that 0.3α0.7 is frequently adopted as the best fitting value in modeling the dynamic properties of viscoelastic materials [1115]. For this reason, the above three values are used for the fractional order. Other system and excitation parameters are set as ω0=1, ε=0.01, λ=7, and S0=1. Further, Δη=0.3, ρ=3.7, and m=220 are adopted in solving Eq. (44).

For comparisons of the analytical results with simulation results, MCS is performed on Eq. (42). The difficulty inherent in simulations of fractional-order systems is that the fractional derivatives depend on the entire past time history of the system response, as can be seen from Eq. (8). As an efficient and accurate numerical technique, this study utilizes a method to approximately convert an equation of motion including fractional-order derivatives into a set of integer-order differential equations [63]. A brief summary of this methodology is provided in Appendix. Once the integer-order differential equations are obtained, standard numerical integration algorithms can be used to generate realizations of the system response. In this study, the fourth-order Runge–Kutta method is employed. The time increment size Δt for the numerical integration is Δt=104. MCS results are computed from 5×106 realizations. Δt and the number of realizations were chosen based on a preliminary study in order to obtain as accurate a solution as possible, since the MCS results are used as substitutes for the true response PDFs of the oscillator. The computer system used in the MCS has an Intel Core i9-10,900 K CPU with 32 GB RAM, and the operating system is Microsoft Windows 10. The source code for the MCS was prepared using the c language, which was compiled with the GNU compiler collection (version 12.3). The simulation method described herein is also adopted in the second and third examples shown later.

In Fig. 1, the PDFs pA(a,t), pX(x,t), and pX˙(x˙,t) of amplitude A(t), displacement X(t), and velocity X˙(t) for the case of α=0.5 are shown. The solid lines indicate the analytical results, whereas the discrete symbols denote the pertinent MCS results. The three times in the figure were chosen to show the early phase of the response, the stationary phase and a phase in between the two, while also taking into account the ease of viewing the figure. From the standpoint of evaluating system reliability, it is crucial to obtain the tails of response PDFs with high accuracy. Therefore, semi-log plots of pX(x,t) and pX˙(x˙,t) are also presented to confirm the agreement between the analytical and MCS results in the tail region. The results of pA(a,t) demonstrate the validity of the assumptions about pA(a,t) used in deriving Eq. (32) from Eq. (31). Slight deviations can be seen in the peaks of the displacement and velocity PDFs. However, the accuracy of the analytical results is overall quite satisfactory, including the tail region. Although not shown in the figure, good agreement was also achieved for the cases of α=0.3 and α=0.7 (some results are provided in Fig. 2). Additionally, the authors also carefully checked the validity for various other values of α in the range 0<α<1, and the analytical method accurately yielded transient response PDF solutions for any α, including the cases of α=0.1 and 0.9.

Fig. 1
Transient and stationary response PDFs of a linear oscillator with a fractional derivative of order 0.5. Solid lines: present analysis. Discrete symbols: MCS. (a) amplitude, (b) displacement, (c) displacement (semi-log plot), (d) velocity, and (e) velocity (semi-log plot).
Fig. 1
Transient and stationary response PDFs of a linear oscillator with a fractional derivative of order 0.5. Solid lines: present analysis. Discrete symbols: MCS. (a) amplitude, (b) displacement, (c) displacement (semi-log plot), (d) velocity, and (e) velocity (semi-log plot).
Close modal
Fig. 2
Comparison of response PDFs between three types of fractional order α. Solid lines: present analysis. Discrete symbols: MCS. (a) displacement and (b) velocity.
Fig. 2
Comparison of response PDFs between three types of fractional order α. Solid lines: present analysis. Discrete symbols: MCS. (a) displacement and (b) velocity.
Close modal

Figure 2 shows a comparison of displacement and velocity PDFs between the three types of fractional order α. The following tendencies can be observed in the response PDFs for larger α:

  • The widths of the PDFs are narrower (i.e., their variance decreases).

  • They reach stationary PDFs faster.

The reason for the above is that as α approaches 1, the contribution of the fractional derivative element to viscous damping increases. The present method accurately predicts such changes in the response PDF depending on the fractional order α.

Finally, the computation time for the present analytical method was 2.7 s, whereas the MCS took about 25 h (the computation times were similar in the second and third examples shown later as well). Thus, the computational time was significantly reduced compared to the corresponding MCS.

6.2 Oscillator With Fractional Derivative Element and Nonlinearities in Both Damping and Stiffness.

The applicability to systems that have not only fractional derivatives but also nonlinearities is illustrated. Let h(X,X˙)=δ1X˙+δ3X˙3 and g(X)=κ3X3+κ5X5, and a fractional nonlinear oscillator with cubic nonlinearity in damping and cubic and quintic nonlinearities in stiffness is considered, whose motion is governed by
X¨+ελD0,tαX+ε(δ1X˙+δ3X˙3)+ω02X+εω02(κ3X3+κ5X5)=ε12W(t)
(46)
ωe(A) and βe(A) for this system are
ωe2(A)=ω02{1+ε(34κ3A2+58κ5A4)}+ελω0αcos(απ2)βe(A)=δ1+34δ3ω02[A2+ε{λω0α2cos(απ2)A2+34κ3A4+58κ5A6}]+λω0α1sin(απ2)
(47)
The governing equations for amplitude CFMs are written in the form
tMpA(γs1,t)=ε2{δ1+λω0α1sin(απ2)}(γs1)MpA(γs1,t)3εδ3ω028{1+ελω0α2cos(απ2)}(γs1)×k=mmMpA(γk1,t)cks(2)9ε2δ3κ3ω0232(γs1)k=mmMpA(γk1,t)cks(4)15ε2δ3κ5ω0264(γs1)k=mmMpA(γk1,t)cks(6)+επS02ω02{1ελω0α2cos(απ2)}(γs1)2×k=mmMpA(γk1,t)cks(2)3ε2πκ3S08ω02γs(γs1)MpA(γs1,t)5ε2πκ5S016ω02(γs21)k=mmMpA(γk1,t)cks(2)
(48)
where s=m,,1,1,m. Equation (48) is solved with Eqs. (35) and (36). Equations (38) and (47) lead to |J| in the joint PDF pXX˙(x,x˙,t) of Eq. (37) as
|J|=4ωe3(a){4ωe4(a)+3εκ3ω02x˙2+5εκ5ω02x˙2a2}a
(49)
For a in Eqs. (37) and (49), the square root of the solution of the following cubic equation for a2 is used
58εκ5a6+ε4(3κ352κ5x2)a4+(C34εκ3x2)a2(Cx2+x˙2ω02)=0
(50)
where
C=1+ελω0α2cos(απ2)
Equation (50) is derived from Eqs. (39) and (47).

The following parameter values are used: ω0=1, ε=0.01, λ=7, α=0.5, δ1=8, δ3=8, κ3=1, κ5=1, and S0=1. The parameters when solving Eq. (48) are given by Δη=0.5, ρ=3.7, and m=220.

Figure 3 displays the analytical results of transient and stationary response PDFs and the corresponding MCS results. There are small errors in the peaks of the displacement and velocity PDFs for the initial short time period, but otherwise the analytical results are in good agreement with the simulation results, including the tail region. The results of the amplitude PDF pA(a,t) show that the assumptions about pA(a,t) are valid for the system containing nonlinearities as well. Further, a comparison of Figs. 1 and 3 reveals that the response PDFs in this example are narrower than those of the linear oscillator with same α. This is because the nonlinear damping and nonlinear stiffness of the system behave to suppress large responses. The analytical result accurately captures such a narrowed PDF shape caused by the system nonlinearities. Thus, the present analytical method is able to provide accurate solutions even for a system with nonlinear damping and nonlinear stiffness as well as a fractional derivative element.

Fig. 3
Transient and stationary response PDFs of a oscillator with a fractional derivative element and nonlinearities in both damping and stiffness. Solid lines: present analysis. Discrete symbols: MCS. (a) amplitude, (b) displacement, (c) displacement (semi-log plot), (d) velocity, and (e) velocity (semi-log plot).
Fig. 3
Transient and stationary response PDFs of a oscillator with a fractional derivative element and nonlinearities in both damping and stiffness. Solid lines: present analysis. Discrete symbols: MCS. (a) amplitude, (b) displacement, (c) displacement (semi-log plot), (d) velocity, and (e) velocity (semi-log plot).
Close modal

6.3 Duffing-van der Pol Oscillator With Fractional Derivative Element.

Finally, an example of a fractional nonlinear oscillator that exhibits limit-cycle behavior is presented in order to further demonstrate the versatility of the present method. Consider a Duffing-van der Pol oscillator with a fractional derivative element, whose equation of motion is described by
X¨+ελD0,tαX+ε(δ1+δ3X2)X˙+ω02X+εω02X3=ε12W(t)
(51)
which corresponds to the case of h(X,X˙)=δ1X˙+δ3X2X,˙ and g(X)=X3 in Eq. (7).
For this system, ωe(A) and βe(A) are given by
ωe2(A)=ω02(1+34εA2)+ελω0αcos(απ2)βe(A)=δ1+14δ3A2+λω0α1sin(απ2)
(52)
The differential equations for amplitude CFMs are
tMpA(γs1,t)=ε2{δ1+λω0α1sin(απ2)}(γs1)MpA(γs1,t)εδ38(γs1)k=mmMpA(γk1,t)cks(2)+επS02ω02{1ελω0α2cos(απ2)}(γs1)2×k=mmMpA(γk1,t)cks(2)3ε2πS08ω02γs(γs1)MpA(γs1,t)
(53)
where s=m,,1,1,m. Equation (53) is solved with Eqs. (35) and (36). Since the nonlinear stiffness term of this system is g(X)=X3, a and |J| in the joint PDF pXX˙(x,x˙,t) of Eq. (37) are given by Eq. (41).

The linear damping coefficient δ1 is set as δ1=5 so that the oscillator (51) takes diffusive limit-cycle motion. The other parameter values are chosen as ω0=1, ε=0.01, λ=1, α=0.5, δ3=1, and S0=1. Further, Δη=0.5, ρ=3.7, and m=220 are used when solving Eq. (53).

The analytical results of response PDFs and the corresponding MCS results for comparison are depicted in Fig. 4. The oscillator behavior approaches diffusive limit-cycle motion with time. Accordingly, the displacement and velocity PDFs change from unimodal shapes to bimodal ones. It can be seen that the present analytical technique estimates such highly non-Gaussian PDFs and the time evolution of their shapes with high accuracy. Furthermore, due to the behavior of diffusive limit cycle, the amplitude PDF also has a unique shape that is different from the PDFs in the previous two examples. Nevertheless, the assumptions about the amplitude PDF are still valid in this example.

Fig. 4
Transient and stationary response PDFs of a Duffing-van der Pol oscillator with a fractional derivative element. Solid lines: present analysis. Discrete symbols: MCS. (a) amplitude, (b) displacement, (c) displacement (semi-log plot), (d) velocity, and (e) velocity (semi-log plot).
Fig. 4
Transient and stationary response PDFs of a Duffing-van der Pol oscillator with a fractional derivative element. Solid lines: present analysis. Discrete symbols: MCS. (a) amplitude, (b) displacement, (c) displacement (semi-log plot), (d) velocity, and (e) velocity (semi-log plot).
Close modal

7 Conclusions

A CFM-based analysis has been performed for obtaining the transient response PDF of nonlinear oscillators with fractional derivative elements under Gaussian white noise. The oscillators have polynomial-type nonlinearities in both stiffness and damping.

First, the equivalent natural frequency and equivalent damping coefficient of the oscillator were determined. The fractional derivative element of the oscillator contributes to both elasticity and viscosity. For this reason, the equivalent coefficients need to be determined by taking account of not only the effects of nonlinearities but also the contributions of the fractional element. Furthermore, the coefficients must be given in polynomial form of response amplitude to utilize the CFM-based approach. This paper proposed formulas to determine the equivalent coefficients that satisfy the above requirements, based on an equivalent linearization technique. Then, the FP equation was derived using stochastic averaging. Further, by the use of the Mellin transform and the relation between CFMs with different real parts of order, the FP equation was converted into a set of linear ordinary differential equations for amplitude CFMs. Solving it combined with a constraint that corresponds to the normalization condition for a PDF, the amplitude CFMs were obtained. Finally, the inverse Mellin transform of the resulting CFMs yielded the amplitude PDF. The joint PDF of displacement and velocity was also attained from the amplitude PDF.

In numerical examples, three linear and nonlinear oscillators with fractional derivatives have been considered to demonstrate the accuracy and versatility of the present analytical technique. It has been shown that the present method accurately determines the shapes of response PDFs that vary depending on both the fractional order and nonlinearities of the oscillator.

The practical importance of acquiring response PDFs is that system reliability can be assessed from the PDF. Therefore, it is highly significant to perform reliability analysis of fractional nonlinear systems using response PDF solutions obtained by the present analytical technique. Such a study would also be an effective and practical way to quantitatively evaluate the accuracy of the present analytical method. This topic will be addressed in our future work.

Acknowledgment

This work was supported by JSPS KAKENHI Grant No. JP21K03944.

Conflict of Interest

The authors have no competing interests to declare that are relevant to the content of this paper.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Appendix

This appendix provides a brief summary of a method to approximately convert an equation of motion including fractional derivative elements into a set of integer-order differential equations, which is used in MCS. Readers interested in more details of the methodology are referred to Ref. [63].

The gamma function Γ(α) is defined as
Γ(α)=0ezzα1dz
(A1)
It is known that Γ(α) satisfies the following identity:
Γ(α)Γ(1α)=πsin(πα)
(A2)
By employing Eqs. (A1) and (A2), the α-order Caputo fractional derivative (8) can be rewritten in the form
D0,tαX(t)=sin(απ)π0t0ezzα1(ts)αX˙(s)dzds
(A3)
Performing a variable transformation z=(ts)y2, we obtain
D0,tαX(t)=2sin(απ)π0y2α1{0te(ts)y2X˙(s)ds}dy
(A4)
Herein, the following variable is introduced:
Φ(y,t)=y2α10te(ts)y2X˙(s)ds
(A5)
Then, Eq. (A4) is written as
D0,tαX(t)=μ0Φ(y,t)dy
(A6)
where μ=2sin(πα)/π. Further, from Eq. (A5), Φ(y,t) is the solution of the following first-order ordinary differential equation:
ddtΦ(y,t)=y2Φ(y,t)+y2α1X˙(t)
(A7)
with zero initial condition Φ(y,0)=0.
In order to obtain an efficient and accurate numerical scheme, the integral in Eq. (A6) is discretized by utilizing the Gauss–Jacobi quadrature. To this aim, the following variable transformation is used for Eq. (A6):
u=1y1+y
(A8)
Then, Eq. (A6) becomes
D0,tαX(t)=μ11(1u)2α1(1+u)12αΦ¯(u,t)du
(A9)
where
Φ¯(u,t)=2(1u)12α(1+u)2α3Φ(1u1+u,t)
(A10)
Next, the Gauss–Jacobi quadrature is applied to Eq. (A9).
D0,tαX(t)=μ11(1u)2α1(1+u)12αΦ¯(u,t)duμi=1nwiΦ¯(ui,t)
(A11)
where ui is the nodal point determined as the ith root of the Jacobi polynomial Pn(2α1,12α)(u), whose definition is
Pn(2α1,12α)(u)=(1)n2nn!(1u)2α1(1+u)12α×dndun[(1u)2α1+n(1+u)12α+n]
(A12)
wi is the weight computed as
wi=2Γ(n+2α)Γ(n+22α)Γ(n+1)(n+1)!Pn(2α1,12α)(ui)Pn+1(2α1,12α)(ui)
(A13)
where Pn(2α1,12α)(ui) denotes the first-order derivative of Pn(2α1,12α)(ui). From Eqs. (A10) and (A7), we have
Φ¯(ui,t)=2(1ui)12α(1+ui)2α3Φ(1ui1+ui,t)
(A14)
ddtΦ(yi,t)=yi2Φ(yi,t)+yi2α1X˙(t),yi=1ui1+ui
(A15)
To summarize the above, Eq. (7) can be approximated as follows:
X¨+ελμi=1n2wi(1ui)12α(1+ui)2α3Φ(yi,t)+εh(X,X˙)+ω02X(t)+εω02g(X)=ε12W(t)
(A16)
where the time evolution of Φ(yi,t) is governed by Eq. (A15). Further, Eq. (A16) is recast in the form
dXdt=X˙
(A17)
dX˙dt=ελμi=1n2wi(1ui)12α(1+ui)2α3Φ(yi,t)εh(X,X˙)ω02X(t)εω02g(X)+ε12W(t)
(A18)
In this manner, the system equation of motion with a fractional derivative term is approximately converted into a set of (n+2) first-order ordinary differential equations consisting of Eqs. (A15), (A17), and (A18). In this study, n was selected as n=10, and thus, 12 simultaneous equations were solved using the fourth-order Runge–Kutta scheme.

References

1.
Lin
,
Y. K.
, and
Cai
,
G. Q.
,
1995
,
Probabilistic Structural Dynamics
,
McGraw-Hill
,
New York
.
2.
Podlubny
,
I.
,
1999
,
Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications
,
Academic Press
,
San Diego, CA
.
3.
Hilfer
,
R.
,
2000
,
Applications Of Fractional Calculus in Physics
,
World Scientific
,
Singapore
.
4.
Sabatier
,
J.
,
Agrawal
,
O. P.
, and
Tenreiro Machado
,
J. A.
,
2000
,
Advances in Fractional Calculus: Theoretical Developments and Applications in Physics and Engineering
,
Springer
,
Dordrecht
.
5.
Mainardi
,
F.
,
2010
,
Fractional Calculus and Waves in Linear Viscoelasticity: An Introduction to Mathematical Models
,
Imperial College Press
,
London
.
6.
Atanackovic
,
T. M.
,
Pilipovic
,
S.
,
Stankovic
,
B.
, and
Zorica
,
D.
,
2010
,
Fractional Calculus With Applications in Mechanics: Vibrations and Diffusion Processes
,
Wiley
,
London
.
7.
Herrmann
,
R.
,
2018
,
Fractional Calculus: An Introduction for Physicists
, 3rd ed.,
World Scientific
,
Singapore
.
8.
Rossikhin
,
Y. A.
, and
Shitikova
,
M. V.
,
1997
, “
Applications of Fractional Calculus to Dynamic Problems of Linear and Nonlinear Hereditary Mechanics of Solids
,”
ASME Appl. Mech. Rev.
,
50
(
1
), pp.
15
67
.
9.
Rossikhin
,
Y. A.
and
Shitikova
,
M. V.
,
2010
, “
Application of Fractional Calculus for Dynamic Problems of Solid Mechanics: Novel Trends and Recent Results
,”
ASME Appl. Mech. Rev.
,
63
(
1
), p.
010801
.
10.
Sun
,
H. G.
,
Zhang
,
Y.
,
Baleanu
,
D.
,
Chen
,
W.
, and
Chen
,
Y. Q.
,
2018
, “
A New Collection of Real World Applications of Fractional Calculus in Science and Engineering
,”
Commun. Nonlinear Sci. Numer. Simul.
,
64
(
1
), pp.
213
231
.
11.
Jones
,
D. I. G.
,
2001
,
Handbook of Viscoelastic Vibration Damping
,
Wiley
,
New York
.
12.
Bagley
,
R. L.
, and
Torvik
,
P. J.
,
1983
, “
Fractional Calculus – A Different Approach to the Analysis of Viscoelastically Damped Structures
,”
AIAA J.
,
21
(
5
), pp.
741
748
.
13.
Torvik
,
P. J.
, and
Bagley
,
R. L.
,
1984
, “
On the Appearance of the Fractional Derivative in the Behavior of Real Materials
,”
ASME J. Appl. Mech.
,
51
(
2
), pp.
294
298
.
14.
Bagley
,
R. L.
, and
Torvik
,
P. J.
,
1985
, “
Fractional Calculus in the Transient Analysis of Viscoelastically Damped Structures
,”
AIAA J.
,
23
(
6
), pp.
918
925
.
15.
Bagley
,
R. L.
, and
Torvik
,
P. J.
,
1986
, “
On the Fractional Calculus Model of Viscoelastic Behavior
,”
J. Rheol.
,
30
(
1
), pp.
133
155
.
16.
Koh
,
C. G.
, and
Kelly
,
J. M.
,
1990
, “
Application of Fractional Derivatives to Seismic Analysis of Base-Isolated Models
,”
Earthq. Eng. Struct. Dyn.
,
19
(
2
), pp.
229
241
.
17.
Sasso
,
M.
,
Palmieri
,
G.
, and
Amodio
,
D.
,
2011
, “
Application of Fractional Derivative Models in Linear Viscoelastic Problems
,”
Mech. Time-Depend. Mater.
,
15
(
1
), pp.
367
387
.
18.
Di Paola
,
M.
,
Pirrotta
,
A.
, and
Valenza
,
A.
,
2011
, “
Visco-Elastic Behavior Through Fractional Calculus: An Easier Method for Best Fitting Experimental Results
,”
Mech. Mater.
,
43
(
12
), pp.
799
806
.
19.
Makris
,
N.
, and
Constantinou
,
M. C.
,
1992
, “
Spring-Viscous Damper Systems for Combined Seismic and Vibration Isolation
,”
Earthq. Eng. Struct. Dyn.
,
21
(
8
), pp.
649
664
.
20.
Lee
,
H. H.
, and
Tsai
,
C.-S.
,
1994
, “
Analytical Model of Viscoelastic Dampers for Seismic Mitigation of Structures
,”
Comput. Struct.
,
50
(
1
), pp.
111
121
.
21.
Shen
,
K. L.
, and
Soong
,
T. T.
,
1995
, “
Modeling of Viscoelastic Dampers for Structural Applications
,”
J. Eng. Mech.
,
121
(
6
), pp.
694
701
.
22.
Rüdinger
,
F.
,
2006
, “
Tuned Mass Damper With Fractional Derivative Damping
,”
Eng. Struct.
,
28
(
13
), pp.
1774
1779
.
23.
Singh
,
M. P.
,
Chang
,
T.-S.
, and
Nandan
,
H.
,
2011
, “
Algorithms for Seismic Analysis of MDOF Systems With Fractional Derivatives
,”
Eng. Struct.
,
33
(
8
), pp.
2371
2381
.
24.
Agrawal
,
O. P.
,
2001
, “
Stochastic Analysis of Dynamic Systems Containing Fractional Derivatives
,”
J. Sound Vib.
,
247
(
5
), pp.
927
938
.
25.
Huang
,
Z. L.
,
Jin
,
X. L.
,
Lim
,
C. W.
, and
Wang
,
Y.
,
2010
, “
Statistical Analysis for Stochastic Systems Including Fractional Derivatives
,”
Nonlinear Dyn.
,
59
(
1
), pp.
339
349
.
26.
Cao
,
Q.
,
Hu
,
S.-L. J.
, and
Li
,
H.
,
2021
, “
Nonstationary Response Statistics of Fractional Oscillators to Evolutionary Stochastic Excitation
,”
Commun. Nonlinear Sci. Numer. Simul.
,
103
(
1
), p.
105962
.
27.
Di Paola
,
M.
,
Failla
,
G.
, and
Pirrotta
,
A.
,
2012
, “
Stationary and Non-Stationary Stochastic Response of Linear Fractional Viscoelastic Systems
,”
Probab. Eng. Mech.
,
28
(
1
), pp.
85
90
.
28.
Spanos
,
P. D.
, and
Zhang
,
W.
,
2022
, “
Nonstationary Stochastic Response Determination of Nonlinear Oscillators Endowed With Fractional Derivatives
,”
Int. J. Non-Linear Mech.
,
146
(
1
), p.
104170
.
29.
Zhang
,
W.
,
Spanos
,
P. D.
, and
Di Matteo
,
A.
,
2023
, “
Nonstationary Stochastic Response of Hysteretic Systems Endowed With Fractional Derivative Elements
,”
ASME J. Appl. Mech.
,
90
(
6
), p.
061011
.
30.
Spanos
,
P. D.
, and
Evangelatos
,
G. I.
,
2010
, “
Response of a Non-Linear System With Restoring Forces Governed by Fractional Derivatives - Time Domain Simulation and Statistical Linearization Solution
,”
Soil Dyn. Earthq. Eng.
,
30
(
9
), pp.
811
821
.
31.
Kong
,
F.
,
Zhang
,
H.
,
Zhang
,
Y.
,
Chao
,
P.
, and
He
,
W.
,
2022
, “
Stationary Response Determination of MDOF Fractional Nonlinear Systems Subjected to Combined Colored Noise and Periodic Excitation
,”
Commun. Nonlinear Sci. Numer. Simul.
,
110
(
1
), p.
106392
.
32.
Kong
,
F.
,
Han
,
R.
, and
Zhang
,
Y.
,
2022
, “
Approximate Stochastic Response of Hysteretic System With Fractional Element and Subjected to Combined Stochastic and Periodic Excitation
,”
Nonlinear Dyn.
,
107
(
1
), pp.
375
390
.
33.
Spanos
,
P. D.
,
Kougioumtzoglou
,
I. A.
,
dos Santos
,
K. R. M
, and
Beck
,
A. T.
,
2018
, “
Stochastic Averaging of Nonlinear Oscillators: Hilbert Transform Perspective
,”
J. Eng. Mech.
,
144
(
2
), p.
04017173
.
34.
Huang
,
Z. L.
, and
Jin
,
X. L.
,
2009
, “
Response and Stability of a SDOF Strongly Nonlinear Stochastic System With Light Damping Modeled by a Fractional Derivative
,”
J. Sound Vib.
,
319
(
3
), pp.
1121
1135
.
35.
Hu
,
F.
,
Chen
,
L. C.
, and
Zhu
,
W. Q.
,
2012
, “
Stationary Response of Strongly Non-Linear Oscillator With Fractional Derivative Damping Under Bounded Noise Excitation
,”
Int. J. Non-Linear Mech.
,
47
(
10
), pp.
1081
1087
.
36.
Yang
,
Y.
,
Xu
,
W.
,
Gu
,
X.
, and
Sun
,
Y.
,
2015
, “
Stochastic Response of a Class of Self-Excited Systems With Caputo-Type Fractional Derivative Driven by Gaussian White Noise
,”
Chaos Solit. Fractals
,
77
(
1
), pp.
190
204
.
37.
Wang
,
D.
,
Xu
,
W.
,
Gu
,
X.
, and
Pei
,
H.
,
2016
, “
Response Analysis of Nonlinear Vibro-Impact System Coupled With Viscoelastic Force Under Colored Noise Excitations
,”
Int. J. Non-Linear Mech.
,
86
(
1
), pp.
55
65
.
38.
Xiao
,
Y.
,
Xu
,
W.
, and
Wang
,
L.
,
2016
, “
Stochastic Responses of Van Der Pol Vibro-Impact System With Fractional Derivative Damping Excited by Gaussian White Noise
,”
Chaos
,
26
(
1
), p.
033110
.
39.
Ning
,
X.
,
Ma
,
Y.
,
Li
,
S.
,
Zhang
,
J.
, and
Li
,
Y.
,
2018
, “
Response of Non-Linear Oscillator Driven by Fractional Derivative Term Under Gaussian White Noise
,”
Chaos Solit. Fractals
,
113
(
1
), pp.
102
107
.
40.
Yang
,
Y.-G.
, and
Xu
,
W.
,
2018
, “
Stochastic Analysis of Monostable Vibration Energy Harvesters With Fractional Derivative Damping Under Gaussian White Noise Excitation
,”
Nonlinear Dyn.
,
94
(
1
), pp.
639
648
.
41.
Sun
,
Y.-H.
,
Yang
,
Y.-G.
,
Zhang
,
Y.
, and
Xu
,
W.
,
2021
, “
Probabilistic Response of a Fractional-Order Hybrid Vibration Energy Harvester Driven by Random Excitation
,”
Chaos
,
31
(
1
), p.
013111
.
42.
Di Matteo
,
A.
,
2023
, “
Response of Nonlinear Oscillators With Fractional Derivative Elements Under Evolutionary Stochastic Excitations: A Path Integral Approach Based on Laplace’s Method of Integration
,”
Probab. Eng. Mech.
,
71
(
1
), p.
103402
.
43.
Fragkoulis
,
V. C.
,
Kougioumtzoglou
,
I. A.
,
Pantelous
,
A. A.
, and
Beer
,
M.
,
2019
, “
Non-Stationary Response Statistics of Nonlinear Oscillators With Fractional Derivative Elements Under Evolutionary Stochastic Excitation
,”
Nonlinear Dyn.
,
97
(
1
), pp.
2291
2303
.
44.
Cottone
,
G.
,
Di Paola
,
M.
, and
Metzler
,
R.
,
2010
, “
Fractional Calculus Approach to the Statistical Characterization of Random Variables and Vectors
,”
Physica A
,
389
(
5
), pp.
909
920
.
45.
Di Paola
,
M.
, and
Pinnola
,
F. P.
,
2012
, “
Riesz Fractional Integrals and Complex Fractional Moments for the Probabilistic Characterization of Random Variables
,”
Probab. Eng. Mech.
,
29
(
1
), pp.
149
156
.
46.
Di Paola
,
M.
,
2014
, “
Fokker Planck Equation Solved in Terms of Complex Fractional Moments
,”
Probab. Eng. Mech.
,
38
(
1
), pp.
70
76
.
47.
Di Matteo
,
A.
,
Di Paola
,
M.
, and
Pirrotta
,
A.
,
2014
, “
Probabilistic Characterization of Nonlinear Systems Under Poisson White Noise Via Complex Fractional Moments
,”
Nonlinear Dyn.
,
77
(
1
), pp.
729
738
.
48.
Alotta
,
G.
, and
Di Paola
,
M.
,
2015
, “
Probabilistic Characterization of Nonlinear Systems Under α-stable White Noise Via Complex Fractional Moments
,”
Physica A
,
420
(
1
), pp.
265
276
.
49.
Itoh
,
D.
,
Tsuchida
,
T.
, and
Kimura
,
K.
,
2018
, “
An Analysis of a Nonlinear System Excited by Combined Gaussian and Poisson White Noises Using Complex Fractional Moments
,”
Theor. and Appl. Mech. Japan
,
64
(
1
), pp.
103
114
.
50.
Jin
,
X.
,
Wang
,
Y.
,
Huang
,
Z.
, and
Di Paola
,
M.
,
2014
, “
Constructing Transient Response Probability Density of Non-Linear System Through Complex Fractional Moments
,”
Int. J. Nonlinear Mech.
,
65
(
1
), pp.
253
259
.
51.
Itoh
,
D.
, and
Tsuchida
,
T.
,
2022
, “
Transient Response Analysis of a System With Nonlinear Stiffness and Nonlinear Damping Excited by Gaussian White Noise Based on Complex Fractional Moments
,”
Acta Mech.
,
233
(
1
), pp.
2781
2796
.
52.
Xie
,
X.
,
Li
,
J.
,
Liu
,
D.
, and
Guo
,
R.
,
2017
, “
Transient Response of Nonlinear Vibro-Impact System Under Gaussian White Noise Excitation Through Complex Fractional Moments
,”
Acta Mech.
,
228
(
1
), pp.
1153
1163
.
53.
Niu
,
L.
,
Xu
,
W.
, and
Guo
,
Q.
,
2021
, “
Transient Response of the Time-Delay System Excited by Gaussian Noise Based on Complex Fractional Moments
,”
Chaos
,
31
(
1
), p.
053111
.
54.
Dalzell
,
J. E.
,
1978
, “
A Note on the Form of Ship Roll Damping
,”
J. Ship Res.
,
22
(
1
), pp.
178
185
.
55.
Muscolino
,
G.
,
Ricciardi
,
G.
, and
Vasta
,
M.
,
1997
, “
Stationary and Non-Stationary Probability Density Function for Non-Linear Oscillators
,”
Int. J. Nonlinear Mech.
,
32
(
6
), pp.
1051
1064
.
56.
Roberts
,
J. B.
, and
Spanos
,
P. D.
,
2003
,
Random Vibration and Statistical Linearization
,
Dover Publications
,
New York
.
57.
Spanos
,
P. D.
,
Di Matteo
,
A.
,
Cheng
,
Y.
,
Pirrotta
,
A.
, and
Li
,
J.
,
2016
, “
Galerkin Scheme-Based Determination of Survival Probability of Oscillators With Fractional Derivative Elements
,”
ASME J. Appl. Mech.
,
83
(
12
), p.
121003
.
58.
Di Matteo
,
A.
,
Spanos
,
P. D.
, and
Pirrotta
,
A.
,
2018
, “
Approximate Survival Probability Determination of Hysteretic Systems With Fractional Derivative Elements
,”
Probab. Eng. Mech.
,
54
(
1
), pp.
138
146
.
59.
Roberts
,
J. B.
, and
Spanos
,
P. D.
,
1986
, “
Stochastic Averaging: An Approximate Method of Solving Random Vibration Problems
,”
Int. J. Nonlinear Mech.
,
21
(
2
), pp.
111
134
.
60.
Stratonovich
,
R. L.
,
1986
,
Topics in the Theory of Random Noise, Vols. 1 and 2
,
Gordon & Breach
,
New York
.
61.
Spanos
,
P.-T. D.
,
1976
, “Linearization Techniques for Non-Linear Dynamical Systems,” Earthquake Engineering Research Laboratory, California Institute of Technology, Report No. EERL 7Q-04.
62.
Iwan
,
W. D.
, and
Spanos
,
P.-T. D.
,
1978
, “
Response Envelope Statistics for Nonlinear Oscillators With Random Excitation
,”
ASME J. Appl. Mech.
,
45
(
1
), pp.
170
174
.
63.
Diethelm
,
K.
,
2008
, “
An Investigation of Some Nonclassical Methods for the Numerical Approximation of Caputo-Type Fractional Derivatives
,”
Numer. Algor.
,
47
(
1
), pp.
361
390
.