Nanoscale piezoelectric vibration energy harvester design

Foruzande, Hamid Reza; Hajnayeb, Ali; Yaghootian, Amin

doi:10.1063/1.4994577

Development of new nanoscale devices has increased the demand for new types of small-scale energy resources such as ambient vibrations energy harvesters. Among the vibration energy harvesters, piezoelectric energy harvesters (PEHs) can be easily miniaturized and fabricated in micro and nano scales. This change in the dimensions of a PEH leads to a change in its governing equations of motion, and consequently, the predicted harvested energy comparing to a macroscale PEH. In this research, effects of small scale dimensions on the nonlinear vibration and harvested voltage of a nanoscale PEH is studied. The PEH is modeled as a cantilever piezoelectric bimorph nanobeam with a tip mass, using the Euler-Bernoulli beam theory in conjunction with Hamilton’s principle. A harmonic base excitation is applied as a model of the ambient vibrations. The nonlocal elasticity theory is used to consider the size effects in the developed model. The derived equations of motion are discretized using the assumed-modes method and solved using the method of multiple scales. Sensitivity analysis for the effect of different parameters of the system in addition to size effects is conducted. The results show the significance of nonlocal elasticity theory in the prediction of system dynamic nonlinear behavior. It is also observed that neglecting the size effects results in lower estimates of the PEH vibration amplitudes. The results pave the way for designing new nanoscale sensors in addition to PEHs.

I. INTRODUCTION

Several energy harvesters for powering nano-scale portable electronic systems have been developed as the demand for self-powered electronic devices grows rapidly.^1,2 One of the most common types of energy harvesters is piezoelectric harvesters. Piezoelectric energy harvesters have been studied and applied in nanoscale extensively because of their simpler structures, scalability and higher energy density in small scales.^3,4 Wang⁴ reviewed the theoretical and experimental aspects of previous research on piezoelectric nanowires. Chen et al.⁵ used piezoelectric nanofibers for energy harvesting. Wang and Wang⁶ studied a nanoscale unimorph PEH. Wang and Wang⁷ and Deng et al.⁸ studied the flexoelectric effects in nanoscale PEHs. After determining the piezoelectric nano-generators as one of the best type of energy harvester in nanoscale, the next step is modeling.

Expensive fabrication processes in nanoscale makes it more necessary to prepare an accurate model during the designing process of a nanoscale device. Because of similarity between the structure and elements of a PEH in nano and macro scales, the basis of modeling process in nanoscale is the same as macro scale. Thereby, before modeling a nanoscale PEH, a review on the techniques of modeling macroscale PEHs is necessary.

The modeling of macro scale PEHs has been investigated extensively in previous studies.^9–13 A typical PEH is comprised of a cantilevered piezoelectric beam with a tip mass, where the clamped end of the beam is attached to the host vibrating structure. Erturk and Inman¹⁴ modeled a unimorph PEH as a continuous system using Euler-Bernoulli beam theory, and compared its vibrations with the simple single-degree-of-freedom (SDOF) models of previous research works. They showed that the SDOF model may predict highly inaccurate vibrations comparing to the beam model. Although a PEH with linear behavior is easily modeled and analyzed, but they have serious shortcomings such as limited performance bandwidth. High values of output voltage in linear PEHs are limited to a narrow band of frequency close to the resonance of the PEH. The simplest solution to low output voltage in any vibration frequency is increasing the amplitude of vibrations. It does not only boost the output voltage of a PEH in all frequencies, but also broadens its bandwidth as a result of nonlinear behavior. Therefore, nonlinearity and size which is also an important issue in modeling, are discussed in the following lines.

There are several studies on nonlinear vibrations of PEHs.^15–17 Sebald et al.¹⁶ experimentally studied the harvested power and bandwidth of a PEH with nonlinearities in stiffness and damping and compared them with the performance of a linear PEH. The results showed significant improvements in the bandwidth and harvested power. The type of stability also affects the performance of a PEH. Cottone et al.¹⁵ theoretically and experimentally studied the nonlinear vibrations of a bistable oscillator as a PEH. The PEH consisted of an axially loaded piezoelectric beam. They observed a superior performance of the proposed PEH under a wideband random excitation. In many studies, an external nonlinear force like the force of a permanent magnet is applied to induce a nonlinearity and use its benefits in the PEH. Tang et al.¹⁷ used magnetic force to investigate influence of both monostable and bistable configurations in a PEH under random excitations. They proved the superiority of these configurations over linear PEHs. In summary, the nonlinearity boosts the performance of a PEH and is desired in energy harvesting.

In addition to nonlinearity, another important factor in modeling of PEHs is the size of the structure. Molecular dynamics simulations and experiments show that the mechanical behaviors of micro and nanoscale structures are significantly affected by size effects.¹⁸ Therefore, the classical continuum theory as a scale-independent theory is unable to precisely model small structures. Recent studies suggest using size-dependent theories (e.g. the modified couple stress theory, the strain gradient theory and etc.).¹⁹ The nonlocal elasticity theory introduced by Eringen²⁰ is one of the size-dependent theories that has been used successfully in the majority of previous research. The theory is also applicable to piezoelectric materials and structures. It has been used to study the vibrations of the nanobeam in a limited number of research works. Ke and Wang²¹ applied the nonlocal elasticity theory to study the linear vibrations of a piezoelectric nanobeam. Differential quadrature (DQ) method was used to obtain the natural frequencies and mode shapes of the nanobeam with different boundary conditions. It was found that the applied continuum theory has an important role in finding more realistic values for the natural frequencies. Ke et al.²² investigated the nonlinear vibrations of a nanobeam by using nonlocal elasticity theory. The DQ method and a direct iterative method were used to obtain the nonlinear natural frequencies of the nanobeam. Because piezoelectric beams in small scales are mostly fabricated in unimorph or bimorph forms, modeling bimorph and unimorph nanobeams in addition to applying appropriate continuum theory is necessary. Nazemizadeh and Bakhtiari-nejad²³ analytically investigated the free vibrations of nano/microbeams using the nonlocal elasticity theory. The effects of nonlocal parameter in addition to different parameters of the system on the natural frequencies and mode shapes of the beam were studied. In conclusion, the resulting error from neglecting the size effects in modeling a nanoscale PEH is significant.

To the best knowledge of authors no study has been published on the dynamic modeling of nanoscale vibratory PEHs. They are used in powering small-scale electronic device and have to be modeled before fabrication. Large amplitude nonlinear vibrations of a nanoscale PEH is investigated in this paper. Nanoscale PEH is considered to be a piezoelectric bimorph cantilever nanobeam with a tip mass. The nanobeam is subjected to harmonic base excitation. It is modeled using Euler-Bernoulli beam theory combined with nonlocal elasticity theory to address the size effects in nanoscale. The nonlinear equations of motion and related boundary conditions are derived employing Hamilton’s principle. Then, the obtained equations are discretized using assumed modes method and solved utilizing the multiple scales method. To validate the obtained approximate solutions, they are compared with a numeric solution, which shows good agreement. The results show the behavior of the nanoscale PEH for different values of system parameters.

II. THEORETICAL MODELING

A bimorph piezoelectric nanobeam of length L and width b has been considered (figure 1). The thickness of the piezoelectric layers is assumed to be h_p while h_e indicates the thickness of the elastic layer. Also, polarity direction of the piezoelectric layers is considered parallel to the z axis. In order to consider direct and inverse effects of piezoelectric material concurrently, electric potential is assumed as a summation of triogonometric and linear terms, according to Wang’s research:²⁴

\begin{matrix} Φ^{(1)} (x, z, t) = - cos (β z_{1}) ϕ (x, t) + \frac{z_{1} V_{0}}{h_{p}} exp (i Ω t), \\ Φ^{(3)} (x, z, t) = - cos (β z_{2}) ϕ (x, t) + \frac{z_{2} V_{0}}{h_{p}} exp (i Ω t), \end{matrix}

(1)

where $Φ^{1} (x, z, t)$ and $Φ^{3} (x, z, t)$ are total electric potential functions of the upper and lower piezoelectric layers, respectively. The superscripts determines the layer number of the beam. Moreover, a numerical value of $β = \frac{π}{h_{p}}$ is known, and ϕ(x, t) indicates the generated voltage due to piezoelectric direct effect. $V_{0}$ is the amplitude of an external voltage and Ω is the frequency of this voltage. The local coordinates of z₁ and z₂ are shown in figure 1. Because there is no external voltage in PEHs, the second term in the right-hand-side of both Equations (1a) and (1b) is not considered in the rest of this study. z₁ and z₂ can be written as the following functions of z:

\begin{matrix} z_{1} = - z - \frac{1}{2} (h_{e} + h_{p}), \\ z_{2} = - z + \frac{1}{2} (h_{e} + h_{p}) . \end{matrix}

(2)

By substituting Equations (2) into Equations (1), the total electric potential functions can be rearranged as:

\begin{matrix} Φ^{(1)} (x, z, t) = - cos [β (- z - \frac{1}{2} (h_{e} + h_{p}))] ϕ (x, t), \\ Φ^{(3)} (x, z, t) = - cos [β (- z + \frac{1}{2} (h_{e} + h_{p}))] ϕ (x, t) . \end{matrix}

