# Stochastic Responses of Lightly Nonlinear Vibroimpact System with Inelastic Impact Subjected to External Poisson White Noise Excitation

Authors: Guidong Yang; Wei Xu; Dongmei Huang; Mengli Hao
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

A procedure for analyzing stationary responses of lightly nonlinear vibroimpact system with inelastic impact subjected to external Poisson white noise excitation is proposed. First, the original vibroimpact system is transformed to a new system without velocity jump in terms of the Zhuravlev nonsmooth coordinate transformation and the Dirac delta function. Second, the averaged generalized Fokker-Planck-Kolmogorov (FPK) equation for transformed system under parametric excitation of Poisson white noise is derived by stochastic averaging method. Third, the averaged generalized FPK equation is solved by using the perturbation technique and inverse transformation of the Zhuravlev nonsmooth coordinate transformation to obtain the approximately stationary solutions for response probability density functions of original vibroimpact system. Last, analytical and numerical results for two typical lightly nonlinear vibroimpact systems are presented to assess the effectiveness of the proposed method. It is found that they are in good agreement and the proposed method is quite effective.


---

## Body

## 1. Introduction

Vibroimpact system, as a class of typical nonsmooth system [1], has attracted much attention in the past few decades due to its existence in various engineering applications [2, 3]. Certain applications or problems involved impacts include vibratory pile drivers, heat exchanger tube wear in nuclear power stations, large gaps in labyrinth seals for rotating shafts of high-pressure pumps, and a moored body with an inextensible mooring line [4–7]. There are two commonly used methods to model impacts. One is the classical or hard impact model [1], and the other is the Hertz contact model or soft impact model [8, 9]. The classical impact model consists of two parts: the ordinary motion described by a differential equation and the impact condition. In the Hertz contact model, the contact force between two elastic bodies is assumed to be proportional to the 3/2 (sometimes with different values) power of the relative displacement between them. By means of the two models, many interesting dynamical behaviors were observed such as chaos motion and bifurcations [10, 11], Hopf bifurcation [12], grazing bifurcation [13, 14], torus bifurcation [15], chattering bifurcation [16, 17], sticking motion [18, 19], and border collision bifurcation [20].As is known to all, random excitations widely exist in vibroimpact systems. Their effects usually can not be neglected and sometimes are even quite large [21]. Therefore, studies on the prediction of stochastic responses of vibroimpact systems are naturally important but difficult due to the presence of random excitations and nonsmooth factors. In recent years, some methods [22–27] have been developed to explore the response probability density function (PDF) which is an important characteristic of a stochastic system. For the classical impact model, Dimentberg et al. [28, 29] studied stochastic response of linear vibroimpact system under Gaussian white noise. Feng et al. [30, 31] explored the response PDFs of Duffing-Van der Pol vibroimpact system under independent and correlated Gaussian white noises, respectively. Zhu [32] obtained the response PDFs of a vibroimpact Duffing system with a randomly varying damping term. For the Hertz contact model, Huang et al. [33] investigated stationary responses of multi-degree-of-freedom vibroimpact system under white noise excitations. Xu et al. [34] analyzed random vibration problems of vibroimpact system with inelastic impact. Although stochastic responses of vibroimpact system have been reported, the above studies are limited to responses of vibroimpact systems subject to Gaussian white noise. In practice, many random excitations are essentially non-Gaussian, for example, those due to sea wave [35] and wind [36]. The Poisson white noise is the most used model for non-Gaussian random excitations [37–40]. However, the research on stochastic responses of vibroimpact system subject to Poisson white noise has been rarely addressed [41] and is far from enough. This paper is devoted to presenting a procedure to predict the stochastic responses of lightly nonlinear vibroimpact system subject to external Poisson white noise excitation. In the following parts, the Zhuravlev nonsmooth coordinate transformation [22], the stochastic averaging method [42], and the perturbation technique [43] are applied in succession to obtain the stationary solutions of response PDFs. Two examples are also presented to assess the effectiveness of the procedure. Comparison between results obtained from the proposed method and those from Monte-Carlo simulation shows that they are in good agreement. Therefore, the proposed procedure is quite effective.


## 2. Problem Statement

