# Stochastic P-Bifurcation of a Bistable Viscoelastic Beam with Fractional Constitutive Relation under Gaussian White Noise

Authors: Yajie Li; Zhiqiang Wu; Guoqi Zhang; Feng Wang
Journal: Discrete Dynamics in Nature and Society
 (2018)

---

## Abstract

In this paper, we study the stochastic P-bifurcation problem for axially moving of a bistable viscoelastic beam with fractional derivatives of high order nonlinear terms under Gaussian white noise excitation. First, using the principle for minimum mean square error, we show that the fractional derivative term is equivalent to a linear combination of the damping force and restoring force, so that the original system can be simplified to an equivalent system. Second, we obtain the stationary Probability Density Function (PDF) of the system’s amplitude by stochastic averaging, and using singularity theory, we find the critical parametric condition for stochastic P-bifurcation of amplitude of the system. Finally, we analyze the types of the stationary PDF curves of the system qualitatively by choosing parameters corresponding to each region within the transition set curve. We verify the theoretical analysis and calculation of the transition set by showing the consistency of the numerical results obtained by Monte Carlo simulation with the analytical results. The method used in this paper directly guides the design of the fractional order viscoelastic material model to adjust the response of the system.


---

## Body

## 1. Introduction

Fractional calculus is a generalization of integer-order calculus, it extends the order of calculus operation from the traditional integer order to the case of noninteger order, and it has a history of more than 300 years as so far. Due to the limitation of the definition of integer-order derivative, it cannot express the memory property of viscoelastic substances. The definition of fractional derivative contains convolution, which can well express the memory effect and show the cumulative effect over time. Compared with the traditional integer-order calculus, fractional calculus has more advantages and is a suitable mathematical tool for describing the memory characteristics [1–11] and in recent years, it has become the powerful mathematical tool in many disciplines, especially in the study of viscoelastic materials.The fractional derivative can accurately describe the constitutive relation of viscoelastic materials with fewer parameters, so the studies of fractional differential equations on the typical mechanical properties and the influences of fractional order parameters on the system are very necessary and have important significance. In recent years, many scholars have done a lot of work and achieved fruitful results in this field: Li and Tang studied the nonlinear parametric vibration of an axially moving string made by rubber-like materials, a new nonlinear fractional mathematical model governing transverse motion of the string is derived based on Newton’s second law, the Euler beam theory, and the Lagrangian strain, and the principal parametric resonance is analytically investigated via applying the direct multiscale method [12]. Liu et al. introduced a transfer entropy and surrogate data algorithm to identify the nonlinearity level of the system by using a numerical solution of nonlinear response of beams, the Galerkin method was applied to discretize the dimensionless differential governing equation of the forced vibration, and then the fourth-order Runge-Kutta method was used to obtain the time history response of the lateral displacement [13]. Liu et al. investigated the stochastic stability of coupled viscoelastic system with nonviscously damping driven by white noise through moment Lyapunov exponents and Lyapunov exponents, obtained the coupled Itô stochastic differential equations of the norm of the response and angles process by using the coordinate transformation, and discussed the effects of various physical quantities of stochastic coupled system on the stochastic stability [14]. Nutting, Gemant and Scott-Blair et al. [15–17] first proposed the fractional derivative models to study the constitutive relation of viscoelastic materials and the research on the viscoelastic materials with fractional derivative is also increasing, and so far, it is still a research hotspot [18–25]. Rodr Guez et al. calculated the correlation function of transverse wave in linear and homogeneous viscoelastic liquid by the Generalized Langevin Equation (GLE) method and the influence of fractional correlation function on the dynamic behavior of the system is analyzed [26]. Bagley and Torvik used fractional calculus to study the dynamic behavior of viscoelastic damping structure and the responses of the system under general load as well as step load are analyzed respectively [27, 28]. Pakdemirli and Ulsoy studied the primary parametric resonance and combined resonance of the axial acceleration rope based on the discrete perturbation method and the multiscale method [29]. Zhang and Zhu analyzed the stability and dynamic response of viscoelastic belt under parametric excitation by the multiscale method [30, 31]. Chen et al. studied the dynamic behavior and steady-state response of axially accelerating viscoelastic beam by the Galerkin method [32–35], derived the differential equation of nonlinear vibration for axially moving viscoelastic rope, and then pointed out that the damping of viscoelastic rope only exists in the nonlinear term [36, 37]. Leung et al. studied the steady-state response of a simply supported viscoelastic column under the axial harmonic excitation based on the fractional derivative constitutive model of cubic nonlinear and derived the generalized Mathieu-Duffing equation with time delay by the Galerkin discrete method, then the bifurcation behavior of the system caused by the order of the fractional derivative is analyzed [38]. Ghayesh and Moradian developed the Kelvin-Voigt viscoelastic model of the axially moving and the tensile belt, and then found the existence of nontrivial limit cycle in this system [39]. Liu et al. studied the dynamic response of an axially moving viscoelastic beam under random disorder periodic excitation, the first order expression of the solution is obtained by the multiscale method, and the stochastic jump phenomenon between the steady-state solutions is carried out [40]. Yang and Fang derived the system equation based on Newton’s second law and the fractional Kelvin constitutive relation and then studied the stability of the axially moving beam under the parametric resonance condition [41]. Leung et al. studied the single mode dynamic characteristics of the nonlinear arch with the fractional derivative, the steady-state solution of the system is obtained based on the residual harmonic homotopy method, and the influence of the parametric variation on the dynamic behaviors of the viscoelastic damping material is analyzed [42]. Galucio et al. obtained the fractional derivative model to describe the viscoelasticity of the system based on the Timoshenko theory and Euler-Bernoulli hypothesis and proposed a finite element formula for analyzing the sandwich beam of viscoelastic material with fractional derivative and the results were verified numerically [43].Due to the complexity of fractional derivative, the analysis method of it becomes more difficult, the study on the vibration characteristics of the parameters can only be qualitatively analyzed, and the critical conditions of the parametric influences cannot be found, which affect the analysis and design of such systems, as well as the stochastic P-bifurcation of bistability for the viscoelastic beam with fractional derivatives of high order nonlinear terms under random noise excitation has not been reported. In view of the above situation, the nonlinear vibration of viscoelastic beam with fractional constitutive relation under Gaussian white noise excitation is taken as an example, the transition set curve of the fractional order system as well as the critical parametric condition for stochastic P-bifurcation of the system is obtained by the singularity theory, and then the types of stationary PDF curves of the system in each region in the parametric plane divided by the transition set are analyzed. By the method of Monte Carlo simulation, the numerical results are compared with the analytical results obtained in this paper, it can be seen that the numerical solutions are in good agreements with the analytical solutions, and thus the correctness of the theoretical analysis in this paper is verified.