(3)

After determining the total electric potential functions, the electric filed components can be obtained as:

\begin{matrix} E_{x}^{(1)} = - \frac{\partial Φ^{(1)}}{\partial x} = cos [β (- z - \frac{1}{2} (h_{e} + h_{p}))] \frac{\partial ϕ}{\partial x}, \\ E_{z}^{(1)} = - \frac{\partial Φ^{(1)}}{\partial z_{1}} = - \frac{\partial Φ^{(1)}}{\partial z} \frac{\partial z}{\partial z_{1}} = - β sin [β (- z - \frac{1}{2} (h_{e} + h_{p}))] ϕ, \\ E_{x}^{(3)} = - \frac{\partial Φ^{(3)}}{\partial x} = cos [β (- z + \frac{1}{2} (h_{e} + h_{p}))] \frac{\partial ϕ}{\partial x}, \\ E_{z}^{(3)} = - \frac{\partial Φ^{(3)}}{\partial z_{2}} = - \frac{\partial Φ^{(3)}}{\partial z} \frac{\partial z}{\partial z_{2}} = - β s i n [β (- z + \frac{1}{2} (h_{e} + h_{p}))] ϕ . \end{matrix}

(4)

Based on Euler-Bernoulli beam theory, the only nonzero component of the strain tensor of an element of a beam is given as:

ε_{x x} {= ε}_{x x}^{°} - z K,

(5)

where $ε_{x x}^{°}$ is strain at neutral axis and K is beam curvature. $ε_{x x}^{°}$ is neglected because of applying shortening assumption. Hence, the strain component of Equation (5) can be rewritten as:²⁵