Without loss of generality, a simple mass-spring-damper model (Figure1) for single-degree-of-freedom lightly nonlinear vibroimpact system with inelastic impact subjected to external Poisson white noise excitation is considered, and the corresponding equation of motion can be expressed as(1a)x¨+ε2fx,x˙+kx=εξt,x<0,(1b)x˙+=-rx˙-,x=0,where ε>0 is a small parameter; ε2fx,x˙ denotes lightly linear or nonlinear damping force; k is the stiffness coefficient; ξt is Poisson white noise; x˙+ and x˙- denote the rebound and impact velocity, respectively; the subscripts “+” and “−” refer to values of velocity just before and after impact; r is the restitution factor whose value is among zero and one. Particularly, when r equals one, the system reduces to an elastic impact one. x=0 denotes that the rigid barrier is located at the system’s equilibrium position. Equation (1a) governs the system’s motion between impacts as x<0. Equation (1b) describes the impact law, which means that the oscillator’s velocity has an instantaneous change when it impacts with the barrier at x=0.Figure 1
Vibroimpact system with a rigid barrier on right side.The Poisson white noise can be treated as the formal derivative of the compound Poisson processCt; that is,(2)ξt=dCtdt,in which(3)Ct=∑i=1NtYiUt-ti,where Nt is the Poisson counting process with a constant impulse arrival rate λ which also means the expected number of event occurrences per second. Yi is a random variable representing the random magnitude of ith impulse arriving at time ti. U· is a unit step function. Assume that Yi are independent identically distributed random variables with zero mean which are independent of the impulse arrival time ti. Ct is a process with independent increments, and the kth correlation function of whose increment process dCt=Ct+dt-Ct is as follows [44]:(4)RkdCt1,dCt2,…,dCtk=λEYkδt2-t1⋯δtk-t1dt1dt2⋯dtk.By the relation between the correlation function and the moment function of random processes, omitting the small quantity that higher thandt, thekth moment function can be derived as(5)EdCk=λEYkdt.Given that the nonsmooth characteristic of vibroimpact system described by (1a) and (1b), the following coordinate transformation introduced by Zhuravlev [22] is used as(6)x=x1=y,x˙=x2=y˙sgn⁡y,sgn⁡y=1,y>00,y=0-1,y<0.The relation between the original variables x, x˙ and the new variables y, y˙ can be simply represented by Figures 2 and 3.Figure 2
Relation between the original variables and the new variables: (a)x and y; (b) x˙ and y˙.
(a) (b)Figure 3
Phase plane for vibroimpact system: (a) original vibroimpact system; (b) the transformed system.
(a) (b)Substituting (6) into (1a) and (1b), the transformation maps the domain x<0 of original phase plane x,x˙ onto the whole phase plane y,y˙ (see Figure 3). Since d/dysgn⁡y=0 when y≠0, this implies also that x¨t=y¨tsgn⁡y. When y=0, the impact/rebound condition for the transformed state variables is found to be y˙+=ry˙- according to Dimentberg and Menyailov [28]. Note that the identity sgn⁡y2=1, and the original equation of motion can be replaced by(7a)y¨+ε2fysgn⁡y,y˙sgn⁡ysgn⁡y+ky=εsgn⁡yξt,t≠t∗,(7b)y˙+=ry˙-,t=t∗,where t∗ represents the time instant of impact corresponding to y=0, that is, yt∗=0. Also, the velocity jump becomes proportional to 1-r instead of 1+r for original one. Using the Dirac delta function, an additional impulsive term can be introduced as follows:(8)y˙--y˙+δt-t∗=1-ry˙y˙δy.According to [28], it is reasonable to treat (7b) as an additional impulse of (7a) as t=t∗. Therefore, (7a), (7b), and (8) can be combined into one new equation as(9)y¨+ε2fysgn⁡y,y˙sgn⁡ysgn⁡y+ky-r-1y˙y˙δy=εsgn⁡yξt.The above transformed equation may permit rigorous analytical study since there contains no velocity jump and the impact condition is taken place by an additional impulsive term.


## 3. Stochastic Responses

### 3.1. Stochastic Averaging

