Stochastic P-bifurcation in a nonlinear viscoelastic beam model with fractional constitutive relation under colored noise excitation

Abstract

In this paper, we study the stochastic P-bifurcation problem for an axially moving bistable viscoelastic beam with fractional derivatives of high-order nonlinear terms under colored noise excitation. Firstly, using the principle for minimal 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. Secondly, we obtain the stationary probability density function of the system amplitude by the stochastic averaging and the singularity theory, we find the critical parametric conditions for stochastic P-bifurcation of the system amplitude. Finally, we analyze different types of the stationary probability density function curves of the system qualitatively by choosing parameters corresponding to each region divided by the transition set curves. 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.

Keywords

Viscoelastic beam fractional derivative colored noise stochastic P-bifurcation transition set Monte Carlo simulation

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 non-integer 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 can’t express the memory property of viscoelastic substances. The definition of fractional derivative contains convolution, which can well express a 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–3 In recent years, it has become the powerful mathematical tool in many disciplines, especially in the study of viscoelastic materials.

The classical integer operators do not have sufficient parameters to handle the different shapes of the hysteresis loops describing the behaviors of viscoelastic materials and structures, while 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 significant. In recent years, many scholars have done a lot of work and achieved fruitful results in this field: Nutting, Gemant, and Scott et al.^4–6 first proposed the fractional derivative models to study the constitutive relation of viscoelastic materials and the research on the viscoelastic material with fractional derivative is also increasing, and so far, it is still a research hotspot.^7–9 Wang et al.¹⁰ produced nanoscale crimped fibers using stuffer box crimping and bubble electrospinning and established the governing equations for nonlinear transverse vibration of an axially moving viscoelastic beam with finite deformation using the Hamiltonian principle, the obtained governing equations can be further used for numerical or analytical study of the crimping mechanism. Li et al.¹¹ solved a paradox in an electrochemical sensor by a fractal modification of the surface coverage model and elucidated a simple solution process to the fractal model. Wang and An adopted He’s fractional derivative which is defined through the variational iteration algorithm to describe a nonlinear vibration in microphysics, and used Ji–Huan He’s amplitude-frequency formulation to solve the fractional Duffing equation.¹² Xin et al.¹³ considered the existence of a positive periodic solution for a kind of high-order p-Laplacian neutral singular Rayleigh equation with variable coefficient based on Mawhin’s continuation theory, and verified the result by an example. Bagley and Torvik used fractional calculus to study the dynamic behavior of viscoelastic damping structure, and analyzed the responses of the system under general load and step load, respectively.^14,15 Zhang and Zhu analyzed the stability and dynamic response of viscoelastic belt under parametric excitation by the multiscale method.^16,17 Chen et al.^18–20 studied the dynamic behavior and steady-state response of axially accelerating viscoelastic beam by the Galerkin method and derived the differential equation of nonlinear vibration for axially moving viscoelastic rope, then pointed out that the damping of viscoelastic rope only exists in the nonlinear term. 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 analyzed the bifurcation behavior of the system caused by the order of the fractional derivative. Liu et al.²² studied the dynamic response of an axially moving viscoelastic beam under random disorder periodic excitation, then obtained the first-order expression of the solution by the multiscale method and carried out the stochastic jump phenomenon between the steady-state solutions. 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.²³ Leung et al.²⁴ studied the single mode dynamic characteristics of the nonlinear arch with the fractional derivative, then obtained the steady-state solution of the system based on the residual harmonic homotopy method and analyzed the influence of the parametric variation on the dynamic behaviors of the viscoelastic damping material. Galucio et al.²⁵ obtained the fractional derivative model to describe the viscoelasticity of the system based on the Timoshenko theory and Euler–Bernoulli hypothesis, then proposed a finite element formula for analyzing the sandwich beam of viscoelastic material with fractional derivative and the results were verified numerically.

Due to the complexity of fractional derivatives, analyzing them is difficult, the vibration characteristics of the parameters can only be qualitatively analyzed, and critical conditions of parametric influences cannot be found. These problems affect the analysis and design of such systems, in part because the stochastic P-bifurcation of bi-stability for the viscoelastic beam with fractional derivatives of high-order nonlinear terms under colored noise excitation has not been reported. Accordingly, we take the nonlinear vibration of a viscoelastic beam with fractional constitutive relation under colored noise excitation as an example, derive the governing equation of the viscoelastic beam and generate the nonlinear fractional oscillator based on the Galerkin discrete procedure, then obtain the transition set curves of the fractional order system and critical parametric conditions for stochastic P-bifurcation of the system by the singularity method, then analyze the different types of stationary probability density function (PDF) curves of the system in each region of the parameter plane. We compare numerical results from a Monte Carlo simulation with the analytical results obtained in this paper; it can be seen that the numerical solutions are in good agreement with the analytical solutions, verifying our theoretical analysis.

Equation of axially moving viscoelastic beam

There are many definitions of fractional derivatives, the following definitions are mainly introduced:

