# Resonance Mechanism of Nonlinear Vibrational Multistable Energy Harvesters under Narrow-Band Stochastic Parametric Excitations

Authors: Dongmei Huang; Shengxi Zhou; Zhichun Yang
Journal: Complexity
 (2019)

---

## Abstract

To improve energy harvesting performance, this paper investigates the resonance mechanism of nonlinear vibrational multistable energy harvesters under narrow-band stochastic parametric excitations. Based on the method of multiple scales, the largest Lyapunov exponent which determines the stability of the trivial steady-state solutions is derived. The first kind modified Bessel function is utilized to derive the solutions of the responses of multistable energy harvesters. Then, the first-order and second-order nontrivial steady-state moments of multistable energy harvesters are considered. To explore the stochastic bifurcation phenomenon between the nontrivial and trivial steady-state solutions, the Fokker–Planck–Kolmogorov equation corresponding to the two-dimensional Itô stochastic differential equations is solved by using the finite difference method. In addition, the mechanism of the stochastic bifurcation of multistable energy harvesters is analyzed for revealing their unique dynamic response characteristics.


---

## Body

## 1. Introduction

Vibration energy harvesting is expected to replace some small traditional chemical batteries and permanently supply power for the low-powered wireless sensors and actuators and also can promote structural health monitoring [1–5]. Recently, vibrational nonlinear energy harvesters have been receiving more and more attention because they have the better energy harvesting performance than traditional linear energy harvesters under time-varying ambient environmental excitations [4–9]. By far, different kinds of high-performance nonlinear energy harvesters have been designed based on geometric nonlinearity, additional nonlinear magnetic force, residual thermal stress, or active/passive control strategies [9–12].In past several years, the novel methods to broaden the bandwidth are achieved by introducing nonlinearity and multioscillator structures, and the nonlinear behaviors may extend the operating bandwidth of the energy harvester, allowing for efficient energy harvesting under random vibrations. Nonlinear energy harvesters with monostable, bistable, tristable, and even quadstable characteristics have been extensively studied to improve the harvesting performance in natural. Cottone et al. [12] designed a bistable energy harvester (BEH) based on additional nonlinear magnetic force and a cantilever beam structure, and they experimentally verified its advantage for harvesting energy from random vibrations. Later on, Litak et al. [13, 14] deeply investigated and verified the energy harvesting improvement of the BEH under random base excitations. Chen and Jiang [15] presented the internal resonance mechanism to enhance vibration energy harvesting. Wang et al. [16] experimentally tested a BEH to harvest energy from human motions, and a better energy harvesting performance over the linear monostable harvester was archived. Then, Zhou et al. [17] designed a tristable energy harvester (TEH) to improve the broadband energy harvesting performance from low-level ambient vibrations. The TEH physically has three potential energy wells, which may have lower potential barriers than the BEH [18]. The advantages of the TEH in broadband vibration energy harvesting over the BEH under low-level harmonic and random base excitations were experimentally and numerically verified [17, 19–21]. For example, Xu et al. [22] revealed the stochastic resonance mechanism of the TEH. Tékam et al. [23] studied the influence of the fractional-order viscoelastic material of a TEH on its energy harvesting performance. The influence mechanism of the system parameters of the TEH on the dynamic responses and the energy harvesting performance was completely provided in [24–26].In recent years, based on the magnetic coupled energy harvesting structures, Zhou et al. [27, 28] developed multistable energy harvesters with four or five stable equilibrium positions, and a better energy harvesting performance over the BEH was verified. For the real application, Younesian and Alam [29] designed broadband multistable wave energy harvesters based on the multistable mechanism, and the designed harvesters were verified to have high-efficiency energy harvesting performance. Gao et al. [30] presented a multistable energy harvester based on the magnetic levitation oscillation, and it has the quadstable configuration. In the railway test, the RMS power output is 440.98 mW, which is large enough to power some small wireless sensor nodes. Under harmonic base excitations, the theoretical solutions of nonlinear vibrational multistable energy harvesters were provided by Huang et al. [31]. These research studies verify the high-efficiency energy harvesting performance of multistable energy harvesters in experiments and real environmental tests.In practice, the operation of the piezoelectric system is highly susceptible to the various random factors from internal and external environment fluctuation [32–39], for example, the devices which are manufactured by the piezoelectric ceramic widely used in meteorological observation and household appliances; the normal operation of these devices are often influenced by the atmospheric turbulence, voltage fluctuation, and even other random disturbance. Thus, the nonlinear energy harvesters, especially in stochastic vibratory environments, have become an unavoidable topic [34–40]. Xu et al. [34] developed the stochastic averaging technique for the energy harvester under Gaussian white noise excitation to evaluate the mean-square electric voltage and mean output power. They also considered the responses of bistable energy harvester under additive and multiplicative white noises [35]. Jiang and Chen [36, 37] analyzed the response of nonlinear energy harvesters subject to external Gaussian white noise and parametric excitations. He and Daqaq [38] illustrated the optimal potential shape on the mean output power through the statistical linearization techniques. Recently, Liu et al. [39] proposed a new quasiconservative stochastic averaging method to discuss the probabilistic response of the nonlinear energy harvesting system subjected to correlated Gaussian colored noise. Zhang and Jin [40] discussed the stationary probability density and signal-to-noise ratio of a piezoelectric energy harvester with correlated additive and multiplicative colored noises. He et al. [41] explored the parameter resonance mechanism of nanocomposite energy harvesters based on Galerkin’s method and the perturbation principle of Poincaré–Lindstedt. They found that a functionally graded reinforcement has a significant influence on the bifurcation buckling, the postbuckling path, the output voltage, and the harvested power.Most previous results studied the energy harvesting performance of the energy harvesters subjected to stochastic external excitation by means of the stochastic averaging technique or the equivalent linearization technique. However, due to the complexity and the multistable characteristics, it is difficult to use these techniques to study multistable energy harvesters under random excitations, and there is only theoretical, numerical, and experimental investigation for multistable energy harvesters under deterministic excitations by far [27–31]. For optimal vibration energy harvesting, it is necessary and of importance to theoretically reveal the response mechanism of multistable energy harvesters subject to random parametric excitations.The novelty of this paper is to theoretically investigate the resonance mechanism of multistable energy harvesters under narrow-band stochastic parametric excitations for improving vibration energy harvesting performance. In Section2, the theoretical analysis framework is presented. In Section 3, based on the method of multiple scales, the largest Lyapunov exponent which determines the stability of the trivial solutions is derived. The first kind modified Bessel function is utilized to solve the theoretical solutions. In Section 4, the nontrivial solutions and their stability are provided. In Section 5, the stochastic bifurcation is analyzed for revealing the response mechanism of multistable energy harvesters by using the finite difference method. Main conclusions are addressed at last.