When the values of1-r approaches to zero, the energy loss of impact is very small, and the system becomes a quasiconservative one as long as the damping coefficients and random excitations are weak. Therefore, the stochastic averaging method can be applied to study the response of the transformed system analytically. For simplicity, 1-r is assumed of order oε2. The transformed equation (9) is equivalent to a pair of first-order equations:(10)dy1=y2dt,dy2=-ε2f~y1,y2-ky1+r-1y2y2δy1+εsgn⁡y1dCt,where(11)f~y1,y2=fy1sgn⁡y1,y2sgn⁡y1sgn⁡y1.The total energy of system is(12)H=12y22+12ky12.Substituting (12) into (10) yields the following Stratonovich stochastic differential equations about displacement y1 and energy H with parametric excitation of Poisson white noise:(13)dy1=∂H∂y2dt,dH=-ε2f~y1,∂H∂y2∂H∂y2dt+r-1∂H∂y22∂H∂y2δy1dt+εsgn⁡y1∂H∂y2∘dCt.The corresponding Itô stochastic differential equations are derived using the converting rule for systems with Poisson white noise excitation proposed by Di Paola and Falsone [37]:(14)dy1=∂H∂y2dt,dH=-ε2f~y1,∂H∂y2∂H∂y2dt+r-1∂H∂y22∂H∂y2δy1dt+εsgn⁡y1∂H∂y2dCt+ε22!sgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2dCt2+ε33!sgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2dCt3+ε44!sgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2dCt4+⋯.Give that(15)sgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2=12∂∂Hsgn⁡y1∂H∂y22=1.Equation (14) can be simplified as(16)dy1=∂H∂y2dt,dH=-ε2f~y1,∂H∂y2∂H∂y2dt+r-1∂H∂y22∂H∂y2δy1dt+εsgn⁡y1∂H∂y2dCt+ε22!dCt2.The generalized FPK equation associated with (16) is(17)∂ρ∂t=-∂∂y1∂h∂y2ρ-∂∂h-ε2f~y1,∂h∂y2∂h∂y2+r-1∂h∂y22∂h∂y2δy1+εsgn⁡y1∂h∂y2EY+ε22λEY2ρ+12!∂2∂h2ε2∂h∂y22λEY2+ε3sgn⁡y13∂h∂y2λEY3+ε44λEY4ρ-13!∂3∂h3ε3sgn⁡y13∂h∂y23λEY3+32ε4∂h∂y22λEY4ρ+14!∂4∂h4ε4∂h∂y24λEY4ρ+Oε5,where ρ=ρy1,h,t.Particularly, if the probability distribution ofY is symmetric, all odd order moments are zero. Thus, we have(18)∂ρ∂t=-∂∂y1∂h∂y2ρ-∂∂h-ε2f~y1,∂h∂y2∂h∂y2+r-1∂h∂y22∂h∂y2δy1+ε22λEY2ρ+12!∂2∂h2ε2∂h∂y22λEY2+ε44λEY4ρ-13!∂3∂h332ε4∂h∂y22λEY4ρ+14!∂4∂h4ε4∂h∂y24λEY4ρ+Oε5.The conditional probability density ρy1∣h,t given by Stratonovich [45] is(19)ρy1∣h,t=k2π2h-ky12forky12<2h,0forky12>2h.By means of the formulation ρy1∣h,t=ρy1,h,t/ρh,t and integrating (18) with respect to y1, the following averaged generalized FPK equation for probability density ρh,t is obtained:(20)∂ρh,t∂t=∂∂hε2Fh+Gh-ε22λEY2ρh,t+12!∂2∂h2ε2λEY2h+ε44λEY4ρh,t-13!∂3∂h33ε42hλEY4ρh,t+14!∂4∂h43ε42h2λEY4ρh,t+Oε5,where(21)Fh=∮k2πf~y1,∂h∂y2dy1(22)Gh=∮-r-1∂h∂y2∂h∂y2δy1k2πdy1=-2kr-1πh.It is worth noting that the impulsive damping term should be averaged over a half period since there have two impacts in each period (see Figure 3).


### 3.2. Approximately Stationary Solution

The averaged generalized FPK (20) is too complicated to be solved analytically. Herein, we seek the approximately stationary solutions. Then, the left-hand side of (20) vanishes and the original equation can be rewritten as(23)dJhdh=0,where(24)Jh=ε2Fhρh+Ghρh+ε2λEY22hdρhdh+ε4λEY48hd2ρhdh2+ε4λEY416d3ρhdh3+oε5.It is easy to know from (23) that Jh is a constant. In another way, the probability current function satisfies Jhh→∞=0. Thus, we only need to solve the following equation to obtain the approximately stationary solutions:(25)ε2Fhρh+Ghρh+ε2λEY22hdρhdh+ε4λEY48hd2ρhdh2+ε4λEY416d3ρhdh3+oε5=0.The above equation can be solved by a perturbation technique [43]. The perturbation solution is assumed to be of the form(26)ρh=ρ0h+ε2ρ1h+ε4ρ2h+⋯,where ρh satisfies the boundary condition:(27)ρh≥0,h∈0,∞,ρ⟶0,h⟶∞.Substituting (26) into (25) yields the following equations:(28)Fhρ0h+Ghε2ρ0h+λEY22hdρ0hdh=0,Fhρ1h+Ghε2ρ1h+λEY22hdρ1hdh=-λEY48hd2ρ0hdh2+h22d3ρ0hdh3,….Generally speaking, the approximately stationary solution of ρh can be obtained analytically and works well with the numerical simulation results up to order of ε4. The higher terms can be neglected since their effects are relatively small.The joint PDF for the variablesy1 and y2 can be further obtained as(29)ρy1,y2=ρhρy1∣hy2=k2πρhh=1/2y22+1/2ky12.Using the inverse transformation of formula (6), the joint PDF for original displacement x1 and velocity x2 is calculated as(30)ρx1,x2=ρy1,y2∂y1,y2∂x1,x2=ρy1,y2+ρ-y1,-y2=2ρy1,y2.Thus, the marginal PDFs for displacement x1 and velocity x2 are easy to be derived, respectively:(31)ρx1=∫-∞+∞ρx1,sds,ρx2=∫-∞0ρs,x2ds.