The Caputo derivative of the function $x (t)$ defined on the interval $[a, b]$ is formulated as

D a C p [x (t)] = \frac{1}{Γ (m - p)} \int_{a}^{t} \frac{x^{(m)} (u)}{{(t - u)}^{1 + p - m}} d u

(1)

where

p

represents the order of the derivative,

m - 1 < p \leq m, m \in N

t \in [a, b]

Γ (m)

is the Euler Gamma function, and

x^{(m)} (t)

is the m order derivative of

x (t)

The Riemann–Liouville derivative of the function $x (t)$ defined on the interval $[a, b]$ is formulated as

D a p [x (t)] = \frac{1}{Γ (m - p)} \frac{d^{m}}{d t^{m}} \int_{a}^{t} \frac{x (u)}{{(t - u)}^{1 + p - m}} d u

(2)

where

p

represents the order of the derivative,

m - 1 < p \leq m, m \in N

t \in [a, b]

Γ (m)

is the Euler Gamma function.

The Grunwald–Letnikov derivative of the function $x (t)$ is formulated as

D^{p} [x (t)]= {\frac{d^{p} [x (t)]}{d t^{p}} |}_{t = k h} = lim_{h \to 0} \frac{1}{h^{p}} \sum_{i = 0}^{k} {(-1)}^{i} (\begin{matrix} p \\ i \end{matrix}) x (t - i h) = h^{- p} \sum_{i = 0}^{k} C_{i}^{(p)} x (t_{k - i})

(3)

where

p

represents the order of the derivative,

C_{i}^{(p)}

is the Bernoulli coefficient, which satisfies the following recursive relation

C_{0}^{(p)} = 1, C_{i}^{(p)} = (1 - \frac{1 + p}{i}) C_{i - 1}^{(p)}, i = 0, 1, …

The fractal derivative of the function $x (t)$ is formulated as^26,27