## 2. Equation of Axially Moving Viscoelastic Beam

There are many definitions of fractional derivatives, and the Riemann-Liouville derivative and Caputo derivative are commonly used. The initial conditions corresponding to the Riemann-Liouville derivative have no physical meanings, however, the initial conditions of the systems described by the Caputo derivative have clear physical meanings and their forms are the same as the initial conditions for the differential equations of integer order. So in this paper, the Caputo-type fractional derivative is adopted as follows:(1)aCDpxt=1Γm-p∫atxmut-u1+p-mdu,where m-1<p≤m,m∈N, t∈[a,b], Γm is the Euler Gamma function, and x(m)(t) is the m order derivative of x(t).For a given physical system, due to the fact that the initial moment of the oscillator ist=0, the following form of the Caputo derivative is often used:(2)0CDpxt=1Γm-p∫0txmut-u1+p-mdu,In this paper, the transverse vibrationy(x,t) of a viscoelastic simply supported beam under lateral excitation F(x,t) as shown in Figure 1 is considered; applying the d’Alambert principle, the governing equation can be written as [42](3)∂Q∂x=Fx,t-ρA∂2y∂t2∂M∂x=Q-N∂y∂x,where ρ is the mass density of the beam, A is the area of cross-section, M is the bending moment, Q is the lateral force, and N is the horizontal force. From (3), we have(4)∂2M∂x2+N∂2y∂x2+ρA∂y2∂t2-Fx,t=0.Assuming the material of the beam obeys a fractional derivative viscoelastic constitutive relation:(5)σx,z,t=E0εx,z,t+E10CDpεx,z,t=E01+η0CDpεx,z,t≜Sεx,z,t,where p is the order of fractional derivative 0CDp[ε(x,z,t)] as is defined in (2), η=E1/E0 is the material modulus ratio, and ε(x,z,t) is axial strain component.Figure 1
Schematic stress diagram of viscoelastic beam.When the deformation of the beam is small, the axial strainε and lateral displacement y(x,t) satisfy the relationship as follows:(6)εx,z,t=-z∂2yx,t∂x2.Substituting (6) into (5) yields(7)σx,z,t=S-z∂2yx,t∂x2.The relationship between bending moment M(x,t) and axial stress σ(x,z,t) of the beam can be expressed as follows:(8)Mx,t=∫-h/2h/2zσx,z,tdz,where h is the thickness of the beam.From (7) and (8), the expression of the bending moment M can be obtained as follows:(9)Mx,t=IE0∂2y∂x2+E1·0CDp∂2y∂x2=IS∂2y∂x2,where I=h3/12.The expression of the horizontal tension is(10)N=E0A2L∫0L∂y∂x2-∂y∂x4dx+E1A2L∫0L0CDp∂y∂x2-∂y∂x4dx=A2LS∫0L∂y∂x2-∂y∂x4dx.Substituting (9) and (10) into system (4), system (4) can be rewritten as(11)IS∂4y∂x4+SA2L∂2y∂x2∫0L0CDp∂y∂x2-∂y∂x4dx+ρA∂2y∂t2-Fx,t=0.The boundary conditions are(12)∂2y∂x2=0,whenx=0,x=L.According to the boundary conditions (12), the solution of system (11) can be expressed as the Fourier series:(13)yx,t=∑n=1∞untsin⁡nπLx.Assume that the initial transverse vibration of the system is y0(x)=0 and the transverse load of the system satisfies the following form:(14)Fx,t=ρA·sin⁡πLxξt,where ξ(t) is the Gaussian white noise, which satisfies E[ξ(t)]=0, E[ξ(t)ξ(t-τ)]=2Dδ(τ), D represents the noise intensity, and δ(τ) is the Dirac function.By the discrete format based on Galerkin method, system (11) can be reduced to the fractional differential equation as follows:(15)u¨+k1u+ηk10CDpu-k3u3-ηk3u0CDpu2+k5u5+ηk5u0CDpu4=ξt,where 0CDtpu is the p  (0<p<1) order Caputo derivative of u(t) as is defined in (2), and(16)k1=πL4IE0ρA,k3=πL4E04ρ,k5=πL63E016ρ.For convenience, system (15) can be represented as follows:(17)u¨+w2u+ηw20CDpu-k3u3-ηk3u0CDpu2+k5u5+ηk5u0CDpu4=ξt,where w=k1.The fractional derivative term has contributions to both damping and restoring forces [44–47], hence, introducing the following equivalent system:(18)u¨+w2u+ηw2cau˙+w2au-k3u3-ηk3u0CDpu2+k5u5+ηk5u0CDpu4=ξt,where c(a), w2(a) are the coefficients of equivalent damping and restoring forces of fractional derivative 0CDpu, respectively.The error between system (17) and (18) is(19)e=ηw2cau˙+w2au-0CDpu.The necessary conditions for minimum mean square error are [48](20)∂Ee2∂ca=0,∂Ee2∂wa=0.Substituting (19) into (20) yields(21)Ecau˙2+w2auu˙-u˙0CDpu=limT→∞⁡1T∫0Tcau˙2+w2auu˙-u˙0CDpudt=0Ecau˙u+w2au2-u0CDpu=limT→∞⁡1T∫0Tcau˙u+w2au2-u0CDpudt=0.Assume that the solution of system (18) has the following form:(22)ut=atcos⁡φtφt=wt+θ,where w is the natural frequency of system (17).Based on (22), we can obtain(23)u˙t=-atwsin⁡φt.Substituting (22) and (23) into (21) yields(24)limT→∞⁡1T∫0Tcau˙2+w2auu˙-u˙0CDpudt=limT→∞⁡1T∫0Tcaa2w2sin2⁡φ-a2ww2asin⁡φcos⁡φ-awsin⁡φ0CDpudt≈caa2w22+1Γ1-plimT→∞⁡1T∫oTawsin⁡φ∫0tx˙t-τταdτdt=caa2w22-1Γ1-plimT→∞⁡1T∫oTa2w2sin⁡φ∫0tsin⁡φcos⁡wτ-cos⁡φsin⁡wτταdt=0(25)limT→∞⁡1T∫0Tcau˙u+w2auu-u0CDpudt=limT→∞⁡1T∫0T-caa2wsin⁡φcos⁡φ+a2w2acos2⁡φ-acos⁡φ0CDpudt≈a2w2a2-1Γ1-plimT→∞⁡1T∫oTacos⁡φ∫0tx˙t-τταdτdt=a2w2a2+1Γ1-plimT→∞⁡1T∫oTa2wcos⁡φ∫0tsin⁡φcos⁡wτ-cos⁡φsin⁡wτταdt=0.To simplify (24) and (25) further, asymptotic integrals are introduced as follows:(26)∫0tcos⁡wττpdτ=wp-1Γ1-psin⁡pπ2+sin⁡wtwtp+owt-p-1∫0tsin⁡wττpdτ=wp-1Γ1-pcos⁡pπ2-cos⁡wtwtp+owt-p-1.Substituting (26) into (24) and (25) and averaging them across φ produce the ultimate forms of c(a) and w2(a) as follows:(27)ca=wp-1sin⁡pπ2,w2a=wpcos⁡pπ2.Therefore, the equivalent system associated with system (18) can be expressed as follows:(28)u¨+w02u+γu˙-k3u3-ηk3u0CDpu2+k5u5+ηk5u0CDpu4=ξt,where w02=w21+ηwpcos⁡pπ/2 and γ=ηwp+1sin⁡(pπ/2).Next, we consider the stochastic P-bifurcation of system (28) which comprises the fractional derivatives of high order nonlinear terms and analyze the influence of parametric variation on the system response.