## 3.1. Stochastic Averaging

When the values of1-r approaches to zero, the energy loss of impact is very small, and the system becomes a quasiconservative one as long as the damping coefficients and random excitations are weak. Therefore, the stochastic averaging method can be applied to study the response of the transformed system analytically. For simplicity, 1-r is assumed of order oε2. The transformed equation (9) is equivalent to a pair of first-order equations:(10)dy1=y2dt,dy2=-ε2f~y1,y2-ky1+r-1y2y2δy1+εsgn⁡y1dCt,where(11)f~y1,y2=fy1sgn⁡y1,y2sgn⁡y1sgn⁡y1.The total energy of system is(12)H=12y22+12ky12.Substituting (12) into (10) yields the following Stratonovich stochastic differential equations about displacement y1 and energy H with parametric excitation of Poisson white noise:(13)dy1=∂H∂y2dt,dH=-ε2f~y1,∂H∂y2∂H∂y2dt+r-1∂H∂y22∂H∂y2δy1dt+εsgn⁡y1∂H∂y2∘dCt.The corresponding Itô stochastic differential equations are derived using the converting rule for systems with Poisson white noise excitation proposed by Di Paola and Falsone [37]:(14)dy1=∂H∂y2dt,dH=-ε2f~y1,∂H∂y2∂H∂y2dt+r-1∂H∂y22∂H∂y2δy1dt+εsgn⁡y1∂H∂y2dCt+ε22!sgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2dCt2+ε33!sgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2dCt3+ε44!sgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2dCt4+⋯.Give that(15)sgn⁡y1∂H∂y2∂∂Hsgn⁡y1∂H∂y2=12∂∂Hsgn⁡y1∂H∂y22=1.Equation (14) can be simplified as(16)dy1=∂H∂y2dt,dH=-ε2f~y1,∂H∂y2∂H∂y2dt+r-1∂H∂y22∂H∂y2δy1dt+εsgn⁡y1∂H∂y2dCt+ε22!dCt2.The generalized FPK equation associated with (16) is(17)∂ρ∂t=-∂∂y1∂h∂y2ρ-∂∂h-ε2f~y1,∂h∂y2∂h∂y2+r-1∂h∂y22∂h∂y2δy1+εsgn⁡y1∂h∂y2EY+ε22λEY2ρ+12!∂2∂h2ε2∂h∂y22λEY2+ε3sgn⁡y13∂h∂y2λEY3+ε44λEY4ρ-13!∂3∂h3ε3sgn⁡y13∂h∂y23λEY3+32ε4∂h∂y22λEY4ρ+14!∂4∂h4ε4∂h∂y24λEY4ρ+Oε5,where ρ=ρy1,h,t.Particularly, if the probability distribution ofY is symmetric, all odd order moments are zero. Thus, we have(18)∂ρ∂t=-∂∂y1∂h∂y2ρ-∂∂h-ε2f~y1,∂h∂y2∂h∂y2+r-1∂h∂y22∂h∂y2δy1+ε22λEY2ρ+12!∂2∂h2ε2∂h∂y22λEY2+ε44λEY4ρ-13!∂3∂h332ε4∂h∂y22λEY4ρ+14!∂4∂h4ε4∂h∂y24λEY4ρ+Oε5.The conditional probability density ρy1∣h,t given by Stratonovich [45] is(19)ρy1∣h,t=k2π2h-ky12forky12<2h,0forky12>2h.By means of the formulation ρy1∣h,t=ρy1,h,t/ρh,t and integrating (18) with respect to y1, the following averaged generalized FPK equation for probability density ρh,t is obtained:(20)∂ρh,t∂t=∂∂hε2Fh+Gh-ε22λEY2ρh,t+12!∂2∂h2ε2λEY2h+ε44λEY4ρh,t-13!∂3∂h33ε42hλEY4ρh,t+14!∂4∂h43ε42h2λEY4ρh,t+Oε5,where(21)Fh=∮k2πf~y1,∂h∂y2dy1(22)Gh=∮-r-1∂h∂y2∂h∂y2δy1k2πdy1=-2kr-1πh.It is worth noting that the impulsive damping term should be averaged over a half period since there have two impacts in each period (see Figure 3).


