Subharmonic Response of Linear Vibroimpact System with Fractional Derivative Damping to a Randomly Disrobed Periodic Excitation

The subharmonic response of single-degree-of-freedom vibroimpact oscillator with fractional derivative damping and one-sided barrier under narrow-band random excitation is investigated. With the help of a special Zhuravlev transformation, the system is reduced to one without impacts, thereby permitting the applications of asymptotic averaging over the period for slowly varying random process. The analytical expression of the response amplitude is obtained in the case without random disorder, while only the approximate analytical expressions for the steady-statemoments of the response amplitude are obtained in the casewith random disorder. The effects of the fractional order derivative term, damping term, random disorder, and the coefficient of restitution and other system parameters on the system response are discussed.Theoretical analyses and numerical simulations show that fractional derivative makes both the system damping and stiffness coefficients increase, such that it changes the system parameters region at which the response amplitude reaches the maximum. The system energy loss in collision is equivalent to increasing the damping coefficient of the system. System response amplitude will increase when the excitation frequency is close to the resonant frequency and will decay rapidly when the excitation frequency gradually deviates from the resonance frequency.


Introduction
Vibroimapct systems, that is, vibrations systems including impact interaction, widely exist in mechanical engineering, ocean engineering, and civil engineering [1,2].Various dynamical phenomena are observed in the vibroimpact systems.For example, under deterministic periodical loadings, the stability, bifurcation, and chaos of the vibroimpact system have been explored by the Poincare map [3,4].In practice, engineering structures are often subjected to random loadings, for example, wind load, earthquakes, ocean waves, and random disturbance or noise.The behavior of vibroimpact systems has been widely studied in the case of both deterministic and random excitations.In the case of random excitation, many researches have been conducted and reported.For example, Feng and He [5] explored the mean response of impact systems by introducing the mean Poincare map.Namachchivaya and Park [6] developed an averaging approach to study the dynamic actions of a vibroimpact system excited by random perturbations.With the help of the nonsmooth transformation proposed by Zhuravlev [7], the stochastic averaging method has been extended to study the stochastic response of the vibroimpact systems in [8][9][10].Numerical simulation method is also used by Li et al. [10] and Iourtchenko and Song [11].Dimentberg and Iourtchenko presented a comprehensive review on random vibration with impact [12], in which the authors have tried to review and summarize the existing methods, results, and literatures available for solving problem of stochastic vibroimpact systems.
On the other hand, fractional calculus is shown to be very suitable for describing the constitutive relationship of materials with frequency-dependent damping behaviors.In this regard, Gemant [13] provided a pioneering work of this field, which first uses fractional calculus to study phenomenological constitutive equations for material behavior.The fractional derivatives have a wide range of applications in high-energy physics, anomalous diffusion, viscoelastic 2 Mathematical Problems in Engineering mechanical constitutive relations of materials, system control, rheology, geophysics, biomedical engineering, economics, and other fields; the research of theory and application has become a hot topic in recent decades [14,15].Rossikhin has provided an excellent review of the research in this field.Response problem under random excitation with fractional derivative damping system also has many researches; some analyses methods, for example, Laplace transform method used by Agrawal [16], Fourier transform method used by Ye et al. [17], residue harmonic balance method proposed by Leung et al. [18], multiple scales method proposed by Liu et al. [19], and widely used stochastic averaging method [20][21][22], have been developed.
However, the study on response of vibroimpact system under random excitation with fractional derivative damping has not yet been seen according to the author's knowledge, due to complexity of the impact, fractional derivative and random noise being considered simultaneously.In this paper, the subharmonic response of single-degree-of-freedom linear vibroimpact oscillator with fractional derivative damping and one-sided barrier, which is slightly offset from the system's static equilibrium position, is investigated.The system is excited by a sinusoidal force with disorder or random phase modulation.The impact considered here is an instantaneous impact with restitution factor .The paper is organized as follows.In Section 2, the Zhuravlev transformation and stochastic averaging method are used to obtain the averaging equations.In Section 3, the averaged equations are solved exactly and the analytical expression of the response amplitude is obtained.In Section 4, the directly numerical simulations verify the analytical result.Conclusions and discussion are presented in Section 5.