## 2. Theoretical Analysis Framework

As introduced above, several different nonlinear vibrational multistable energy harvesters were designed by far. Taking the piezoelectric pentastable energy harvester as an example to demonstrate the structure of multistable energy harvesters, the structural schematic is shown in Figure1. The piezoelectric materials (PZTs) attached at the clamped end of the cantilever beam are used to convert vibration energy into usable electric energy, which may supply power for wireless sensors and small portable electromechanical devices. Because of the nonlinear magnetic force produced by the interaction between the tip magnet and external magnets, the pentastable energy harvester has nine equilibrium positions. 1, 3, 5, 7, and 9 are stable equilibrium positions, and intercalary 2, 4, 6, and 8 are unstable equilibrium positions. Figure 2 describes the equivalent schematic diagram of the general piezoelectric multistable energy harvester which is subjected to the base excitation xbt. In principle, the harvester may exhibit nonlinear monostable, bistable, tristable, or other multistable characteristics which mainly depend on the property of the nonlinear spring.Figure 1
Structural diagram of the piezoelectric vibrational pentastable energy harvester.Figure 2
Equivalent schematic diagram of the general piezoelectric vibrational multistable energy harvester.For the cantilever beam-based nonlinear energy harvesters, the theoretical model can be obtained based on Euler–Bernoulli beam theory, piezoelectric theory, Kirchhoff’s law, and experimental identification [17, 42]. Under the stochastic parametric excitation, the mechanical governing equation and the electrical governing equation of multistable energy harvester shown in Figures 1 and 2 can be described as follows [17, 42]:(1)Mx¨t+Cx˙t+Frt−θ0υt=Fξtxt,(2)Cpυ˙t+υtRl−1+θ0x˙t=0,where xt is the tip displacement of the harvester relative to the base. M, C, and Cp are the equivalent mass, the equivalent damping, and the equivalent capacitance, respectively. Rl and θ0 are the load resistance and the equivalent electromechanical coupling coefficient, respectively. υt is the output voltage. ξt is the stochastic excitation, and F is the density of the stochastic excitation (F>0). Fr is the equivalent nonlinear restoring force (it should be noted that its opposite number −Fr is the real restoring force in physics and experiments), which is the vector sum of the linear restoring force and the nonlinear magnetic force. Its general expression can be expressed as a polynomial, as follows [17, 42]:(3)Fr=m0+m1x+m2x2+⋯+mnxn,where m0, m1, m2, …, mn are the polynomial coefficients. For the pentastable energy harvester, Fr=ω02/x12x22x32x42xx2−x12x2−x22x2−x32x2−x42, where ω0 can be considered as a coefficient. It should be noted that x = 0 is the middle stable equilibrium position of the harvester.The stochastic excitationξt is assumed as an ergodic narrow-band stochastic process with zero mean. It is governed by the following equation advanced by Wedig [43]:(4)ξt=cosΩt+γ^Wt,where Ω is the center frequency, Wt is the standard Wiener process, and γ^ is the bandwidth of the stochastic excitation. The corresponding power spectrum Sξϖ can be written as follows [43]:(5)Sξϖ=12γ^2Ω2+ϖ2+γ^4/4Ω2−ϖ2+γ^4/42+ϖ2γ^4.Under the limiting procedureγ^⟶∞, it results in the uniformly distributed power spectrum of white noise. However, when γ^⟶0, the power spectrum Sξϖ vanishes in the entire frequency range except at the singular frequency ϖ=±Ω, where Sξϖ goes to infinity. The power spectrum in equation (5) is same as the power spectral density of the narrow-band filtered Gaussian white noise. In this paper, we focus on the case that γ^ is small; thus, ξt is bound to a narrow-band random process.Then, for the convenience of the following analysis, the nondimensional governing equations of equations (1) and (2) are provided, as follows:(6)x¨t+cx˙t+ω2xt+a1x3t+a2x5t+a3x7t+a4x9t−θυt=fξtxt,(7)υ˙t+λυt+gx˙t=0,where c=C/M,θ=θ0/M,f=F/M,λ=1/RlCp,g=θ0/Cp, ω2=ω02/M, a1=−ω2x12x22x32+x12x22x42+x12x32x42+x22x32x42/x12x22x32x42, a2=ω2x12x22+x12x32+x12x42+x22x32+x22x42+x32x42/x12x22x32x42, and a3=−ω2x12+x22+x32+x42/x12x22x32x42; a4=ω2/x12x22x32x42.The method of multiple scales [44] is applied to find the second-order analytically approximate solutions of equations (6) and (7). Firstly, the expressions of the solutions can be expressed as(8)xt,ε=x0T0,T1+εx1T0,T1+⋯,(9)υτ,ε=υ0T0,T1+ευ1T0,T1+⋯,where Tn=εnt and ε is a small nondimensional bookkeeping parameter.Denoting the differential operators byD0=∂/∂T0, D1=∂/∂T1, and D2=∂/∂T2 and using the chain rule, the following equations can be obtained:(10)∂∂t=D0+εD1+ε2D2+⋯,(11)∂2∂t2=D02+2εD0D1+ε2D12+2D0D2+⋯.To use the method of multiple scales, the system parameters could be expressed asc=εc^, a1=εa^1, a2=εa^2, a3=εa^3, a4=εa^4, and f=εf^. Then, substituting equations (8)–(11) into equations (6) and (7) and collecting the terms with the identical powers of ε, we have the following:ε0:(12)D02x0T0,T1+ω2x0T0,T1=0,(13)D0υ0T0,T1+λυ0T0,T1=−gD0x0T0,T1.ε1:(14)D02x1T0,T1+ω2x1T0,T1=−2D0D1x0T0,T1−c^D0x0T0,T1−a^1x03T0,T1−a^2x05T0,T1−a^3x07T0,T1−a^4x09T0,T1+θ^υ0T0,T1+f^cosΩT0+γWT1x0T0,T1,where the statistical property of the standard Wiener process γ^Wt=γ^Wεt/ε=γWT1 is used.According to equations (12)–(14), the first-order solutions can be written as follows:(15)x0T0,T1=AT1cosωT0+φT1,(16)υ0T0,T1=C0T1e−λT0+gωAT1λ2+ω2λsinωT0+φT1−ωcosωT0+φT1,where A=AT1 and φ=φT1 are the functions of the slow time scale T1.Denoting the excitation frequency asΩ=2ω+εσ (σ is the detuning parameter), we introduce a new variable ηT1=σT1+γWT1−2φT1. Substituting equations (15) and (16) into equation (14), the following equation can be obtained:(17)D02x1T0,T1+ω2x1T0,T1=2ωA′T1sinωT0+φ+2ωAT1φ′cosωT0+φ+c^AT1ωsinωT0+φ−a^1A3T1cos3ωT0+φ−a^2A5T1cos5ωT0+φ−a^3A7T1cos7ωT0+φ−a^4A9T1cos9ωT0+φ+θ^C1T1e−λT0+gωAT1λ2+ω2λsinωT0+φ−ωcosωT0+φ+f^2AT1cosηcosωT0+φ−f^2AT1sinηsinωT0+φ+f^2AT1cosηcos3ωT0+3φ−f^2AT1sinηsin3ωT0+3φ.Eliminating the secular terms in equation (17), we can derive(18)A′T1=−θ^gλAT12λ2+ω2−c^2AT1+f^AT14ωsinη,(19)AT1η′T1=σAT1−θ^gωAT1λ2+ω2−3a^14ωA3T1−5a^28ωA5T1−35a^364ωA7T1−63a^4128ωA9T1+f^AT12ωcosη+γAT1W′T1.Equations (18) and (19) are the first-order equations which govern the amplitude and the phase of the response displacement and the output voltage of the multistable energy harvester. The first-order uniform expansion of the solutions of equations (6) and (7) is given by(20)x=2AεtcosΩ2t−ηεt2+Oε.


