1 Introduction

Vibration energy harvesters (VEHs) excited by random excitations have flourished as a focal research topic in recent years. Such sources of powers are currently finding applications in different fields of technology, such as wireless sensors, data transmitters, vivo biomedical implants, and health monitoring of structures and machines^[1-5].

Many studies have been devoted to exploiting different random vibratory energy harvesting systems. Cottone et al.^[6] found that the bistable oscillators can outperform the linear ones under the Gaussian white noise excitation. Daqaq^[7] derived an approximate expression for the mean power under the exponentially correlated noise and demonstrated the existence of an optimal potential shape maximizing the output power. Green et al.^[8] demonstrated that Duffing-type nonlinearities can reduce the size of electromagnetic energy harvesting devices. Daqaq^[9] used the method of moment differential equations to calculate response statistics. Masana and Daqaq^[10] investigated the influence of stiffness-type nonlinearities on the transduction of the buckling piezoelectric beam under the band-limited noise. He and Daqaq^[11-12] employed the statistical linearization technique and the finite element method to investigate how the shape of the potential energy function influences the mean steady-state approximate output power. Xu et al.^[13] introduced the stochastic averaging of energy envelope for Duffing-type VEHs. Kumar et al.^[14] used the finite element method to solve the Fokker-Planck-Kolmogorov (FPK) equation of the associated bistable energy harvester. Jin et al.^[15] introduced the equivalent nonlinearization technique to derive a semi-analytical solution of the corresponding nonlinear VEHs. Jiang and Chen^[16-17] used the standard stochastic averaging method and the generalized stochastic averaging method to solve the response of nonlinear random vibratory energy harvesting systems. Liu et al. introduced VEHs with fractional-order nonlinear properties^[18] and colored noise excitation^[19], respectively. Xiao and Jin^[20] proposed the energy harvester under the correlated white noise excitation, and computed the response with the generalized harmonic function method. Yang and Xu^[21] exploited the VEHs with fractional derivative damping, and analyzed the dynamical response with the stochastic averaging method. As seen from the above description, the response statistical characteristics of random energy harvesting systems can be solved from the stationary FPK equation. However, to solve the probability density function (PDF) from the transition FPK equation has not been reported. To address the lack of studies in this aspect, the present work develops the path integration technique to solve the transient FPK equation.

The PDF evolution is governed by the FPK equation. Thus, the FPK equation provides a powerful tool for addressing the response characteristics of the nonlinear stochastic system. However, the FPK equation cannot be exactly solved. Many efforts have been devoted to developing different approximate or numerical approaches. Among various approaches^[22-24], the path integration method has been receiving much attention, which can be traced back to Wehner and Wolfer^[25], and also developed by Crandall et al.^[26], Hsu and Chiu^[27-28], Sun and Hsu^[29], Naess and Johnson^[30], Iourtchenkoa et al.^[31] and Feng et al.^[32]. Recently, a Gauss-Legendre integration scheme to improve the computational efficiency and accuracy of the path integration method was developed ^[33-34]. Subsequently, this approach was extended to the Duffing-Rayleigh oscillator^[35], structural dynamic systems^[36], and ship roll motion^[37]. So far, to the authors' best knowledge, there is no path integration analysis on energy harvesting. To address the lack of studies in this aspect, the present work develops the path integration technique to determine the response of nonlinear energy harvesters under the Gaussian white noise excitation.

The paper is organized as follows. Section 2 reviews the governing equation of motion for vibratory energy harvesting systems. Section 3 establishes the PDF solution of vibratory energy harvesting systems under Gaussian noises excitations with the path integration method. Sections 4-6 investigate the stationary PDF solution of three different types of nonlinear VEHs. Section 7 ends the paper with concluding remarks.

2 Vibration energy harvesting system

The electromechanical coupling equations of motion for vibration energy harvesting systems can be expressed as

(1)

(2)

where x represents the displacement of the mass, c is the damping coefficient, χ is the linear electromechanical coupling coefficient, _b is the base acceleration, C_p is the capacitance, L is the inductance, U(x) represents the potential function symmetries and asymmetries, and z is the electric quantity representing the induced voltage in piezoelectric harvesters and the induced current in the electromagnetic ones. These are measured across an equivalent resistive load R.

The non-dimensional electromechanical coupling equations can be described as^[8-9]

(3)

(4)

where , and ξ(t)=-_b is the base acceleration. ξ(t) is the Gaussian white noise with a zero mean and autocorrelation function

(5)

in which E[·] denotes the expected value, D is the intensity of the excitation, and δ(·) is the Dirac function.

It is worth noting that Eqs. (3)-(4) are universal since g(X) represents the restoring force symmetries and asymmetries which include cubic nonlinearity, quadratic-cubic nonlinearity, and cubic-quintic nonlinearity.

The transition FPK equation corresponding to Eqs. (3) and (4) is governed as

(6)

where x=[x, y, z]=[X, , Z], f₁ = y, f₂ = -(2ζy + g(x) + κz), and f₃=-αz+y.

3 Path integration method for energy harvesting systems

For the PDF p(x, y, z, t) at any given moment, it can be calculated through the integration from the transition probability density q(x, y, z, t) and the initial distribution as

(7)

where Ω is the integral interval of x, y, and z.

For a long-time duration into a sequence of appropriate short-time intervals from the initial distribution as follows:

(8)

Next, the Gauss-Legendre quadrature integration scheme is used to calculate Eq. (8). First, the domain along the x-direction is divided into two Gauss points as

(9)

along the y-direction it can be discretized into

(10)

and along the z-direction it can be written as

(11)

In this case, the integral (8) can be discretized into the following composite Gauss-Legendre quadrature form:

(12)

where Δ_x, Δ_y, and Δ_z are the lengths of sub-intervals along the x-, y-, and z-directions, respectively.

