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.

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 [15]. 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 [49]. 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 [912].

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, 1921]. 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 [2426].

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 [3239], 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 [3440]. 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 [2731]. 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 Section 2, 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 Figure 1. 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 . In principle, the harvester may exhibit nonlinear monostable, bistable, tristable, or other multistable characteristics which mainly depend on the property of the nonlinear spring.

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]:where is the tip displacement of the harvester relative to the base. M, C, and are the equivalent mass, the equivalent damping, and the equivalent capacitance, respectively. and are the load resistance and the equivalent electromechanical coupling coefficient, respectively. is the output voltage. is the stochastic excitation, and is the density of the stochastic excitation (). is the equivalent nonlinear restoring force (it should be noted that its opposite number 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]:where , , , …, are the polynomial coefficients. For the pentastable energy harvester, , where 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 is assumed as an ergodic narrow-band stochastic process with zero mean. It is governed by the following equation advanced by Wedig [43]:where is the center frequency, is the standard Wiener process, and is the bandwidth of the stochastic excitation. The corresponding power spectrum can be written as follows [43]:

Under the limiting procedure , it results in the uniformly distributed power spectrum of white noise. However, when , the power spectrum vanishes in the entire frequency range except at the singular frequency , where 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, 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:where , , , , and ; .

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 aswhere and is a small nondimensional bookkeeping parameter.

Denoting the differential operators by , , and and using the chain rule, the following equations can be obtained:

To use the method of multiple scales, the system parameters could be expressed as , , , , , and . Then, substituting equations (8)–(11) into equations (6) and (7) and collecting the terms with the identical powers of , we have the following:ε0:ε1:where the statistical property of the standard Wiener process is used.

According to equations (12)–(14), the first-order solutions can be written as follows:where and are the functions of the slow time scale .

Denoting the excitation frequency as ( is the detuning parameter), we introduce a new variable . Substituting equations (15) and (16) into equation (14), the following equation can be obtained:

Eliminating the secular terms in equation (17), we can derive

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

3. Trivial Solutions and Their Stability

Obviously, equations (18) and (19) have a trivial solution of , 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 can be written as follows:

Let , and according to Itô’s formula, equations (21) and (22) can be written as the following Itô stochastic differential equations:

According to equations (23) and (24), the steady-state probability density is governed by the following Fokker–Planck–Kolmogorov (FPK) equation:where and .

Using both the periodicity condition and the normality condition [45], the solution of equations (18) and (19) can be derived, as follows:where and 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:where 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 , , , , and . 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 . Corresponding isohypse curves of 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 ( and ) 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 ( and ).

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:

For the case of , the nontrivial steady-state solutions can be denoted as and . Combining the conditions of , , and , the following relations can be derived based on equations (28) and (29):

When , the stability of the nontrivial steady-state solutions can be examined by using the perturbation terms, as follows:where and are determined by equations (30) and (31). and 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 and , as follows:

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:

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:

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:

For the steady-state case, is satisfied, which corresponds the term in equation (16) satisfied with the conditions bounded and ; then, the amplitude of the response voltage could be obtained:

The first-order and second-order steady-state moments of the voltage amplitude are

For the case of , , , and , according to the expressions below equation (7), the corresponding nonlinear coefficients are , , , and . 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).

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 of (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 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 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 and in Figure 6(a) and 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 (response property III). With the increase in , the continuous curve moves to the right side of the transverse axis (the cases of and shown in Figure 6(b); response property IV). The threshold for appearing this kind of phenomenon is (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 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 (response property VII). Therefore, the response curve separates into two independent continuous curves, which can be clearly observed from the case of shown in Figure 6(c) (response property VIII).

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 () 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 shown in Figure 6(d); response property X). Increasing leads to the appearance of the new phenomenon, such as in the case of (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 and ; 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 ). 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 and 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 (response property XVI). Only one part remains as shown in Figure 6(e) for the case of (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 Figure 7. 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) ( and ). The peak point O at the right side turns inward, which is shown in Figure 7(a) (). Furthermore, the lower branch from the original case with a single solution degenerates to the curve with two solutions for the case of , 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 ( at point Q). Finally, the peak value in the middle part attaches with the vertical axis, which appears at (point R), as shown in Figure 7(b).

Then, the response curve separates into two parts ( 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 . Therefore, the lower branch finally disappears in Figure 7(c) for the critical case of . For , 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 in Figure 7(c). Finally, only the higher branch remains for the cases of , , and 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 , , , , , , , , , , , and . In Figure 8(b), the noise intensity is , 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.

Then, in Figure 10, with the variation in excitation , the bifurcation diagram of the output voltage with , , 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.

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 into and , 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.

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 . The statistics of the responses of the multistable energy harvester can be obtained by solving the following FPK equation of equations (28) and (29):where is the transition probability density.

The two functions in equation (40) are defined as follows:

The boundary condition is the natural boundary condition, while the initial distribution is assumed to be Gaussian with a probability density, as follows [49]:where , , , and .

In the following analysis, the system parameters are set as , , , , , , and . The values of other nonlinear coefficients are same with those in Section 4. The numerical simulation is obtained in a reduced state space . It is verified that the joint probability density can be stationary when with the time step size of . The integral time can be defined as . 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 .

It is found that the joint probability density concentrates at the nontrivial steady-state solution branch when , as shown in Figure 13(a). The joint probability density apparently jumps downward to the trivial steady-state solution branch for , 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 in Figure 13(c) and 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 (, where and 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 and , 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, degenerates into , 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.

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.

Data Availability

All data were obtained from the methods and the equations used in this paper. Therefore, any person can obtain the data and repeat the study once the paper is published.

Conflicts of Interest

The authors declare that there are no conflicts of interest concerning the publication of this article.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant nos. 11702201 and 11802237), the China Postdoctoral Science Foundation (Grant no. 2018M641012), the Fundamental Research Funds for the Central Universities (Grant no. G2018KY0306), and the 111 Project (no. BP0719007).