Piezoelectric buckled beams for random vibration energy harvesting

The energy harvesting system that we consider is based on a piezoelectric beam clamped on both ends on a base excited vertically by a shaker. One of the two clamps is fixed, while the other can be moved through a micrometric stage.

Figure 1 shows a diagram of the piezoelectric buckled bridge. In our case, it consists of a support of a thin steel beam with one or two bonded layers of piezoelectric material of width b. These are assumed to be polarized only along the z axis. An additional inertial mass can be arranged on the centre of the bridge in order to vary the dynamical response of the system. Two electrodes are placed on both surfaces of the beam and collect the current which then flows into a resistive load. The clamping support is assumed to vibrate only vertically. The total beam deflection is equal to w(x,t). As a first approximation, the main contribution to the electrical field along this direction is due to the electromechanical coupling related to the horizontal strain of the beam. Therefore, the constitutive equations for the piezoelectric effect are as follows [27]:

$$\begin{eqnarray} \begin{array}{@{}rc@{}} \displaystyle &\displaystyle {T}_{x}={c}_{x x}^{\mathrm{E}}{S}_{x}-{e}_{z x}{E}_{z},\\ \displaystyle &\displaystyle {D}_{z}={e}_{z x}{S}_{x}+{\varepsilon }_{z z}^{\mathrm{S}}{E}_{z}, \end{array} \end{eqnarray} \tag{ 1 }$$

where T_x and S_x represent mechanical stress and strain, E_z and D_z are respectively the electric field component and electric displacement along the z axis, ${c}_{x x}^{\mathrm{E}}$ is the elastic stiffness of the piezoceramic layer at constant electric field, and e_zx and ${\varepsilon }_{z z}^{\mathrm{S}}$ are the piezoelectric coupling constant and electric permittivity at zero strain. The structure is initially compressed by moving one end towards the other so as to fit in the distance L. By gradually changing the compression distance ΔL we explore the system electrodynamics from the monostable to the bistable regime.

Zoom In Zoom Out Reset image size

The analytical model of the system is derived from the first-order composite plate theory according to Reddy [28]. The piezoelectric beam is supposed to be homogeneous with constant moment of inertia and stiffness over the whole structure. In order to derive the motion equations of the transverse deflection, we will first proceed to calculate the Lagrangian functional. When the clamping base oscillates with amplitude y with respect to the ground, the total kinetic energy of the beam is the sum of the energy of the steel shim and that related to n piezoelectric layers:

$$\begin{eqnarray} \displaystyle \fl K=\frac{1}{2}{\rho }_{\mathrm{s}}{A}_{\mathrm{s}}\int \nolimits _{0}^{L}{\left(\frac{\partial w(x,t)}{\partial t}+\frac{\mathrm{d}y}{\mathrm{d}t}\right)}^{2}\hspace{0.167em} \mathrm{d}x&&\displaystyle \\ \displaystyle +~\frac{n}{2}{\rho }_{\mathrm{p}}{A}_{\mathrm{p}}\int \nolimits _{0}^{{L}_{\mathrm{p}}}{\left(\frac{\partial w(x,t)}{\partial t}+\frac{\mathrm{d}y}{\mathrm{d}t}\right)}^{2}\hspace{0.167em} \mathrm{d}x&&\displaystyle \\ \displaystyle +~\frac{1}{2}{M}_{0}{\left(\frac{\partial w(L/2,t)}{\partial t}+\frac{\mathrm{d}y}{\mathrm{d}t}\right)}^{2}.&&\displaystyle \end{eqnarray} \tag{ 2 }$$

A symmetrical bimorph piezoelectric beam is assumed for simplicity with L_p = L, whereas w(x,t) is the mid-plane deflection along the z axis, and ρ_s and ρ_p represent the steel shim and piezoelectric material density with sectional area, respectively equal to A_s and A_p. Actually, in the experimental device the piezoelectric layer is not continuous along the x axis. A small gap is left in the middle of the beam in order to have room for an additional mass M₀ to be added. The last term accounts for its kinetic energy. A comprehensive treatment of this case was made in [29]. Moreover, the electrodes are connected in parallel only at one of the two sides of the beam, with length L_e < L. Therefore, in the following, we will not consider a self-cancelling electromechanical coupling coefficient. The strain along the x axis, associated with the mid-plane displacement field (u,w), respectively along the x and z axes, is calculated using von Kàrmàn nonlinear relations [28]:

$$\begin{eqnarray} \displaystyle \fl {S}_{x}={\varepsilon }_{x x}^{0}+z{\varepsilon }_{x x}^{1},\tqs \text{with~ }{\varepsilon }_{x x}^{0}=\frac{\partial u}{\partial x}+\frac{1}{2}{\left(\frac{\partial w}{\partial x}\right)}^{2},&&\displaystyle \\ \displaystyle {\varepsilon }_{x x}^{1}=-\frac{{\partial }^{2}w}{\partial {x}^{2}},&&\displaystyle \end{eqnarray} \tag{ 3 }$$

where the first term of the strain is related to the axial in-plane forces and the second term is the curvature due to the bending moment. Since the electromechanical coupling is present, the total potential energy of the composite laminate includes the bending energy, the work done against the induced electric field and the work done by the external compressive forces. Consequently, the total potential energy is expressed as

$$\begin{eqnarray} &&\Pi =\frac{1}{2}\int \nolimits _{V}{T}_{x}{S}_{x}\hspace{0.167em} \mathrm{d}V-\frac{1}{2}\int \nolimits _{V}{D}_{z}{E}_{z}\hspace{0.167em} \mathrm{d}V-W \end{eqnarray} \tag{ 4 }$$

where V is the entire domain of the composite beam and W is the work done by the external compressive force P applied at the movable tip which is given by

$$\begin{eqnarray} &&W=P \Delta L=P\left({d}_{\mathrm{b}}+\frac{1}{2}\int \nolimits _{0}^{L}\mathop{(w ^{\prime})}\nolimits ^{2}\hspace{0.167em} \mathrm{d}x\right) \end{eqnarray} \tag{ 5 }$$

where w' stands for a derivation with respect to x. The last integral above is nonzero only in the post-buckling case. The compressive load is always taken as larger than the critical load so that P > P_cr in order to realize the buckling. d_b represents the contraction length from side wall pressure corresponding to the critical load P_cr. Manipulation of equation (4) is simplified by introducing the following relations:

$$\begin{eqnarray} \displaystyle (A,B,D)=b\sum _{k=1}^{N}\int _{{z}_{k}}^{{z}_{k+1}}{Q}^{(k)}(1,z,{z}^{2})\hspace{0.167em} \mathrm{d}z,&&\displaystyle \end{eqnarray} \tag{ 6a }$$

$$\begin{eqnarray} \displaystyle {N}_{\mathrm{p}}=b\sum _{k=1}^{N}\int _{{z}_{k}}^{{z}_{k+1}}e_{z x}^{(k)}E_{z}^{(k)}\hspace{0.167em} \mathrm{d}z,&&\displaystyle \end{eqnarray} \tag{ 6b }$$

$$\begin{eqnarray} \displaystyle {M}_{\mathrm{p}}=b\sum _{k=1}^{N}\int _{{z}_{k}}^{{z}_{k+1}}e_{z x}^{(k)}E_{z}^{(k)}z \hspace{0.167em} \mathrm{d}z,&&\displaystyle \end{eqnarray} \tag{ 6c }$$

where Q^(k) is the first element of the in-plane kth layer stress-reduced stiffness matrix, which corresponds for isotropic material to the Young's modulus. N_p and M_p are respectively the in-plane force and bending moment due to the piezoelectric coupling. The position of the kth layer is indicated by z_k. By substituting equations (1) and relations (3) in (4) and using (6) yields

$$\begin{eqnarray} \displaystyle \fl \prod =\frac{1}{2}\int _{0}^{L}A\mathop{(\varepsilon _{x x}^{0})}\nolimits ^{2}\hspace{0.167em} \mathrm{d}x+\int _{0}^{L}B\varepsilon _{x x}^{0}\varepsilon _{x x}^{1}\hspace{0.167em} \mathrm{d}x&&\displaystyle \\ \displaystyle +~\frac{1}{2}\int \nolimits _{0}^{L}D\mathop{({\varepsilon }_{x x}^{1})}\nolimits ^{2}\hspace{0.167em} \mathrm{d}x-\int \nolimits _{0}^{L}{N}_{\mathrm{p}}{\varepsilon }_{x x}^{0}\hspace{0.167em} \mathrm{d}x&&\displaystyle \\ \displaystyle -~\int \nolimits _{0}^{L}{M}_{\mathrm{p}}{\varepsilon }_{x x}^{1}\hspace{0.167em} \mathrm{d}x-b{h}_{\mathrm{p}}\int \nolimits _{0}^{L}{\varepsilon }_{z z}^{\mathrm{S}}{E}_{z}^{2}\hspace{0.167em} \mathrm{d}x-W.&&\displaystyle \end{eqnarray} \tag{ 7 }$$