System Description and Averaging Method
Considering a single-degree-of-freedom linear vibroimpact with fractional derivative damping oscillator to random excitations, where dot indicates differentiation with respect to time ; and ÿ , ẏ , and  are acceleration, velocity, and displacement of the oscillator, respectively. is the damping coefficient, Ω is the natural frequency, ℎ ≥ 0 represents the distance from the system's static equilibrium position to the single rigid barrier, and 0 <  ≤ 1 is the restitution factor to be a known parameter of impact losses, whereas subscripts "minus" and "plus" refer to values of response velocity just before and after the instantaneous impact.Thus ẏ + and ẏ − are actually rebound and impact velocities of the mass, respectively.They have the same magnitude whenever  = 1; therefore, this special case is that of elastic impacts, whereas in case  < 1, some impact losses are observed.   ( > 0) is the viscoelastic damping which has the characteristics of memory, in which    denotes the fractional derivative.For fractional derivative [14,15], there are many kinds of definitions.From a physical point of view Caputo's fractional derivative may be more suitable than a Riemann-Liouville one, since it appears for the viscoelastic component of the oscillator.However, from a mathematical point of view the two definitions coalesce each other for the cases, the steady-state solutions and quiescent systems at  = 0, and the Riemann-Liouville definition is more easy to do mathematical derivation and therefore is chosen in this investigation, which can be written as () is a random process governed by the following equation: where  > 0 and Ω 1 > 0 are the amplitude and frequency of the random excitation, respectively, and () is a stationary Gaussian white noise of unit intensity, which describes random temporal deviations of the excitation frequency from its expected or mean Ω 1 .() is a boundary noise since |()| ≤ .Compared with the classical white noise, boundary noise is an ideal model of the actual noise and has attracted many scholars' attention.The process () has the following power spectral density [23]: Obviously, if  → 0, the power spectrum   () vanishes in the entire frequency range except at the singular frequency  = ±Ω 1 , where   () goes to infinity, and thus this process will be assumed to be narrow-band in this case.In this paper, we restrict our attention to the case when  is sufficiently small so that () is bound to a narrow-band random process.Following Zhuravlev [7], the nonsmooth transformation of state variables is introduced as follows: where sgn  is the signum function and Obviously, this transformation makes the transformed velocity ẋ continuous at the impact instants (i.e.,  = 0) in the special case of elastic impact (i.e.,  = 1), thereby reducing the problem to one without velocity jumps.However, this is not the case with a general vibroimpact system with impact losses; the jump of the transformed velocity ẋ becomes proportional to 1 −  instead of 1 +  for the jump of original velocity ẏ .This jump may be included in the transformed differential equation of motion by using the Dirac delta function ().Since ( * ) = 0 at the impact instant  * and ( −  * ) = | ẋ |(), the impulsive term can be obtained as The transformed equation of motion can be written by substituting ( 5) into (1) as Thus, the original impact system (1) is reduced to the "common" vibration system (8) without impact.The term (1 − ) ẋ | ẋ |() on the right hand side of ( 8) describes the impact losses of system, which can be regarded as an impulsive damping term.For simplicity it is assumed that system (1) approximates elastic impact system (i.e.,  ≈ 1), and all coefficients , , , ℎ, 1 −  in the right hand side of ( 8) are small and proportional to a small parameter, such that the transformed equation ( 8) permits rigorous analytical study by the asymptotic method of averaging over the period.Moreover, only subharmonic resonant responses will be considered; that is, frequency Ω 1 of the random excitation is near the subharmonic resonant responses 2Ω, Ω 1 ≈ 2Ω, where  is an arbitrary positive integer.The detuning parameter  is defined according to  = Ω 1 − 2Ω, and  is also assumed to be small and proportional to a small parameter.Then the response of ( 8) can be approximated, represented as By introducing a new slowly varying phase shift () = () − 2Φ(), (8) can be transformed into the following pair of firstorder equations: Under the foregoing assumption that , , , , , ℎ, (1−) are small parameters, it can be known from ( 10) that Φ ≈ Ω, Ȧ ≈ 0, and φ ≈ 0, then Φ is a fast varying random process with respect to time , and  and  are two slowly varying random processes.By averaging over the fast state variable Φ [24], the shortened equations of the slowly varying random processes  and  can be obtained.One now discusses the averaging of the first formula in (10).The averaging over the first three terms in the right hand side of the first formula in ( 10) is The averaging of fourth term is The averaging of fifth term is By using the subsection integral method (13) can be transformed into the following equation: Since Φ is a fast varying random process with respect to time , Φ ≈ Ω, we can obtain the approximate relation if  is small; that is, Substituting ( 15) into ( 14), one obtains To simplify ( 16), the following approximation of integrals [25] are proposed: Substituting ( 17) into ( 16), ( 16) can be simplified as follows: Then the averaging of the first formula in (10) can be obtained by adding (11), (12), and (18).The averaging of the second formula in (10) can also be obtained by using similar derivation.This leads to a pair of "shortened" equations: where The Krylov-Bogoliubov average method [24] is used in the derivation from (10) to (19), and similar skill is used by Dimentberg et al. [26] and the authors of this paper [27] in the studying of the response of the vibroimpact system.Equation ( 19) becomes a Markovian one, where the presence of the fractional operator in (1) destroys the Markovianity.In fact, all coefficients , , , ℎ, 1 −  in the right hand side of ( 8) are assumed to be small and proportional to a small parameter, so (1) can be taken approximately as a Markovian one.However, when these coefficients are not small, much effort should be done.
It can be seen from ( 19) and ( 20) that the difference between elastic impact ( = 1) and inelastic impact ( < 1) is that inelastic impact increases the damping of the system in the term ((1 − )/)Ω.The effects of the fractional derivative    are embodied in two aspects: first it increases the damping of the system in term (/2)Ω −1 sin(/2) > 0 (0 <  ≤ 1), and the second one is that it decreases the detuning parameter from  = Ω 1 − 2Ω to  =  − Ω   cos(/2); therefore it increases natural frequency of the system in term (1/2)Ω   cos(/2) according to the definition of .Namely, the fractional order derivative term increases both damping and stiffness of the system, which is consistent with the results on the fractional derivative action obtained by Chen et al. [21] in studying the response of a Duffing oscillator with fractional derivative by using the stochastic averaging method.In fact, consider two extreme cases  = 0, 1, and in these cases the fractional derivative    changes to the traditional derivative  0  = ,  1  =  ẏ according to [14].Therefore the fractional derivative    increases the damping of the system in term /2 = (/2)Ω −1 sin(/2) > 0,  = 1 in the case  = 1, while increasing stiffness of the system from Ω 2 to Ω 2 + , that is, increasing the natural frequency of the system from Ω to √Ω 2 +  ≈ Ω + /2Ω = Ω+(1/2)Ω −1  cos(/2),  = 0 in the case  = 0.In general, the fractional derivative will influence the system damping more when  increases and tends to 1 and will influence the system stiffness more when  decreases and tends to 0.

