Modeling and Experimental Validation of Bi-Stable Beam Based Piezoelectric Energy Harvester

One Degree-of-Freedom Electromechanical Lumped Parameter Model

The energy harvester we are modeling has a clamped-clamped beam structure of a stack of thin films including structural layer, seed layer, piezoelectric layer and passivation layer (Figure 1(a)). A heavy proof mass is concentrated at the middle of the beam to capture the external vibration and excites the whole beam to vibrate out-of-plane. Piezoelectric elements work in d₃₃ mode with top interdigitated electrodes, coupling the mechanical response with electrical domain. The multi-layer beam is designed to buckle by incorporating the controlled compressive residual stress of the micro-fabricated thin films. Statically, the beam either buckles up or down (two equilibriums), and the dynamics becomes complex when the system is continuously excited in post-buckling regime. To capture the essential dynamics phenomena of the complex system, the beam vibration mode has been assumed so that a one degree-of-freedom model can be constructed. The nonhomogeneous beam structure has been taken into account by considering the different thickness and material properties of each layer. Furthermore, residual stress of each layer is built in as part of the stiffness of the beam to induce buckling. The electrical and mechanical domains are coupled with linear and nonlinear coupling, so that the generated electrical power can be calculated.

Figure 1:

Schematics of the doubly clamped beam based energy harvester. (a) The piezoelectric element is in d₃₃ mode with top interdigitated electrodes. The piezoelectric element and electrode span symmetrically from –b*L to –a*L and a*L to b*L on the substrate. (b) The beam has n layers of different material thin films with various thicknesses.

The lumped parameter model is formulated by Lagrange’s method. In classical mechanics, the Lagrangian is defined as,

where KE is the kinetic energy of the system, PE is the potential of the system. In this energy harvester as in many other vibration energy harvesters, the proof mass is much heavier than the beam’s distributed mass, so that the kinetic energy of the system is simplified as that of the center-concentrated proof mass,

where m is the proof mass, w˙ is the time derivative of the beam center displacement, i. e. the velocity of the proof mass. To find out the thermodynamic potential of the system including the piezoelectric material, we start by considering the electrical enthalpy volume density:

and piezoelectric constitutive equations in d₃₃ mode (Standards Committee of the IEEE Ultrasonics and Ferroelectrics 1988),

where T3, S3, D3, E3 are the stress, strain, electric displacement and electric field in 3-direction respectively; c33E, e33 and ε33S are the elastic modulus, piezoelectric constant, and permittivity of the piezoelectric material; the superscripts E and S denote the parameters are at constant electric field and strain respectively. Substitute T3 and D3 into [3], and add the strain energy contributed by the residual stress T0, ∫0S3T0ds=T0S3:

The residual stress in thin films micro-fabricated on a substrate causes a change in the radius of curvature of the substrate, which can be measured. The stress can then be calculated from the Stoney’s formula:

where E_sub, υ_sub, h_sub are the Young’s modulus, Poisson’s ratio and thickness of the substrate, R is the substrate radius of curvature, and t_film is the thickness of the film. The Lagrangian of the system can now be evaluated by integrating the enthalpy density over the beam’s volume layer by layer,

where v_i is the volume of i-th layer and n is the total number of layers.

The strains developed in the beam need to be evaluated before carrying out the integrations in [9]. The total strain S_T developed in the beam has two components: bending strain, which changes linearly across the beam thickness, and axial strain due to large deflection,

where l is the beam length. The strain is calculated from the neutral axis of the beam (Figure 1(b)). It should be stressed that the beam composes multilayers with various elastic properties, the neutral axis therefore does not coincide with the mid-plane. The general formula for calculating the position of the neutral axis (distance from the bottom surface) of a n-layer beam is,

The beam vibrates up and down in z-axis, and by assuming that it vibrates predominantly in one mode, simplification can be made when evaluating the lumped parameters. The first buckling mode of the beam is adopted, which satisfies the boundary conditions of clamped-clamped beam and has been verified as the vibration mode shape at the largest deflection in experiment, by measuring the deflection on different spots along the beam using laser vibrometer (Figure 2). The deflection of the beam can then be separated into time and space,

Figure 2:

The comparison between measured maximum deflection versus the assumed sinusoidal function in [9]. The deflection was measured from single point laser vibrometer at various points on the vibrating beam. The assumed shape agrees well with the beam vibration shape.

where w(t) is the deflection of the beam center varying with time.

Lagrange equations are,

where ξi is the i-th independent generalized coordinate, QiForce and QiDissipation are the generalized external force and the generalized dissipative force. We choose the deflection of the mid-point of the beam w and the output voltage V as the generalized coordinates. The Lagrange equation with respect to the first coordinate w is then,

Evaluating the integrations in [8] and substituting into [13], the governing equation of the mechanical domain can be obtained and written in a compact form,

where k_L, k_N, b, c_L, c_N and F are the linear stiffness and nonlinear stiffness of the beam, the mechanical damping coefficient, the linear and nonlinear electromechanical coupling, and the external excitation force respectively. These lumped parameters are functions of the device dimensions and material properties, which are useful for device design,