The terms containing B vanish in the case of a symmetric laminate, as is supposed here. Note that stress-stiffening effects result from external compressive forces and residual piezoelectric stress ${N}_{\mathrm{r}}=A{\varepsilon }_{x x}^{0}-{N}_{\mathrm{p}}$. For slender beams, the longitudinal motion can also be neglected compared to the transversal one. As a result ${\varepsilon }_{x x}^{0}$ can be approximated as follows:

$$\begin{eqnarray} &&{\varepsilon }_{x x}^{0}=u ^{\prime}+\frac{1}{2}{~}{w ^{\prime}}^{2}=\frac{1}{2 L}\int \nolimits _{\mathrm{L}}{~}{w ^{\prime}}^{2}\hspace{0.167em} \mathrm{d}x. \end{eqnarray} \tag{ 8 }$$

Evaluating expressions (6), considering two piezoelectric layers, gives

$$\begin{eqnarray} \displaystyle \fl A={E}_{\mathrm{s}}{A}_{\mathrm{s}}+2{E}_{\mathrm{p}}{A}_{\mathrm{p}},&&\displaystyle \end{eqnarray} \tag{ 9a }$$

$$\begin{eqnarray} \displaystyle \fl D=\frac{b{h}_{\mathrm{s}}^{3}{E}_{\mathrm{s}}}{12}+\frac{2}{12}b{h}_{\mathrm{p}}(4{h}_{\mathrm{p}}^{2}+6{h}_{\mathrm{p}}{h}_{\mathrm{s}}+~3{h}_{\mathrm{s}}^{2}){E}_{\mathrm{p}}&&\displaystyle \\ \displaystyle ={E}_{\mathrm{s}}{I}_{\mathrm{s}}+2{E}_{\mathrm{p}}{I}_{\mathrm{p}},&&\displaystyle \end{eqnarray} \tag{ 9b }$$

$$\begin{eqnarray} \displaystyle \fl {N}_{\mathrm{p}}=-b{e}_{z x}\dot {\lambda },&&\displaystyle \end{eqnarray} \tag{ 9c }$$

$$\begin{eqnarray} \displaystyle \fl {M}_{\mathrm{p}}=-b{e}_{z x}({h}_{\mathrm{s}}+{h}_{\mathrm{p}})\dot {\lambda }/2&&\displaystyle \end{eqnarray} \tag{ 9d }$$

where I_s and I_p are the second moment of inertia of the steel shim and piezoelectric layers, while the electric field is expressed in terms of a time derivative of the flux linkage as ${E}_{z}=-\dot {\lambda }/2{h}_{\mathrm{p}}$ in case active layers are connected in series. The Lagrangian functional of the buckled piezoelectric beam results in

$$\begin{eqnarray} &&\mathfrak{L}(\dot {w},w ^{\prime},w^{\prime\prime},\dot {\lambda })=K-\Pi . \end{eqnarray} \tag{ 10 }$$

Before proceeding with the calculus of the Euler–Lagrange governing equations, it is worth making some assumptions on the trial functions of vibration mode shapes. We are interested in the dynamic response of the system around the first buckling mode, which is supposed to be dominant. The total deflection with respect to the clamping base can be defined as

$$\begin{eqnarray} &&w(x,t)={w}_{1}(x)+v(x,t) \end{eqnarray} \tag{ 11 }$$

where w₁ = h₀φ₁(x) is the initial equilibrium shape of the beam and h₀ = w(L/2,0) being the height of the central point at t = 0, whereas v(x,t) is the time-dependent deflection around the initial shape. According to the Galerkin discretization, v(x,t) can be expanded into a superposition of N orthonormal basis functions φ_i(x) as

$$\begin{eqnarray} &&v(x,t)=\sum _{i=1}^{N}{r}_{i}(t){\varphi }_{i}(x) \end{eqnarray} \tag{ 12 }$$

where r_i(t) are the time-dependent expansion coefficients of the modal configurations. Hence, by replacing (12) in equation (11) and keeping only the first mode, the total deflection results in

$$\begin{eqnarray} &&w(x,t)={h}_{0}{\varphi }_{1}(x)+{r}_{1}(t){\varphi }_{1}(x) \end{eqnarray} \tag{ 13 }$$

where r₁(t) is the time-dependent displacement of the middle point of the clamped–clamped beam around the initial buckled shape. Discontinuity of the piezoelectric beam can be treated by using stepwise functions with the continuity and compatibility conditions at the boundaries of each uniform beam section. The boundary conditions of a clamped–clamped beam are w(0,t) = 0, w(L,0) = 0, ∂w(0,t)/∂x = 0 and ∂w(L,t)/∂x = 0 and mode shape functions are taken from prior works [28, 30, 31].