ε_{x x} = - z K = - z (W^{″} - W^{'} U^{″} - W^{″} U^{'} - W^{″} {W^{'}}^{2}),

(6)

where U(x, t) and W(x, t) are the axial and transverse displacements of the nanobeam at point x and time t, respectively.

III. NONLOCAL ELASTICITY THEORY

Based on Nonlocal elasticity theory, the constitutive equations of piezoelectric materials in absence of body forces are expressed as:²¹

\begin{matrix} σ_{i j} - {(e_{0} a)}^{2} \nabla^{2} σ_{i j} = {C_{i j k l} ε}_{k l} - e_{k i j} E_{k}, \\ D_{i} - {(e_{0} a)}^{2} \nabla^{2} D_{i} = {e_{i k l} ε}_{k l} + 𝜖_{i k} E_{k}, \end{matrix}

(7)

where $σ_{i j}$ ⁠, $ε_{k l}$ ⁠, D_i and E_k are stress, strain, electric displacement and electric filed components, respectively. Also, C_ijkl, e_kij and $𝜖_{i k}$ are elastic, piezoelectric and dielectric constants, respectively. Since the thickness to the length ratio and the width to the length ration of the studied beam are assumed to be small, Equation (7) can be rewritten in a reduced form as:

σ_{x x} - {(e_{0} a)}^{2} \frac{\partial^{2} σ_{x x}}{\partial x^{2}} = {C_{11} ε}_{x x} - e_{31} E_{z},

(8)

D_{x} - {(e_{0} a)}^{2} \frac{\partial^{2} D_{x}}{\partial x^{2}} = 𝜖_{11} E_{x},

(9)

D_{z} - {(e_{0} a)}^{2} \frac{\partial^{2} D_{z}}{\partial z^{2}} = {e_{31} ε}_{x x} + 𝜖_{33} E_{z} .

(10)

IV. DERIVING GOVERNING EQUATIONS

The strain energy of a bimorph piezoelectric nanobeam can be expressed as:

\begin{matrix} Π_{s} = & \frac{1}{2} \int_{0}^{b} \int_{0}^{l} \int_{- \frac{h_{e}}{2} - h_{p}}^{- \frac{h_{e}}{2}} (σ_{x x}^{(1)} ε_{x x} - D_{x}^{(1)} E_{x}^{(1)} - D_{z}^{(1)} E_{z}^{(1)}) d z d x d y \\ + \frac{1}{2} \int_{0}^{b} \int_{0}^{l} \int_{\frac{h_{e}}{2}}^{\frac{h_{e}}{2} + h_{p}} (σ_{x x}^{(3)} ε_{x x} - D_{x}^{(3)} E_{x}^{(3)} - D_{z}^{(3)} E_{z}^{(3)}) d z d x d y \\ + \frac{1}{2} \int_{0}^{b} \int_{0}^{l} \int_{- \frac{h_{e}}{2}}^{\frac{h_{e}}{2}} (σ_{x x}^{(2)} ε_{x x}) d z d x d y . \end{matrix}

(11)

By defining the bending moments as below:

\begin{matrix} M_{1} = \int_{- \frac{h_{e}}{2} - h_{p}}^{- \frac{h_{e}}{2}} b (σ_{x x}^{(1)} z) d z, \\ M_{2} = \int_{- \frac{h_{e}}{2}}^{\frac{h_{e}}{2}} b (σ_{x x}^{(2)} z) d z, \\ M_{3} = \int_{\frac{h_{e}}{2}}^{\frac{h_{e}}{2} + h_{p}} b (σ_{x x}^{(3)} z) d z, \end{matrix}

(12)

variation of the strain energy can be given as:

\begin{matrix} δ Π_{s} = & \int_{0}^{l} (- M_{1} δ k - M_{2} δ k - M_{3} δ k) d x \\ + \int_{0}^{b} \int_{0}^{l} \int_{- \frac{h_{e}}{2} - h_{p}}^{- \frac{h_{e}}{2}} (- D_{x}^{(1)} δ E_{x}^{(1)} + D_{z}^{(1)} δ E_{z}^{(1)}) d z d x d y \\ + \int_{0}^{b} \int_{0}^{l} \int_{\frac{h_{e}}{2}}^{\frac{h_{e}}{2} + h_{p}} (- D_{x}^{(3)} δ E_{x}^{(3)} + D_{z}^{(3)} δ E_{z}^{(3)}) d z d x d y . \end{matrix}

(13)

The kinetic energy of the nanobeam is defined as:

\begin{matrix} Π_{k} = \frac{1}{2} \int_{0}^{l} (m_{1} + m_{2} + m_{3} + \frac{m_{tip}}{l_{2}} H (x - l_{2})) [{(\frac{\partial W}{\partial t})}^{2} + {(\frac{\partial U}{\partial t})}^{2}] d x, \\ m_{1} = m_{3} = ρ_{p} b h_{p}, m_{2} = ρ_{e} b h_{e}, \end{matrix}

(14)

where ρ_p, ρ_e and m_i are density of piezoelectric layers, density of elastic layer and mass per unit length of the i-th layer of the nanobeam, respectively. Also, m_tip and l₂ are the quantity and location of the tip mass, respectively. In order to add a tip mass, Heaviside function dented by H(x) is used. Therefore, the variation of kinetic energy is:

δ Π_{k} = \frac{1}{2} \int_{0}^{l} (m_{1} + m_{2} + m_{3} + \frac{m_{tip}}{l_{2}} H (x - l_{2})) δ [{(\frac{\partial W}{\partial t})}^{2} + {(\frac{\partial U}{\partial t})}^{2}] d x .

(15)

The work done by external forces i.e. viscous damping force is expressed as:

Π_{f} = \int_{0}^{l} C \frac{\partial W}{\partial t} W d x,

(16)

Where C is the viscous damping coefficient.The work done by the external forces is given as:

δ Π_{f} = \int_{0}^{l} C \frac{\partial W}{\partial t} δ W d x .

(17)

To derive equations of motion and related boundary conditions, Hamilton’s principle is employed. Hamilton’s principle can be expressed as:

δ H = \int_{0}^{t} (δ Π_{s} + δ Π_{f} - δ Π_{k}) d t - \frac{1}{2} \int_{0}^{t} [\int_{0}^{l} λ δ (1 - {(1 + U^{'})}^{2} - {W^{'}}^{2})] d x = 0 .

(18)

The last term in Equation (18) is used to satisfy the inextensionality condition of the beam. λ is a Lagrange multiplier and would be calculated later. By substituting Equations (13), (15) and (16) in Equation (18) and then equating coefficients of δU, δW and δϕ by zero, the equations of motion are obtained:

\begin{matrix} \begin{matrix} {(λ (1 + U^{'}))}^{'} = {(M_{1} W^{'})}^{''} + {(M_{2} W^{'})}^{''} + {(M_{3} W^{'})}^{''} - {(M_{1} W^{″})}^{'} - {(M_{2} W^{″})}^{'} - {(M_{3} W^{″})}^{'} \\ + (m_{1} + m_{2} + m_{3} + \frac{m_{tip}}{l_{2}} H (x - l_{2})) Ü, \end{matrix} \end{matrix}

(19)

\begin{matrix} \begin{matrix} - {M_{1}}^{″} - {M_{2}}^{″} - {M_{3}}^{″} - {(M_{1} U^{″})}^{'} - {(M_{2} U^{″})}^{'} - {(M_{3} U^{″})}^{'} + {(M_{1} U^{'})}^{''} + {(M_{2} U^{'})}^{''} + {(M_{3} U^{'})}^{''} \\ - 2 {(W^{″} W^{'} M_{1})}^{'} - 2 {(W^{″} W^{'} M_{2})}^{'} - 2 {(W^{″} W^{'} M_{3})}^{'} + {(M_{1} {W^{'}}^{2})}^{''} \\ + {(M_{2} {W^{'}}^{2})}^{''} + {(M_{3} {W^{'}}^{2})}^{''} + (m_{1} + m_{2} + m_{3} + \frac{m_{tip}}{l_{2}} H (x - l_{2})) Ẅ \\ + C \frac{\partial W}{\partial t} - {(λ W^{'})}^{'} = 0, \end{matrix} \end{matrix}

(20)

\begin{matrix} b \int_{- \frac{h_{e}}{2} - h_{p}}^{- \frac{h_{e}}{2}} [D_{z}^{(1)} β sin [β (- z - \frac{1}{2} (h_{e} + h_{p}))] + \frac{\partial D_{x}^{(1)}}{\partial x} cos [β (- z - \frac{1}{2} (h_{e} + h_{p}))]] d z \\ + b \int_{\frac{h_{e}}{2}}^{\frac{h_{e}}{2} + h_{p}} [D_{z}^{(3)} β sin [β (- z + \frac{1}{2} (h_{e} + h_{p}))]) \\ (+ \frac{\partial D_{x}^{(3)}}{\partial x} cos [β (- z + \frac{1}{2} (h_{e} + h_{p}))]] d z = 0 . \end{matrix}

(21)

Also boundary conditions at x = 0 and x = l are given as:

\begin{matrix} (M_{1} + M_{2} + M_{3}) W^{″} - {[(M_{1} + M_{2} + M_{3}) W^{'}]}^{'} + λ (1 + U^{'}) = 0, o r U = 0, \\ (M_{1} + M_{2} + M_{3}) W^{'} = 0, o r U = 0 \\ (M_{1}^{'} + M_{2}^{'} + M_{3}^{'}) + (M_{1} + M_{2} + M_{3}) U^{″} - {[(M_{1} + M_{2} + M_{3}) U^{'}]}^{'} \\ + 2 (M_{1} + M_{2} + M_{3}) W^{″} W^{'} - {[(M_{1} + M_{2} + M_{3}) {W^{'}}^{2}]}^{'} + λ W^{'} = 0, \\ o r W = 0, \end{matrix}

- (M_{1} + M_{2} + M_{3}) + (M_{1} + M_{2} + M_{3}) U^{'} + (M_{1} + M_{2} + M_{3}) W^{' 2} = 0, o r W^{'} = 0,

\begin{matrix} \int_{\frac{h_{e}}{2}}^{\frac{h_{e}}{2} + h_{p}} D_{x}^{(3)} cos [β (- z + \frac{1}{2} (h_{e} + h_{p}))] d z + \int_{- \frac{h_{e}}{2} - h_{p}}^{- \frac{h_{e}}{2}} D_{x}^{(1)} cos [β (- z - \frac{1}{2} (h_{e} + h_{p}))] d z, \\ o r ϕ = 0 . \end{matrix}

(22)

It can be proved that the strain at neutral axis of the nanobeam is obtained as:²⁵

ε_{x x}^{°} = \sqrt{{(1 + U^{'})}^{2} + {W^{'}}^{2}} - 1 .

(23)

Using Taylor expansion, Equation (23) can be approximated as:

ε_{x x}^{°} = U^{'} + \frac{1}{2} W^{' 2} .

(24)

Shortening effect reduces the number of independent motion parameters from three to two by relating the axial displacement to the transverse displacement as the following equation shows:²⁶

ε_{x x}^{°} = 0 \Rightarrow U^{'} = - \frac{1}{2} W^{' 2},

(25)

and also:

U^{″} = - W^{'} W^{″}, Ü = - \frac{1}{2} \frac{\partial^{2}}{\partial t^{2}} (\int_{0}^{x} {W^{'}}^{2} d x) .

(26)

λ is calculated by introducing Equations (25) and (26) in Equation (19):

\begin{matrix} λ = \frac{1}{(1 - \frac{1}{2} {w^{'}}^{2})} [) {(M_{1} w^{'})}^{'} + {(M_{2} w^{'})}^{'} + {(M_{3} w^{'})}^{'} - M_{1} w^{″} - M_{2} w^{″} - M_{3} w^{″} \\ - \frac{1}{2} \int_{l}^{x} (m_{1} + m_{2} + m_{3} + \frac{m_{tip}}{l_{2}} H (x - l_{2})) \frac{\partial^{2}}{\partial t^{2}} (\int_{0}^{x} {w^{'}}^{2} d x) d x d x]) . \end{matrix}

(27)

The denominator of Equation (27) can be substituted by its Taylor expansion as below:

\frac{1}{(1 - \frac{1}{2} {W^{'}}^{2})} \approx (1 + \frac{1}{2} {W^{'}}^{2} + \frac{1}{4} {W^{'}}^{4}) .

(28)

Because of keeping the nonlinear terms up to the third order in this study, only the first term of the expansion is kept and the other two terms are neglected. By substituting Equation (28) into Equation (27), the Lagrange multiplier can be obtained as:

\begin{matrix} λ = ({M^{'}}_{1} + {M^{'}}_{2} + {M^{'}}_{3}) w^{'} \\ - \frac{1}{2} \int_{l}^{x} (m_{1} + m_{2} + m_{3} + \frac{m_{tip}}{l_{2}} H (x - l_{2})) \frac{\partial^{2}}{\partial t^{2}} (\int_{0}^{x} {w^{'}}^{2} d x) d x d x . \end{matrix}

(29)

By submitting Equations (25), (26) and (29) in Equation (20), Equation (30) is obtained:

\begin{matrix} - ({M^{″}}_{1} + {M^{″}}_{2} + {M^{″}}_{3}) - \frac{1}{2} ({M^{″}}_{1} + {M^{″}}_{2} + {M^{″}}_{3}) {W^{'}}^{2} - ({M^{'}}_{1} + {M^{'}}_{2} + {M^{'}}_{3}) W^{'} W^{″} + C \dot{W} \\ + (m_{1} + m_{2} + m_{3} + \frac{m_{tip}}{l_{2}} H (x - l_{2})) Ẅ \\ + {[W^{'} \int_{l}^{x} (m_{1} + m_{2} + m_{3} + \frac{m_{tip}}{l_{2}} H (x - l_{2})) \int_{0}^{x} ({\dot{W^{'}}}^{2} + W^{'} Ẅ^{'}) d x d x]}^{'} = 0 . \end{matrix}

(30)

In order to achieve an explicit expression for the summation of second derivatives of the bending moments, ${M^{″}}_{1} + {M^{″}}_{2} + {M^{″}}_{3}$ ⁠, higher order terms of Equation (30) are ignored:

\begin{matrix} {M^{″}}_{1} + {M^{″}}_{2} + {M^{″}}_{3} \\ = C \dot{W} + (m_{1} + m_{2} + m_{3} + \frac{m_{tip}}{l_{2}} H (x - l_{2})) Ẅ \\ + {[W^{'} \int_{l}^{x} (m_{1} + m_{2} + m_{3} + \frac{m_{tip}}{l_{2}} H (x - l_{2})) \int_{0}^{x} ({\dot{W^{'}}}^{2} + W^{'} Ẅ^{'}) d x d x]}^{'} . \end{matrix}

(31)

On other hand, by multiplying Equation (8) by z and then integrating on the thickness of each layer one obtains the following equations:

M_{1} - {(e_{0} a)}^{2} {M^{″}}_{1} = F_{31}^{(1)} ϕ - D_{11}^{(1)} (W^{″} + \frac{1}{2} W^{″} {W^{'}}^{2}),

(32-a)

M_{2} - {(e_{0} a)}^{2} {M^{″}}_{2} = - D_{11}^{(2)} (W^{″} + \frac{1}{2} W^{″} {W^{'}}^{2}),

(32-b)

M_{3} - {(e_{0} a)}^{2} {M^{″}}_{3} = F_{31}^{(3)} ϕ - D_{11}^{(3)} (W^{″} + \frac{1}{2} W^{″} {W^{'}}^{2}),

(32-c)

where

\begin{matrix} F_{31}^{(1)} & = & \int_{- \frac{h_{e}}{2} - h_{p}}^{- \frac{h_{e}}{2}} b e_{31} β z sin [β (- z - \frac{1}{2} (h_{e} + h_{p}))] d z, \\ F_{31}^{(3)} & = & \int_{\frac{h_{e}}{2}}^{\frac{h_{e}}{2} + h_{p}} b e_{31} β z sin [β (- z + \frac{1}{2} (h_{e} + h_{p}))] d z, \end{matrix}

(33)

\begin{matrix} D_{11}^{(1)} = \int_{- \frac{h_{e}}{2} - h_{p}}^{- \frac{h_{e}}{2}} b C_{11}^{p} z^{2} d z, \\ D_{11}^{(2)} = \int_{- \frac{h_{e}}{2}}^{\frac{h_{e}}{2}} b C_{11}^{e} z^{2} d z, \\ D_{11}^{(3)} = \int_{\frac{h_{e}}{2}}^{\frac{h_{e}}{2} + h_{p}} b C_{11}^{p} z^{2} d z . \end{matrix}

(34)

From Equation (9), the following equations can be obtained:

\begin{matrix} b \int_{- \frac{h_{e}}{2} - h_{p}}^{- \frac{h_{e}}{2}} (D_{x}^{(1)} - {(e_{0} a)}^{2} \frac{\partial^{2} D_{x}^{(1)}}{\partial x^{2}}) cos [β (- z - \frac{1}{2} (h_{e} + h_{p}))] = X_{11}^{(1)} \frac{\partial ϕ}{\partial x}, \\ b \int_{\frac{h_{e}}{2}}^{\frac{h_{e}}{2} + h_{p}} (D_{x}^{(3)} - {(e_{0} a)}^{2} \frac{\partial^{2} D_{x}^{(3)}}{\partial x^{2}}) cos [β (- z + \frac{1}{2} (h_{e} + h_{p}))] = X_{11}^{(3)} \frac{\partial ϕ}{\partial x}, \end{matrix}

(35)

where

\begin{matrix} X_{11}^{(1)} = \int_{- \frac{h_{e}}{2} - h_{p}}^{- \frac{h_{e}}{2}} b \in_{11} {\{cos [β (- z - \frac{1}{2} (h_{e} + h_{p}))]\}}^{2} d z, \\ X_{11}^{(3)} = \int_{\frac{h_{e}}{2}}^{\frac{h_{e}}{2} + h_{p}} b \in_{11} {\{cos [β (- z + \frac{1}{2} (h_{e} + h_{p}))]\}}^{2} d z . \end{matrix}

(36)

Similarly, following equations can be obtained using Equation (10):

\begin{matrix} b \int_{- \frac{h_{e}}{2} - h_{p}}^{- \frac{h_{e}}{2}} [D_{z}^{(1)} - {(e_{0} a)}^{2} \frac{\partial^{2} D_{z}^{(1)}}{\partial z^{2}}] β sin [β (- z - \frac{1}{2} (h_{e} + h_{p}))] = - F_{31}^{(1)} \frac{\partial^{2} W}{\partial x^{2}} - X_{33}^{(1)} ϕ, \\ b \int_{\frac{h_{e}}{2}}^{\frac{h_{e}}{2} + h_{p}} [D_{z}^{(3)} - {(e_{0} a)}^{2} \frac{\partial^{2} D_{z}^{(3)}}{\partial z^{2}}] β sin [β (- z + \frac{1}{2} (h_{e} + h_{p}))] = - F_{31}^{(3)} \frac{\partial^{2} W}{\partial x^{2}} - X_{33}^{(3)} ϕ, \end{matrix}

(37)

where

\begin{matrix} X_{33}^{(1)} = \int_{- \frac{h_{e}}{2} - h_{p}}^{- \frac{h_{e}}{2}} b \in_{33} {\{β sin [β (- z - \frac{1}{2} (h_{e} + h_{p}))]\}}^{2} d z, \\ X_{33}^{(3)} = \int_{\frac{h_{e}}{2}}^{\frac{h_{e}}{2} + h_{p}} b \in_{33} {\{β sin [β (- z + \frac{1}{2} (h_{e} + h_{p}))]\}}^{2} d z . \end{matrix}

(38)

Summing Equations (32-a), (32-b) and (32-c) yields:

\begin{matrix} M_{1} + M_{2} + M_{3} & = & {(e_{0} a)}^{2} ({M^{″}}_{1} + {M^{″}}_{2} + {M^{″}}_{3}) + ϕ (F_{31}^{(1)} + F_{31}^{(3)}) \\ - (W^{″} + \frac{1}{2} W^{″} {W^{'}}^{2}) (D_{11}^{(1)} + D_{11}^{(2)} + D_{11}^{(3)}) . \end{matrix}

(39)

By inserting Equation (31) in (39), the total bending moment is given:

\begin{matrix} M_{T} & = & F_{31} ϕ - D_{11} (W^{″} + \frac{1}{2} W^{″} {W^{'}}^{2}) + {(e_{0} a)}^{2} (C \dot{W} + m_{T} Ẅ) \\ + {(e_{0} a)}^{2} {(W^{'} \int_{l}^{x} m_{T} \int_{0}^{x} ({\dot{W^{'}}}^{2} + W^{'} Ẅ^{'}) d x d x)}^{'}, \end{matrix}

(40)

where

\begin{matrix} M_{T} = M_{1} + M_{2} + M_{3}, \\ m_{T} = m_{t o t} + \frac{m_{t i p}}{l_{2}} H (x - l_{2}), m_{tot} = m_{1} + m_{2} + m_{3}, \\ F_{31} = F_{31}^{(1)} + F_{31}^{(3)}, \\ D_{11} = D_{11}^{(1)} + D_{11}^{(2)} + D_{11}^{(3)} \end{matrix}

(41)

Eventually, the first equation governing on the transverse vibrations of the nanobeam is obtained by submitting Equation (40) in (30) and keeping terms up to the third order as:

\begin{matrix} D_{11} (W^{''''} + {W^{″}}^{3} + 4 W^{'} W^{″} W^{‴} + W^{''''} {W^{'}}^{2}) - (1 + \frac{1}{2} {W^{'}}^{2}) F_{31} ϕ^{″} - W^{'} W^{″} F_{31} ϕ^{'} \\ - W^{'} W^{″} {(e_{0} a)}^{2} {(C \dot{W} + m_{T} Ẅ)}^{'} + C \dot{W} + m_{T} Ẅ \\ + {(W^{'} \int_{l}^{x} m_{T} \int_{0}^{x} ({\dot{W^{'}}}^{2} + W^{'} Ẅ^{'}) d x d x)}^{'} \\ - {(e_{0} a)}^{2} {(W^{'} \int_{l}^{x} m_{T} \int_{0}^{x} ({\dot{W^{'}}}^{2} + W^{'} Ẅ^{'}) d x d x)}^{'''} \\ - (1 + \frac{1}{2} {W^{'}}^{2}) {(e_{0} a)}^{2} {(C \dot{W} + m_{T} Ẅ)}^{''} = 0 . \end{matrix}

(42)

After deriving the first governing equation, the electrical charge equation has to be derived. For that purpose, Equations (33) and (34) are introduced into (21):

X_{11} \frac{\partial^{2} ϕ}{\partial x^{2}} - X_{33} ϕ - F_{31} \frac{\partial^{2} W}{\partial x^{2}} = 0,

(43)

where

\begin{matrix} X_{11} = X_{11}^{(1)} + X_{11}^{(3)}, \\ X_{33} = X_{33}^{(1)} + X_{33}^{(3)}, \end{matrix}

(44)

Because the clamped end of the nanobeam has a translational motion, i.e. W_(0,t) = W_base(t), the absolute transverse vibrations of the nanobeamcan be determined as below:

W (x, t) = W_{rel} (x, t) + W_{base} (t),

(45)

where W_rel is the transverse displacement relative to the base and W_base is the harmonic motion of the base that is assumed in the following form:

W_{base} = A_{b} sin (ω_{b} t),

(46)

where ω_b and A_b are frequency and amplitude of the base vibrations, respectively. By applying Equation (46), Equation (45) can be rewritten as:

W = W_{rel} + A_{b} sin (ω_{b} t) .

(47)

Finally, equations of motion of a nanobeam subjected to a harmonic base excitation is obtained by submitting Equation (47) in Equations (42) and (43):

\begin{matrix} D_{11} (W_{rel}^{''''} + {W_{rel}^{″}}^{3} + 4 W_{rel}^{'} W_{rel}^{″} W_{rel}^{‴} + W_{rel}^{''''} {W_{rel}^{'}}^{2}) - (1 + \frac{1}{2} {W_{rel}^{'}}^{2}) F_{31} ϕ^{″} \\ - W_{rel}^{'} W_{rel}^{″} F_{31} ϕ^{'} - (1 + \frac{1}{2} {W_{rel}^{'}}^{2}) {(e_{0} a)}^{2} {(C {\dot{W}}_{rel} + m_{T} Ẅ_{rel})}^{''} \\ - {(e_{0} a)}^{2} {(W_{rel}^{'} \int_{l}^{x} m_{T} \int_{0}^{x} ({\dot{W}}_{rel}^{'}^{2} + W_{rel}^{'} {Ẅ^{'}}_{rel}) d x d x)}^{'''} \\ - W_{rel}^{'} W_{rel}^{″} {(e_{0} a)}^{2} {(C {\dot{W}}_{rel} + m_{T} Ẅ_{rel})}^{'} + C {\dot{W}}_{rel} + m_{T} Ẅ_{rel} \\ + {(W_{rel}^{'} \int_{l}^{x} m_{T} \int_{0}^{x} ({\dot{W}}_{rel}^{'}^{2} + W_{rel}^{'} Ẅ_{rel}^{'}) d x d x)}^{'} \\ = m_{T} A_{b} {ω_{b}}^{2} s i n (ω_{b} t) - C A_{b} ω_{b} c o s (ω_{b} t) \end{matrix}

(48)

X_{11} \frac{\partial^{2} ϕ}{\partial x^{2}} - X_{33} ϕ - F_{31} \frac{\partial^{2} W_{rel}}{\partial x^{2}} = 0 .

(49)

By defining the following dimensionless parameters:

\begin{matrix} ξ = \frac{x}{L}, w = \frac{W_{rel}}{h_{e}}, φ = \frac{ϕ}{ϕ_{0}}, ϕ_{0} = \sqrt{\frac{A_{11}}{X_{11}^{3}}}, A_{11} = \int_{- \frac{h_{e}}{2}}^{\frac{h_{e}}{2}} b C_{11}^{e} d z, \\ τ = \frac{t}{t_{0}}, η = \frac{h_{e}}{L}, μ = \frac{e_{0} a}{L}, {\bar{ω}}_{b} = ω_{b} t_{0}, \\ t_{0} = L \sqrt{\frac{m_{2}}{A_{11}}}, {\bar{F}}_{31}^{(1)} = \frac{F_{31}^{(1)} ϕ_{0}}{A_{11} h_{e}}, {\bar{F}}_{31}^{(3)} = \frac{F_{31}^{(3)} ϕ_{0}}{A_{11} h_{e}}, {\bar{F}}_{31} = {\bar{F}}_{31}^{(1)} + {\bar{F}}_{31}^{(3)}, \\ {\bar{D}}_{11}^{(1)} = \frac{D_{11}^{(1)}}{L^{2} A_{11}}, {\bar{D}}_{11}^{(2)} = \frac{D_{11}^{(2)}}{L^{2} A_{11}}, {\bar{D}}_{11}^{(3)} = \frac{D_{11}^{(3)}}{L^{2} A_{11}}, {\bar{D}}_{11} = {\bar{D}}_{11}^{(1)} + {\bar{D}}_{11}^{(2)} + {\bar{D}}_{11}^{(3)}, \\ {\bar{X}}_{11}^{(1)} = \frac{X_{11}^{(1)} ϕ_{0}^{2}}{A_{11} {h_{e}}^{2}}, {\bar{X}}_{11}^{(3)} = \frac{X_{11}^{(3)} ϕ_{0}^{2}}{A_{11} {h_{e}}^{2}}, {\bar{X}}_{11} = {\bar{X}}_{11}^{(1)} + {\bar{X}}_{11}^{(3)}, {\bar{X}}_{33}^{(1)} = \frac{X_{33}^{(1)} ϕ_{0}^{2} L^{2}}{A_{11} {h_{e}}^{2}}, \\ {\bar{X}}_{33}^{(3)} = \frac{X_{33}^{(3)} ϕ_{0}^{2} L^{2}}{A_{11} {h_{e}}^{2}}, {\bar{X}}_{33} = {\bar{X}}_{33}^{(1)} + {\bar{X}}_{33}^{(3)}, \bar{C} = \frac{C L^{2}}{t_{0} A_{11}} = \frac{C L}{\sqrt{A_{11} m_{2}}} . \end{matrix}

(50)

The nondimensional form of Equations (48) and (49) are given as:

\begin{matrix} {\bar{D}}_{11} w^{''''} + {\bar{D}}_{11} η^{2} ({w^{″}}^{3} + 4 w^{'} w^{″} w^{‴} + w^{''''} w^{' 2}) - {\bar{F}}_{31} (1 + \frac{η^{2}}{2} {w^{'}}^{2}) φ^{″} - {\bar{F}}_{31} η^{2} w^{'} w^{″} φ^{'} \\ - μ^{2} (1 + \frac{η^{2}}{2} {w^{'}}^{2}) {(\bar{C} \dot{w} + \frac{m_{T}}{m_{2}} ẅ)}^{''} \\ - μ^{2} η^{2} {[w^{'} \int_{1}^{ξ} \frac{m_{T}}{m_{2}} \int_{0}^{ξ} ({\dot{w^{'}}}^{2} + w^{'} ẅ^{'}) d ξ d ξ]}^{'''} - μ^{2} η^{2} w^{'} w^{″} {(\bar{C} \dot{w} + \frac{m_{T}}{m_{2}} ẅ)}^{'} \\ + \bar{C} \dot{w} + \frac{m_{T}}{m_{2}} ẅ + η^{2} {[w^{'} \int_{1}^{ξ} \frac{m_{T}}{m_{2}} \int_{0}^{ξ} ({\dot{w^{'}}}^{2} + w^{'} ẅ^{'}) d ξ d ξ]}^{'} \\ = \frac{m_{T}}{m_{2}} \frac{A_{b}}{h_{e}} {\bar{ω}}_{b}^{2} sin ({\bar{ω}}_{b} τ) - \bar{C} \frac{A_{b}}{h_{e}} {\bar{ω}}_{b} cos ({\bar{ω}}_{b} τ) \end{matrix}

(51)

{\bar{X}}_{11} \frac{\partial^{2} φ}{\partial ξ^{2}} - {\bar{X}}_{33} φ - {\bar{F}}_{31} \frac{\partial^{2} w}{\partial ξ^{2}} = 0 .

(52)

Where ${\bar{ω}}_{b}$ is the dimensionless frequency of base excitation.

V. DISCRETIZATION

In order to discretize the obtained partial differential equations (Equations (51) and (52)) and reducing them to ordinary differential equations, assumed mode method is applied. In this method, the transverse displacement w and the generated voltage φ are expressed as:

w (ξ, τ) = \sum_{i = 1}^{n} q_{i} (τ) P_{i} (ξ),

(53)

φ (ξ, τ) = \sum_{i = 1}^{n} s_{i} (τ) r_{i} (ξ),

(54)

where n is number of applied modes in these two expansions. $P_{i} (ξ)$ and $r_{i} (ξ)$ are test functions of transverse displacements and genereated voltage, respectively. The applied test functions have to satisfy kinematic boundary conditions. $q_{i} (τ)$ and $s_{i} (τ)$ are generalized coordinates that have to be found. $P_{i} (ξ)$ are considered as linear mode shapes of a cantilever beam:²⁷

P_{i} (ξ) = cosh (ψ_{i} ξ) - cos (ψ_{i} ξ) - \frac{sinh (ψ_{i}) - sin (ψ_{i})}{cosh (ψ_{i}) + cos (ψ_{i})} [sinh (ψ_{i} ξ) - sin (ψ_{i} ξ)],

(55)

where $ψ_{i}$ are roots of the following characteristic equation:

cosh (ψ_{i}) cos (ψ_{i}) = - 1 .

(56)

In order to satisfy the electrical charge boundary conditions, $r_{i} (ξ)$ is assumed as below:

r_{i} (ξ) = cos [(\frac{2 i - 1}{2}) π ξ] .

(57)

In this research, only one mode has been employed. Therefore, subscripts in above equations will be removed in the following equations. Now by introducing Equations (53) and (54) into Equation (52) then multiplying the resulting equation by $P (ξ)$ and integrating along the length of the beam, s is obtained with respect to q:

s = q \frac{\int_{0}^{1} {\bar{F}}_{31} P P^{″} d ξ}{\int_{0}^{1} {\bar{X}}_{11} P r^{″} d ξ - \int_{0}^{1} P {\bar{X}}_{33} r d ξ} .

(58)

Then by substituting Equations (53), (54) and (58) in (51) then multiplying by $P (ξ)$ and integrating along the length of the beam, an ordinary differential equation is obtained:

α_{1} \ddot{q} + α_{2} q + α_{3} q^{3} + α_{4} q {\dot{q}}^{2} + α_{5} q^{2} \ddot{q} + α_{6} \dot{q} + α_{7} q^{2} \dot{q} = α_{8} {\bar{ω}}_{b}^{2} sin ({\bar{ω}}_{b} τ) + α_{9} {\bar{ω}}_{b} cos ({\bar{ω}}_{b} τ),

(59)

where

α_{1} = \int_{0}^{1} \frac{m_{T}}{m_{2}} P^{2} d ξ - μ^{2} \int_{0}^{1} \frac{m_{T}}{m_{2}} P P^{″} d ξ - μ^{2} \frac{m_{tip}}{m_{2} l_{2}} [\int_{0}^{1} P^{2} H^{″} (ξ - \frac{l_{2}}{l}) d ξ + \int_{0}^{1} 2 P P^{'} H^{'} (ξ - \frac{l_{2}}{l}) d ξ],

α_{2} = {\bar{D}}_{11} \int_{0}^{1} P P^{″ ″} d ξ - \frac{\int_{0}^{1} {\bar{F}}_{31} P P^{″} d ξ \int_{0}^{1} {\bar{F}}_{31} P r^{″} d ξ}{\int_{0}^{1} {\bar{X}}_{11} P r^{″} d ξ - \int_{0}^{1} {\bar{X}}_{33} P r d ξ},

α_{3} = {\bar{D}}_{11} η^{2} \int_{0}^{1} ({P^{″}}^{3} P + 4 P P^{'} P^{″} P^{‴} + P P^{''''} {P^{'}}^{2}) d ξ - \frac{1}{2} η^{2} \frac{\int_{0}^{1} {\bar{F}}_{31} P P^{″} d ξ \int_{0}^{1} {P^{'}}^{2} {\bar{F}}_{31} P r^{″} d ξ}{\int_{0}^{1} {\bar{X}}_{11} P r^{″} d ξ - \int_{0}^{1} {\bar{X}}_{33} P r d ξ} - η^{2} \frac{\int_{0}^{1} {\bar{F}}_{31} P P^{″} d ξ \int_{0}^{1} {\bar{F}}_{31} P P^{'} P^{″} r^{'} d ξ}{\int_{0}^{1} {\bar{X}}_{11} P r^{″} d ξ - \int_{0}^{1} {\bar{X}}_{33} P r d ξ},

(60)

α_{4} = η^{2} \int_{0}^{1} \{P {[P^{'} \int_{1}^{ξ} \frac{m_{T}}{m_{2}} \int_{0}^{ξ} {P^{'}}^{2} d ξ d ξ]}^{'}\} d ξ - μ^{2} η^{2} \int_{0}^{1} \{P {[P^{'} \int_{1}^{ξ} \frac{m_{T}}{m_{2}} \int_{0}^{ξ} {P^{'}}^{2} d ξ d ξ]}^{'''}\} d ξ,

α_{5} = α_{4} - μ^{2} η^{2} [\int_{0}^{1} \frac{m_{T}}{m_{2}} P {P^{'}}^{2} P^{″} d ξ + \int_{0}^{1} \frac{m_{tip}}{m_{2} l_{2}} P^{'} P^{2} P^{″} d ξ] - \frac{η^{2}}{2} μ^{2} \{\int_{0}^{1} \frac{m_{T}}{m_{2}} P {P^{'}}^{2} P^{″} d ξ) (+ \frac{m_{tip}}{m_{2} l_{2}} [\int_{0}^{1} P^{2} {P^{'}}^{2} H^{″} (ξ - \frac{l_{2}}{l}) d ξ + \int_{0}^{1} 2 {P^{'}}^{2} P P^{'} H^{'} (ξ - \frac{l_{2}}{l})]\},

α_{6} = \bar{C} \int_{0}^{1} P^{2} d ξ - μ^{2} \bar{C} \int_{0}^{1} P P^{″} d ξ,

α_{7} = - \frac{3}{2} μ^{2} η^{2} \bar{C} \int_{0}^{1} P {P^{'}}^{2} P^{″} d ξ,

α_{8} = \int_{0}^{1} P \frac{m_{T}}{m_{2}} \frac{A_{b}}{h_{e}} d ξ,

α_{9} = - \int_{0}^{1} P \bar{C} \frac{A_{b}}{h_{e}} d ξ .

Equation (59) can be rewritten as:

\ddot{q} + β_{1} q + β_{2} q^{3} + β_{3} q {\dot{q}}^{2} + β_{4} q^{2} \ddot{q} + β_{5} \dot{q} + β_{6} q^{2} \dot{q} = β_{7} {\bar{ω}}_{b}^{2} sin ({\bar{ω}}_{b} τ) + β_{8} {\bar{ω}}_{b} cos ({\bar{ω}}_{b} τ) .

(61)

VI. METHOD OF MULTIPLE SCALES

In order to plot the frequency-response and force-response diagrams of the nanoscale PEH, multiple scale method has been employed. In order to apply this method, nonlinear terms have to be balanced with the terms of damping and excitation by considering appropriate orders of ε for them. For this purpose, Equation (61) is rewritten as:

\ddot{q} + β_{1} q + β_{2} q^{3} + β_{3} q {\dot{q}}^{2} + β_{4} q^{2} \ddot{q} + ε^{2} β_{5} \dot{q} + β_{6} q^{2} \dot{q} = ε^{3} β_{7} {\bar{ω}}_{b}^{2} sin ({\bar{ω}}_{b} τ) + ε^{3} β_{8} {\bar{ω}}_{b} cos ({\bar{ω}}_{b} τ) .

(62)

q can be expanded as

q = ε q_{1} (T_{0}, T_{2}) + ε^{3} q_{3} (T_{0}, T_{2}),

(63)

where T₀ and T₂ are time scales that are defined as:

T_{0} = τ, T_{2} = ε^{2} τ .

(64)

It is worth mentioning that Equation (63) excluded q₂ because there is no quadratic nonlinearity in Equation (62). Regarding to the defined time scales, time derivatives are rearranged as:

\frac{d}{d τ} = D_{0} + ε^{2} D_{2}, \frac{d^{2}}{d τ^{2}} = {D_{0}}^{2} + 2 ε^{2} D_{0} D_{2} .

(65)

The dimensionless frequency of the base excitation is considered to be close to the dimensionless linear natural frequency of the nanobeam:

{\bar{ω}}_{b} = \bar{ω} + ε^{2} σ,

(66)

where $\bar{ω}$ is dimensionless linear natural frequency and σ is detuning parameter. Substituting Equations (63), (65) and (66) into (62) yields:

\begin{matrix} ({D_{0}}^{2} + 2 ε^{2} D_{0} D_{2}) (ε q_{1} + ε^{3} q_{3}) + β_{1} (ε q_{1} + ε^{3} q_{3}) + β_{2} {(ε q_{1} + ε^{3} q_{3})}^{3} \\ + β_{3} (ε q_{1} + ε^{3} q_{3}) {[(D_{0} + ε^{2} D_{2}) (ε q_{1} + ε^{3} q_{3})]}^{2} \\ + β_{4} {(ε q_{1} + ε^{3} q_{3})}^{2} ({D_{0}}^{2} + 2 ε^{2} D_{0} D_{2}) (ε q_{1} + ε^{3} q_{3}) \\ + ε^{2} β_{5} (D_{0} + ε^{2} D_{2}) (ε q_{1} + ε^{3} q_{3}) \\ + β_{6} {(ε q_{1} + ε^{3} q_{3})}^{2} (D_{0} + ε^{2} D_{2}) (ε q_{1} + ε^{3} q_{3}) \\ - ε^{3} β_{7} {(\bar{ω} + ε^{2} σ)}^{2} sin ((\bar{ω} + ε^{2} σ) T_{0}) - ε^{3} β_{8} (\bar{ω} + ε^{2} σ) cos ((\bar{ω} + ε^{2} σ) T_{0}) \\ = 0 . \end{matrix}

(67)

By equating coefficients of $ε^{1}$ and $ε^{3}$ to zero, the following two equations are acquired:

ε^{1} : D_{0}^{2} q_{1} + β_{1} q_{1} = 0, β_{1} = {\bar{ω}}^{2},

(68)

\begin{matrix} ε^{3} : D_{0}^{2} q_{3} + β_{1} q_{3} = - 2 D_{0} D_{2} q_{1} - β_{2} q_{1}^{3} - β_{3} q_{1} {(D_{0} q_{1})}^{2} - β_{4} q_{1}^{2} D_{0}^{2} q_{1} - β_{5} D_{0} q_{1} - β_{6} q_{1}^{2} D_{0} q_{1} \\ + β_{7} {\bar{ω}}^{2} sin ((\bar{ω} + ε^{2} σ) T_{0}) + β_{8} \bar{ω} cos ((\bar{ω} + ε^{2} σ) T_{0}) . \end{matrix}

(69)

By solving linear Equation (68), $q_{1}$ is given as:

q_{1} = A (T_{2}) exp (i \bar{ω} T_{0}) + Ā (T_{2}) exp (- i \bar{ω} T_{0}) .

(70)

By introducing Equation (70) into (69), then setting coefficients of secular terms, i.e. $exp (i \bar{ω} T_{0})$ ⁠, to zero, solvability condition is obtained as:

\begin{matrix} - 2 i \bar{ω} A^{'} - 3 β_{2} A^{2} Ā - β_{3} A^{2} Ā {\bar{ω}}^{2} + 3 β_{4} A^{2} Ā {\bar{ω}}^{2} - i β_{5} A \bar{ω} - i β_{6} A^{2} Ā \bar{ω} \\ + \frac{1}{2} β_{8} \bar{ω} exp (i σ T_{2}) - \frac{1}{2} i β_{7} {\bar{ω}}^{2} exp (i σ T_{2}) = 0 . \end{matrix}

(71)

Writing A in polar form as

A = \frac{1}{2} P (T_{2}) exp (i γ (T_{2}))

(72)

and then inserting it in Equation (71) yields:

\begin{matrix} - 2 i \bar{ω} (\frac{1}{2} P^{'} exp (i γ) + \frac{1}{2} P i γ^{'} exp (i γ)) - 3 β_{2} (\frac{1}{4} P^{2} exp (i 2 γ)) (\frac{1}{2} P exp (- i γ)) \\ - β_{3} {\bar{ω}}^{2} (\frac{1}{4} P^{2} exp (i 2 γ)) (\frac{1}{2} P exp (- i γ)) + 3 β_{4} {\bar{ω}}^{2} (\frac{1}{4} P^{2} e^{i 2 γ}) (\frac{1}{2} P exp (- i γ)) \\ - i β_{5} \bar{ω} (\frac{1}{2} P exp (i γ)) - i \bar{ω} β_{6} (\frac{1}{4} P^{2} exp (i 2 γ)) (\frac{1}{2} P exp (- i γ)) \\ + \frac{1}{2} β_{8} \bar{ω} exp (i σ T_{2}) - \frac{1}{2} i β_{7} {\bar{ω}}^{2} exp (i σ T_{2}) = 0 . \end{matrix}

(73)

Collecting the coefficients of $exp (i γ)$ in Equation (73) results in:

\begin{matrix} - i \bar{ω} (P^{'} + i P γ^{'}) - \frac{3}{8} β_{2} P^{3} - \frac{1}{8} β_{3} {\bar{ω}}^{2} P^{3} + \frac{3}{8} β_{4} {\bar{ω}}^{2} P^{3} - \frac{1}{2} P \bar{ω} β_{5} i - \frac{1}{8} i β_{6} \bar{ω} P^{3} \\ - \frac{1}{2} i β_{7} {\bar{ω}}^{2} exp (i σ T_{2} - i γ) + \frac{1}{2} β_{8} \bar{ω} exp (i σ T_{2} - i γ) = 0 . \end{matrix}

(74)

By defining

δ = σ T_{2} - γ .

(75)

and employing Euler’s formula,

exp (i δ) = cos (δ) + i sin (δ),

(76)

Equation (74) can be rewritten as:

\begin{matrix} - i \bar{ω} (P^{'} + i P (σ - δ^{'})) - \frac{3}{8} β_{2} P^{3} - \frac{1}{8} β_{3} {\bar{ω}}^{2} P^{3} + \frac{3}{8} β_{4} {\bar{ω}}^{2} P^{3} - \frac{1}{2} P \bar{ω} β_{5} i - \frac{1}{8} i β_{6} \bar{ω} P^{3} \\ - \frac{1}{2} i β_{7} {\bar{ω}}^{2} (cos (δ) + i sin (δ)) + \frac{1}{2} β_{8} \bar{ω} (cos (δ) + i sin (δ)) = 0 . \end{matrix}

(77)

By separating real and imaginary parts, modulating equations are obtained as:

P^{'} = - \frac{1}{2} P β_{5} - \frac{1}{8} β_{6} P^{3} - \frac{1}{2} β_{7} \bar{ω} cos (δ) + \frac{1}{2} β_{8} sin (δ),

(78)

P δ^{'} = P σ - \frac{3}{8 \bar{ω}} β_{2} P^{3} - \frac{1}{8} β_{3} \bar{ω} P^{3} + \frac{3}{8} β_{4} \bar{ω} P^{3} - \frac{1}{2} β_{7} \bar{ω} sin (δ) - \frac{1}{2} β_{8} cos (δ) .

(79)

In order to achieve a steady state solution, the time derivatives have to be set to zero:

P^{'} = P δ^{'} = 0,

(80)

By inserting Equation (80) in Equations (78) and (79), the frequency-response equation is obtained:

- \frac{1}{4} {β_{7}}^{2} {\bar{ω}}^{2} - \frac{1}{4} {β_{8}}^{2} + {(- \frac{1}{2} P β_{5} - \frac{1}{8} β_{6} P^{3})}^{2} + {(P σ - \frac{3}{8 \bar{ω}} β_{2} P^{3} - \frac{1}{8} β_{3} \bar{ω} P^{3} + \frac{3}{8} β_{4} \bar{ω} P^{3})}^{2} = 0,

(81)

where

σ = {\bar{ω}}_{b} - \bar{ω} .

(82)

VII. RESULTS AND DISCUSSION

In this section, the results of a numerical study on the proposed model are presented. Material and geometric properties of the studied system is listed in Table I. First, the frequency-response curve of the nanoscale PEH, obtained by a multiple scale method, is compared with a numeric method. Then, frequency-response curves are investigated for different values of system parameters.

figure 2 shows comparison between frequency-response that is obtained by the multiple-scale method and also a numeric method. The applied numeric method is the forth-order Runge-Kutta. In that figure, a good agreement between the two applied solution methods can be seen. Furthermore, stable and unstable branches are represented in figure 2.

The importance of considering the size effects in the model can also be manifested. Frequency-response curves for different values of scale factor in absence of the tip mass is depicted in figure 3. An interesting characteristic can be seen: For small values of the scale factor (i.e. lower than μ ≈ 0.165) the nonlinear behavior of the nanobeam is of hardening type. For larger values of the scale factor (larger than μ ≈ 0.165), nonlinear behavior is softening. Therefore, ignoring nonlocal effects results in an incorrect predictions of a nanoscale PEH behavior. In addition, it is observed that increasing the scale factor yields increases in both maximum vibration amplitude and generated voltage of the PEH.

figure 4 displays the frequency-response curves for different values of the scale factor and tip mass. It is observed that in presence of a tip mass, nanobeam demonstrates softening behavior and this behavior is intensified by increasing the scale factor. Increasing the scale factor also increases the maximum amplitude of the generated voltage. It should be noted that increasing the value of the tip mass reduces the values of nonlinear natural frequencies more than linear natural frequencies of the PEH. Therefore, the frequency-response curve bends to higher values of the detuning parameter.

The frequency-response curves for different values of the tip mass are demonstrated in figure 5. It can be seen that increasing the value of the tip mass bends the frequency-response curve away from σ = 0 to the left and increases the amplitude of the generated voltage in a broad band of frequency.

After investigating the frequency–response curves, the effects of scale factor and excitation frequency on the force response curves are studied. A Force-response curve determines the multiple-valued regions for a known harmonic excitation. Multiple-valued region is a range of values of excitation amplitude that multiple values of response amplitude are observed for each excitation amplitude.

The force-response curves for different values of scale factor and detuning parameter are plotted in figures 6 and 7, respectively. According to figure 6, by increasing the value of scale factor, the jumps occur in higher values of the excitation amplitude. In figure 7, it can be observed that for some values of the detuning parameter, there is no multiple-value region.

VIII. CONCLUSION

With the growing application of PEH in small electronic devices, preparing a more accurate model for this type of energy harvesters is an essential task to maximize the harvested energy. For this purpose, the nonlinear vibrations and generated voltage of a nanoscale PEH based on nonlocal elasticity theory was studied in this research. Nanoscale PEH was modeled as a piezoelectric bimorph cantilever nanobeam using Euler-Bernoulli beam theory. The nanobeam was assumed to be under a harmonic base excitation. In order to consider the size effects in the prepared model, the nanlocal elasticity theory as a strong nonclassic theory was applied. Nonlinear equations of motion were derived and then solved using Hamilton’s principle and a multiple scales method, respectively. A numerical method was also used to validate the obtained analytical solution. After solving the equations of motion, Frequency-response and Force-response curves were plotted. Also, effects of scale factor and tip mass on the generated voltage were investigated. Results showed the significant effect of the scale factor in the nonlinear behavior of the nanoscale PEH. Therefore neglecting nonlocal effects leads to an incorrect prediction of behavior of nanoscale PEH. Moreover, the results showed an increase in the generated voltage and amplitude of vibrations by increasing the scale factor and also nanobeam tip mass. To sum up, this paper provided a more accurate model for nanoscale piezoelectric energy harvesters, and the presented modeling procedure could pave the way for modeling other types of nanoscale vibration energy harvesters.

ACKNOWLEDGMENTS

The authors would like to thank Shahid Chamran University of Ahvaz for its financial support.

REFERENCES

1.

F. R.

Fan

,

W.

Tang

, and

Z. L.

Wang

, “

Flexible nanogenerators for energy harvesting and self-powered electronics

,”

Advanced Materials

28

(

22

),

4283

–

4305

(

2016

).

https://doi.org/10.1002/adma.201504299

Google Scholar

Crossref

PubMed

2.

T.

Fan

and

L.

Yang

, “

Surface effect on nano piezoelectric energy harvester based on flexural mode

,”

Polymer Composites

(

2016

).

https://doi.org/10.1002/pc.24020

Google Scholar

3.

A.

Toprak

and

O.

Tigli

, “

Piezoelectric energy harvesting: State-of-the-art and challenges

,”

Applied Physics Reviews

1

(

3

),

031104

(

2014

).

https://doi.org/10.1063/1.4896166

Google Scholar

Crossref

4.

X.

Wang

, “

Piezoelectric nanogenerators—Harvesting ambient mechanical energy at the nanometer scale

,”

Nano Energy

1

(

1

),

13

–

24

(

2012

).

https://doi.org/10.1016/j.nanoen.2011.09.001

Google Scholar

Crossref

5.

X.

Chen

et al, “

1.6 V nanogenerator for mechanical energy harvesting using PZT nanofibers

,”

Nano Letters

10

(

6

),

2133

–

2137

(

2010

).

https://doi.org/10.1021/nl100812k

Google Scholar

Crossref

PubMed

6.

K.

Wang

and

B.

Wang

, “

An analytical model for nanoscale unimorph piezoelectric energy harvesters with flexoelectric effect

,”

Composite Structures

153

,

253

–

261

(

2016

).

https://doi.org/10.1016/j.compstruct.2016.05.104

Google Scholar

Crossref

7.

K.

Wang

and

B.

Wang

, “

Non-linear flexoelectricity in energy harvesting

,”

International Journal of Engineering Science

116

,

88

–

103

(

2017

).

https://doi.org/10.1016/j.ijengsci.2017.02.010

Google Scholar

Crossref

8.

Q.

Deng

et al, “

Nanoscale flexoelectric energy harvesting

,”

International Journal of Solids and Structures

51

(

18

),

3218

–

3225

(

2014

).

https://doi.org/10.1016/j.ijsolstr.2014.05.018

Google Scholar

Crossref

9.

O.

Aldraihem

and

A.

Baz

, “

Energy harvester with a dynamic magnifier

,”

Journal of Intelligent Material Systems and Structures

22

(

6

),

521

–

530

(

2011

).

https://doi.org/10.1177/1045389x11402706

Google Scholar

Crossref

10.

S.

Ali

,

M.

Friswell

, and

S.

Adhikari

, “

Piezoelectric energy harvesting with parametric uncertainty

,”

Smart Materials and Structures

19

(

10

),

105010

(

2010

).

https://doi.org/10.1088/0964-1726/19/10/105010

Google Scholar

Crossref

11.

M.

Fakhzan

and

A. G.

Muthalif

, “

Harvesting vibration energy using piezoelectric material: Modeling, simulation and experimental verifications

,”

Mechatronics

23

(

1

),

61

–

66

(

2013

).

https://doi.org/10.1016/j.mechatronics.2012.10.009

Google Scholar

Crossref

12.

H. S.

Kim

,

J.-H.

Kim

, and

J.

Kim

, “

A review of piezoelectric energy harvesting based on vibration

,”

International Journal of Precision Engineering and Manufacturing

12

(

6

),

1129

–

1141

(

2011

).

https://doi.org/10.1007/s12541-011-0151-3

Google Scholar

Crossref

13.

G.

Wang

, “

Analysis of bimorph piezoelectric beam energy harvesters using Timoshenko and Euler–Bernoulli beam theory

,”

Journal of Intelligent Material Systems and Structures

24

(

2

),

226

–

239

(

2013

).

https://doi.org/10.1177/1045389x12461080

Google Scholar

Crossref

14.

A.

Erturk

and

D. J.

Inman

, “

On mechanical modeling of cantilevered piezoelectric vibration energy harvesters

,”

Journal of Intelligent Material Systems and Structures

19

(

11

),

1311

–

1325

(

2008

).

https://doi.org/10.1177/1045389x07085639

Google Scholar

Crossref

15.

F.

Cottone

et al, “

Piezoelectric buckled beams for random vibration energy harvesting

,”

Smart Materials and Structures

21

(

3

),

035021

(

2012

).

https://doi.org/10.1088/0964-1726/21/3/035021

Google Scholar

Crossref

16.

G.

Sebald

et al, “

Experimental Duffing oscillator for broadband piezoelectric energy harvesting

,”

Smart Materials and Structures

20

(

10

),

102001

(

2011

).

https://doi.org/10.1088/0964-1726/20/10/102001

Google Scholar

Crossref

17.

L.

Tang

and

Y.

Yang

, “

A nonlinear piezoelectric energy harvester with magnetic oscillator

,”

Applied Physics Letters

101

(

9

),

094102

(

2012

).

https://doi.org/10.1063/1.4748794

Google Scholar

Crossref

18.

M.

Eltaher

,

S. A.

Emam

, and

F.

Mahmoud

, “

Free vibration analysis of functionally graded size-dependent nanobeams

,”

Applied Mathematics and Computation

218

(

14

),

7406

–

7420

(

2012

).

https://doi.org/10.1016/j.amc.2011.12.090

Google Scholar

Crossref

19.

S.

Pradhan

and

T.

Murmu

, “

Application of nonlocal elasticity and DQM in the flapwise bending vibration of a rotating nanocantilever

,”

Physica E: Low-Dimensional Systems and Nanostructures

42

(

7

),

1944

–

1949

(

2010

).

https://doi.org/10.1016/j.physe.2010.03.004

Google Scholar

Crossref

20.

A. C.

Eringen

and

D.

Edelen

, “

On nonlocal elasticity

,”

International Journal of Engineering Science

10

(

3

),

233

–

248

(

1972

).

https://doi.org/10.1016/0020-7225(72)90039-0

Google Scholar

Crossref

21.

L.-L.

Ke

and

Y.-S.

Wang

, “

Thermoelectric-mechanical vibration of piezoelectric nanobeams based on the nonlocal theory

,”

Smart Materials and Structures

21

(

2

),

025018

(

2012

).

https://doi.org/10.1088/0964-1726/21/2/025018

Google Scholar

Crossref

22.

L.-L.

Ke

,

Y.-S.

Wang

, and

Z.-D.

Wang

, “

Nonlinear vibration of the piezoelectric nanobeams based on the nonlocal theory

,”

Composite Structures

94

(

6

),

2038

–

2047

(

2012

).

https://doi.org/10.1016/j.compstruct.2012.01.023

Google Scholar

Crossref

23.

M.

Nazemizadeh

and

F.

Bakhtiari-Nejad

, “

Size-dependent free vibration of nano/microbeams with piezo-layered actuators

,”

Micro & Nano Letters

10

(

2

),

93

–

98

(

2015

).

https://doi.org/10.1049/mnl.2014.0317

Google Scholar

Crossref

24.

Q.

Wang

, “

On buckling of column structures with a pair of piezoelectric layers

,”

Engineering Structures

24

(

2

),

199

–

205

(

2002

).

https://doi.org/10.1016/s0141-0296(01)00088-8

Google Scholar

Crossref

25.

S. N.

Mahmoodi

and

N.

Jalili

, “

Non-linear vibrations and frequency response analysis of piezoelectrically driven microcantilevers

,”

International Journal of Non-Linear Mechanics

42

(

4

),

577

–

587

(

2007

).

https://doi.org/10.1016/j.ijnonlinmec.2007.01.019

Google Scholar

Crossref

26.

A. H.

Nayfeh

and

P. F.

Pai

,

Linear and nonlinear structural mechanics

.

2008

:

John Wiley & Sons

.

Google Scholar

27.

D. J.

Inman

and

R. C.

Singh

,

Engineering vibration

. Vol. 3.

2001

:

Prentice Hall Upper Saddle River

.

Google Scholar

2017

Author(s)

Parameter	Unit	Numerical value
$C_{11}^{p}$	GPa	132
$C_{11}^{e}$	GPa	169
e₃₁	Cm⁻²	-4.1
$𝜖_{11}$	CV⁻¹m⁻¹	5.841 × 10⁻⁹
$𝜖_{33}$	CV⁻¹m⁻¹	7.124 × 10⁻⁹
$ρ_{e}$	Kg m⁻³	2330
ρ_p	kg m⁻³	7500
L	m	600 ×10⁻⁹
l₁	m	0.6 × 10⁻⁹
b	m	30 × 10⁻⁹
h_p	m	5 × 10⁻⁹
h_e	m	5 × 10⁻⁹

Parameter	Unit	Numerical value
$C_{11}^{p}$	GPa	132
$C_{11}^{e}$	GPa	169
e₃₁	Cm⁻²	-4.1
$𝜖_{11}$	CV⁻¹m⁻¹	5.841 × 10⁻⁹
$𝜖_{33}$	CV⁻¹m⁻¹	7.124 × 10⁻⁹
$ρ_{e}$	Kg m⁻³	2330
ρ_p	kg m⁻³	7500
L	m	600 ×10⁻⁹
l₁	m	0.6 × 10⁻⁹
b	m	30 × 10⁻⁹
h_p	m	5 × 10⁻⁹
h_e	m	5 × 10⁻⁹

Nanoscale piezoelectric vibration energy harvester design

I. INTRODUCTION

II. THEORETICAL MODELING

III. NONLOCAL ELASTICITY THEORY

IV. DERIVING GOVERNING EQUATIONS

V. DISCRETIZATION

VI. METHOD OF MULTIPLE SCALES

VII. RESULTS AND DISCUSSION

VIII. CONCLUSION

ACKNOWLEDGMENTS

REFERENCES

Contents

Figures & Tables

Supplementary Material

References

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Nanoscale piezoelectric vibration energy harvester design

I. INTRODUCTION

II. THEORETICAL MODELING

III. NONLOCAL ELASTICITY THEORY

IV. DERIVING GOVERNING EQUATIONS

V. DISCRETIZATION

VI. METHOD OF MULTIPLE SCALES

VII. RESULTS AND DISCUSSION

VIII. CONCLUSION

ACKNOWLEDGMENTS

REFERENCES

Contents

Figures & Tables

Supplementary Material

References

Related

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only