## 3. Trivial Solutions and Their Stability

Obviously, equations (18) and (19) have a trivial solution of A=0, which corresponds to the steady-state response. Hence, the trivial steady-state solutions and their stability should be analyzed firstly. The induced linearization equations of equations (18) and (19) at the point 0,0 can be written as follows:(21)A′T1=−θ^gλA2λ2+ω2−c^2A+f^A4ωsinη,(22)Aη′T1=σA−θ^gωAλ2+ω2+f^A2ωcosη+γAW′T1.Letρ=lnA, and according to Itô’s formula, equations (21) and (22) can be written as the following Itô stochastic differential equations:(23)dρ=−θ^gλ2λ2+ω2−c^2+f^4ωsinηdT1,(24)dη=σ−θ^gωλ2+ω2+f^2ωcosηdT1+γdW.According to equations (23) and (24), the steady-state probability density pη is governed by the following Fokker–Planck–Kolmogorov (FPK) equation:(25)d2pdη2−ddησ¯−f¯cosηp=0,where σ¯=2σ−θ^gω/λ2+ω2/γ2 and f¯=−2f^/2ωγ2.Using both the periodicity condition and the normality condition [45], the solution of equations (18) and (19) can be derived, as follows:(26)pη=expσ¯η+π−f¯sinη4π2Iiσ¯f¯I−iσ¯f¯∫ηη+2πexp−σ¯y+f¯sinydy,where Iiσ¯f¯ and I−iσ¯f¯ are the first kind modified Bessel functions. η is any real or complex number. According to the multiplicative ergodic theorem [46], the Lyapunov exponent of the trivial steady-state solutions of equations (21) and (22) can be written as follows:(27)λ0=limT1⟶∞1T1lnAT1A0=limT1⟶∞1T1ρT1−ρ0=−θ^gλ2λ2+ω2−c^2+f^4ωlimT1⟶∞1T1∫0T1sinηθdθ=−θ^gλ2λ2+ω2−c^2+f^4ωEsinη=−θ^gλ2λ2+ω2−c^2+f^4ω∫02πsinηpηdη=−θ^gλ2λ2+ω2−c^2+f^8ωI1−iσ¯−f¯I−iσ¯−f¯+I1+iσ¯−f¯Iiσ¯−f¯,where E⋅ denotes the mathematical expectation.Since the response mechanism of nonlinear systems are very complex [42, 44, 47], in order to verify the aforementioned theoretical results, the parameters of equations (6) and (7) are chosen as θ^=0.5, g=1.0, λ=1.0, ω=1.0, and c^=0.1. Figure 3 shows the variations of the largest Lyapunov exponent determined by equation (27). Specifically, Figure 3(a) depicts the three-dimensional plot of the largest Lyapunov exponent λ0. Corresponding isohypse curves of λ0 are shown in Figure 3(b). Obviously, there exist two different ranges of the solutions for the largest Lyapunov exponents. Near the resonance area, the largest Lyapunov exponent increases to the maximum value in the center of the instability region. Outside this area, there exists a complete plane where the Lyapunov exponents are zero. To verify the results shown in Figure 3, the time-domain results in Figure 4(a) show that the trivial solution is unstable at point A (f^=1.5 and σ=0.5) in Figure 4(b). However, the results shown in Figure 4(b) indicate a contrary phenomenon that the trivial solution is stable at point B (f^=0.5 and σ=0.5).Figure 3
Largest Lyapunov exponents: (a) mesh surface; (b) isohypse curves.
(a)
(b)Figure 4
Time-domain response displacement of the multistable energy harvester: (a)f^=1.5; (b) f^=0.5.
(a)
(b)