In the following, only the first mode will be considered. Therefore, injecting (13) into the Lagrangian functional (10) gives

$$\begin{eqnarray} \displaystyle \fl \mathfrak{L}(q,\dot {q},\dot {\lambda })=\frac{1}{2}m({\dot {q}}^{2}+{\dot {y}}^{2})+\eta \dot {q}\dot {y}-\frac{1}{4}{k}_{3}{q}^{4}&&\displaystyle \\ \displaystyle -~\frac{1}{2}({k}_{2}+{k}_{1}\dot {\lambda }){q}^{2}+{k}_{0}\dot {\lambda }q+\frac{1}{4}{C}_{\mathrm{p}}{\dot {\lambda }}^{2}-P{d}_{\mathrm{b}},&&\displaystyle \end{eqnarray} \tag{ 14 }$$

where the superscript dot denotes the time derivative and the sixth term is rewritten in terms of equivalent capacitance ${C}_{\mathrm{p}}=b{\varepsilon }_{z z}^{\mathrm{s}}{L}_{\mathrm{p}}/{h}_{\mathrm{p}}$ of the piezoelectric beam, qhile m, η, k₀, k₁, k₂ and k₃ are the effective mass, inertia force factor, piezoelectric coupling factor, in-plane piezoelectric force factor, linear and nonlinear stiffness of the beam, respectively. These parameters are evaluated as follows:

$$\begin{eqnarray} \displaystyle m=L({\rho }_{\mathrm{s}}{A}_{\mathrm{s}}+2{\rho }_{\mathrm{p}}{A}_{\mathrm{p}})+{M}_{0},&&\displaystyle \end{eqnarray} \tag{ 15a }$$

$$\begin{eqnarray} \displaystyle \eta =({\rho }_{\mathrm{s}}{A}_{\mathrm{s}}+2{\rho }_{\mathrm{p}}{A}_{\mathrm{p}})\int \nolimits _{0}^{L}\psi \hspace{0.167em} \mathrm{d}x+{M}_{0}\psi (L/2)&&\displaystyle \end{eqnarray} \tag{ 15b }$$

$$\begin{eqnarray} \displaystyle {k}_{0}=\frac{1}{2}b{e}_{z x}({h}_{\mathrm{s}}+{h}_{\mathrm{p}})\int \nolimits _{0}^{L}\psi ^{\prime\prime}\hspace{0.167em} \mathrm{d}x,&&\displaystyle \end{eqnarray} \tag{ 15c }$$

$$\begin{eqnarray} \displaystyle {k}_{1}=b{e}_{z x}\int \nolimits _{0}^{L}\mathop{(\psi ^{\prime})}\nolimits ^{2}\hspace{0.167em} \mathrm{d}x,&&\displaystyle \end{eqnarray} \tag{ 15d }$$

$$\begin{eqnarray} \displaystyle {k}_{2}=D\int \nolimits _{0}^{L}\mathop{(\psi ^{\prime\prime})}\nolimits ^{2}\hspace{0.167em} \mathrm{d}x-P\int \nolimits _{0}^{L}\mathop{(\psi ^{\prime})}\nolimits ^{2}\hspace{0.167em} \mathrm{d}x,&&\displaystyle \end{eqnarray} \tag{ 15e }$$

$$\begin{eqnarray} \displaystyle {k}_{3}=\frac{A}{2 L}{\left(\int \nolimits _{0}^{L}\mathop{(\psi ^{\prime})}\nolimits ^{2}\hspace{0.167em} \mathrm{d}x\right)}^{2},&&\displaystyle \end{eqnarray} \tag{ 15f }$$

in which the time-dependent first mode is expressed as (h₀ + r₁(t))φ₁ = q(t)ψ, where the initial buckling shape function results in ψ(x) = (1 − cos(2πx/L))/2.

At this stage, single-degree-of-freedom motion equations for vertical deflection of the beam midpoint can be evaluated by applying Euler–Lagrange equations to (14):

$$\begin{eqnarray} \begin{array}{@{}rc@{}} \displaystyle &\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial \mathfrak{L}}{\partial \dot {q}}\right)-\frac{\partial \mathfrak{L}}{\partial q}=F(t), \\ \displaystyle &\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}\left(\frac{\partial \mathfrak{L}}{\partial \dot {\lambda }}\right)-\frac{\partial \mathfrak{L}}{\partial \lambda }=I(t) \end{array} \end{eqnarray} \tag{ 16 }$$