\frac{d x (t)}{d t^{p}} = lim_{t^{'} - t} \frac{x (t) - x (t^{'})}{t^{p} - {t^{'}}^{p}}

(4)

where

p

represents the fractal in time and the order of the fractal derivative.

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 is defined in equation (1).

For a given physical system, due to the initial moment of the oscillator is $t = 0$ , so the following form of the Caputo derivative is often used

D 0 C^{p} [x (t)] = \frac{1}{Γ (m - p)} \int_{0}^{t} \frac{x^{(m)} (u)}{{(t - u)}^{1 + p - m}} d u

(5)

In this paper, the transverse vibration $y (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²⁴

{\begin{array}{l} \frac{\partial Q}{\partial x} = F (x, t) - ρ A \frac{\partial^{2} y}{\partial t^{2}} \\ \frac{\partial M}{\partial x} = Q - N \frac{\partial y}{\partial x} \end{array}

(6)

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 equation (6), we have

\frac{\partial^{2} M}{\partial x^{2}} + N \frac{\partial^{2} y}{\partial x^{2}} + ρ A \frac{\partial y^{2}}{\partial t^{2}} - F (x, t) = 0

(7)

Figure 1. Schematic stress diagram of viscoelastic beam.

Assuming the material of the beam obeys a viscoelastic constitutive relation in the form of fractional derivative

σ (x, z, t) = E_{0} ε (x, z, t) + E_{1} D 0 C^{p} [ε (x, z, t)] = E_{0} (1 + η D 0 C^{p}) ε (x, z, t) ≜ S ε (x, z, t)

(8)

where

η = E_{1} / E_{0}

is the material modulus ratio,

ε (x, z, t)

is axial strain component.

When the deformation of the beam is small, the axial strain $ε$ and lateral displacement $y (x, t)$ satisfy the relationship as follows

ε (x, z, t) = - z \frac{\partial^{2} y (x, t)}{\partial x^{2}}

(9)

Substituting equation (9) into equation (8) yields

σ (x, z, t) = S [- z \frac{\partial^{2} y (x, t)}{\partial x^{2}}]

(10)

The relationship between bending moment $M (x, t)$ and axial stress $σ (x, z, t)$ of the beam can be expressed as follows

M (x, t) = \int_{- h / 2}^{h / 2} z σ (x, z, t) d z

(11)

where

h

is the thickness of the beam.

From equations (10) and (11), the expression of the bending moment $M$ can be obtained as follows

M (x, t) = I [E_{0} (\frac{\partial^{2} y}{\partial x^{2}}) + E_{1} \cdot D 0 C^{p} (\frac{\partial^{2} y}{\partial x^{2}})] = I S (\frac{\partial^{2} y}{\partial x^{2}})

(12)

where

I = \frac{h^{3}}{12}

The expression of the horizontal tension is

\begin{matrix} N = \frac{E_{0} A}{2 L} \int_{0}^{L} [{(\frac{\partial y}{\partial x})}^{2} - {(\frac{\partial y}{\partial x})}^{4}] d x + \frac{E_{1} A}{2 L} \int_{0}^{L} D 0 C^{p} [{(\frac{\partial y}{\partial x})}^{2} - {(\frac{\partial y}{\partial x})}^{4}] d x \\ = \frac{A}{2 L} S \int_{0}^{L} [{(\frac{\partial y}{\partial x})}^{2} - {(\frac{\partial y}{\partial x})}^{4}] d x \end{matrix}

(13)

Substituting equations (12) and (13) into system (7), the system (7) can be rewritten as

I S (\frac{\partial^{4} y}{\partial x^{4}}) + S \frac{A}{2 L} \frac{\partial^{2} y}{\partial x^{2}} \int_{0}^{L} D 0 C^{p} [{(\frac{\partial y}{\partial x})}^{2} - {(\frac{\partial y}{\partial x})}^{4}] d x + ρ A \frac{\partial^{2} y}{\partial t^{2}} - F (x, t) = 0

(14)

The boundary conditions are

\frac{\partial^{2} y}{\partial x^{2}} = 0, when x = 0, x = L

(15)

According to the boundary conditions (15), the solution of system (14) can be expressed as the Fourier series

y (x, t) = \sum_{n = 1}^{\infty} u_{n} (t) \sin (\frac{n π}{L} x)

(16)

Assuming that the initial transverse vibration of the system is $y_{0} (x) = 0$ and the transverse load of the system satisfies the following form

F (x, t) = ρ A \cdot \sin (\frac{π}{L} x) ξ (t)

(17)

where

ξ (t)

is the Gaussian colored noise with zero mean, and its correlation function and power spectral density are as follows

{\begin{array}{l} R (t, t + h) = \frac{D}{τ} e^{- \frac{| h |}{τ}} \\ S (w) = \frac{D}{π (1 + τ^{2} w^{2})} \end{array}

(18)

where $D$ represents the noise intensity and $τ$ is the correlation time.

By the discrete format based on Galerkin method, the system (14) can be reduced to the ordinary differential equation as follows

\ddot{u} + k_{1} u + η k_{1} D 0 C^{p} u - k_{3} u^{3} - η k_{3} u D 0 C^{p} (u^{2}) + k_{5} u^{5} + η k_{5} u D 0 C^{p} (u^{4}) = ξ (t)

(19)

where

D 0 C^{p} u

is the

p

(

0 < p < 1

) order Caputo derivative of

u (t)

as is defined in equation (5), and

k_{1} = {(\frac{π}{L})}^{4} \frac{I E_{0}}{ρ A}, k_{3} = {(\frac{π}{L})}^{4} \frac{E_{0}}{4 ρ}, k_{5} = {(\frac{π}{L})}^{6} \frac{3 E_{0}}{16 ρ}

(20)

For convenience, the system (19) can be represented as follows

\ddot{u} + w^{2} u + η w^{2} D 0 C^{p} u - k_{3} u^{3} - η k_{3} u D 0 C^{p} (u^{2}) + k_{5} u^{5} + η k_{5} u D 0 C^{p} (u^{4}) = ξ (t)

(21)

where

w = \sqrt{k_{1}}

The fractional derivative term has contributions to both damping and restoring forces,^28,29 hence, introducing the following equivalent system

\ddot{u} + w^{2} u + η w^{2} [c (a) \dot{u} + w^{2} (a) u] - k_{3} u^{3} - η k_{3} u D 0 C^{p} (u^{2}) + k_{5} u^{5} + η k_{5} u D 0 C^{p} (u^{4}) = ξ (t)

(22)

where

c (a)

and

w^{2} (a)

are the coefficients of equivalent damping and restoring forces of the fractional derivative

D 0 C^{p} u

, respectively.

The error between systems (21) and (22) is

e = η w^{2} [c (a) \dot{u} + w^{2} (a) u - D 0 C^{p} u]

(23)

The necessary conditions for minimal mean square error are³⁰

\partial E [e^{2}] / \partial (c (a)) = 0, \partial E [e^{2}] / \partial (w (a)) = 0

(24)

Substituting equation (23) into equation (24) yields

{\begin{array}{l} E [c (a) {\dot{u}}^{2} + w^{2} (a) u \dot{u} - \dot{u} D 0 C^{p} u] = \lim_{T \to \infty} \frac{1}{T} \int_{0}^{T} (c (a) {\dot{u}}^{2} + w^{2} (a) u \dot{u} - \dot{u} D 0 C^{p} u) d t = 0 \\ E [c (a) \dot{u} u + w^{2} (a) u^{2} - u D 0 C^{p} u] = \lim_{T \to \infty} \frac{1}{T} \int_{0}^{T} (c (a) \dot{u} u + w^{2} (a) u^{2} - u D 0 C^{p} u) d t = 0 \end{array}

(25)

Assuming the solution of system (22) has the following form

u (t) = a (t) \cos φ (t)

(26)

where

φ (t) = w t + θ

Based on equation (26), we can obtain

\dot{u} (t) = - a (t) w sin φ (t)

(27)

Substituting equations (26) and (27) into equation (25) yields

\begin{array}{l} \lim_{T \to \infty} \frac{1}{T} \int_{0}^{T} (c (a) {\dot{u}}^{2} + w^{2} (a) u \dot{u} - \dot{u} D 0 C^{p} u) d t \\ = \lim_{T \to \infty} \frac{1}{T} \int_{0}^{T} (c (a) a^{2} w^{2} \sin^{2} φ - a^{2} w w^{2} (a) \sin φ \cos φ - a w \sin φ D 0 C^{p} u) d t \\ \approx \frac{c (a) a^{2} w^{2}}{2} + \frac{1}{Γ (1 - p)} \lim_{T \to \infty} \frac{1}{T} \int_{o}^{T} [(a w s i n φ) \int_{0}^{t} \frac{\dot{x} (t - τ)}{τ^{α}} d τ] d t \\ = \frac{c (a) a^{2} w^{2}}{2} - \frac{1}{Γ (1 - p)} \lim_{T \to \infty} \frac{1}{T} \int_{o}^{T} (a^{2} w^{2} \sin φ \int_{0}^{t} \frac{\sin φ cos (w τ) - \cos φ \sin (w τ)}{τ^{α}}) d t = 0 \end{array}

(28)

\begin{array}{l} \lim_{T \to \infty} \frac{1}{T} \int_{0}^{T} (c (a) \dot{u} u + w^{2} (a) u u - u D 0 C^{p} u) d t \\ = \lim_{T \to \infty} \frac{1}{T} \int_{0}^{T} (- c (a) a^{2} w \sin φ \cos φ + a^{2} w^{2} (a) \cos^{2} φ - a \cos φ D 0 C^{p} u) d t \\ \approx \frac{a^{2} w^{2} (a)}{2} - \frac{1}{Γ (1 - p)} \lim_{T \to \infty} \frac{1}{T} \int_{o}^{T} [(a \cos φ) \int_{0}^{t} \frac{\dot{x} (t - τ)}{τ^{α}} d τ] d t \\ = \frac{a^{2} w^{2} (a)}{2} + \frac{1}{Γ (1 - p)} \lim_{T \to \infty} \frac{1}{T} \int_{o}^{T} (a^{2} w \cos φ \int_{0}^{t} \frac{\sin φ cos (w τ) - \cos φ \sin (w τ)}{τ^{α}}) d t = 0 \end{array}

(29)

To simplify equations (28) and (29) further, asymptotic integrals are introduced as follows

{\begin{array}{l} \int_{0}^{t} \frac{\cos (w τ)}{τ^{p}} d τ = w^{p - 1} (Γ (1 - p) \sin (\frac{p π}{2}) + \frac{\sin (w t)}{{(w t)}^{p}} + o ({(w t)}^{- p - 1})) \\ \int_{0}^{t} \frac{\sin (w τ)}{τ^{p}} d τ = w^{p - 1} (Γ (1 - p) \cos (\frac{p π}{2}) - \frac{\cos (w t)}{{(w t)}^{p}} + o ({(w t)}^{- p - 1})) \end{array}

(30)

Substituting equation (30) into equations (28) and (29), and averaging them across $φ$ produces the ultimate forms of $c (a)$ and $w^{2} (a)$ as follows

c (a) = w^{p - 1} \sin (\frac{p π}{2}), w^{2} (a) = w^{p} \cos (\frac{p π}{2})

(31)

Therefore, the equivalent system associated with system (22) can be expressed as follows

\ddot{u} + w_{0}^{2} u + γ \dot{u} - k_{3} u^{3} - η k_{3} u D 0 C^{p} (u^{2}) + k_{5} u^{5} + η k_{5} u D 0 C p (u^{4}) = ξ (t)

(32)

where

w_{0}^{2} = w^{2} [1 + η w^{p} \cos (\frac{p π}{2})]

γ = η w^{p + 1} \sin (\frac{p π}{2})

Next, we consider the stochastic P-bifurcation of system (32) which comprises the fractional derivatives of high-order nonlinear terms, and analyze the influence of parametric variation on the system response.

The stationary PDF of amplitude

For the system (32), the material modulus ratio is given as $η = 0.5$ , coefficients of nonlinear terms are given as $k_{3} = 7.8$ , $k_{5} = 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 (32) at $D = 0$ .

As can be seen from Figure 2, 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 where $0 \leq p \leq 0.6817$ (equilibrium, as shown in Figure 3(a)); there are two attractors where $0.6818 \leq p \leq 1$ (equilibrium and limit cycle, as shown in Figure 3(b)).

Figure 3.

Phase diagrams of system (32) at $D = 0$ . (a) Equilibrium and (b) Equilibrium and limit cycle.

In order to obtain the stationary PDF of the amplitude of system (32), the following transformation is introduced³¹

{\begin{array}{l} u (t) = a (t) \cos Φ (t) \\ \dot{u} (t) = - a (t) w_{0} \sin Φ (t) \\ Φ (t) = w_{0} t + θ (t) \end{array}

(33)

where

a (t)

and

θ (t)

represent the amplitude process and the phase process of the system response, respectively, and they are all random processes.

Substituting equation (33) into equation (32), we can obtain

{\begin{array}{l} \frac{d a}{d t} = F_{11} (a, θ) + G_{11} (a, θ) ξ (t) \\ \frac{d θ}{d t} = F_{21} (a, θ) + G_{21} (a, θ) ξ (t) \end{array}

(34)

in which

{\begin{array}{l} F_{11} = - \frac{Δ \sin Φ}{w_{0}}, F_{21} = - \frac{Δ \cos Φ}{a w_{0}} \\ G_{11} = - \frac{\sin Φ}{w_{0}}, G_{21} = - \frac{\cos Φ}{a w_{0}} \end{array}

(35)

and

\begin{array}{l} Δ = η w^{p + 1} a w_{0} \sin (\frac{p π}{2}) \sin Φ - k_{3} a^{3} \cos^{3} Φ - η k_{3} \frac{a^{3}}{2} (2 w_{0})^{p} \cos Φ \cos (2 Φ + \frac{p π}{2}) \\ - k_{5} a^{5} \cos^{5} Φ - η k_{5} a^{5} \cos Φ [\frac{1}{8} {(4 w_{0})}^{p} \cos (4 Φ + \frac{p π}{2}) + \frac{1}{2} {(2 w_{0})}^{p} \cos (2 Φ + \frac{p π}{2})] \end{array}

(36)

Equation (34) 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

{\begin{array}{l} d a = m_{1} (a) d t + σ_{1} (a) d B (t) \\ d θ = m_{2} (a) d t + σ_{2} (a) d B (t) \end{array}

(37)

where

B (t)

is a standard Wiener process and

{\begin{array}{l} \begin{array}{l} m_{1} (a) = G_{11} + \frac{1}{2} \int_{- \infty}^{\infty} [\frac{\partial G_{11} (a, θ)}{\partial a} G_{11} (a, θ, t + h) + \frac{\partial G_{11} (a, θ)}{\partial θ} G_{21} (a, θ, t + h)] R (h) d h \\ = G_{11} + \frac{D \cos^{2} (Φ)}{a w^{2}^{0} (1 + w^{2}^{0} τ^{2})} \end{array} \\ \begin{array}{l} σ_{1}^{2} (a) = \int_{- \infty}^{\infty} G_{11} (a, θ) G_{11} (a, θ, t + h) R (h) d h \\ = \frac{2 D \sin^{2} (Φ)}{w^{2}^{0} (1 + w^{2}^{0} τ^{2})} \\ m_{2} (a) = G_{21} + \frac{1}{2} \int_{- \infty}^{\infty} [\frac{\partial G_{21} (a, θ)}{\partial a} G_{11} (a, θ, t + h) + \frac{\partial G_{21} (a, θ)}{\partial θ} G_{21} (a, θ, t + h)] R (h) d h \\ = G G_{21} - \frac{D \sin (2 Φ)}{a^{2} w^{2}^{0} (1 + w^{2}^{0} τ^{2})} \\ σ_{2}^{2} (a) = \int_{- \infty}^{\infty} G_{21} (a, θ) G_{21} (a, θ, t + h) R (h) d h \\ = \frac{2 D \cos^{2} (Φ)}{a^{2} w^{2}^{0} (1 + w^{2}^{0} τ^{2})} \end{array} \end{array}

(38)

By the stochastic averaging method,³² averaging equation (37) regarding $Φ$ , the following averaged Itô equation can be obtained as follows

{\begin{array}{l} d a = \bar{m_{1}} (a) d t + \bar{σ_{1}} (a) d B (t) \\ d θ = \bar{m_{2}} (a) d t + \bar{σ_{2}} (a) d B (t) \end{array}

(39)

where

{\begin{array}{l} \bar{m_{1}} (a) = - \frac{η k_{5} 2^{p} w_{0}^{p - 1} a^{5} \sin (p π / 2)}{8} - \frac{η k_{3} 2^{p} w_{0}^{p - 1} a^{3} \sin (p π / 2)}{8} - \frac{η w^{p + 1} a \sin (p π / 2)}{2} + \frac{D}{2 a w^{2}^{0} (1 + w^{2}^{0} τ^{2})} \\ {\bar{σ_{1}}}^{2} (a) = \frac{D}{w^{2}^{0} (1 + w^{2}^{0} τ^{2})} \\ \bar{m_{2}} (a) = \frac{k_{5} (2 η \cdot 2^{p} w_{0}^{p} \cos (p π / 2) + 5) a^{4}}{16 w_{0}} + \frac{k_{3} (η \cdot 2^{p} w_{0}^{p} \cos (p π / 2) + 3) a^{2}}{8 w_{0}} \\ {\bar{σ_{2}}}^{2} (a) = \frac{D}{a^{2} w^{2}^{0} (1 + w^{2}^{0} τ^{2})} \end{array}

(40)

Equations (39) and (40) show that $d a$ does not depend on $θ$ , the averaged Itô equation for $a (t)$ 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 wrote as

\frac{\partial p}{\partial t} = - \frac{\partial}{\partial a} [\bar{m_{1}} (a) p] + \frac{1}{2} \frac{\partial^{2}}{\partial a^{2}} [{\bar{σ_{1}}}^{2} (a) p]

(41)

The boundary conditions are

{\begin{array}{l} p = c c \in (- \infty, + \infty) when a =0 \\ p \to 0, \partial p / \partial a \to 0 when a \to \infty \end{array}

(42)

Thus, based on these boundary conditions (42), we can obtain the stationary PDF of amplitude $a$ as follows

p (a) = \frac{C}{{\bar{σ}}_{1}^{2} (a)} \exp [\int_{0}^{a} \frac{2 \bar{m_{1}} (u)}{{\bar{σ_{1}}}^{2} (u)} d u]

(43)

where C is a normalized constant, which satisfies

C = {[\int_{0}^{\infty} (\frac{1}{{\bar{σ}}_{1}^{2} (a)} \exp [\int_{0}^{a} \frac{2 \bar{m_{1}} (u)}{{\bar{σ_{1}}}^{2} (u)} d u]) d a]}^{- 1}

(44)

Substituting equation (40) into equation (43), the explicit expression for the stationary PDF of amplitude $a$ can be obtained as follows

p (a) = \frac{C a w_{0}^{2} (1 + w_{0}^{2} τ^{2})}{D} \exp [- \frac{η a^{2} (1 + w_{0}^{2} τ^{2}) \sin (p π / 2) (k_{5} 2^{p + 1} a^{4} w_{0}^{p + 1} + 3 k_{3} 2^{p} a^{2} w_{0}^{p + 1} + 24 w_{0}^{2} w^{p + 1})}{48 D}]

(45)

where

w_{0}^{2} = w^{2} [1 + η w^{p} \cos (\frac{p π}{2})]

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 conditions 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 write $p (a)$ as follows

p (a) = C R (a, D, w, η, p, τ, k_{3}, k_{5}) \exp [Q (a, D, w, η, p, τ, k_{3}, k_{5})]

(46)

in which

{\begin{array}{l} R (a, D, w, η, p, τ, k_{3}, k_{5}) = \frac{a w_{0}^{2} (1 + w_{0}^{2} τ^{2})}{D} \\ Q (a, D, w, η, p, τ, k_{3}, k_{5}) = - \frac{η a^{2} (1 + w_{0}^{2} τ^{2}) \sin (p π / 2) (k_{5} 2^{p + 1} a^{4} w_{0}^{p + 1} + 3 k_{3} 2^{p} a^{2} w_{0}^{p + 1} + 24 w_{0}^{2} w^{p + 1})}{48 D} \end{array}

(47)

According to the singularity theory,³³ the stationary PDF needs to satisfy the following two conditions

\frac{\partial p (a)}{\partial a} = 0, \frac{\partial^{2} p (a)}{\partial a^{2}} = 0

(48)

Substituting equation (46) into equation (48), we can obtain the following condition^34,35

H = {R^{'} + R Q^{'} = 0, R^{″} + 2 R^{'} Q^{'} + R Q^{″} + R {Q^{'}}^{2} = 0}

(49)

where

H

represents the condition for the changes of the number of peaks in the PDF curve.

Substituting equation (47) into equation (49), we can get the critical parametric conditions for stochastic P-bifurcation of the system as follows

D = \frac{1}{8} η \sin (p π / 2) (1 + w_{0}^{2} τ^{2}) a^{2} (k_{5} 2^{p + 1} w_{0}^{p + 1} a^{4} + k_{3} 2^{p + 1} w_{0}^{p + 1} a^{2} + 8 w_{0}^{2} w^{p + 1})

(50)

where

w_{0}^{2} = w^{2} [1 + η w^{p} \cos (\frac{p π}{2})]

And amplitude $a$ satisfies

3 k_{5} 2^{p} w_{0}^{p + 1} a^{4} + 2 k_{3} 2^{p} w_{0}^{p + 1} a^{2} + 4 w_{0}^{2} w^{p + 1} = 0

(51)

It can be seen from equations (50) and (51), the fractional order $p$ , the correlation time $τ$ and the noise intensity $D$ can each cause the stochastic P-bifurcation behavior of the system, the stochastic P-bifurcation of the system is analyzed by taking ( $p, D$ ) and ( $τ, D$ ) as the unfolding parameters separately below.

Taking ( $p, D$ ) as the unfolding parameters

Taking the parameters as $η = 0.5$ , $k_{3} = 7.8$ , $k_{5} = 5.9$ , $w = 1$ , $τ = 0.5$ , according to equations (50) and (51), 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.

Transition set curves (taking $p$ and $D$ as the unfolding parameters).

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 curves, 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.

Stationary PDF curves (taking $p$ and $D$ as the unfolding parameters). (a) Parameter ( $p, D$ ) in region 1 of Figure 4 and (b) Parameter ( $p, D$ ) in region 2 of Figure 4.

As can be seen from Figure 4, 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 a distinct peak near the origin, as shown in Figure 5(a); in region 2, the PDF curve still 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).

Taking ( $τ, D$ ) as the unfolding parameters

Taking the parameters as $η = 0.5$ , $k_{3} = 7.8$ , $k_{5} = 5.9$ , $w = 1$ , $p = 0.5$ , according to equations (50) and (51), the transition set for stochastic P-bifurcation of the system with the unfolding parameters $τ$ and $D$ can be obtained, as shown in Figure 6.

Figure 6.

Transition set curves (taking $τ$ and $D$ as the unfolding parameters).

As can be seen that the unfolding parameter $τ - D$ plane is divided into three subregions by the transition set curves; for convenience, each region in Figure 6 is marked with a number. We analyze the characteristics of stationary PDF curves at a point ( $τ, D$ ) in each of the three subregions of Figure 6, and show the corresponding results in Figure 7.

Figure 7.

Stationary PDF curves (taking $τ$ and $D$ as the unfolding parameters). (a) Parameter ( $τ, D$ ) in region 1 of Figure 6. (b) Parameter ( $τ, D$ ) in region 2 of Figure 6, and (c) Parameter ( $τ, D$ ) in region 3 of Figure 6.

As can be seen from Figure 6, the parameter region where the PDF curve appears bimodal is surrounded by two lines. When the parameter ( $τ, D$ ) is taken in region 1, the PDF curve has a distinct peak far away from the origin, as shown in Figure 7(a); in region 2, the PDF curve still has a distinct peak farther away from the origin, but the probability is not zero near the origin, as shown in Figure 7(b), both the equilibrium and limit cycle exist in the system simultaneously; in region 3, the PDF curve has a distinct peak near the origin, and the system has a small vibration near the equilibrium, as shown in Figure 7(c).

The above results show that the stationary PDF curves of the system amplitude in any two adjacent regions in Figures 4 and 6 are qualitatively different. Regardless of the exact values of the unfolding parameters, if they cross any line in these figures, the system will demonstrate stochastic P-bifurcation behavior. Therefore, the transition set curves are just the critical parametric conditions for the stochastic P-bifurcation of the system. The analytic results shown in Figures 5 and 7 are in good agreement with the numerical results obtained by Monte Carlo simulation, further verifying the theoretical analysis .

Conclusions

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 colored noise excitation. According to the minimal 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 conditions for stochastic P-bifurcation of the system using singularity theory, based on this, the system response can be maintained at the 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 fractional order $p$ , the correlation time $τ$ and the noise intensity $D$ can each cause stochastic P-bifurcation of the system and the number of peaks in the stationary PDF curves of system amplitude can change from 2 to 1 by selecting the appropriate unfolding parameters ( $p, D$ ) and ( $τ, 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.

Footnotes

Acknowledgements

The authors would like to thank the editors and reviewers for their conscientious reading of this article and their numerous comments for improvement which were extremely useful and helpful in modifying the article.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The authors gratefully acknowledge the financial support from the National 973 Project of China (grant no. 2014CB046805) and the Natural Science Foundation of China (grant no. 11372211).

References

Tan

MY.

Plane surface suddenly set in motion in a viscoelastic fluid with fractional Maxwell model. Acta Mech Sin 2002; 18: 342–349.

Yildirim

Gülkanat

Analytical approach to fractional Zakharov–Kuznetsov equations by He’s homotopy perturbation method. Commun Theor Phys 2010; 53: 1005–1010.

Roul

Analytical approach for nonlinear partial differential equations of fractional order. Commun Theor Phys 2013; 60: 269–277.

Nutting

PG.

A new general law deformation. J Franklin Inst 1921; 191: 679–685.

Gemant

Method of analyzing experimental results obtained from elasto-viscous bodies. Physics 1936; 7: 311–317.

Scott

Blair

An application of the theory of quasiproperties to the treatment of anomalous strain-stress relations. Philos Mag 1949; 40: 80–94.

Wang

Cao

Jing

, et al. Dynamic characteristics and stability of axially moving viscoelastic plate with piezoelectric layer. J Low Freq Noise Vib Active Control 2014; 33: 341–356.

Ghayesh

Kafiabad

Reid

Sub- and super-critical nonlinear dynamics of a harmonically excited axially moving beam. Int J Solids Struct 2012; 49: 227–243.

Ghayesh

MH.

Stability and bifurcations of an axially moving beam with an intermediate spring support. Nonlinear Dyn 2012; 69: 193–210.

10.

Wang

Liu

Kong

, et al. Nonlinear vibration mechanism for fabrication of crimped nanofibers with bubble electrospinning and stuffer box crimping method. Text Res J 2016; 87: 1–5.

11.

Tian

, et al. A fractal modification of the surface coverage model for an electrochemical arsenic sensor. Electrochim Acta 2018. DOI: 10.1016/j.electacta.2018.11.042.

12.

Wang

JY.

Amplitude-frequency relationship to a fractional Duffing oscillator arising in microphysics and tsunami motion.

J Low Freq Noise Vib Active Control. Epub ahead of print. DOI: 10.1177/1461348418795813.

13.

Xin

Yao

Cheng

ZB.

Positive periodic solution for higher-order p -Laplacian neutral singular Rayleigh equation with variable coefficient. Bound Value Probl 2017; 2017: 153.

14.

Bagley

Torvik

Fractional calculus – a different approach to the analysis of viscoelastically damped structures. AIAA J 2012; 21: 741–748.

15.

Bagley

Torvik

PJ.

Fractional calculus in the transient analysis of viscoelastically damped structures. AIAA J 2012; 23: 918–925.

16.

Zhang

JW.

Nonlinear vibration of parametrically excited moving belts, part I: dynamic response. J Appl Mech 1999; 66: 396–402.

17.

Zhang

JW.

Nonlinear vibration of parametrically excited moving belts, part II: stability analysis. J Appl Mech 1999; 66: 403–409.

18.

Chen

Yang

Cheng

CJ.

Dynamic stability of an axially accelerating viscoelastic beam. Eur J Mech 2004; 23: 659–666.

19.

Chen

Yang

XD.

Steady-state response of axially moving viscoelastic beams with pulsating speed: comparison of two nonlinear models. Int Solids Struct 2005; 42: 37–50.

20.

Chen

Yang

XD.

Transverse nonlinear dynamics of axially accelerating viscoelastic beams based on 4-term Galerkin truncation. Chaos Soliton Fract 2006; 27: 748–757.

21.

Leung

AYT

Yang

Chen

JY.

Parametric bifurcation of a viscoelastic column subject to axial harmonic force and time-delayed control. Comput Struct 2014; 136: 47–55.

22.

Liu

Dynamic responses of axially moving viscoelastic beam under a randomly disordered periodic excitation. J Sound Vib 2012; 331: 4045–4056.

23.

Yang

Fang

Stability in parametric resonance of an axially moving beam constituted by fractional order material. Arch Appl Mech 2012; 82: 1763–1770.

24.

Leung

AYT

Yang

Zhu

, et al. Steady state response of fractionally damped nonlinear viscoelastic arches by residue harmonic homotopy. Comput Struct 2013; 121: 10–21.

25.

Galucio

Deü

Ohayon

Finite element formulation of viscoelastic sandwich beams using fractional derivative operators. Comput Struct 2004; 33: 282–291.

26.

Cai

Chen

WX.

Characterizing the creep of viscoelastic materials by fractal derivative models. Int J Nonlin Mech 2016; 87: 58–63.

27.

JH.

Fractal calculus and its geometrical explanation. Results Phys 2018; 10: 272–276.

28.

Chen

Wang

, et al. Stationary response of Duffing oscillator with hardening stiffness and fractional derivative. Int J Nonlin Mech 2013; 48: 44–50.

29.

Chen

Zhuang

, et al. First-passage failure of single-degree-of-freedom nonlinear oscillators with fractional derivative. J Vib Control 2013; 19: 2154–2163.

30.

Chen

Zhu

WQ.

Stochastic response of fractional-order van der Pol oscillator. Theor Appl Mech Lett 2014; 4: 68–72.

31.

Spanos

Zeldin

BA.

Random vibration of systems with frequency-dependent parameters or fractional derivatives. J Eng Mech 1997; 123: 290–292.

32.

Zhu

WQ.

Random vibration. Beijing: Science Press, 1992 (in Chinese).

33.

Ling

FH.

Catastrophe theory and its applications. Shanghai: Shang Hai Jiao Tong University Press, 1987 (in Chinese).

34.

Hao

Stochastic P-bifurcations in tri-stable van der Pol-Duffing oscillator with multiplicative colored noise. Acta Phys Sin 2015; 64. DOI:10.7498/aps.64.060501 (in Chinese).

35.

Zhang

Zhao

JF.

Stochastic bifurcations of generalized Duffing-van der Pol system with fractional derivative under colored noise. Chin Phys B 2017; 26: 62–69.

Stochastic P-bifurcation in a nonlinear viscoelastic beam model with fractional constitutive relation under colored noise excitation

Abstract

Keywords

Introduction

Equation of axially moving viscoelastic beam

The stationary PDF of amplitude

Stochastic P-bifurcation

Taking ( p , D ) as the unfolding parameters

Taking ( τ , D ) as the unfolding parameters

Conclusions

Footnotes

Acknowledgements

Declaration of conflicting interests

Funding

References

Taking ( $p, D$ ) as the unfolding parameters

Taking ( $τ, D$ ) as the unfolding parameters