Approximate Fokker–Planck–Kolmo-gorov equation analysis for asymmetric multistable energy harvesters excited by white noise

The typical structure for multistable piezoelectric energy harvesters is illustrated in figure 1. It consists of a cantilever beam, two piezoelectric layers pasted at the root, and tip magnets coupled with two rotational external magnets. The multistable energy harvesters can be achieved by appropriately adjusting the parameters d, h, and ϑ, as well as the polarity of magnets. Due to the difficulties in modulation, it is impossible to attain a multistable harvester with a completely symmetric potential function.

Zoom In Zoom Out Reset image size

By applying the Hamilton principle and Euler–Bernoulli beam theory, the dynamic model of the multistable energy harvesters is expressed as follows:

$$\begin{equation}\begin{cases}m{r}^{{\prime\prime}}\left(\bar{t}\right)+c{r}^{\prime }\left(\bar{t}\right)+\bar{f}\hspace{-2.5pt}\left(r\right)-\theta \bar{v}\left(\bar{t}\right)=F\left(\bar{t}\right)\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \left(1\mathrm{a}\right)\\ {C}_{p}{\bar{v}}^{\prime }\left(\bar{t}\right)+\frac{\bar{v}\left(\bar{t}\right)}{R}+\theta {r}^{\prime }\left(\bar{t}\right)=0\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \left(1\mathrm{b}\right)\end{cases}\end{equation}$$

in which m, c, $\bar{f}$, θ, C_p and R are the equivalent mass, equivalent damping, nonlinear restoring force, equivalent electromechanical coupling coefficient, equivalent capacitance, and the load resistance respectively. Also, r, $\bar{v}$ and F are the tip displacement of the cantilever, the voltage output, and the external excitation respectively. Further, ${\left(\;\right)}^{\prime }$ represents the derivate with respect to time $\bar{t}$. By introducing the variables expressed by $x=\frac{r}{{l}_{0}}$, $t=\bar{t}{\omega }_{n}$ and $v=\frac{{C}_{p}\bar{v}}{\theta {l}_{0}}$, the nondimensional model of the system is

$$\begin{equation}\begin{cases}\ddot {x}+2\xi \dot {x}+\bar{u}\left(x\right)-{\kappa }^{2}v=W\left(t\right)\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \enspace \left(2\mathrm{a}\right)\\ \dot {v}+\alpha v+\dot {x}=0\quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \enspace \left(2\mathrm{b}\right)\end{cases}\end{equation}$$

where l₀ is a length scale, and ω_n is the resonance frequency. In this equation, x, v, and W respectively represent the dimensionless displacement, voltage, and the Gaussian white noise excitation, ξ is the damping ratio, α is the time constant ratio, κ² is dimensionless electromechanical coupling coefficient, and $\bar{u}\left(x\right)$ is the dimensionless nonlinear restoring force. Also, the upper dot is the derivate to dimensionless time t. Furthermore, W(t) is the Gaussian white noise.

To characterize the dynamic behavior of multistable harvester under Gaussian white noise, the FPK equation governing equation (2) will be approximately solved by the method of detailed balance. Based on equation (2), the solution in global could be assumed as [55, 70]

$$\begin{equation}X=x\left(t\right)=A\enspace \mathrm{cos}\enspace \phi ,\qquad Y=\dot {x}\left(t\right)=-A{\omega }_{0}\enspace \mathrm{sin}\enspace \phi ,\qquad \phi \left(t\right)={\omega }_{0}t+\varphi .\end{equation} \tag{ 3 }$$

By combining the assumed solution with equation (2b), the response of voltage in global can be expressed by

$$\begin{equation}v\hspace{-2.5pt}\left(t\right)\approx -\left[\frac{\alpha }{{\alpha }^{2}+{\omega }_{0}^{2}}\dot {x}+\frac{{\omega }_{0}^{2}}{{\alpha }^{2}+{\omega }_{0}^{2}}x\right].\end{equation} \tag{ 4 }$$

Substituting equation (4) to equation (2a), one can get the following approximate equation characterizing the model:

$$\begin{equation}\ddot {x}+{\Lambda}\dot {x}+u\left(x\right)=W\left(t\right)\end{equation} \tag{ 5 }$$