where F(t) and I(t) respectively represent the generalized force and current within the energy harvesting circuit. Equivalently, the current will be related to effective electrical load R_L by Ohm's law $I(t)=-\dot {\lambda }/{R}_{\mathrm{L}}$. Hence, computing (16) gives two second-order nonlinear differential equations that link the motion of the midpoint to the output voltage across the resistive load:

$$\begin{eqnarray} \displaystyle m\ddot{q}+c\dot {q}+{k}_{3}{q}^{3}+({k}_{2}+{k}_{1}\dot {\lambda })q-{k}_{0}\dot {\lambda }=-\eta \ddot{y},&&\displaystyle \end{eqnarray} \tag{ 17a }$$

$$\begin{eqnarray} \displaystyle \frac{1}{2}{C}_{\mathrm{p}}\ddot{\lambda }+\frac{\dot {\lambda }}{{R}_{\mathrm{L}}}={k}_{1}q\dot {q}-{k}_{0}\dot {q}.&&\displaystyle \end{eqnarray} \tag{ 17b }$$

The viscous damping is accounted for by the second term in the first equation. Here, it can be noted that time-dependent variables are not completely separated. In the case k₁ was zero, the equations reduce to those of the nonlinear piezoelectric cantilever seen in [22, 17, 19]. In addition, the conservative force is a Duffing-like type, but the linear stiffness parameter is a function of the output voltage. As a continuous piezoelectric layer is assumed over the whole length of the beam, the stiffness parameter k₀ should vanish. However, in our experimental model, the electrodes which capture the charge are of length L_e < L and are connected only from one side of the beam. Thus, the parameters k₀ and k₁ are evaluated from zero to L_e. Accordingly, all the coefficients (15) are then calculated as

$$\begin{eqnarray} \displaystyle {k}_{0}=b{e}_{z x}({h}_{\mathrm{s}}+{h}_{\mathrm{p}})\frac{\pi }{2 L}\sin \left(\frac{4 \pi {L}_{\mathrm{e}}}{L}\right),&&\displaystyle \end{eqnarray} \tag{ 18a }$$

$$\begin{eqnarray} \displaystyle {k}_{1}=b{e}_{z x}\frac{\pi }{8{L}^{2}}\left(4 \pi {L}_{\mathrm{e}}-L\sin \left(\frac{4 \pi {L}_{\mathrm{e}}}{L}\right)\right),&&\displaystyle \end{eqnarray} \tag{ 18b }$$

$$\begin{eqnarray} \displaystyle {k}_{2}=D\left(\frac{2{\pi }^{4}}{{L}^{3}}\right)-P\left(\frac{{\pi }^{2}}{2 L}\right),&&\displaystyle \end{eqnarray} \tag{ 18c }$$

$$\begin{eqnarray} \displaystyle {k}_{3}=\frac{A{\pi }^{4}}{8{L}^{3}}.&&\displaystyle \end{eqnarray} \tag{ 18d }$$

In this paper, our intent is to solve the dynamical equations with stochastic base excitation for an increasing contraction length ΔL. For this purpose, the compression load P, which is unknown in the experiment, can be expressed in terms of ΔL, using the following relationships:

$$\begin{eqnarray} \displaystyle P={P}_{\mathrm{cr}}+\frac{A}{2 L}\int \nolimits _{0}^{L}\mathop{({h}_{0}\psi ^{\prime})}\nolimits ^{2}\hspace{0.167em} \mathrm{d}x&&\displaystyle \end{eqnarray} \tag{ 19a }$$

$$\begin{eqnarray} \displaystyle {P}_{\mathrm{cr}}=4{\pi }^{2}D/{L}^{2},&&\displaystyle \end{eqnarray} \tag{ 19b }$$

$$\begin{eqnarray} \displaystyle {h}_{0}^{2}=4 L(\Delta L-{d}_{\mathrm{b}})/{\pi }^{2},&&\displaystyle \end{eqnarray} \tag{ 19c }$$

$$\begin{eqnarray} \displaystyle {d}_{\mathrm{b}}={P}_{\mathrm{cr}}L/A,&&\displaystyle \end{eqnarray} \tag{ 19d }$$

and, in turn, ΔL can be expressed as initial midpoint height h₀ and buckling displacement d_b. This relation is obviously valid only for the post-buckling situation. The critical load for such a clamped–clamped beam is retrieved by Reddy [28]. In the following sections we will make use of numerical methods to solve the motion equations in the post-buckling configuration (17) in order to compare them with experimental results.

In this section, we consider the case of an analytical system similar to the experimental device, in order to investigate the main dynamical behaviour under random noise excitation. However, for the real system, some corrections to geometrical parameters will be taken into account, as described in section 3. Table 1 summarizes the main model parameters used for simulations and experiments. As a first instance, we can consider the situation at open circuit. Thus, the effective potential relative to the midpoint vertical coordinate q reduces to a quartic function given by