The short time transition probability density is

(13)

where , and x^{I - 1} = [x_i^I-1, y_j^I-1, z_k^I-1].

Moreover, an initial distribution is selected as

(14)

where m_x, m_y, and m_z are means, and σ_x, σ_y, and σ_z are standard deviations.

4 Application to the Duffing-type VEH

The Duffing-type VEH is an important model of the nonlinear VEH, which had been widely reported^[1]. Moreover, the non-dimensional electromechanical coupling equations can be described as

(15)

(16)

The transition probability in Eq. (13) is assumed to be Gaussian, and depends only on the mean and variances of x, y, and z. The mean and variances for the system can be derived from the non-dimensional electromechanical coupling equations (15) and (16), which yields

(17)

where m_ijk = E[xⁱy^jz^k].

These moment equations are not closed because the system is nonlinear. The moments of higher than the second-order can be truncated by using the Gaussian closure procedure,

(18)

Substitute Eq. (18) into Eq. (17) and obtain the closed moment equations, i.e.,

(19)

Then, Eq. (19) can be solved by the fourth-order Runge-Kutta (RK4) algorithm. By substituting the solution of Eq. (19) into Eq. (13), the path integration PDFs can been calculated by using the initial distribution (14) and Eq. (12). The numerical implementation of the path integration method has been described above, and the main steps for the numerical implementation are concluded and given in Fig. 1.

The PDFs and logarithmic PDFs of the displacement, velocity, and electricity obtained with the path integration method are compared with the Monte Carlo simulation (MCS) results of the systems (15) and (16) shown in Fig. 2. The system parameters are set as ζ=0.8, δ=0.5, κ=0.05, α=0.05, and D=0.01.

Figure 2 shows that the PDFs and logarithmic PDFs of displacement, velocity and electricity obtained with the path integration method are close to the MCS results. The logarithmic PDFs which show the effectiveness of the solution procedure can illustrate the accuracy of the tail region. The validity and accuracy of the path integration method have been verified by the MCS.

Figure 3 describes the joint PDFs of the displacement and velocity obtained with the path integration method and the MCS results of the systems (15) and (16). It is demonstrated that the joint PDF which is based on the path integration method agrees very well with the simulation result.

5 Application to the quadratic-cubic VEH

The quadratic-cubic VEH is an asymmetric system which includes quadratic and cubic nonlinearity in the restoring force^[8]. One of the mathematical models for the quadratic-cubic VEH can be expressed as follows:

(20)

(21)

The mean and variances for the system can be derived from the non-dimensional electromechanical coupling equations (20) and (21), which yields

(22)

where m_ijk = E[xⁱy^jz^k].

Substitute Eq. (18) into Eq. (22) and obtain closed moment equations by the RK4 algorithm. The system parameters are selected as ζ=0.2, δ=0.5, κ=0.5, α=0.5, λ=1, and D=0.01. Therefore, the path integration PDFs can be rendered by using the initial distribution (14) and Eq. (12).

Figure 4 depicts the PDFs and logarithmic PDFs of the displacement, velocity, and electricity for the quadratic-cubic VEH, and Fig. 5 shows the joint PDFs of displacement and velocity obtained with the path integration method and the MCS results of the systems (20) and (21). These results show that the approximate PDFs and logarithmic PDFs of the path integration method coincide with the MCS ones.

6 Application to cubic-quintic VEH

A cubic-quintic Duffing-type VEH was designed by Cao et al.^[38], and they discussed the influence of the potential well depth on the harvesting performance under the harmonic excitation. Jiang and Chen^[17] investigated the stochastic response of the cubic-quintic Duffing-type VEH under the Gaussian white noise excitation by using the stochastic averaging method. The non-dimensional electromechanical coupling equations can be written as

(23)

(24)

The transition probability in Eq. (13) is assumed to be Gaussian, and depends only on the mean and variances of x, y, and z. The mean and variances for the system can be derived from the non-dimensional electromechanical coupling equations (23) and (24), which yields

(25)

where m_ijk = E[xⁱy^jz^k].

These moment equations are not closed because the system is nonlinear. The moments of higher than the second-order can be truncated by using the Gaussian closure procedure. By substituting Eq. (18) into Eq. (25) and obtaining closed moment equations, the path integration PDFs can be calculated with the closed moment. The system parameters are selected as ζ=0.4, δ=5, κ=0.05, α=0.05, β=1, and D=0.01.

For the cubic-quintic VEH, the PDFs and logarithmic PDFs of the displacement, velocity, and electricity obtained with the path integration method are compared with the MCS results in Fig. 6, and the joint PDFs of displacement and velocity are presented in Fig. 7. The PDFs and logarithmic PDFs obtained from the path integration agree well with the MCS ones.

7 Conclusions

In this paper, we investigate the response of the vibration-based energy harvester under the Gaussian white noise. To achieve this goal, the path integration solution of the transition FPK equation of energy harvesting systems is used to obtain the response statistics of the harvester. The path integration process is determined by using the Gauss-Legendre integration scheme, and the short-time transition PDF is introduced by the short-time Gaussian approximation. The procedure is applied to three different types of nonlinear VEHs under Gaussian white excitations. It is observed that, both of the path integration solution and the direct Monte Carlo numerical integration of the stochastic differential equations provide close predictions for the marginal PDFs of the displacement, velocity, and electricity, and the joint PDFs of the displacement and velocity.

The investigation yields the following conclusions.

(Ⅰ) The information of how probability distribution evolves over time is formulated numerically by the values of transition PDF at the discrete Gaussian quadrature points.

(Ⅱ) The PDF from the transition FPK equation of energy harvesting systems is obtained.

(Ⅲ) The path integration outcomes are qualitatively and quantitatively supported by the MCS.