## 3.2. Approximately Stationary Solution

The averaged generalized FPK (20) is too complicated to be solved analytically. Herein, we seek the approximately stationary solutions. Then, the left-hand side of (20) vanishes and the original equation can be rewritten as(23)dJhdh=0,where(24)Jh=ε2Fhρh+Ghρh+ε2λEY22hdρhdh+ε4λEY48hd2ρhdh2+ε4λEY416d3ρhdh3+oε5.It is easy to know from (23) that Jh is a constant. In another way, the probability current function satisfies Jhh→∞=0. Thus, we only need to solve the following equation to obtain the approximately stationary solutions:(25)ε2Fhρh+Ghρh+ε2λEY22hdρhdh+ε4λEY48hd2ρhdh2+ε4λEY416d3ρhdh3+oε5=0.The above equation can be solved by a perturbation technique [43]. The perturbation solution is assumed to be of the form(26)ρh=ρ0h+ε2ρ1h+ε4ρ2h+⋯,where ρh satisfies the boundary condition:(27)ρh≥0,h∈0,∞,ρ⟶0,h⟶∞.Substituting (26) into (25) yields the following equations:(28)Fhρ0h+Ghε2ρ0h+λEY22hdρ0hdh=0,Fhρ1h+Ghε2ρ1h+λEY22hdρ1hdh=-λEY48hd2ρ0hdh2+h22d3ρ0hdh3,….Generally speaking, the approximately stationary solution of ρh can be obtained analytically and works well with the numerical simulation results up to order of ε4. The higher terms can be neglected since their effects are relatively small.The joint PDF for the variablesy1 and y2 can be further obtained as(29)ρy1,y2=ρhρy1∣hy2=k2πρhh=1/2y22+1/2ky12.Using the inverse transformation of formula (6), the joint PDF for original displacement x1 and velocity x2 is calculated as(30)ρx1,x2=ρy1,y2∂y1,y2∂x1,x2=ρy1,y2+ρ-y1,-y2=2ρy1,y2.Thus, the marginal PDFs for displacement x1 and velocity x2 are easy to be derived, respectively:(31)ρx1=∫-∞+∞ρx1,sds,ρx2=∫-∞0ρs,x2ds.


## 4. Illustrative Examples

In this section, two examples are shown to illustrate the proposed procedure. In order to compare the analytical results and numerical results, the impulses magnitudesY of Poisson white noise are assumed to be Gaussian distributed with zero mean. Thus, an interesting result can be used to simplify the practical computation.(32)EY4=3EY22.


### 4.1. Rayleigh Oscillator