$$\begin{eqnarray} &&\Pi (q)=\frac{1}{4}{k}_{3}{q}^{4}+\frac{1}{2}{k}_{2}{q}^{2}+{P}_{\mathrm{cr}}{d}_{\mathrm{b}}. \end{eqnarray} \tag{ 20 }$$

When the beam is compressed by increasing ΔL or, equivalently, the lateral load P, the stiffness parameter k₂ becomes negative. As a consequence, the potential energy varies from being a nearly flatted parabola to a double-well function. So that, as expected, the systems behaves as a Duffing oscillator for k₂ ≥ 0, which corresponds to a small compression (P ≤ P_cr) or even a tension P < 0, and as a Duffing bistable oscillator for k₂ < 0 (P > P_cr). Clearly, a supercritical bifurcation occurs when P − P_cr (or equivalently k₂) changes sign, as can be seen in figure 2.

Zoom In Zoom Out Reset image size

Table 1.  Model parameters used for numerical and experimental tests.

Steel shim		Piezoelectric layer
Parameter, symbol	Value	Parameter, symbol	Value
Material	Steel	Material	PZT film
Length, L0	55 mm	Length of layers, Lp	55 mm
Width, b	11 mm	Length of electrodes, Le	21 mm
Thickness, hs	0.1 mm	Width, bp = b	10 mm
Density, ρs	7850 kg m−3	Thickness, hp	0.08 mm
Young's modulus, Es	2.03 × 1011 Pa	Density, ρp	4000 kg m−3
Damping ratio, ζ	0.018 (unb)	Young's modulus, Ep	4 × 1010 Pa
	0.025 (buck ΔL/L = 0.4%)	Coupling coefficient, d31	−10 pC N−1
		Electrical permittivity, ${\varepsilon }_{z z}^{\mathrm{S}}$	100ε0
		Vacuum permittivity, ε0	8.854 × 10−12 F m−1

For k₂ ≥ 0, the potential has an absolute minimum at q = 0 around which the system oscillates with a natural frequency ${\omega }_{0}^{2}=\mathop{(1/m)\Pi ^{\prime\prime}\vert }\nolimits _{q=0}={k}_{2}/m$ in the undamped case. Whereas, for k₂ < 0, the quartic potential has two local minima corresponding to stability positions q₋ =− (|k₂|/k₃)^−1/2 and q₊ = (|k₂|/k₃)^−1/2. The barrier height h_b between the two wells is equivalent to ${k}_{2}^{2}/2{k}_{3}$. Depending on it and on the vibration noise strength, the system tends to oscillate around either minima and, from time to time, it jumps across the barrier following Kramer's transition rate [32].

Actually, in the dynamical regime, the stiffness coefficient k₂ depends on the output voltage. Therefore, the potential energy varies in time. Moreover, in the real experiment, the beam has only one piezoelectric layer, which means that the term B in expression equation (7) is nonzero, so making the potential energy asymmetrical. Such an asymmetry, though, can be considered as not important for a qualitative comparison of the power performance.

The equipment used for experimental tests is shown in figure 3(a). The piezoelectric beam is clamped on both its ends as shown in figure 3(b). One of the two clamps is kept fixed, while the other one was manually moved by a micrometric stage with spatial resolution of 0.01 mm. The piezoelectric sample was manufactured by Ferrari's research group at the University of Brescia (Italy) by screen-printing low-curing-temperature lead zirconate titanate (PZT) films on steel lamina. The main geometric and electrical properties are listed in table 1. However, some differences in the real sample with respect to the numerical model must be specified. The experimental device is made of a steel shim upon which four piezoelectric layers are printed using a PZT milled powder (Piezokeramica APC-856) and low temperature polymeric vehicle-binder mixed in a ratio of 2:1 wt. They are polarized in opposing directions under a 4 MV m⁻¹ electric field perpendicular to the lamina. The total piezoelectric film thickness is 0.08 mm. Hence, the testing system results in a unimorph asymmetrical beam. Measured electrical permittivity and density, after the printing process, are respectively equal to ${\varepsilon }_{z z}^{\mathrm{S}}=100$ (±10) and ρ_p = 4000 kg m⁻³. The piezoelectric film is also discontinuous in the middle, where a 10 mm gap is depleted to leave room for an additional mass. Thus, the actual length of each piezoelectric strip is 21 mm. The effective capacitance of each piezoelectric side was measured with an RLC meter resulting in around C_p = 1.15 ± 0.05 nF, in very good agreement with the calculated value of 1.16 nF. Two electrodes were screen-printed and fired at low temperature using a conductive silver-based polymeric ink. Electrical contacts are clamped on only one side of the bridge connecting a variable resistance (ranging from 10 Ω to 1 MΩ). Two laser displacement sensors (Keyence LK-501) are used in differential mode to detect the position of the beam midpoint relative to the vibrating base. An accelerometer (PCB Piezotronics 338A35) is also used in determining the base acceleration level in order to minimize inaccuracies. Input vibrations are generated by using a signal generator (Agilent 33220A) and a National Instruments DAQ USB card connected to the power amplifier of the electromagnetic shaker (TIRA). In addition, the input vibration spectrum is monitored by a 1 GHz oscilloscope (Le Croy LC574A). The voltage drop across the resistance together with all other measured signals are captured by a computer through a NI DAQ card with sampling frequency of 2 kHz.