The Steady-State Moments of the Response
In this section the steady-state response of system (19) will be discussed.Consider first the steady-state response for the case without random disorder as  = 0, and ( 19) becomes The steady-state solutions of ( 21) can be found by putting  =  0 ,  =  0 and Ȧ = 0, φ = 0, which leads to the following result: Squaring and adding (22) yields the frequency-response equation: There are two steady-state responses corresponding to upper and lower signs in (23), and the stability of these steady-state responses can be examined by introducing some perturbation terms as where  0 and  0 are governed by ( 22) and ( 23) and  1 and  1 are perturbation terms.Substituting ( 24) into ( 21) and neglecting the nonlinear terms, one obtains the following linearization of the modulation equations ( 21) at  0 and  0 : Substituting ( 22) into (25), one obtains The eigenvalues of the coefficient matrix of system (25) are From (20) it can be known that  > 0, and therefore the necessary and sufficient condition of the stability of the steady-state solutions  0 and  0 is that the real parts of the eigenvalues  1,2 are less than zero; that is,  ≥ 0, or Calculation shows that the steady-state response corresponding to upper sign in (23) satisfies the stability condition (28) and thus is stable, while the response corresponding to lower sign is unstable.Next, one determines the steady-state response of system (19) in the stochastic case as  > 0. Introducing another new pair of state variables Equation ( 19) can be transformed into In general, the noise () in the right hand side of ( 30) is taken as physical noise; then (30) should be transformed to Ito ones by adding Eugene and Moshe [28] correction terms to the Ito-type stochastic differential equations as follows: where () is a unit Wiener process.
An exact analytical study of system (31) seems impossible due to nonlinear nature.Thus, approximate solutions of the second-order moments of the subharmonic response are proposed.For the case  = 0, (31) can be written as Whilst the response moments of any order can be predicted easily, only the mean square amplitude [ 2 ] = [ 2 + V 2 ] will be considered here, where  denotes the mathematics expectation.For steady-state responses [] and [V], one has []/ = 0 and [V]/ = 0. Taking mathematical expectation on both sides of (32), one obtains the following equations of [] and [V]: Equations ( 33) have the following solutions: For the second-order moment of the steady-state response , one has [ 2 ]/ = 0. From (32), applying the Ito rule [29] and expectation operator yields Substituting (34) into (35), one obtains The more general case ℎ > 0 will be taken into account in the following discussion, and the approximate expression of the second order will be proposed.Substituting two terms  2 + V 2 in (31) by its steady-state moment  2 , (31) can be reduced to the following equations: where . Equations (37) can be solved in the same way as (32), and one obtains Equation ( 38) has the following solutions: It can be seen clearly that with ℎ → 0, that is,  → 0, the squared quantity ( * ) 2 approaches the exact mean square amplitude as governed by (36) and with  → 0 the squared quantity ( * ) 2 approaches the exact one as governed by (23).Similar technique has been used in [27] in solving the steady-state response.However, more efforts have been done in the derivation of the fractional derivatives.