where ${\Lambda}=2\xi +\frac{{\kappa }^{2}\alpha }{{\alpha }^{2}+{\omega }_{0}^{2}}$ and $u\left(x\right)=\bar{u}\left(x\right)+\frac{{\kappa }^{2}{\omega }_{0}^{2}}{{\alpha }^{2}+{\omega }_{0}^{2}}x$. For equation (5), it can be converted into the following Ito stochastic differential equation

$$\begin{equation}\begin{gathered}\mathrm{d}{x}_{1}={x}_{2}\enspace \mathrm{d}t\\ \mathrm{d}{x}_{2}=-\left[{\Lambda}{x}_{2}+u\left({x}_{1}\right)\right]\mathrm{d}t+\sqrt{2D}\enspace \mathrm{d}B\left(t\right)\end{gathered}\end{equation} \tag{ 6 }$$

where $B\left(t\right)$ is a Brownian motion process and D is the noise intensity. The corresponding FPK equation of equation (6) is expressed as

$$\begin{equation}\frac{\partial p}{\partial t}=-\frac{\partial }{\partial {x}_{1}}\left({x}_{2}p\right)+\frac{\partial }{\partial {x}_{2}}\left[\left(\mathbf{\Lambda }{x}_{2}+u\left({x}_{1}\right)\right)p\right]+D\frac{{\partial }^{2}}{\partial {{x}_{2}}^{2}}p\end{equation} \tag{ 7 }$$

where p is the stationary probability density function with the boundary conditions $p\left(-\infty ,t\right)=p\left(\infty ,t\right)=0$. The stationary solution of equation (7) can be supposed as

$$\begin{equation}p\left(\boldsymbol{x}\right)=C\enspace \mathrm{exp}\left[-\phi \left(\boldsymbol{x}\right)\right].\end{equation} \tag{ 8 }$$

Based on the principle of detailed balance, the first and second-order derivative moment of equation (6) are

$$\begin{equation}{a}_{1}={x}_{2},\qquad {a}_{2}=-{\Lambda}{x}_{2}-u\left({x}_{1}\right),\qquad {b}_{11}={b}_{12}={b}_{21}=0,\qquad {b}_{22}=2D\end{equation} \tag{ 9 }$$

in which the reversible and irreversible part of the first-order derivative moment are respectively

$$\begin{equation}{{a}_{1}}^{R}={x}_{2},\qquad {{a}_{1}}^{I}=0,\qquad {{a}_{2}}^{R}=-u\left({x}_{1}\right),\qquad {{a}_{2}}^{I}=-{\Lambda}{x}_{2}.\end{equation} \tag{ 10 }$$

On account of detailed balance, the following equations explaining $\phi \left(\boldsymbol{x}\right)$ are obtained

$$\begin{equation}D\frac{\partial \phi }{\partial {x}_{2}}={\Lambda}{x}_{2}\end{equation} \tag{ 11 }$$

$$\begin{equation}{x}_{2}\frac{\partial \phi }{\partial {x}_{1}}=u\left({x}_{1}\right)\frac{\partial \phi }{\partial {x}_{2}}.\end{equation} \tag{ 12 }$$

By solving equations (11)–(12), the expression of $\phi \left(\boldsymbol{x}\right)$ is

$$\begin{equation}\phi \left({x}_{1},{x}_{2}\right)=\frac{{\Lambda}}{D}\left[{\int }_{0}^{{x}_{1}}u\left(z\right)\mathrm{d}z+\frac{1}{2}{{x}_{2}}^{2}\right].\end{equation} \tag{ 13 }$$

Thus, the stationary probability density function is achieved as

$$\begin{equation}p\left({x}_{1},{x}_{2}\right)=C\enspace \mathrm{exp}\left\{-\frac{{\Lambda}}{D}\left[{\int }_{0}^{{x}_{1}}u\left(z\right)\mathrm{d}z+\frac{1}{2}{{x}_{2}}^{2}\right]\right\}\end{equation} \tag{ 14 }$$

where ${C}^{-1}={\int }_{-\infty }^{\infty }\mathrm{exp}\left\{-\frac{{\Lambda}}{D}\left[{\int }_{0}^{{x}_{1}}u\left(z\right)\mathrm{d}z+\frac{1}{2}{{x}_{2}}^{2}\right]\right\}\mathrm{d}{x}_{1}\enspace \mathrm{d}{x}_{2}$. Particularly, equation (14) could be expressed as follows