Zoom In Zoom Out Reset image size

In the experimental system, a small mass of 1.47 g was added at the centre of the beam in order to improve dynamic flexibility. For the theoretical model, this variation implies an additional lumped mass to be included in the kinetic energy term. However, deviations from the theoretical model were expected due to the asymmetry of the structure, parasite capacitances of the electrodes and a not constant damping. Natural frequencies were determined by first exciting the beam with white noise at open circuit. The uncompressed configuration shows a resonance f₀ = 215 Hz, while the predicted value is f₀ = 220 Hz. The damping ratio was estimated by using both the decay rate and bandwidth method measuring the resonances at open and short circuit, under repeating impulsive vibration input. It resulted in ζ = 0.018 for the unbuckled configuration up to ζ = 0.025 in the buckled configuration at ΔL/L = 0.4%.

The experimental tests were mainly performed under exponentially correlated Gaussian noise vibration with zero mean, autocorrelation time of τ = 0.001 s and root-mean-squared (rms) acceleration from σ = 1 to 3g (where g stands for the gravitational acceleration: 9.81 m s⁻²).

Figure 4 shows the measured central point displacement of the beam and the output voltage compared to the prebuckled configuration with ΔL = 0.07 mm, at σ = 3g base acceleration. Such high vibration amplitude was required to obtain a good bistability of the oscillator. This is due to the relatively high damping. The plots of the left (a), (c) and right (b), (d) columns in figure 4 refer respectively to 1 MΩ and 25 kΩ of electrical load. The latter value was chosen as a reference value for the comparison, because it corresponds to the optimal electrical load for the maximum power in the unbuckled case. The enhancement in terms of rms displacement and output voltage is significant. The bistable configuration shows indeed an amplification factor of 8.46 times (for the load of 25 kΩ) up to 16.22 times (for 1 MΩ) in the rms voltage.

The bistable behaviour is evident in figure 4(a), where two stable minima are located around ±0.9 mm. The theoretical bimorph model gives slightly larger values q_± =± 1.2 mm. Nevertheless, the expected qualitative behaviour is not significantly influenced, as shown in the following.

Zoom In Zoom Out Reset image size

In effect, in figure 4(b), the deflection distribution clearly appears asymmetrical. This is due to the unimorph structure of the experimental laminate, which determines the asymmetry of the potential well. It can be noted that output voltage strongly decreases with decreasing load resistance, while the displacement is only slightly affected. This means that the shunt effect does not affect very much the bistable dynamics of the oscillator.

Figures 5(a) and (b) illustrate the comparison of the rms voltage and the average electrical power respectively versus load variation between linear and nonlinear configurations. In particular, the unbuckled beam is compared to two buckled configurations corresponding to ΔL = 0.05 and 0.07 mm of compression lengths (equivalent to 0.09% and 0.12% of the unbuckled beam length of 55 mm). In this first test, the increasing trends are qualitatively similar for both linear and nonlinear regimes but with a remarkable difference in the amplitude. The maximum rms voltage read from the linear beam is 0.018 V, whereas, it corresponds to 0.095 V and 0.29 V respectively for the first and second buckling cases. In figure 5(b), the average electrical power delivered to the load was calculated by the formula ${P}_{\mathrm{avg}}={V}_{\mathrm{rms}}^{2}/{R}_{\mathrm{L}}$. What is more important to stress is that the energy harvester in a nonlinear regime is capable of generating one order of magnitude more power than its resonant counterpart, at the same 3g base acceleration. In addition, the optimal load increases from 25 to 150 kΩ in the first buckled example, while for higher compression the power curve keeps increasing for the entire interval. Therefore, the nonlinear nature of buckled systems widens the impedance bandwidth so as to make this configuration more suitable for different operations like battery charging, sensing, radio frequency transmission and sleep mode.