In this example, a typical Rayleigh vibroimpact system under external Poisson white noise excitation is considered:(33a)x¨+ε2x˙2-1x˙+kx=εξt,x<0,(33b)x˙+=-rx˙-,x=0.In this case, fx,x˙=x˙2-1x˙. The following transformed equation without velocity jump is obtained by means of the Zhuravlev transformation and Dirac delta function.(34)y¨+ε2y˙2-1y˙+ky-r-1y˙y˙δy=εsgn⁡yξt.The corresponding Itô stochastic differential equation is derived as(35)dy1=∂H∂y2dt,dH=-ε2∂H∂y22-1∂H∂y22dt+r-1∂H∂y22∂H∂y2δy1dt+εsgn⁡y1∂H∂y2adCt+ε22dCt2.According to (20) and (21), Fh is calculated as follows in the averaged generalized FPK equation for probability density ρh,t:(36)Fh=∮-k2πβ∂h∂y22-1∂h∂y2dy1=32h2-h.Using the perturbation technique, the approximately stationary solution for energy PDF is obtained:(37)ρh=C~e-N1+ε2Dh+Eh2+Fh3+Gh4+Mh5,where(38)A=λEY2,B=λEY4,C=-1-2r-1kπε2,D=3B4A2-BC2A3,E=BC32A4-21BC8A3,F=3BC22A4-15B8A3,G=27BC16A4,M=27B40A4,N=2A3h24+Ch.Furthermore, the joint PDF for original displacement x1 and velocity x2 can be calculated by substituting (37) and (38) into (29) and (30) as(39)ρx1,x2=kπρhh=1/2y22+1/2ky12.To verify the effectiveness of the proposed procedure, analytical results and Monte-Carlo simulation results are shown in Figure4 for PDFs of system energy with k=1, r=0.99, λ=1, and EY2=1 and two different values of small parameter ε (ε=0.1 and 0.2). The solid line is the analytical results obtained by the proposed method while • represents the Monte-Carlo simulation results based on the original system (33a) and (33b). It is obvious that the two results agree well with each other. In order to illustrate the effects of average arrival rate of external Poisson white noise on the accuracy of the approximately stationary solution, Figure 5 shows PDFs of system energy analytically and numerically with different values of λ (λ=1 and 2) when ε is fixed as 0.1. It is seen that the proposed method works well.Figure 4
Stationary PDFs of total energy obtained numerically and analytically for different values of small parameterε.Figure 5
Stationary PDFs of total energy obtained numerically and analytically for different values of average arrival rateλ.To investigate the effects of restitution coefficient on the stationary response of vibroimpact system, numerical and analytical results are shown in Figures6–8 with the parameters k=1, ε=0.1, λ=1, and EY2=1 and different values of r. It is seen from Figure 6 that the stationary PDFs of total energy are shown for r = 0.97, 0.98, 0.99, and 0.995. The PDFs of total energy present lower peaks with the increase of restitution coefficient. It is found that the values of peak are at the system’s static equilibrium position when r = 0.97, 0.98 while they do not stay there when r = 0.99, 0.995. This may mean that the steady state has a qualitative change with the increase of restitution coefficient. To show it more clearly, Figures 7 and 8 provide analytical and numerical results of joint PDFs of the displacement x1 and velocity x2 for r=0.98 and r=0.995, respectively. It is seen that the analytical results agree well with the numerical results. In Figure 7, the joint PDFs of the displacement x1 and velocity x2 have a singular peak which appears at the static equilibrium position. When r increases to r=0.995, the joint PDFs present the shape of crater in Figure 8. This illuminates the occurrence of P-bifurcation according to the concept of stochastic bifurcation. In other words, the restitution coefficient r can induce the occurrence of P-bifurcation.Figure 6
Stationary PDFs of total energy obtained numerically and analytically for different values of restitution coefficientr.Figure 7
Joint PDFs of the displacementx1 and velocity x2 for r=0.98: (a) Monte-Carlo results; (b) analytical results.
(a) (b)Figure 8
Joint PDFs of the displacementx1 and velocity x2 for r=0.995: (a) Monte-Carlo results; (b) analytical results.
(a) (b)


### 4.2. Van der Pol Oscillator

As an example, a typical Van der Pol vibroimpact system under external Poisson white noise excitation is considered:(40a)x¨-ε2β1-x2x˙+kx=εξt,x<0,(40b)x˙+=-rx˙-,x=0.Based on (9), the transformed system without velocity is obtained as(41)y¨-ε2β1-y2y˙+ky-r-1y˙y˙δy=εsgn⁡yξt.The Itô stochastic differential equation about displacement y1 and energy H is derived as(42)dy1=∂H∂y2dt,dH=ε2β1-y12∂H∂y22dt+r-1∂H∂y22∂H∂y2δy1dt+εsgn⁡y1∂H∂y2dCt+ε22dCt2.In the averaged generalized FPK (20), Fh here is calculated as(43)Fh=∮-k2πβ1-y12∂h∂y2dy1=β2kh2-βh.Thus, ρh is obtained analytically by solving the corresponding averaged generalized FPK equation as(44)ρh=C~e-N1+ε2Dh+Eh2+Fh3+Gh4+Mh5,where(45)A=λEY2,B=λEY4,C=-β-2r-1kπε2,D=Bβ4A2k-BC2A3,E=BC32A4-7BCβ8A3k,F=BC2β2A4k-5Bβ224A3k2,G=3BCβ216A4k2,M=Bβ340A4k3,N=2Aβh24k+Ch.Furthermore, the PDFs ρx1,x2, ρx1, and ρx2 can be calculated following (29)–(31), respectively.To demonstrate the effectiveness of the proposed method, numerical and analytical results of response PDFs are provided in Figures9–13. The effects of damping coefficient on the accuracy of the energy PDFs are shown in Figure 9 with ε=0.1, r=0.99, λ=1, and EY2=1 and two different values of β (β=1.0 and 2.0). Stationary PDFs of total energy obtained numerically and analytically for different values of restitution coefficient r are presented in Figure 10. It is clear from Figures 9 and 10 that the analytical results give powerful validation with Monte-Carlo simulation results. The joint PDFs of the displacement x1 and velocity x2 for different impulse intensities are shown numerically and analytically in Figures 11 and 12. The corresponding marginal PDFs ρx1 and ρx2 are provided in Figures 13(a) and 13(b), respectively. It is seen that the results obtained by using the proposed method and those from Monte-Carlo simulation are in good agreement. Additionally, the increase of impulse intensity only effects the values of peaks of stationary response PDFs, and there is no occurrence of P-bifurcation.Figure 9
Stationary PDFs of total energy obtained numerically and analytically for different values ofβ.Figure 10
Stationary PDFs of total energy obtained numerically and analytically for different values of restitution coefficientr.Figure 11
Joint PDFs of the displacementx1 and velocity x2 for EY2=1.0: (a) Monte-Carlo results; (b) analytical results.
(a) (b)Figure 12
Joint PDFs of the displacementx1 and velocity x2 for EY2=2.0: (a) Monte-Carlo results; (b) analytical results.
(a) (b)Figure 13
Stationary PDFs obtained numerically and analytically for different values of random impulse intensityEY2. -, analytical results; ·, Monte-Carlo results: (a) PDFs of displacement x1; (b) PDFs of velocity x2.
(a) (b)