where W, H, are the width and thickness, a and b denotes the span of the electrodes on the beam, since they do not cover the whole beam (Figure 1(a)), and g is the gap between two electrode fingers; the subscript p denotes the variable is associated to the piezoelectric layer. It should be noted that the linear stiffness has two parts: the first part is from bending of the beam and the second comes from the residual stress. More particularly, when the residual stress is negative (compressive) and large enough, the linear stiffness k_L will be negative, so that [14] becomes a characteristic bi-stable Duffing equation.

The second Lagrange equation with respect to the coordinate V is,

Taking time derivative of the equation gives the governing equation for the electrical domain,

where IL=cLw˙ and IN=cNww˙ are two parts of the electrical current generated by piezoelectric element through coupling, and L and N denote linear and nonlinear coupling respectively; the induced voltages on the electrical port are written in separate parts V_L and V_N, due to the fact that they come from two parts of the current respectively and have different frequencies due to different coupling, and this differentiation makes the assumptions on their function simple (next section). C0 is the internal capacitance of the piezoelectric element and is calculated as

Analytical Solution of the Frequency Response

The nonlinear governing equations [14] and [20] could be numerically integrated to obtain the solution in time domain, but analytical solutions provide more insights on the dynamic behavior; more importantly, the explicit relations between system parameters and the performance are significant for device design. Therefore, the heuristic harmonic balance method is adopted to approximate the frequency response analytically.

Bi-stable oscillator has characteristic double-well potential (Xu and Kim 2015). Depending on if the oscillator has enough energy to overcome the energy barrier, it crosses the well and has the so called inter-well oscillation, or stays in one well and oscillates intra-well. To differentiate the two modes of oscillations, we assumed that when the beam oscillates around the buckled configuration (intra-well), the beam is oscillating around a non-zero bias position δ from the flat position; while the bias δ is zero, the beam oscillates around the flat position and the oscillation becomes symmetric on both sides of the flat position (inter-well). Consider the intra-well case first: plug in the assumed functions of w to [14], and balance the terms with δ but without sinusoidal functions, we find δ, which is a function of the deflection amplitude,

By assuming the external force is sinusoidal and continuous, from trigonometric relations, it is easily found that the frequency of the induced electrical current from nonlinear coupling doubles that of the current from linear coupling. Physically, it is due to the developed stretching strain has half the cycle of the bending strain. The induced electrical current and voltage are then written in different parts with different frequencies shown in Table 1. The subscript 0 denotes amplitude. Unknown phase differences have been added to the sinusoidal functions and will be solved. Writing [20] into two separate equations, we can solve for V_L and V_N separately:

Multiply [23] and [24] by sinωt and cosωt, integrate from 0 to 2π and set both to zeros, we obtain six equations; and multiply [24] by sin2ωt and cos2ωt, integrate from 0 to π and set both to zeros, we obtain another four equations. Solving the ten equations, the amplitudes of voltages and the phase constants can be found,

It is interesting to note that the voltage due to nonlinear coupling is a function of the deflection amplitude squared, while the voltage due to linear coupling is proportional to the deflection amplitude. This indicates that when the deflection is beyond some point, the nonlinear response will dominate the total response. The deflection amplitude is not known yet and will be solved from the mechanical domain equation [14].

Substitute assumed functions listed in Table 1 to [14], multiply by sinωt, integrate from 0 to 2π, set to zero, and multiply cosωt, integrate from 0 to 2π, set to zero, we obtain two equations:

Take the square of both sides of [30] and [31] and add them together, one equation with one unknown w₀ is formed,

where f₁, f₂, f₃, f₄, f₅ are functions of the known system’s parameters:

The inter-well oscillation is symmetric with respect to the flat unbuckled position, and hence there is no bias in w(t). Solve [14], [23] and [24] in the same way as solving for the intra-well case, but with δ=0, we obtain the same equation as [32], but with different coefficients:

[32] is an algebraic equation and is straightforward to solve. Since the voltages are functions of w₀, the generated power can be calculated by assuming the harvester is connected to a resistor:

The w₀ in [32] are in even order and hence the third order equation gives three solutions. By solving the equation for intra-well and inter-well cases, we will obtain two sets of solutions, which correspond to two modes of oscillations (Figure 3). The jumps and tilt of the two curves show clearly the two modes of oscillations have characteristic softening and stiffening frequency responses. The softening response tilts to lower frequencies, which can be used for low-frequency energy harvesting at low g’s. When the input excitation is strong enough, the beam can snap between two buckled positions and the inter-well oscillation is triggered. The stiffening response has shifted up low-frequency response, which has large amplitude and increases the power output at low frequencies. The analytical solutions have been compared to the numerical integration results at 0.5 g and 3 g (Figure 4). The governing equations have been integrated numerically from 0 to 1 s with initial conditions as the system has the lowest energy (buckled position, zero velocity). The numerical amplitude has been extracted from the maximum and minimum after half a second (assumed steady). The two modes of oscillations and power enhancement have also been verified by experiment.

Figure 3:

Simulated frequency response of bi-stable oscillator. The blue and red lines are the solutions of intra-well and inter-well oscillations respectively. The simulation is based on the meso-scale model parameters shown in Figure 5 with acceleration amplitude of 1 g.

Figure 4:

Comparison between analytical solutions and numerical integration solutions. (a) and (b) compares the deflection and power in frequency domain at low excitation amplitude of 0.5 g so that the harvester oscillates intra-well. (c) and (d) compares the deflection and power in frequency domain at 3 g so that the harvester oscillates inter-well.