## 3. The Stationary PDF of Amplitude

For the system (28), the material modulus ratio is given as η=0.5, coefficients of nonlinear terms are given as k3=7.8, k5=5.9, respectively, and nature frequency is given as w=1. For the convenience to discuss the parametric influence, the bifurcation diagram of amplitude of the limit cycle along with variation of the fractional order p is shown in Figure 2 when D=0.Figure 2
Bifurcation diagram of amplitude of system (28) at D=0.As can be seen from Figure2, the solution corresponding to the solid line is almost completely coincided with the numerical solution, which proves the correctness and accuracy of the approximate analytical result of the deterministic system; at the same time, it shows that the solution corresponding to the solid line is stable and the solution corresponding to the dotted line is unstable. And it also can be seen that there is 1 attractor in the system when p changes in the interval [0,0.6817]: equilibrium, as shown in Figure 3(a); there are 2 attractors when p changes in the interval [0.6818, 1]: equilibrium and limit cycle, as shown in Figure 3(b).Figure 3
Phase diagrams of system (28) at D=0.
(a)
Equilibrium (b)
Equilibrium and limit cycleIn order to obtain the stationary PDF of the amplitude of system (28), the following transformation is introduced:(29)ut=atcos⁡Φtu˙t=-atw0sin⁡ΦtΦt=w0t+θt,where w0 is the natural frequency of the equivalent system (28), a(t) and θ(t) represent the amplitude process and the phase process of the system response, respectively, and they are all random processes.Substituting (29) into (28), we can obtain(30)dadt=F11a,θ+G11a,θξtdθdt=F21a,θ+G21a,θξt,in which(31)F11=-Δsin⁡Φw0F21=-Δcos⁡Φaw0G11=-sin⁡Φw0G21=-cos⁡Φaw0,and(32)Δ=ηwp+1aw0sin⁡pπ2sin⁡Φ-k3a3cos3⁡Φ-ηk3a322w0pcos⁡Φcos⁡2Φ+pπ2-k5a5cos5⁡Φ-ηk5a5cos⁡Φ184w0pcos⁡4Φ+pπ2+122w0pcos⁡2Φ+pπ2.Equation (30) can be regarded as a Stratonovich stochastic differential equation, by adding the corresponding Wong-Zakai correction term: it can be transformed into the following Itô stochastic differential equation:(33)da=F11a,θ+F12a,θdt+2DG11a,θdBtdθ=F21a,θ+F22a,θdt+2DG21a,θdBt,where B(t) is a standard Wiener process and(34)F12a,θ=D∂G11∂aG11+D∂G11∂θG21F22a,θ=D∂G21∂aG11+D∂G21∂θG21.By the stochastic averaging method, averaging (33) regarding Φ, the following averaged Itô equation can be obtained as follows:(35)da=m1adt+σ1adBtdθ=m2adt+σ2adBt,where(36)m1a=-ηk52pw0p-1a5sin⁡pπ/28-ηk32pw0p-1a3sin⁡pπ/28-ηwp+1asin⁡pπ/22+D2aw02σ12a=Dw02m2a=k52η·2pw0pcos⁡pπ/2+5a416w0+k3η·2pw0pcos⁡pπ/2+3a28w0σ22a=Da2w02.Equations (35) and (36) show that the averaged Itô equation for at is independent of θ(t) and the random process a(t) is a one-dimensional diffusion process. The corresponding FPK equation associated with a(t) can be written as(37)∂p∂t=-∂∂am1ap+12∂2∂a2σ12ap.The boundary conditions are(38)p=c,c∈-∞,+∞whena=0p→0,∂p∂a→0whena→∞.Thus, based on these boundary conditions (38), we can obtain the stationary PDF of amplitude a as follows:(39)pa=Cσ12aexp⁡∫0a2m1uσ12udu,where C is a normalized constant, which satisfies(40)C=∫0∞1σ12aexp⁡∫0a2m1uσ12ududa-1.Substituting (36) into (39), the explicit expression for the stationary PDF of amplitude a can be obtained as follows:(41)pa=Caw02Dexp⁡-ηa2sin⁡pπ/2k52p+1a4w0p+1+3k32pa2w0p+1+24w02wp+148D,where w02=w2[1+ηwpcos⁡(pπ/2)].