## 4.1. Rayleigh Oscillator

In this example, a typical Rayleigh vibroimpact system under external Poisson white noise excitation is considered:(33a)x¨+ε2x˙2-1x˙+kx=εξt,x<0,(33b)x˙+=-rx˙-,x=0.In this case, fx,x˙=x˙2-1x˙. The following transformed equation without velocity jump is obtained by means of the Zhuravlev transformation and Dirac delta function.(34)y¨+ε2y˙2-1y˙+ky-r-1y˙y˙δy=εsgn⁡yξt.The corresponding Itô stochastic differential equation is derived as(35)dy1=∂H∂y2dt,dH=-ε2∂H∂y22-1∂H∂y22dt+r-1∂H∂y22∂H∂y2δy1dt+εsgn⁡y1∂H∂y2adCt+ε22dCt2.According to (20) and (21), Fh is calculated as follows in the averaged generalized FPK equation for probability density ρh,t:(36)Fh=∮-k2πβ∂h∂y22-1∂h∂y2dy1=32h2-h.Using the perturbation technique, the approximately stationary solution for energy PDF is obtained:(37)ρh=C~e-N1+ε2Dh+Eh2+Fh3+Gh4+Mh5,where(38)A=λEY2,B=λEY4,C=-1-2r-1kπε2,D=3B4A2-BC2A3,E=BC32A4-21BC8A3,F=3BC22A4-15B8A3,G=27BC16A4,M=27B40A4,N=2A3h24+Ch.Furthermore, the joint PDF for original displacement x1 and velocity x2 can be calculated by substituting (37) and (38) into (29) and (30) as(39)ρx1,x2=kπρhh=1/2y22+1/2ky12.To verify the effectiveness of the proposed procedure, analytical results and Monte-Carlo simulation results are shown in Figure4 for PDFs of system energy with k=1, r=0.99, λ=1, and EY2=1 and two different values of small parameter ε (ε=0.1 and 0.2). The solid line is the analytical results obtained by the proposed method while • represents the Monte-Carlo simulation results based on the original system (33a) and (33b). It is obvious that the two results agree well with each other. In order to illustrate the effects of average arrival rate of external Poisson white noise on the accuracy of the approximately stationary solution, Figure 5 shows PDFs of system energy analytically and numerically with different values of λ (λ=1 and 2) when ε is fixed as 0.1. It is seen that the proposed method works well.Figure 4
Stationary PDFs of total energy obtained numerically and analytically for different values of small parameterε.Figure 5
Stationary PDFs of total energy obtained numerically and analytically for different values of average arrival rateλ.To investigate the effects of restitution coefficient on the stationary response of vibroimpact system, numerical and analytical results are shown in Figures6–8 with the parameters k=1, ε=0.1, λ=1, and EY2=1 and different values of r. It is seen from Figure 6 that the stationary PDFs of total energy are shown for r = 0.97, 0.98, 0.99, and 0.995. The PDFs of total energy present lower peaks with the increase of restitution coefficient. It is found that the values of peak are at the system’s static equilibrium position when r = 0.97, 0.98 while they do not stay there when r = 0.99, 0.995. This may mean that the steady state has a qualitative change with the increase of restitution coefficient. To show it more clearly, Figures 7 and 8 provide analytical and numerical results of joint PDFs of the displacement x1 and velocity x2 for r=0.98 and r=0.995, respectively. It is seen that the analytical results agree well with the numerical results. In Figure 7, the joint PDFs of the displacement x1 and velocity x2 have a singular peak which appears at the static equilibrium position. When r increases to r=0.995, the joint PDFs present the shape of crater in Figure 8. This illuminates the occurrence of P-bifurcation according to the concept of stochastic bifurcation. In other words, the restitution coefficient r can induce the occurrence of P-bifurcation.Figure 6
Stationary PDFs of total energy obtained numerically and analytically for different values of restitution coefficientr.Figure 7
Joint PDFs of the displacementx1 and velocity x2 for r=0.98: (a) Monte-Carlo results; (b) analytical results.
(a) (b)Figure 8
Joint PDFs of the displacementx1 and velocity x2 for r=0.995: (a) Monte-Carlo results; (b) analytical results.
(a) (b)


