Response Analysis of Nonlinear Viscoelastic Energy Harvester with Bounded Noise Excitation

Abstract

Energy harvesting has become a popular topic in recent years. A number of studies have been conducted in the field of vibration energy harvesting system (VEHS). However, few studies have concentrated on viscoelastic energy harvesters driven by bounded noise excitation. In this paper, the stochastic response of a viscoelastic energy harvester subjected to bounded noise is discussed. Approximate solutions of the system were derived by utilizing the method of multiple scales, and the expressions of the mean square voltage (MSV) and mean output power (MOP) were obtained. The relation between the detuning frequency and first-order steady moment was first revealed. The effectiveness of the approach was verified by a good agreement between theoretical results and numerical results. Furthermore, the variations in the detuning frequency can result in the stochastic jump phenomenon, and stochastic bifurcation is induced with the changes in the viscoelastic parameter and detuning frequency. Finally, the impacts of system parameters on the MSV and the MOP were also analyzed.

1. Introduction

Vibration accompanied with the production of energy is a common phenomenon in nature involving human motion, bridge vibration, airfoil vibration, etc. For the sake of powering remote sensors sustainably [1,2], the vibration energy harvesting system (VEHS) has drawn much attention in recent years [3,4,5,6]. Meanwhile, the stochastic dynamical behavior of the VEHS with random excitation was investigated by utilizing different methods. For energy harvesters subjected to Gaussian white noise (GWN), Daqaq [7] indicated that the time constant ratio of a nonlinear energy harvester has an important effect on the harvesting performance by using the Fokker–Plank–Kolmogorov equation method. Mokem et al. [8] applied a stochastic averaging method to investigate the probabilistic responses of the sandwiched buckled beam, and analyzed the impacts of noise intensity on MSV. Jiang and Chen [9] adopted an equivalent linearization technique to obtain the equivalent linear system. Jin et al. [10] studied the stochastic responses of the nonlinear vibration energy harvester via employing an equivalent nonlinearization technique. Xu et al. [11] adopted stochastic averaging of an energy envelope to derive the analytical solutions of yjr mean-square voltage and mean output power of the nonlinear energy harvester.

However, GWN is too ideal to exist in actual environmental excitations. Colored noise has a wider application in a real environment. For colored noise, modified stochastic averaging, the multiple scales method and other methods have been confirmed to be better research techniques [12,13,14,15,16,17,18]. For example, Barton et al. [13] studied a nonlinear energy harvester by applying the stochastic averaging method. Osório and Daqaq [14] used the expansion of Jacobi–Anger to investigate the performance of the piezoelectric VEHS influenced by sinusoidal frequency variation. Bobryk et al. [15] conducted a study on the piezoelectric energy harvester regarding the mechanical oscillator, and found that colored noise plays a dramatic impact on energy harvesting. Liu et al. [16] presented the quasi-conservative stochastic averaging method to analyze the responses of a nonlinear energy harvester subject to colored noise. Zhang et al. [17] explored the stochastic responses of tri-stable energy harvesters by using the stochastic averaging of an energy envelope. Zhang et al. [19] analyzed the stochastic responses of a tri-stable energy harvester under harmonic excitation by using the method of Chebyshev polynomial approximation and showed that changes in the nonlinear coefficient have an influence on the output voltage. Mei et al. applied the perturbation method to study a tri-stable energy harvester in rotational motion to improve the harvesting performance and demonstrated that, at different rotational speeds, the existence of coefficient $K_c$ has different effects on the performance [20]. Zhang et al. [21] investigated a tri-stable energy harvester subjected to dual-frequency harmonic excitations via the harmonic balance method and showed that the harvesting performance increases with vibration resonance. In addition, based on the method of multiple scales, Huang et al. [18] analyzed the response process of a nonlinear multi-stable VEHS under narrow-band stochastic parametric excitations and discussed the stochastic bifurcation phenomenon. Jin and Zhang [22] explored the effects of time delay and feedback gain on the energy harvesting performance and demonstrated that negative feedback gain has a positive influence on the mean output power.