## 4. Nontrivial Solutions and Their Stability

In this section, the nontrivial steady-state solutions of equations (18) and (19) will be derived. To obtain the first-order and the second-order steady-state moments of the displacement amplitude, the Itô stochastic differential equations of equations (18) and (19) can be derived, as follows:(28)dA=−θ^gλA2λ2+ω2−c^2A+f^A4ωsinηdT1,(29)dη=σ−θ^gωAλ2+ω2−3a^14A2−5a^28A4−35a^364A6−63a^4128A8+f^2ωcosηdT1+γdW.For the case ofγ=0, the nontrivial steady-state solutions can be denoted as A=A0 and η=η0. Combining the conditions of A′=0, η′=0, and A≠0, the following relations can be derived based on equations (28) and (29):(30)−θ^gλ2λ2+ω2−c^2+f^4ωsinη0=0,(31)σ−θ^gωA0λ2+ω2−3a^14A02−5a^28A04−35a^364A06−63a^4128A08+f^2ωcosη0=0.Whenγ≠0, the stability of the nontrivial steady-state solutions can be examined by using the perturbation terms, as follows:(32)A=A0+A1,(33)η=η0+η1,where A0 and η0 are determined by equations (30) and (31). A1 and η1 are the perturbation terms.Substituting equations (32) and (33) into equations (28) and (29), neglecting the nonlinear terms, we can obtain the linearization modulation of equations (28) and (29) at A0 and η0, as follows:(34)dA1=f^4ωA0cosη0η1dT1,(35)dη1=−32a^1A0−52a^2A03−10532a^3A05−6316a^4A07A1dT1−f^2ωsinη0η1dT1+γdW.Since equations (34) and (35) are the linear Itô stochastic differential equations, the property can be obtained by combining the moment method [48], as follows:(36)dEA1dT1=dEη1dT1=dEA12dT1=dEη12dT1=0.The first-order and second-order steady-state moments of the nontrivial steady-state solutions of equations (34) and (35) can be derived, as follows:(37)EA1=Eη1=0,EA1η1=0,Eη12=γ2ωf^sinη0,EA12=−σ+θ^gωA0/λ2+ω2+3a^1/4A02+5a^2/8A04+35a^3/64A06+63a^4/128A08−3/2a^1−5/2a^2A02−105/32a^3A04−63/16a^4A06Eη122.Taking the expectation on both sides of equations (32) and (33), the first-order and second-order steady-state moments of the nontrivial steady-state solutions of equations (18) and (19) can be obtained, as follows: (38)EA=EA0+A1=A0,(39)EA2=A02+EA12=A02+−σ+θ^gωA0/λ2+ω2+3a^1/4A02+5a^2/8A04+35a^3/64A06+63a^4/128A08−3/2a^1−5/2a^2A02−105/32a^3A04−63/16a^4A06⋅Eη122.For the steady-state case,T0=t⟶∞ is satisfied, which corresponds the term C0T1e−λT0 in equation (16) satisfied C0T1e−λT0⟶0 with the conditions C0T1 bounded and λ>0; then, the amplitude of the response voltage could be obtained:(40)V=gωλ2+ω2AT1.The first-order and second-order steady-state moments of the voltage amplitude are(41)EV=gωλ2+ω2A0,(42)EV2=g2ω2λ2+ω2EA2.For the case ofx1=1.8, x2=3.3, x3=4.5, and x4=5.4, according to the expressions below equation (7), the corresponding nonlinear coefficients are a1=−0.4841456090, a2=0.06354506813, a3=−0.003049732679, and a4=0.00004799705192. Figures 5(a) and 5(b) indicate the first-order and second-order steady-state moments of the displacement amplitude determined by equations (38) and (39). The corresponding steady-state moments of the voltage amplitude are plotted in Figures 5(c) and 5(d). Due to the dependence on the steady-state moments of the displacement amplitude in expressions (41) and (42), the similar phenomenon could be observed in Figures 5(a)–5(d).Figure 5
Steady-state moments of displacement and voltage of the multistable energy harvester with different high-order stiffness coefficients (γ=0.1). (a, b) First-order and second-order steady-state moments of the displacement amplitude; (c, d) first-order and second-order steady-state moments of the voltage amplitude.
(a)
(b)
(c)
(d)It can be found that the steady-state moments change along with the appearance of the multivalued response branches, which are induced by the existence of the high-order nonlinear terms. As a comparison, another case ofa4=0.0001 (keeping other parameters being same) is also shown in Figure 5. The results showed that the response of the multistable energy harvester is very sensitive to the high-order stiffness coefficients a4 which also have a strong influence on output voltage.The steady-state moments of the voltage amplitude are proportional to the displacement amplitude; thus, in the following part, only the dynamic properties of the displacement amplitude are given, which also reflects the characteristics of the corresponding voltage amplitude.For a different detuning parameterσ, the first-order steady-state moments with the variation in θ are plotted in Figure 6. The detailed response properties are presented in Table 1. When σ is small (the case of σ=−1.51 shown in Figure 6(a)), the response curve is made up of one semicircle, which is the response property I in Table 1. The larger value of σ leads to the other semicircle appearing in the region with the lower amplitude (the cases of σ=−1.507 and σ=−1.5 in Figure 6(a) and σ=−1.4 in Figure 6(b)); thus, the response curve includes two semicircles (response property II). They combine to one continuous curve, and the critical value of σ is approximate to −1.3 (response property III). With the increase in σ, the continuous curve moves to the right side of the transverse axis (the cases of σ=−1.1 and σ=−0.9 shown in Figure 6(b); response property IV). The threshold for appearing this kind of phenomenon is σ≈−0.732 (response property V), which means that the existence of the other two smaller semicircles could be found in the response curve. The response curve could also consist of three branches (the general case of σ=−0.6 shown in Figure 6(b); response property VI). Furthermore, the new two semicircles become larger gradually, and they finally combine to one continuous curve at σ=−0.525 (response property VII). Therefore, the response curve separates into two independent continuous curves, which can be clearly observed from the case of σ=−0.3 shown in Figure 6(c) (response property VIII).Figure 6
Steady-state moments of the multistable energy harvester under different detuning parameters versusθ.
(a)
(b)
(c)
(d)
(e)Table 1
Response characteristics shown in Figure6.
Response characteristicsθ⟶EA
Interval range off
I
[−1.512, −1.507)
II
[−1.507, −1.3)
III (critical case)
σ
≈
−
1
.
3
IV
(−1.3, −0.732)
V (critical case)
σ
≈
−
0
.
732
VI
(−0.732, −0.525)
VII (critical case)
σ
≈
−
0
.
525
VIII
(−0.525, −0.2666)
IX (critical case)
σ
≈
−
0
.
2666
X
(−0.2666, 0.01)
XI
[0.01, 0.063)
XII
[0.063, 0.387)
XIII
[0.387, 0.4522)
XIV (critical case)
σ
≈
0
.
4522
XV
(0.4522, 0.4582)
XVI
σ
≈
0.4582
XVII
(0.4582, 401)Along with the increase inσ, it can be observed that the left continuous curve moves to the right side quickly although the continuous curve on the right side also shifts to this side. At last, they attach to each other at point H (σ≈−0.2666) as response property IX shown in Table 1. Due to the further increase in σ, the curve degenerates the higher branch and the lower branch (the cases of σ≈−0.2 shown in Figure 6(d); response property X). Increasing σ leads to the appearance of the new phenomenon, such as in the case of σ≈0.0 (response property XI). Then, for the larger σ, the lower branch moves to the origin point and disappears gradually, as shown in Figure 6(e) (the cases of σ=0.1 and σ=0.3; response property XII). Furthermore, this phenomenon is observed in the remaining branch, as shown in Figure 6(e). At the beginning, the response curve is shown by response property XIII (the case of σ=0.45). When σ is increased to 0.4522, the independent branch degenerates to one continuous curve and one closed circle (response property XIV). It locates at the bottom of the curve, which is attached to each other at point K. Naturally, increasing the value of σ makes the response curve separate into two parts: the continuous curve and the closed circle (the cases of σ=0.455 and σ=0.458 plotted in Figure 6(e); response property XV). Due to the increase in σ, the closed circle becomes smaller and smaller. Finally, the critical case appears at σ≈0.4582 (response property XVI). Only one part remains as shown in Figure 6(e) for the case of σ=0.48 (response property XVII). It can be found that the electromechanical coupling coefficient θ plays an important role in the steady-state responses of the multistable energy harvester.The effect of the damping coefficient on the first-order steady-state moment is shown in Figure7. It can be found that the different phenomena also appear along with the change in the damping coefficient under different detuning parameters. When the detuning parameter σ is small, the response curve is one independent semicircle, as shown in Figure 7(a) (σ=−0.7 and σ=−0.5). The peak point O at the right side turns inward, which is shown in Figure 7(a) (σ=−0.35). Furthermore, the lower branch from the original case with a single solution degenerates to the curve with two solutions for the case of σ=−0.268, which is tangent to the vertical axis at point P. This part shifts to the right side of the lateral axis, while the peak value in the middle of the response curve moves to the left side of the lateral axis (σ=−0.2 at point Q). Finally, the peak value in the middle part attaches with the vertical axis, which appears at σ=−0.156 (point R), as shown in Figure 7(b).Figure 7
Steady-state moments of the multistable energy harvester under different detuning parameters versus the damping coefficient.
(a)
(b)
(c)Then, the response curve separates into two parts (σ=0.1 in Figure 7(b)). The higher branch turns inward, while the lower branch moves to the left side gradually, and the critical value can be found in Figure 7(b) for the case of σ=0.355. Therefore, the lower branch finally disappears in Figure 7(c) for the critical case of σ≈0.356. For σ≈0.473, the remaining higher branch turns inward to the vertical axis (tangent to the vertical axis at point S) and the response curve degenerates into two parts. The semicircle in the lower branch becomes smaller and smaller for the case of σ=0.495 in Figure 7(c). Finally, only the higher branch remains for the cases of σ=0.6, σ=0.8, and σ=1.5 in Figure 7(c).As a conclusion, the damping coefficient is found to be one unavoidable factor to induce the complex dynamic properties of the multistable energy harvester.In order to discuss the influence mechanism of stochastic noise, the bifurcation diagrams of the output voltage of the energy harvester with the variation of the excitation frequencyΩ are plotted in Figure 8 under the system parameters g=1.0, λ=1.0, ω0=1.0, θ=0.5, c=0.1, f=2.0, σ=0.5, γ=0.1, x1=1.8, x2=3.3, x3=4.5, and x4=5.4. In Figure 8(b), the noise intensity is γ=0.1, while the deterministic case is plotted in Figure 8(a). As can be seen that the existence of the stochastic noise has the obvious influence on the output voltage, the regular bifurcation behavior disappears; thus, the random factor could not be ignored in the energy harvesters. Two typical examples are shown in Figure 9, and the existence of the stochastic noise could induce the instability of the output voltage, which leads to the difference in the spectrum of the output voltage.Figure 8
Bifurcation diagram of the output voltage with the variation inΩ: (a) γ=0.0; (b) γ=0.1.
(a)
(b)Figure 9
Response of the energy harvester forΩ=1.8: (a) time-domain output voltage for γ=0.0; (b) spectrum of the time-domain output voltage for γ=0.0; (c) time-domain output voltage for γ=0.1; (d) spectrum of the time-domain output voltage for γ=0.1.
(a)
(b)
(c)
(d)Then, in Figure10, with the variation in excitation f, the bifurcation diagram of the output voltage with c=0.5, Ω=3.0, and other parameters is the same with Figure 8. The fluctuation of the voltage is obvious in Figure 10(b), and the existence of the noise results in the change of the bifurcation characteristics of the energy harvesters. One typical example is shown in Figure 11, and it can be seen that the response curve becomes disordered in Figure 11(c), which is very different from the case without the existence of noise. The corresponding spectrum of output voltage shown in Figures 11(b) and 11(d) also makes a big difference.Figure 10
Bifurcation diagram of the output voltage with the variation inf: (a) γ=0.0; (b) γ=0.1.
(a)
(b)Figure 11
Response of the energy harvester forf=2.22: (a) time-domain output voltage for γ=0.0; (b) spectrum of the time-domain output voltage for γ=0.0; (c) time-domain output voltage for γ=0.1; (d) spectrum of the time-domain output voltage for γ=0.1.
(a)
(b)
(c)
(d)Finally, in order to consider the influence of the number of potential wells, the restoring force of the tristable and monostable energy harvesters are degraded intoFr=ω02/x12x22xx2−x12x2−x22 and Fr=ω02x, and the bifurcation diagrams of the tristable and monostable energy harvesters are shown in Figure 12, which are different from the results shown in Figure 8. The monostable energy harvesters keep the steady-state behavior; thus, the existence of noise does not lead to the fluctuation of the voltage.Figure 12
Bifurcation diagram of the output voltage with the variation inΩ: (a) γ=0.0 for the tristable system; (b) γ=0.1 for the tristable system; (c) γ=0.0 for the monostable system; (d) γ=0.1 for the monostable system.
(a)
(b)
(c)
(d)