## 4.2. Van der Pol Oscillator

As an example, a typical Van der Pol vibroimpact system under external Poisson white noise excitation is considered:(40a)x¨-ε2β1-x2x˙+kx=εξt,x<0,(40b)x˙+=-rx˙-,x=0.Based on (9), the transformed system without velocity is obtained as(41)y¨-ε2β1-y2y˙+ky-r-1y˙y˙δy=εsgn⁡yξt.The Itô stochastic differential equation about displacement y1 and energy H is derived as(42)dy1=∂H∂y2dt,dH=ε2β1-y12∂H∂y22dt+r-1∂H∂y22∂H∂y2δy1dt+εsgn⁡y1∂H∂y2dCt+ε22dCt2.In the averaged generalized FPK (20), Fh here is calculated as(43)Fh=∮-k2πβ1-y12∂h∂y2dy1=β2kh2-βh.Thus, ρh is obtained analytically by solving the corresponding averaged generalized FPK equation as(44)ρh=C~e-N1+ε2Dh+Eh2+Fh3+Gh4+Mh5,where(45)A=λEY2,B=λEY4,C=-β-2r-1kπε2,D=Bβ4A2k-BC2A3,E=BC32A4-7BCβ8A3k,F=BC2β2A4k-5Bβ224A3k2,G=3BCβ216A4k2,M=Bβ340A4k3,N=2Aβh24k+Ch.Furthermore, the PDFs ρx1,x2, ρx1, and ρx2 can be calculated following (29)–(31), respectively.To demonstrate the effectiveness of the proposed method, numerical and analytical results of response PDFs are provided in Figures9–13. The effects of damping coefficient on the accuracy of the energy PDFs are shown in Figure 9 with ε=0.1, r=0.99, λ=1, and EY2=1 and two different values of β (β=1.0 and 2.0). Stationary PDFs of total energy obtained numerically and analytically for different values of restitution coefficient r are presented in Figure 10. It is clear from Figures 9 and 10 that the analytical results give powerful validation with Monte-Carlo simulation results. The joint PDFs of the displacement x1 and velocity x2 for different impulse intensities are shown numerically and analytically in Figures 11 and 12. The corresponding marginal PDFs ρx1 and ρx2 are provided in Figures 13(a) and 13(b), respectively. It is seen that the results obtained by using the proposed method and those from Monte-Carlo simulation are in good agreement. Additionally, the increase of impulse intensity only effects the values of peaks of stationary response PDFs, and there is no occurrence of P-bifurcation.Figure 9
Stationary PDFs of total energy obtained numerically and analytically for different values ofβ.Figure 10
Stationary PDFs of total energy obtained numerically and analytically for different values of restitution coefficientr.Figure 11
Joint PDFs of the displacementx1 and velocity x2 for EY2=1.0: (a) Monte-Carlo results; (b) analytical results.
(a) (b)Figure 12
Joint PDFs of the displacementx1 and velocity x2 for EY2=2.0: (a) Monte-Carlo results; (b) analytical results.
(a) (b)Figure 13
Stationary PDFs obtained numerically and analytically for different values of random impulse intensityEY2. -, analytical results; ·, Monte-Carlo results: (a) PDFs of displacement x1; (b) PDFs of velocity x2.
(a) (b)


## 5. Conclusions

The paper is devoted to presenting a solution procedure for predicting stationary responses of single-degree-of-freedom lightly nonlinear vibroimpact systems with inelastic impact subjected to external Poisson white noise excitations. The classical model with instantaneous impacts is used and the rigid barrier is located at the static equilibrium of the system. The solution procedure consists of three steps. First, the Zhuravlev nonsmooth coordinate transformation and the Dirac delta function are adopted to convert the original vibroimpact system to a new system without velocity jump. Second, the Itô stochastic differential equations with correction terms are obtained and the averaged generalized Fokker-Planck-Kolmogorov (FPK) equation is derived by stochastic averaging method. Third, the perturbation method is utilized to solve the averaged generalized FPK equation and the inverse transformation of the Zhuravlev nonsmooth coordinate transformation is used to establish the PDFs for original vibroimpact system. Furthermore, two examples are presented to show the effectiveness of the proposed method. It is found that the analytical results obtained by the proposed method and the Monte-Carlo simulation results agree well. Therefore, the proposed procedure is quite promising and can be extended to more vibroimpact systems subject to Poisson white noise excitation.



---

*Source: 3627195-2018-04-24.xml*