Recently, viscoelastic materials have been of great concern because of their application in various domains involving chemistry, physics, engineering, etc. [23,24,25] Nevertheless, influences of the viscoelastic term on the VEHS have not been discussed in the above articles. Different from conventional elastic materials, viscoelastic materials can not only show the attributes of spring, such as elasticity, but can also dissipate energy, such as damping [26]. Accompanied with the noise existing in nature, the dynamic behavior of viscoelastic systems with random excitation was investigated. Xu et al. [27] explored a viscoelastic system under narrow-band noise by using the multiple scales method. Fan and Shen [28] extended the multiple scales method and applied the approach to explore the dynamics of a nonlinear system with viscoelastic force. Zhu and Cai [29] analyzed stochastic responses of a viscoelastic system driven by broadband noise through a quasi-conservative stochastic averaging approach. Zhao et al. [30] considered the viscoelastic system under both additive GWN and multiplicative stochastic excitation, and discussed the stochastic P-bifurcation caused by parameters. Gu et al. [31] developed an approximate analytical procedure to explore the dynamic behavior of nonlinear viscoelastic systems under bounded noise excitation. Guo et al. [32] proposed a stochastic averaging approach to derive the steady expressions of amplitude for a nonlinear viscoelastic energy harvester and analyzed the effects of system parameters on the MOP.

As far as we are aware, few authors have paid attention to the dynamics of the viscoelastic VEHS under bounded noise. In this paper, we considered a nonlinear energy harvester with an integral viscoelastic term and bounded noise excitation, and utilized the multi-scale method to obtain the expressions of the mean voltage and mean output power. The effects of the viscoelastic parameter and noise intensity on system responses are discussed.

The framework of this paper is as follows. In Section 2, the analytical expressions of a viscoelastic energy harvester are derived by the multiple scales method. In Section 3, responses of the system are discussed, involving the stochastic jump phenomenon, stochastic bifurcation of the first-order steady-state moment and the influences of parameters on MSV and MOP. Meanwhile, the effectiveness of the theoretical method is verified via agreements between the theoretical results and the Monte Carlo simulation results. Section 4 provides concluding remarks.

2. Solving Procedure and Analytical Solutions

A model that symbolizes the dynamics of a family of piezoelectric VEHSs was investigated. The model can be seen as a base-excited viscoelastic-spring–damper system coupled to a capacitive energy harvesting circuit [7,33], as shown in Figure 1.

The electromechanically coupling equations of motion with viscoelastic force are the following: equationsection

$$\begin{array}{cc}\hfill & M\ddot{\overline{X}}(\tau )+{\overline{k}}_{1}H(\overline{\beta },\tau ,\overline{X}(\tau ))+{\overline{k}}_2\dot{\overline{X}}(\tau )+{\overline{\delta }}_1\overline{X}(\tau )+{\overline{\delta }}_2{\overline{X}}^3(\tau )-\overline{\zeta }\overline{V}(\tau )=-M{\ddot{\overline{X}}}_b(\tau ),\hfill \\ \hfill & C_p\dot{\overline{V}}(\tau )+\frac{\overline{V}(\tau )}{R}=-\overline{\lambda }\dot{\overline{X}}(\tau ),\hfill \end{array}$$

(1)

in which, $\overline{X}$ is the displacement and the dot is a derivative about time $\tau$. $\dot{\overline{X}}$ and ${\ddot{\overline{X}}}_b$ denote the velocity and the base acceleration of the mass M, respectively. $\overline{V}$ is the electric voltage measured across the equivalent resistive load R. $C_p$ is the effective capacitance of the piezoelectric element. ${\overline{k}}_1$, $\beta$, ${\overline{k}}_2$, ${\overline{\delta }}_1$, ${\overline{\delta }}_2$, $\overline{\zeta }$ and $\overline{\lambda }$ represent the viscoelastic coefficient, the viscoelastic parameter, the damping force coefficient, the linear stiffness coefficient, the nonlinear stiffness coefficient, the linear electromechanical coupling coefficient and the piezoelectric coupling coefficient in the electrical circuit, respectively. H is the viscoelastic term taking the form of [34]

$$\begin{array}{c}\hfill H(\overline{\beta },\tau ,\overline{X}(\tau ))={\int }_0^{\tau }\frac{1}{\overline{\beta }}\exp (-\frac{\tau -q}{\overline{\beta }})\overline{X}(q)dq.\end{array}$$