Monte-Carlo Simulations Results
In this section, the analytical results will be shown and compared with the directly numerical results.All the directly numerical simulations using Monte-Carlo method are based on the original system dominated by (1) and can give powerful validation with analytical results.For the method of numerical simulation, readers can refer to Zhu [29], Xu et al. [30], and Shinozuka and Jan [31].In this paper, the power spectrum of () is taken as For numerical simulation it is more convenient to use the pseudorandom signal given by [29]: where   's are independent and uniformly distributed in (0, 2] and  is a larger integer number.For the simulation method of the fractional derivative, readers can refer to Kilbas et al. [14] and Samko et al. [15].For a fractional derivative of the Riemann-Liouville definition, the first-order difference is used: The time step is taken as Δ = 0.01 here.Monte-Carlo simulations are focused on the first-order subharmonics ( = 1, Ω ≈ 2), although the higher-order subharmonics (Ω ≈ 2,  = 2, 3, 4, . ..) simulations should be of the same importance.In the numerical simulation, the parameters in system (1) are chosen as  = 3.5,  = 0.95, ℎ = 0.1,  = 1, Ω = 1.The governing equation ( 1) is numerically integrated by the first-order difference method between impacts, which is valid until the first encounter with the barriers, that is, until the equality  = −ℎ is satisfied.The impact condition ẏ + = − ẏ − is then imposed, using the numerical solution ẏ − .This results in the rebound velocity ẏ + , thereby providing the initial values for the next step numerical calculation.The numerical results are shown in Figures 1-6.
The case  = 0 will be considered firstly.The variations of the steady-state response  0 with Ω 1 are shown in Figure 1 in case as  = 0.2, 0.3,  = 0.1,  = 0.6, and the theoretical results of ( 23) are also shown in Figure 1.From ( 5), one has 2 , and therefore the mean square response amplitude will be calculated as  2 0 =  2 * = ⟨( + ℎ) 2 ⟩ + ⟨( ẏ ) 2 ⟩ in numerical simulation, where angular brackets denote common time averaging for the response sample.It can be seen from Figure 1 that the response predicted by the averaging method is accurate in comparing the numerical simulations.It is clearly seen from Figure 1 that the response amplitude will decrease when the damping  increases.The response amplitude will increase when the excitation frequency Ω 1 reaches the resonant frequency 2Ω = 2.However, the peak response amplitude is not reached at Ω 1 = 2Ω = 2 but has some deviation on Ω 1 ≈ 2.04 in Figure 1.From the theoretical analysis in Section 2 one knows that the fractional derivative increases the natural frequency of the system from Ω = 1 to Ω + (1/2)Ω   cos(/2) ≈ 1.03, and then the peak response amplitude should be reached at Ω 1 = 2(Ω + (1/2)Ω   cos(/2)) ≈ 2.06, which is in good agreement with the numerical result (Ω 1 ≈ 2.04) as shown in Figure 1.The response amplitude will decrease strongly when Ω 1 departs from the resonant frequency 2Ω, and the accuracy of the analytical solution is seen to be reduced a little  in the case of large detuning comparing with the numerical solution; this may be partly due to some inaccuracy of the Krylov-Bogoliubov averaging method at large detuning.
The response time history of system (1) and the phase plot are shown in Figure 2 in the case  = 0.3,  = 0.1,  = 0.6,  = 0, Ω 1 = 2.1, where () = ẏ () denotes the velocity of the mass.Clearly, the response is a period one while the phase trajectory is a limit cycle.Now the effect of the random disorder () on the response amplitude  0 of the system is considered.The variations of the steady-state response  0 with Ω 1 are shown in Figure 3 in case as  = 0.2, 0.3,  = 0.1,  = 0.6,  = 0.25, and the theoretical results of (39) are also shown in Figure 3.
From Figure 3, high accuracy of the analytical method can also be claimed for the case with random disorder.Once again, strong reduction of the peak response amplitude due to large damping and large detuning can be seen, and the response amplitude will increase when the damping  decreases or Ω 1 is close to resonant frequency 2Ω = 2.The fractional derivative also makes the peak response amplitude reaching at Ω 1 ≈ 2.04 rather than Ω 1 = 2. Noting that the parameters of system (1) corresponding to Figure 1 ( = 0) and Figure 3 ( = 0.25) are the same except for , a rather drastic reduction of peak response amplitudes due to random disorder () in the excitation can be seen from Figure 3 in comparison with Figure 1.Such phenomena can also be illustrated in Figure 4, where the parameters of system (1) are  = 0.3,  = 0.1,  = 0.6,  = 0.25, Ω 1 = 2.1, which are the same corresponding to Figure 2 ( = 0) except for .
In comparison with Figure 2, it can be seen from Figure 4 that random disorder () will reduce the peak response amplitude and change the steady-state response of system (1) from a periodic solution to a quasiperiodic one such that it changes the phase trajectory from a limit cycle to a diffused limit cycle.Further numerical simulations show that     the width of the diffused limit cycle will be large when the intensity of the random disorder increases.
Finally the effect of the fractional derivative    on the response of the system will be discussed.The variations of the steady-state response  0 with  are shown in Figure 5 for a representative case as  = 0.3,  = 0.25, Ω 1 = 2.1,  = 0.1, 0.2, and the theoretical results given by (39) are also shown in Figure 5 for comparison.It can be seen from Figure 5 that the response amplitude will decrease when  or  increases.
However, the effect of fractional derivative    on the response of the system is more complex.The variations of the steady-state response  0 with  are shown in Figure 6 for another representative case as  = 0.3,  = 0.25, Ω 1 = 1.9,  = 0.1, 0.2, and the theoretical results given by (39) are also shown in Figure 6 for comparison.It can be seen from Figure 6 that the response amplitude will decrease when  increases, which is consistent with the conclusion of Figure 5.But the response of the system will decrease first and then increase when  increases, which is not consistent with the conclusion of Figure 5. Noting that the parameters of system (1) corresponding to Figure 5 (Ω 1 = 2.1) and Figure 6 (Ω 1 = 1.9) are the same except for Ω 1 , it can be known that the effects of  on the system response are not the same for different case.In fact from the theoretical analysis in Section 2 one knows that the fractional derivative increases both the natural frequency and damping coefficient of the system and therefore effects the response of the system in a more complicated way.High precision of the analytical expression given by (39) can also be seen in Figures 5 and  6, in comparison with the numerical results.Small system parameters are required in the theoretical derivation above; however numerical simulations show that the response values predicated by ( 23), (36), and (39) have high precision even if the system parameters are appropriately large; for example,  = 3.5 is not small in the numerical simulation here; of course the high precision of the theoretical solution will also be observed when  is small.

Conclusions
In this paper, the subharmonic response of single-degree-offreedom linear vibroimpact oscillator with fractional derivative damping and one-sided barrier under narrow-band random excitation is investigated.The main advantage of this paper is that the analytical expressions for the steady-state moments of the response amplitude are obtained, such that the effects of the fractional order derivative term, damping term, random disorder, and the coefficient of restitution and other system parameters on the system response can be discussed theoretically.The basic conclusion of the analysis may be that the fractional derivative makes both the system damping and stiffness coefficients increase, such that it changes the system parameters region at which the response amplitude reaches the maximum, thus effecting the system response in a more complicated way; the system energy loss in collision is equivalent to increasing the damping coefficient of the system; the system response amplitude will increase when the excitation frequency is close to the resonant frequency, and the system response will decay rapidly when the excitation frequency gradually deviates from the resonance frequency; the steady-state solution may change from a limit cycle to a diffused limit cycle when intensity of the random disorder increases; numerical simulations show that the proposed method is effective.