Zoom In Zoom Out Reset image size

Figures 6(a) and (b) show the standard deviation of the midpoint deflection q(t) = w(L/2,t) over a time interval of ΔT = 20 s versus compression ratio percentage ΔL/L (%). All the experimental quantity are normalized by a constant and expressed in arbitrary units to be comparable with numerical results. Each graph shows three curves relative to different rms noise accelerations, equal to σ = 1–3g. In this test, the standard deviation of the central point deflection is chosen as an indicator to measure the dynamic dispersion around the average equilibrium position in the unbuckled case, which in our case corresponds approximately to zero. In fact, under Gaussian noise this represents a reliable parameter to quantify the strength of the output signal. The same approach has been used for the output voltage in place of the rms not to take into account accidental bias or dc artefacts. As the beam is compressed, the vibration amplitude of the central point increases, achieving a maximum. Close after this critical point, it suddenly falls down to smaller vibration amplitudes which become about half of the unbuckled system. This phenomenon is in good agreement with the bistable dynamics discussed in section 2.2 and previously explained by the authors [17]. Once the barrier height is too high to be climbed within the chosen time interval, the buckled oscillator is confined to vibrate around one of the two minima of the quartic potential well. Therefore, the vibration amplitude is reduced as q ∝ σ/ω_± which is due to the stiffening effect, wherein the intrawell resonance ${\omega }_{\pm }^{2}=(1/m)\Pi ^{\prime\prime}\vert _{q=0}$ increases. The enhancement of beam deflections with respect to the value at ΔL = 0 mm goes from a factor of 1.5 to 6.7 for the experimental device and from 1.4 to 2.5 for the simulated system. From a qualitative point of view, these ranges can be considered in good agreement. On the other hand, the quantitative mismatch is due to the approximations on the theoretical model which does not fit exactly the real system. Besides, the shaker response was nonlinear at high accelerations.

Zoom In Zoom Out Reset image size

Figures 6(c) and (d) describe the average electrical power versus the compression length ratio respectively for the experimental and numerical models. The average power is calculated here with the formula ${P}_{\mathrm{avg}}={V}_{\mathrm{std}}^{2}/{R}_{\mathrm{L}}$, where R_L = 1 MΩ is the resistive load and V_std is the standard deviation of the voltage through it. These results correspond to the aforementioned deflection dynamics. When the beam is compressed the average power increases, with maximum power amplification factors of 2.9, 13.1 and 3.7 at σ = 1 and 3g rms of base acceleration respectively in the experimental case. The numerical model instead gives gain factors of 1.6, 3.9 and 1.9. It is interesting to note that there exists an optimal value of acceleration corresponding to the maximum amplification factor: in this case the intermediate acceleration of 2g. This effect is also predicted by the simulations. Apart from the beam asymmetry, in the first place, the general discrepancies with the theoretical model can be accounted for by the electrode contact leakage of the testing device and, in the second place, to the unknown experimental value of the piezoelectric film elastic stiffness. This parameter was indeed taken equal to 4 × 10¹⁰ Pa by the typical Young's modulus range for thin piezoelectric films [33].

The effective stress–charge piezoelectric coupling was estimated by using the same method of Roundy and Wright [34]. It was found to be equal to e_zx =− 0.003 C m⁻² which is two orders of magnitude smaller than the theoretical value e_zx =− 0.4 C m⁻². In addition, the damping ratio was assumed to be constant, whereas in general it depends on the nonlinear tuning force, as also indicated in previous works [35].

Despite discrepancies, the trends are increasing in both cases and they do not exhibit critical breakdown points as in the deflection plots. A possible interpretation of such a phenomenology can be given as follows. Different from piezoelectric coupled equations of previous works [17, 36], the governing equations (17) include the voltage–velocity coupled terms with k₁: fifth on the left-hand side of the first equation and first on the right-hand side of the second equation. Being k₁/k₀ ∼ 10⁴ m⁻¹, the voltage–velocity term that is due to the in-plane piezoelectric force outweighs the piezoelectric coupling term due to the bending. This can be seen in figure 7 where numerical simulations were performed comparing the system response with zero in-plane piezoelectric force (k₁ = 0) and with nonzero force (k₁ < 0). As this term is proportional to the velocity, when the system passes over the critical bifurcation point, it seem to be responsible for keeping the output voltage high as intrawell vibrations occur at higher frequency. So the absence of a breakdown in the output voltage is qualitatively predicted by the model. Nevertheless, it must be noticed that the power gain results are more important in the experimental data than for simulations.

Zoom In Zoom Out Reset image size