## 5. Stochastic Bifurcation

For understanding the stochastic bifurcation phenomenon which is induced by stochastic excitations, the finite difference method [32] is used in this section. Apparently, equations (28) and (29) are two-dimensional Itô differential equations when γ≠0. The statistics of the responses of the multistable energy harvester can be obtained by solving the following FPK equation of equations (28) and (29):(43)∂p∂T1=−∂∂Au1A,ηp−∂∂ηu2A,ηp+γ^22∂2∂η2p,where pA,η,T1 is the transition probability density.The two functions in equation (40) are defined as follows:(44)u1A,η=−θ^gλA2λ2+ω2−c^2A+f^A4ωsinη,u2A,η=σ−θ^gωAλ2+ω2−3a^14A2−5a^28A4−35a^364A6−63a^4128A8+f^2ωcosη.The boundary condition is the natural boundary condition, while the initial distribution is assumed to be Gaussian with a probability density, as follows [49]:(45)pA0,η0=12πs1s2exp−A0−μ122s12−η0−μ222s22,where μ1=Aup+Abelow/2, μ2=ηup+ηbelow/2, s1=0.1, and s2=0.2.In the following analysis, the system parameters are set asθ^=0.5, g=1.0, λ=1.0, ω=1.0, σ=0.05, c^=0.1, and f^=0.5. The values of other nonlinear coefficients are same with those in Section 4. The numerical simulation is obtained in a reduced state space 0.0,2.0×−5.0,10.0. It is verified that the joint probability density can be stationary when T≥15 with the time step size of Δt=0.01. The integral time can be defined as T=20. In addition, the upwind scheme, which is accurate enough to discrete, is adopted to deal with the convective term, while the central difference scheme is applied on the diffusion term. Figure 13 shows a series of the joint probability densities of the amplitude and the phase obtained from the FPK equation solved by using the finite difference method with different bandwidths of γ.Figure 13
Joint probability of the multistable energy harvester: (a)γ=0.1; (b) γ=0.2; (c) γ=0.3; (d) γ=0.4.
(a)
(b)
(c)
(d)It is found that the joint probability density concentrates at the nontrivial steady-state solution branch whenγ=0.1, as shown in Figure 13(a). The joint probability density apparently jumps downward to the trivial steady-state solution branch for γ=0.2, as shown in Figure 13(b). Due to the increase in γ, the joint probability density continues to jump downward to the trivial steady-state solution branch for γ=0.3 in Figure 13(c) and γ=0.4 in Figure 13(d). The results in Figure 13 indicate that, for a lower noise intensity γ, the most probable motion of the multistable energy harvester is around the higher branch (nontrivial branch) of the amplitude response curve. However, the most probable motion of the multistable energy harvester is around the lower branch (trivial branch) of the amplitude response curve for a higher noise intensity γ. Therefore, the change in γ will induce the phenomenon of stochastic bifurcation. The direction of the jump moves from the higher branch to the lower branch along with the increase in γ.In order to deeply investigate the probability distribution, three different energy harvesters are studied and the corresponding potential energy functions (Ux=∫−x0x0Frdx, where −x0 and x0 are the integral boundaries) are depicted in Figure 14. In detail, the equivalent nonlinear restoring forces of the tristable energy harvester and the monostable energy harvester are Fr=ω02/x12x22xx2−x12x2−x22 and Fr=ω02x, respectively. The contour surfaces for the three different energy harvesters are compared. It is found that the contour surface of the tristable energy harvester in Figures 15(c) and 15(d) is very similar to that of the pentastable energy harvester in Figures 15(a) and 15(b). For the monostable energy harvester, u2 degenerates into u2A,η=σ−θ^gωA/λ2+ω2+f^/2ωcosη, as shown in Figures 15(e) and 15(f). Its response characteristic is very different with the other two harvesters. In other words, the number of stable equilibrium positions leads to variation in the distribution of the trivial and nontrivial steady-state solutions. Therefore, the appearance of the stochastic bifurcation can be partly determined by the number of stable equilibrium points. This conclusion can be referred for the optimization design of energy harvesters for different external excitation conditions.Figure 14
Potential energy functions: (a) pentastable energy harvester; (b) tristable energy harvester; (c) monostable energy harvester.
(a)
(b)
(c)Figure 15
Contour surfaces: (a, b) pentastable energy harvester; (c, d) tristable energy harvester; (e, f) monostable energy harvester.
(a)
(b)
(c)
(d)
(e)
(f)