(2)

$-{\ddot{\overline{X}}}_b$ is considered as a bounded stochastic process that is governed by $\xi (\tau )=\hat{F}\cos [\mathrm{\Omega }\tau +\overline{\gamma }W(\tau )]$. Here, $W(\tau )$ is the standard Winner process. $\hat{F}$ denotes the deterministic amplitude of the bounded noise with the frequency $\mathrm{\Omega }$ and the noise intensity $\overline{\gamma }$.

Firstly, we introduce the non-dimensional parameters

$$\begin{array}{cc}\hfill & t=\omega \tau ,x=\frac{\overline{X}}{l_c},x_b=\frac{{\overline{X}}_b}{l_c},{\hat{k}}_1=\frac{{\overline{k}}_1}{{\overline{\delta }}_1},\omega =\sqrt{\frac{{\overline{\delta }}_1}{M}},{\hat{k}}_2=\frac{{\overline{k}}_2}{\sqrt{M}{\overline{\delta }}_1},\hfill \\ \hfill & {\hat{\delta }}_2=\frac{{\overline{\delta }}_2{l}_c^2}{{\overline{\delta }}_1},V=\frac{C_p\overline{V}}{\overline{\zeta }{l}_c},\hat{\zeta }=\frac{{\overline{\zeta }}^2}{{\overline{\delta }}_1{C}_p},\mu =\frac{1}{\omega R{C}_p},\lambda =\frac{\overline{\lambda }}{\overline{\zeta }}.\hfill \end{array}$$

(3)

Then, applying the transform in the viscoelastic force and bounded noise, we obtain the following equations:

$$\begin{array}{cc}\hfill H(\overline{\beta },\tau ,\frac{\overline{X}(\tau )}{l_c})& ={\int }_0^{\tau }\frac{1}{\overline{\beta }}\exp (-\frac{\tau -q}{\overline{\beta }})\frac{\overline{X}(q)}{l_c}dq={\int }_0^t\frac{1}{\overline{\beta }}\exp [-\frac{1}{\overline{\beta }}(\frac{t}{\omega })-\frac{s}{\omega }]x(s)d(\frac{s}{\omega })\hfill \\ \hfill & ={\int }_0^t\frac{1}{\overline{\beta }\omega }\exp (-\frac{t-s}{\overline{\beta }\omega })x(s)ds\triangleq H(\beta ,t,x(t)),\hfill \\ \hfill \overline{\gamma }W(\tau )& =\overline{\gamma }W(\frac{t}{\omega })=\frac{\overline{\gamma }}{\sqrt{\omega }}W(t)\triangleq \gamma W(t),\hfill \end{array}$$

(4)

where $\beta =\overline{\beta }\omega$ and $\gamma =\overline{\gamma }/\sqrt{\omega }$.

Subsequently, the dimensionless formula of motion can be written as

$$\begin{array}{cc}\hfill & \ddot{x}+{\hat{k}}_{1}H(\beta ,t,x)+{\hat{k}}_2\dot{x}+x+{\hat{\delta }}_2{x}^3-\hat{\zeta }V=\hat{F}\cos [\frac{\mathrm{\Omega }t}{\omega }+\gamma W(t)],\hfill \\ \hfill & \dot{V}+\mu V=-\lambda \dot{x},\hfill \end{array}$$

(5)

where $l_c$, $\omega$ and $\mu$ denote the length scale, the natural frequency of the system and the dimensionless time constant ratio, respectively, and the others correspond to the dimensional form of Equation (1).

Following the method of multiple scales, a small parameter $\epsilon$ will be introduced to acquire the parameters: $\epsilon \alpha ={\hat{k}}_1$, $\epsilon k_2={\hat{k}}_2$, $\epsilon {\delta }_2={\hat{\delta }}_2$, $\epsilon \theta =\hat{\zeta }$ and $\epsilon F=\hat{F}$. Then, Equation (5) can be transformed into the following form:

$$\begin{array}{cc}\hfill & \ddot{x}+\epsilon \alpha H(\beta ,t,x)+\epsilon k_2\dot{x}+x+\epsilon {\delta }_2{x}^3-\epsilon \theta V=\epsilon F\cos [\frac{\mathrm{\Omega }t}{\omega }+\gamma W(t)],\hfill \\ \hfill & \dot{V}+\mu V=-\lambda \dot{x}.\hfill \end{array}$$

(6)

Moreover, we introduce the time scales as $T_0=t$ and $T_1=\epsilon t$ and expand the solutions x and V of system (6), i.e.,

$$\begin{array}{cc}\hfill & x(t)=x_0(T_0,T_1)+\epsilon x_1(T_0,T_1)+O({\epsilon }^2),\hfill \\ \hfill & V(t)=V_0(T_0,T_1)+\epsilon V_1(T_0,T_1)+O({\epsilon }^2).\hfill \end{array}$$

(7)

For convenience, with the introduction of partial derivative operators $D_0=\partial /\partial T_0$ and $D_1=\partial /\partial T_1$, we can transform the ordinary derivative into a partial derivative expansion:

$$\left\{\begin{array}{c}\frac{d}{dt}=D_0+\epsilon D_1+O\left({\epsilon }^2\right)\\ \frac{d^2}{d{t}^2}=D_0^2+2\epsilon D_0{D}_1+O\left({\epsilon }^2\right)\end{array}\right..$$

(8)

Substituting Equations (7) and (8) into Equation (6), and collecting the terms according to the uniform order of $\epsilon$, we immediately obtain the formulas at first order as

$$\begin{array}{c}{D}_0^2{x}_0+x_0=0,\hfill \\ D_0{V}_0+\mu V_0=-\lambda D_0{x}_0,\hfill \end{array}$$

(9)

and at second order as

$$\begin{array}{c}{D}_0^2{x}_1+x_1=-2D_0{D}_1{x}_0-\alpha H\left(\beta ,T_0,x_0\right)-k_2{D}_0{x}_0-{\delta }_2{x}_0^3+\theta V_0+F\cos [\frac{\mathrm{\Omega }t}{\omega }+\gamma W(t)],\hfill \\ D_0{V}_1+\mu V_1=-D_1{V}_0-\lambda D_1{x}_0-\lambda D_0{x}_1,\hfill \end{array}$$

(10)

The general solutions of Equation (9) can be derived as

$$\begin{array}{c}{x}_0\left(T_0,T_1\right)=A_0\left(T_1\right)e^{i{T}_0}+\overline{A_0\left(T_1\right)}{e}^{-i{T}_0},\hfill \\ V_0\left(T_0,T_1\right)=B_0\left(T_1\right)e^{i{T}_0}+\overline{B_0\left(T_1\right)}{e}^{-i{T}_0},\hfill \end{array}$$

(11)

where $A_0$ is the amplitude and $\overline{A_0}$ is the complex conjugate of $A_0$. Substituting Equation (11) into the second expression of Equation (9), we find that the amplitude $B_0(T_1)$ can be expressed as

$$B_0(T_1)=-\frac{i\lambda }{\mu +i}{A}_0(T_1).$$

(12)

The viscoelastic force is

$$\begin{array}{cc}\hfill H\left(\beta ,T_0,x_0\right)& =\frac{1}{\beta }{\int }_0^{T_0}\exp (-\frac{T_0-s}{\beta })x_0(s)ds\hfill \\ \hfill & =\frac{1-i\beta }{1+{\beta }^2}{A}_0\left(T_1\right)\exp (i{T}_0)-\frac{1-i\beta }{1+{\beta }^2}{A}_0\left(T_1\right)\exp (-\frac{T_0}{\beta })+cc,\hfill \end{array}$$

(13)

where $cc$ represents the conjugate force.

Substituting Equations (11)–(13) into the first expression of Equation (10) and vanishing the secular terms, we have

$$\begin{array}{c}\hfill -2i{A}_0^{\prime }-\frac{1-i\beta }{1+{\beta }^2}\alpha A_0-\frac{i\lambda \theta }{\mu +i}{A}_0-i{k}_2{A}_0-3{\delta }_2{A}_0^2\overline{A_0}+\frac{F}{2}\exp \left[i\left(\frac{\mathrm{\Omega }{T}_0}{\omega }-T_0+\gamma W(T_1)\right)\right]=0,\end{array}$$

(14)

where $A_0^{\prime }$ is the derivative with respect to $T_1$.

For the sake of analyzing the primary resonance, the detuning frequency $\sigma$ is introduced to measure the excitation frequency $\mathrm{\Omega }$ as

$$\mathrm{\Omega }=\omega (1+\epsilon \sigma ).$$

(15)

According to Equations (14) and (15), we obtain

$$-2i{A}_0^{\prime }-\frac{1-i\beta }{1+{\beta }^2}\alpha A_0-\frac{i\lambda \theta }{\mu +i}{A}_0-i{k}_2{A}_0-3{\delta }_2{A}_0^2\overline{A_0}+\frac{F}{2}\exp \left[i\left(\sigma T_1+\gamma W(T_1)\right)\right]=0.$$

(16)

To calculate the function $A_0$, we assume that

$$A_0(T_1)=\frac{a(T_1)}{2}\exp \left[i\varphi (T_1)\right].$$

(17)

Substituting Equation (17) into Equation (16) and simplifying the equations to obtain the expressions of the real and imaginary parts, one obtains the following differential equations:

$$\begin{array}{c}{D}_{1}a=\frac{\alpha \beta }{2(1+{\beta }^2)}a-\frac{\lambda \theta \mu }{2\left({\mu }^2+1\right)}a-\frac{k_{2}a}{2}+\frac{F}{2}\sin \eta ,\hfill \\ D_1\eta =\sigma -\frac{\alpha }{2(1+{\beta }^2)}-\frac{\lambda \theta }{2\left({\mu }^2+1\right)}-\frac{3{\delta }_2{a}^2}{8}+\frac{F}{2a}\cos \eta +\gamma W^{\prime }\left(T_1\right),\hfill \end{array}$$

(18)

where $\eta (T_1)=\sigma T_1+\gamma W(t)-\varphi (T_1)$.

Combining Equations (11), (12), (17) and (18), the first-order approximate solutions of Equation (6) will be derived as follows:

$$x(t)=a(\epsilon t)\cos (t+\varphi (\epsilon t))+O(\epsilon ),$$

(19)

$$V(t)=\frac{\lambda }{\sqrt{{\mu }^2+1}}a(\epsilon t)\cos \left(t+\varphi (\epsilon t)+\arctan \mu \right)+O(\epsilon ).$$

(20)

Steady-State Solution

When $\gamma =0$, the noise becomes a harmonic excitation. According to the conditions $D_{1}a=0$, $D_1\eta =0$, the steady-state solutions $a=a_0$, $\eta ={\eta }_0$ can be obtained as follows:

$$\begin{array}{c}\left[\frac{\alpha \beta }{2{(1+\beta )}^2}-\frac{\lambda \theta \mu }{2\left({\mu }^2+1\right)}-\frac{k_2}{2}\right]a_0+\frac{F}{2}\sin {\eta }_0=0,\hfill \\ \left[\sigma -\frac{\alpha }{2(1+{\beta }^2)}-\frac{\lambda \theta }{2\left({\mu }^2+1\right)}\right]a_0-\frac{3{\delta }_2}{8}{a}_0^3+\frac{F}{2}\cos {\eta }_0=0.\hfill \end{array}$$

(21)

When $\gamma \ne 0$, the solutions of Equation (18) have the hypothesis of the following form:

$$\begin{array}{c}a=a_0+a_1,\hfill \\ \eta ={\eta }_0+{\eta }_1,\hfill \end{array}$$

(22)

where $a_0,{\eta }_0$ are satisfying Equation (21) and $a_1,{\eta }_1$ are the small perturbation terms. First-order Taylor expansion is performed on the equations, substituting Equation (22) into Equation (18). Vanishing the high-order small items, we obtain the linearization of Equation (18) at $(a_0,{\eta }_0)$ as

$$\begin{array}{c}d{a}_1=(M_{11}{a}_1+M_{12}{\eta }_1)d{T}_1,\hfill \\ d{\eta }_1=(M_{21}{a}_1+M_{22}{\eta }_1)d{T}_1+\gamma dW(T_1),\hfill \end{array}$$

(23)

where

$$\begin{array}{cc}\hfill M_{11}=& \frac{\alpha \beta }{2(1+{\beta }^2)}-\frac{\lambda \theta \mu }{2\left({\mu }^2+1\right)}-\frac{k_2}{2},\hfill \\ \hfill M_{12}=& -\left[\sigma -\frac{\alpha }{2(1+{\beta }^2)}\right]a_0+\frac{\lambda \theta }{2\left({\mu }^2+1\right)}{a}_0+\frac{3{\delta }_2}{8}{a}_0^3,\hfill \\ \hfill M_{21}=& \frac{\sigma }{a_0}-\frac{\alpha }{2(1+{\beta }^2)a_0}-\frac{\lambda \theta }{2\left({\mu }^2+1\right)a_0}-\frac{9{\delta }_2}{8}{a}_0,\hfill \\ \hfill M_{22}=& \frac{\alpha \beta }{2(1+{\beta }^2)}-\frac{\lambda \theta \mu }{2\left({\mu }^2+1\right)}-\frac{k_2}{2}.\hfill \end{array}$$

Considering the steady-state situation and applying the method of moment, we have

$$\frac{dE(a_1)}{d{T}_1}=\frac{dE({\eta }_1)}{d{T}_1}=\frac{dE(a_1^2)}{d{T}_1}=\frac{dE(a_1{\eta }_1)}{d{T}_1}=\frac{dE({\eta }_1^2)}{d{T}_1}=0.$$

According to Equation (23), we obtain that the second-order moments can be solved by the following equations:

$$\begin{array}{c}\frac{dE\left(a_1^2\right)}{d{T}_1}=2M_{11}E\left(a_1^2\right)+2M_{12}E\left(a_1{\eta }_1\right)=0,\hfill \\ \frac{dE\left(a_1{\eta }_1\right)}{d{T}_1}=M_{21}E\left(a_1^2\right)+\left(M_{11}+M_{22}\right)E\left(a_1{\eta }_1\right)+M_{12}E\left({\eta }_1^2\right)=0,\hfill \\ \frac{dE\left({\eta }_1^2\right)}{d{T}_1}=2M_{21}E\left(a_1{\eta }_1\right)+2M_{22}E\left({\eta }_1^2\right)+{\gamma }^2=0,\hfill \end{array}$$

(24)

where $E(•)$ represents the mathematical expectation.

The first-order steady-state moments are derived as

$$E(a_1)=0,E({\eta }_1)=0,E(a)=a_0,E(\eta )={\eta }_0.$$

(25)

By applying Cramer’s rule, the second-order moments can be derived as

$$\begin{array}{c}E\left(a_1^2\right)=\frac{-2{\gamma }^2{M}_{12}^2}{\det },\hfill \\ E\left({\eta }_1^2\right)=\frac{-M_{12}{\gamma }^2+2M_{11}{M}_{21}E\left(a_1^2\right)}{2M_{12}{M}_{22}},\hfill \\ E\left(a^2\right)=a_0^2+E\left(a_1^2\right),\hfill \\ E\left({\eta }^2\right)={\eta }_0^2+E\left({\eta }_1^2\right),\hfill \end{array}$$

(26)

where $\det =4\left[M_{11}^2{M}_{22}+M_{11}{M}_{22}^2-M_{11}{M}_{12}{M}_{21}-M_{12}{M}_{21}{M}_{22}\right]$. Then, according to Equation (20), we obtain the expressions of the mean and mean square voltage as

$$\begin{array}{c}E(V)=\frac{\lambda }{\sqrt{{\mu }^2+1}}E(a),\hfill \\ E(V^2)=\frac{{\lambda }^2}{{\mu }^2+1}E(a^2).\hfill \end{array}$$

(27)

Since the output power can be expressed as $P=\mu \theta V^2$, the mean output power is calculated as

$$E(P)=\mu \theta E(V^2).$$

(28)

3. Response Analysis of System

3.1. Numerical Simulation Results

This subsection verifies the effectiveness of the devised method by using the Monte Carlo approach. Some system parameters in Equation (6) are fixed as $\epsilon =0.1$, ${\omega }^2=1.5$, $k_2=1.0$, $F=5.0$, $\alpha =2.0$, $\beta =2.0$, $\lambda =0.50$, ${\delta }_2=1.0$, $\mu =1.50$, unless otherwise mentioned. The other parameters will be specified in subsequent analyses. Based on the Hurwitz criterion, Figure 2 reveals the distribution of the first-order steady moment $E(a)$ with the parameters $\theta =0.5$, $\lambda =0.50$, ${\delta }_2=1.0$, $\mu =1.50$, $\gamma =0.10$. From Figure 2, the solutions with a blue mark are stable and the solutions with a red mark are non-stable. Figure 3 shows the comparison of numerical and theoretical results, which proves the validity of the approach of multiple scales in this paper with the parameters $\theta =0.50$, $\gamma =0.10$.

3.2. Stochastic Jump and Stochastic Bifurcation

The variations in the system regarding different detuning frequencies in time history are presented in Figure 4. One can see from Figure 4 that, as the detuning frequency $\sigma$ increases, the displacement x, velocity $\dot{x}$ and voltage V increase slowly and then decrease rapidly in $\sigma \in (2.0,2.5)$. This result also indicates the phenomenon of the stochastic jump.

From Figure 2 and Figure 3, we know that the number of the first-order steady-state moment $E(a)$ can change from one to three with an increase in the detuning frequency. By considering the viscoelastic parameter $\beta$ and the detuning frequency $\sigma$, we fix the parameters $\gamma =0.10$, $\theta =0.50$ and show the variation diagram in the $\sigma -\beta$ parameter plane. One can see from Figure 5a that the $E(a)$ can transform one steady solution into two steady solutions from area I to area II. If $\sigma <2.25$, $E(a)$ has always one steady-state solution. When $\sigma >2.25$, the number of solutions may change. The results show that $\beta$ and $\sigma$ can induce stochastic bifurcation. Figure 5b shows the variations in $E(a)$ with different $\beta$ and $\sigma$.

3.3. Impacts of System Parameters on the Mean Square Voltage and the Mean Output Power

In this subsection, we focus on the impacts of noise intensity $\gamma$, linear electromechanical coupling coefficient $\theta$, piezoelectric coupling coefficient $\lambda$ and nonlinear stiffness coefficient ${\delta }_2$ on the MSV and MOP. Figure 6 describes the relevance between the MSV and the noise intensity. With the increase in $\gamma$ and $\lambda$, the MSV keeps increasing. However, the MSV decreases when $\theta$ and ${\delta }_2$ increase.

In Figure 7, the MOP keeps increasing as the noise intensity $\gamma$, linear electromechanical coupling coefficient $\theta$, piezoelectric coupling coefficient $\lambda$ and time constant ratio $\mu$ increase. The results indicate that the above parameters have a positive impact on the mean output power.

4. Conclusions

This paper conducted a response analysis of an energy harvester with a viscoelastic element driven by bounded noise excitation. We simultaneously considered the effects of viscoelasticity and bounded noise on the system, and derived the analytical expressions of the steady-state response of the VEHS by using the multiple scales method.

First, the effectiveness of the method was verified by a good agreement between the analytical results and numerical results. Additionally, the stochastic jump phenomenon occurs in the range $(2.0,2.5)$ of the detuning frequency. The bifurcation diagram was given in the $\sigma -\beta$ parameter plane. The results demonstrated that the viscoelastic parameter and detuning frequency can induce stochastic bifurcation. The theoretical analysis in this paper showed that, for the same detuning frequency, there are possibly multiple stable steady-state average locomotions, which results in multiple solutions in the resonant region. Finally, the influences of parameters $\gamma$, $\theta$, $\lambda$, $\mu$ and ${\delta }_2$ on the mean square voltage (MSV) and the mean output power (MOP) were analyzed. Through this paper, the following conclusions are obtained. The MSV decreases as $\theta$ and ${\delta }_2$ increase. However, the MSV increases as the noise intensity and piezoelectric coupling coefficient increase. The MOP increases as $\gamma$, $\theta$, $\mu$ and $\lambda$ increase.