$$\begin{equation}p\left(x,\dot {x}\right)=p\left(x\right)\cdot p\left(\dot {x}\right)={A}_{1}\cdot \enspace \mathrm{exp}\left[-\frac{{\Lambda}}{D}{\int }_{0}^{x}u\left(z\right)\mathrm{d}z\right]\cdot {A}_{2}\cdot \enspace \mathrm{exp}\left(-\frac{{\Lambda}}{D}\frac{1}{2}{\dot {x}}^{2}\right)\end{equation} \tag{ 15 }$$

where ${{A}_{1}}^{-1}={\int }_{-\infty }^{\infty }\mathrm{exp}\left[-\frac{{\Lambda}}{D}{\int }_{0}^{x}u\left(z\right)\mathrm{d}z\right]\mathrm{d}x$ and ${{A}_{2}}^{-1}={\int }_{-\infty }^{\infty }\mathrm{exp}\left(-\frac{{\Lambda}}{D}\frac{1}{2}{\dot {x}}^{2}\right)\mathrm{d}\dot {x}$.

Although the voltage response has a relationship with the displacement response, it depends on the displacement in local potential wells. In order to characterize the power output of the system, the voltage response in local can be expressed as

$$\begin{equation}{v}_{i}^{{\ast}}\left(t\right)\approx -\left[\frac{\alpha }{{\alpha }^{2}+{\omega }_{0}^{2}}\dot {x}+\frac{{\omega }_{0}^{2}}{{\alpha }^{2}+{\omega }_{0}^{2}}\left(x-{x}_{i}^{{\ast}}\right)\right]\end{equation} \tag{ 16 }$$

where i represents the ith response condition and ${x}_{i}^{{\ast}}$ is the mean value of displacement response in the ith condition. For monostable harvesters, i equals to 1. While i is 3 for the bistable system representing two conditions of intrawell motion and one condition for inter-well oscillation. In a similar fashion, i is 4 for tristable energy harvesting systems. Then the mean square value of voltage is characterized by the following expression:

$$\begin{equation}E\left[{v}^{2}\right]=\sum\limits _{i=1}^{N}{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }{v}_{i}^{{\ast}}{\left(t\right)}^{2}p\left(x,\dot {x}\right)\mathrm{d}x\enspace \mathrm{d}\dot {x}\end{equation} \tag{ 17 }$$