## 6. Conclusions

This paper theoretically investigates the resonance mechanism of nonlinear vibrational multistable energy harvesters under narrow-band stochastic parametric excitation. The largest Lyapunov exponent of the trivial steady-state solutions of multistable energy harvesters is deduced to analyze its stability, and the first-order and second-order steady-state moments of the nontrivial steady-state solutions are derived. It is found that the electromechanical coupling coefficient and damping coefficient greatly influence the dynamic characteristic of the steady-state moments. The stochastic bifurcation phenomenon between the nontrivial trivial steady-state solutions is revealed based on the FPK equation corresponding to the two-dimensional Itô stochastic differential equations. According to the numerical simulation results, the stationary joint probability density concentrates at the nontrivial steady-state solution branch when the intensity of the stochastic excitation is small. The probability of the trivial steady-state solution becomes larger along with the increase in the intensity of the stochastic excitation. In addition, the number of stable equilibrium points is found to influence the distribution of the trivial and nontrivial steady-state solutions. This influences the appearance of the stochastic bifurcation phenomenon and further influences the energy harvesting performance. Overall, the resonance mechanism of multistable energy harvesters under narrow-band stochastic parametric excitations has been revealed. In the future, the research framework in this paper can be referred for the optimization design of nonlinear vibrational multistable energy harvesters for different external excitation conditions.



---

*Source: 1050143-2019-12-14.xml*