## 4. Stochastic P-Bifurcation

Stochastic P-bifurcation refers to the changes of the number of peaks in the PDF curve; in order to obtain the critical parametric condition for stochastic P-bifurcation, the influences of the parameters for stochastic P-bifurcation of the system are analyzed by using the singularity theory below.For the sake of convenience, we can writep(a) as follows:(42)pa=CRa,D,w,η,p,k3,k5exp⁡Qa,D,w,η,p,k3,k5,in which(43)Ra,D,w,η,p,k3,k5=aw02DQa,D,w,η,p,k3,k5=-ηa2sin⁡pπ/2k52p+1a4w0p+1+3k32pa2w0p+1+24w02wp+148D.According to the singularity theory, the stationary PDF needs to satisfy the following two conditions:(44)∂pa∂a=0,∂2pa∂a2=0.Substituting (42) into (44), we can obtain the following condition [49]:(45)H=R′+RQ′=0,R′′+2R′Q′+RQ′′+RQ′2=0,where H represents the critical condition for the changes of the number of peaks in the PDF curve.Substituting (43) into (45), we can get the critical parametric condition for stochastic P-bifurcation of the system as follows:(46)D=18ηk52p+1w0p+1sin⁡pπ2a6+18ηk32p+1w0p+1sin⁡pπ2a4+ηw02wp+1sin⁡pπ2a2,where w02=w2[1+ηwpcos⁡(pπ/2)].And amplitudea satisfies(47)3k52pw0p+1a4+2k32pw0p+1a2+4w02wp+1=0.Taking the parameters asη=0.5, k3=7.8, k5=5.9, and w=1, according to (46) and (47), the transition set for stochastic P-bifurcation of the system with the unfolding parameters p and D can be obtained, as shown in Figure 4.Figure 4
The transition set curve of system (28).According to the singularity theory, the topological structures of the stationary PDF curves of different points (p,D) in the same region are qualitatively the same. Taking a point (p,D) in each region, all varieties of the stationary PDF curves which are qualitatively different could be obtained. It can be seen that the unfolding parameter p-D plane is divided into two subregions by the transition set curve; for convenience, each region in Figure 4 is marked with a number.Taking a given point (p,D) in each of the two subregions of Figure 4, the characteristics of stationary PDF curves are analyzed, and the corresponding results are shown in Figure 5.Figure 5
The stationary PDF curves for the amplitude of system (28).
(a)
parameter (p,D) in region 1 of Figure 4 (b)
parameter (p,D) in region 2 of Figure 4As can be seen from Figure4, the parametric region where the PDF curve appears bimodal is surrounded by an approximately triangular region. When the parameter (p,D) is taken in the region 1, the PDF curve has only a distinct peak, as shown in Figure 5(a); in region 2, the PDF curve has a distinct peak near the origin, but the probability is obviously not zero far away from the origin; there are both the equilibrium and limit cycle in the system simultaneously, as shown in Figure 5(b).The analysis results above show that the stationary PDF curves of the system amplitude in any two adjacent regions in Figure4 are qualitatively different. No matter the values of the unfolding parameters cross any line in the figure, the system will occur stochastic P-bifurcation behaviors, so the transition set curve is just the critical parametric condition for the stochastic P-bifurcation of the system, and the analytic results in Figure 5 are in good agreement with the numerical results by Monte Carlo simulation, which further verify the correctness of the theoretical analysis.


## 5. Conclusion

In this paper, we studied the stochastic P-bifurcation for axially moving of a bistable viscoelastic beam model with fractional derivatives of high order nonlinear terms under Gaussian white noise excitation. According to the minimum mean square error principle, we transformed the original system into an equivalent and simplified system and obtained the stationary PDF of the system amplitude using stochastic averaging. In addition, we obtained the critical parametric condition for stochastic P-bifurcation of the system using singularity theory; based on this, the system response can be maintained at the monostability or small amplitude near the equilibrium by selecting the appropriate unfolding parameters, providing theoretical guidance for system design in practical engineering, and avoiding the instability and damage caused by the large amplitude vibration or nonlinear jump phenomenon of the system. Finally, the numerical results by Monte Carlo simulation of the original system also verify the theoretical results obtained in this paper. We conclude that the orderp of fractional derivative and the noise intensity D can both cause stochastic P-bifurcation of the system and the number of peaks in the stationary PDF curve of system amplitude can change from 2 to 1 by selecting the appropriate unfolding parameters (p,D); it also shows that the method used in this paper is feasible to analyze the stochastic P-bifurcation behaviors of viscoelastic material system with fractional constitutive relation.



---

*Source: 6935095-2018-12-02.xml*