In the investigation into the probability distribution of linear harvester, the dimensionless nonlinear restoring force is $\bar{u}\left(x\right)=x$. For the numerical simulation, the sampling frequency is set to be 5 Hz and the data length is 2 × 10⁶. Figures 2(a) and (b) illustrate the numerical and theoretical marginal PDF of displacement and velocity response of the linear harvester under white noise. Besides, the comparison of numerical and theoretical joint PDF of displacement and velocity is shown in figure 3. Results indicate that the response of the linear harvester mainly distributes around the original point, and numerical simulations are in good agreement with the theoretical calculation.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To study the average output power of the linear harvester under white noise, figure 4 compares the numerical and theoretical power under excitations with various intensities. It is seen that the power increase with an increase in the noise intensity of the excitation. Furthermore, figure 5 illustrates the variation of output power with other parameters such as dimensionless electromechanical coupling coefficient, time constant ratio, and damping ratio. Results demonstrate that the power decreases with an increase in the values of electromechanical coupling coefficient and damping ratio. While the power firstly increases and then decreases after a certain value of the time constant ratio, indicating that there is an optimum time constant ratio to maximize the output performance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In the investigation into monostable harvesters, the potential function expressed by $\bar{u}\left(x\right)=x+{x}^{3}$ is firstly considered for the symmetric harvester. Under excitation with an intensity of 0.04, the marginal and joint PDF of displacement and velocity are depicted in figures 6 and 7. It is viewed that numerical results match well with the theoretical calculation. Compared with the linear harvester, there is a larger probability for the monostable one to oscillate near the stable point and this is due to the change in the potential energy function. While for the response probability distribution of velocity, it is not affected by the shape of potential energy functions. For the average output power shown in figure 10(a), numerical simulation agrees well with the theoretical analysis and the variation trend is similar to that of the linear harvester. The only point to be noted is that the introduction of cubic coefficients to the system decreases the output performance.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For an asymmetric monostable harvester, a quadratic coefficient in the nonlinear restoring force is considered with a value of 1.9. In this condition, the marginal PDF of displacement and velocity are shown in figure 8, while the joint distribution is figured in figure 9. Numerical and theoretical results all demonstrate that the probability distribution of velocity is the same as before discussed, that is to say, the potential function does not affect the probability distribution of velocity. While for the probability distribution of displacement, it indicates that the system has a larger probability to oscillate on the side with a smaller nonlinear restoring force. In figure 10(b), the relationship between average power and noise intensity for the asymmetric MEH is presented. With the appearance of asymmetry in a monostable system, the power output increases and the error between numerical and theoretical results shows an increasing tendency.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In terms of nonlinear bistable harvesters, the nonlinear restoring force is considered with the expression of $\bar{u}\left(x\right)=-x+\beta {x}^{2}+{x}^{3}$. For symmetric bistable harvester with β = 0, the comparison between numerical and theoretical marginal and joint PDF of displacement and velocity are explained in figures 11 and 12. From the results, it is seen that the distribution of displacement exhibits obvious bimodal characteristics because of the two potential wells in the bistable system and the oscillator has a larger probability to vibrate in the two potential wells. While for the distribution of velocity, it is not affected by the bistable potential.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When an asymmetric potential is introduced with β = 0.15, the potential energy function has a deeper potential well on the left and a shallower on the right. In this case, the probability distribution is illustrated in figures 13 and 14. Again, the characteristic of velocity response is the same as before discussed and not influenced by potential shape. For the response of displacement, the distribution has the characteristic of asymmetric bimodal and a larger probability is obtained in the deeper well. As the quadratic coefficient further increases to 0.25 and 0.35, the theoretical joint PDF are shown in figure 15. Results demonstrate that the increase in the asymmetric degree increases the probability to oscillate in the deeper potential well, while decrease that in the shallower one. In figures 16(a) and (b), the relationship between average power and noise intensity for the symmetric and asymmetric (β = 0.35) BEH are respectively presented. The output power increases with an increase in the noise intensity and numerical results are in good agreement with the theoretical ones. By careful observation, it is viewed that the power output of the asymmetric (β = 0.35) BEH is a little lower than the symmetric BEH, and less interwell motion results in a small error between numerical and theoretical results.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Regarding the tristable harvester, the nonlinear restoring force is expressed by the polynomial form of $\bar{u}\left(x\right)=\lambda x+\beta {x}^{2}+\delta {x}^{3}+\gamma {x}^{4}+\eta {x}^{5}$. For the reason that the dynamic behavior of the tristable system is very complex, the sampling frequency here is set to be 5 Hz and the data length is 2 × 10⁷. When λ = 1, β = 0, δ = −2.5, γ = 0 and η = 1 for symmetric harvester, the potential energy function is shown in figure 17 with a blue line. Under Gaussian white noise excitation, the PDFs are shown in figures 18 and 19 for marginal and joint ones. Similar to the results before, the probability distribution of velocity is not affected. For the displacement response characteristics, the probability density function is symmetric with respect to x = 0, and a smaller probability is indicated in the shallower potential well in the middle.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When the quartic coefficient in nonlinear restoring force equals 0.10, 0.15, and 0.20, the potential energy functions are shown in figure 17(b). With an increase in the value of the quartic coefficient, the potential well on the left becomes much deeper while the right one becomes much shallower, and the potential well depth in the middle is not affected. Under white noise excitation with an intensity of 0.04, the marginal and joint probability distribution of the tristable system with γ = 0.10 is demonstrated in figures 20 and 21. The system has a larger probability to oscillate in the deeper potential well and a smaller one in the shallower potential well. This phenomenon is further illustrated by the joint probability distribution of the tristable system with γ equaling to 0.15 and 0.20 shown in figure 22. Figure 23 shows the variation of power output with noise intensity for the symmetric and asymmetric (γ = 0.20) TEH. Although the theoretical result predicts the numerical simulation in a certain degree, the error between numerical and theoretical results is larger than that for the monostable and bistable systems. The reason for the phenomenon is that equations (16) and (17) approximately viewed the interwell motion as periodic oscillation. While under random excitations, the interwell motion is always chaotic oscillation. Therefore, approximately viewing the interwell motion as periodic oscillation increases the error between numerical and theoretical results especially for tristable harvesters that possess large-amplitude periodic oscillation crossing the potential wells.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size