Bistable energy harvesting enhancement with an auxiliary linear oscillator

A diagram of the system under investigation is provided in figure 1. An auxiliary linear oscillator of mass m₂ is attached to an electromagnetic induction bistable energy harvester of mass m₁ which is subject to base acceleration $\ddot{z}(\tau )=-W\cos \Omega \tau$. The lumped-parameter modeling approach follows that of many researchers who have shown that the primary mechanical responses of a variety of continuously distributed bistable energy harvester systems may be accurately modeled as 1DOF devices, including beams buckled via axial loads [24], magnetic repulsion [26], or magnetic attraction [27, 28]. In practice, such assumptions may be violated for bistable harvester systems if higher-order buckling modes are induced either as a consequence of the specific interaction between the excitation and the bistable device or due to the location of attachment of the auxiliary linear oscillator on the bistable harvester. In this work, such scenarios are assumed to be avoided by design. Moreover, the lumped modeling method is particularly relevant in light of the present experimentation which utilizes lumped inertial masses for the bistable and inertial bodies. In the modeling, the bistable electromagnetic transduction mechanism is assumed to have a low inductance and its output shunted by a resistive load R, where the relative motion x(τ) between the base z(τ) and the bistable harvester mass m₁ induces a flow of current i(τ) in the harvesting circuit. The double-well restoring potential of the bistable device is

$$\begin{eqnarray} &&U(x)=-\frac{1}{2}{k}_{1}{x}^{2}+\frac{1}{4}{k}_{3}{x}^{4} \end{eqnarray} \tag{ 1 }$$

where position x = 0 is defined as the central, unstable equilibrium position of the double well, as illustrated in figure 1. The governing equations for the system may then be expressed as

$$\begin{eqnarray} &&{m}_{1}(\ddot{x}+\ddot{z})+{d}_{1}\dot {x}-{k}_{1}x+{k}_{3}{x}^{3}-{d}_{2}\dot {y}-{k}_{2}y+\psi \mathrm{i}=0 \end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray} &&{m}_{2}(\ddot{x}+\ddot{y}+\ddot{z})+{d}_{2}\dot {y}+{k}_{2}y=0 \end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray} &&\mathrm{i}R=\psi \dot {x} \end{eqnarray} \tag{ 4 }$$

where y is the relative displacement between the bistable inertial mass and the linear oscillator mass and ψ is an electromechanical coupling coefficient. For energy harvesters having an electrical time constant much smaller than the mechanical natural period, it is well known that the effects of external circuitry upon the mechanical response of the harvester are principally additional dissipative influences [1, 29–33], and this assumption is employed in constructing equation (4) from Kirchhoff's law. Following substitution and nondimensionalization, the governing equations may be represented as follows:

$$\begin{eqnarray} &&\fl {x}^{\prime\prime}+{\gamma }_{1}{x}^{\prime}-x+\beta {x}^{3}-\mu f{\gamma }_{2}{y}^{\prime}-\mu {f}^{2}y=p\cos \omega t \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} &&\fl {y}^{\prime\prime}+(1+\mu )f{\gamma }_{2}{y}^{\prime}+(1+\mu ){f}^{2}y-{\gamma }_{1}{x}^{\prime}+x-\beta {x}^{3}=0 \end{eqnarray} \tag{ 6 }$$

where the following constants are defined:

$$\begin{eqnarray} &&\fl {\omega }_{1}^{2}={k}_{1}/{m}_{1};\tqs {\omega }_{2}^{2}={k}_{2}/{m}_{2};\tqs \mu ={m}_{2}/{m}_{1};\nonumber\\ &&f={\omega }_{2}/{\omega }_{1};\tqs \omega =\Omega /{\omega }_{1};\tqs p=W/{\omega }_{1}^{2} \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray} &&\fl {\gamma }_{1}=({d}_{1}+{\varepsilon }^{2})/{m}_{1}{\omega }_{1};\tqs \gamma ={d}_{1}/{m}_{1}{\omega }_{1};\nonumber\\ &&{\varepsilon }^{2}={\psi }^{2}/R;\tqs {\gamma }_{2}={d}_{2}/{m}_{2}{\omega }_{2};\tqs \beta ={k}_{3}/{k}_{1};\nonumber\\ &&t={\omega }_{1}\tau \end{eqnarray} \tag{ 8 }$$

and ()' indicates a derivative with respect to nondimensional time t.

The harmonic balance method is employed to solve the governing equations (5) and (6). For most basic response prediction, one-term solutions may be assumed of the form

$$\begin{eqnarray} x(t)&=&{c}_{1}(t)+{a}_{1}(t)\sin \omega t+{b}_{1}(t)\cos \omega t \end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray} y(t)&=&{a}_{2}(t)\sin \omega t+{b}_{2}(t)\cos \omega t \end{eqnarray} \tag{ 10 }$$

where c₁ is required due to the possibility that the bistable oscillator may vibrate around a non-zero stable equilibrium. The attached linear oscillator does not have a constant term in the expansion due to the notation that y represents a relative coordinate. The fidelity of results obtained via the harmonic balance method is limited by the number of harmonic terms used in the Fourier series expansion. The approach has been historically shown to be qualitatively and quantitatively accurate for both low-energy intrawell and interwell snap through dynamics [34] and has recently been revisited as an efficient analytical tool for bistable energy harvesting systems [25, 35, 36]. Finally, when applying the results of a steady-state analysis towards a realistic application, one is advised to evaluate the nearness of operating conditions to potential bifurcations which indicate the possibility of co-existing responses including aperiodic dynamics; a variety of useful analytical tools are available for these assessments [37–39].

We substitute equations (9) and (10) into (5) and (6), eliminate higher-order terms, and group the constant, sinωt, and cosωt terms to yield five governing equations for the slow-varying coefficients c₁,a₁,b₁, a₂, and b₂,

$$\begin{eqnarray} &&\fl -{\gamma }_{1}{c}_{1}^{\prime}={\Lambda }_{c}{c}_{1} \end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray} &&\fl -{\gamma }_{1}{a}_{1}^{\prime}+2 \omega {b}_{1}^{\prime}+\mu f{\gamma }_{2}{a}_{2}^{\prime}\nonumber\\ &&=(\Lambda -{\omega }^{2}){a}_{1}-\omega {\gamma }_{1}{b}_{1}-\mu {f}^{2}{a}_{2}+\mu f \omega {\gamma }_{2}{b}_{2} \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray} &&\fl -2 \omega {a}_{1}^{\prime}-{\gamma }_{1}{b}_{1}^{\prime}+\mu f{\gamma }_{2}{b}_{2}^{\prime}\nonumber\\ &&=\omega {\gamma }_{1}{a}_{1}+(\Lambda -{\omega }^{2}){b}_{1}-\mu f \omega {\gamma }_{2}{a}_{2}-\mu {f}^{2}{b}_{2}-p \end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray} &&\fl {\gamma }_{1}{a}_{1}^{\prime}-(1+\mu )f{\gamma }_{2}{a}_{2}^{\prime}+2 \omega {b}_{2}^{\prime}\nonumber\\ &&=-\Lambda {a}_{1}+\omega {\gamma }_{1}{b}_{1}+\Sigma {a}_{2}-\sigma {b}_{2} \end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray} &&\fl {\gamma }_{1}{b}_{1}^{\prime}-2 \omega {a}_{2}^{\prime}-(1+\mu )f{\gamma }_{2}{b}_{2}^{\prime}\nonumber\\ &&=-\omega {\gamma }_{1}{a}_{1}-\Lambda {b}_{1}+\sigma {a}_{2}+\Sigma {b}_{2} \end{eqnarray} \tag{ 15 }$$

where the following are defined:

$$\begin{eqnarray} {\Lambda }_{c}&=&-1+\beta ({c}_{1}^{2}+\frac{3}{2}{r}_{1}^{2}) \end{eqnarray} \tag{ 16 }$$

$$\begin{eqnarray} \Lambda &=&-1+\beta (3{c}_{1}^{2}+\frac{3}{4}{r}_{1}^{2}) \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} \Sigma &=&(1+\mu ){f}^{2}-{\omega }^{2} \end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray} \sigma &=&(1+\mu )f \omega {\gamma }_{2} \end{eqnarray} \tag{ 19 }$$

$$\begin{eqnarray} {r}_{1}^{2}&=&{a}_{1}^{2}+{b}_{1}^{2} \end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray} {r}_{2}^{2}&=&{a}_{2}^{2}+{b}_{2}^{2}. \end{eqnarray} \tag{ 21 }$$

The steady-state responses of the system are determined by solving the coupled equations (11)–(15). As such, equations (14) and (15) may be solved explicitly in terms of a₁ and b₁ to determine constants a₂ and b₂,

$$\begin{eqnarray} {a}_{2}&=&\frac{\Sigma \Lambda +\sigma \omega {\gamma }_{1}}{{\Sigma }^{2}+{\sigma }^{2}}{a}_{1}+\frac{\sigma \Lambda -\Sigma \omega {\gamma }_{1}}{{\Sigma }^{2}+{\sigma }^{2}}{b}_{1} \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray} {b}_{2}&=&-\frac{\sigma \Lambda -\Sigma \omega {\gamma }_{1}}{{\Sigma }^{2}+{\sigma }^{2}}{a}_{1}+\frac{\Sigma \Lambda +\sigma \omega {\gamma }_{1}}{{\Sigma }^{2}+{\sigma }^{2}}{b}_{1}. \end{eqnarray} \tag{ 23 }$$

Equations (22) and (23) are substituted into (12) and (13) which are thereafter squared and summed to yield one equation for the bistable harvester response amplitude squared, ${r}_{1}^{2}$,

$$\begin{eqnarray} &&{p}^{2}={r}_{1}^{2}[(\kappa \Lambda +\delta \omega {\gamma }_{1}-{\omega }^{2})^{2}+(\delta \Lambda -\kappa \omega {\gamma }_{1})^{2}] \end{eqnarray} \tag{ 24 }$$

given the following terms:

$$\begin{eqnarray} \delta &=&\omega {\gamma }_{2}\Sigma \Theta -f \sigma \Theta \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} \kappa &=&1-f \Sigma \Theta -\omega {\gamma }_{2}\sigma \Theta \end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray} \Theta &=&\frac{\mu f}{{\Sigma }^{2}+{\sigma }^{2}}. \end{eqnarray} \tag{ 27 }$$

Given that the term Λ contains unknowns c₁ and ${r}_{1}^{2}$, equation (11) is first solved under steady-state conditions and one choice of the parameter c₁ is selected, depending on whether one is interested in intrawell, low-energy oscillations of the bistable oscillator (c₁ ≠ 0) or interwell, high-energy vibration (c₁ = 0). Then, c₁ is substituted into Λ such that equation (24) is expressed fully in terms of a cubic polynomial of the parameter ${r}_{1}^{2}$. The roots of equation (24) are then determined and the final coefficients are computed explicitly from

$$\begin{eqnarray} {a}_{1}&=&\frac{-\delta \Lambda +\kappa \omega {\gamma }_{1}}{(\kappa \Lambda +\delta \omega {\gamma }_{1}-{\omega }^{2})^{2}+(\delta \Lambda -\kappa \omega {\gamma }_{1})^{2}}p \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} {b}_{1}&=&\frac{\kappa \Lambda +\delta \omega {\gamma }_{1}-{\omega }^{2}}{(\kappa \Lambda +\delta \omega {\gamma }_{1}-{\omega }^{2})^{2}+(\delta \Lambda -\kappa \omega {\gamma }_{1})^{2}}p. \end{eqnarray} \tag{ 29 }$$

Having computed all five of the coefficients, the stability of the solutions ${r}_{1}^{2}$ may be determined via a perturbation method, for example by Jacobian analysis [40].

This work also focuses on the phase relationships amongst the dynamic elements. Thus, the system response may be expressed alternatively from equations (9) and (10) as follows:

$$\begin{eqnarray} x&=&{c}_{1}+{r}_{1}\cos \left (\omega t-{\phi }_{1}\right ) \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} y&=&{r}_{2}\cos \left (\omega t-{\phi }_{2}\right ) \end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray} \tan {\phi }_{1}&=&{a}_{1}/{b}_{1};\tqs \tan {\phi }_{2}={a}_{2}/{b}_{2} \end{eqnarray} \tag{ 32 }$$

where ϕ₁ and ϕ₂ represent phase lags between the input excitation and the bistable harvester and linear oscillator displacements, respectively.

The advantage of the one-term solution is the direct analytic expression, equation (24), used to determine the responses. However, the limited number of terms assumed in the Fourier expansion constrict the wider applicability of the results to realistic observations of nonlinear systems coupled to other linear or nonlinear DOFs which are known to more readily exhibit multi-harmonic dynamics [20–22]. As a result, the assumed solution form is now expanded with terms proportional to the fundamental and third harmonics, which is the result of the cubic nonlinearity assumption [22] and recent numerical and experimental observations of the third superharmonic in bistable energy harvesting studies [23–25].

$$\begin{eqnarray} x(t)&=&{c}_{1}(t)+{a}_{1}(t)\sin \omega t+{b}_{1}(t)\cos \omega t\nonumber\\ &&+~{a}_{3}(t)\sin 3 \omega t+{b}_{3}(t)\cos 3 \omega t \end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray} y(t)&=&{g}_{1}(t)\sin \omega t+{h}_{1}(t)\cos \omega t+{g}_{3}(t)\sin 3 \omega t\nonumber\\ &&+~{h}_{3}(t)\cos 3 \omega t. \end{eqnarray} \tag{ 34 }$$

Equations (33) and (34) are substituted into (5) and (6) and the coefficients of constant, sinωt,cosωt, sin3ωt, and cos3ωt terms are collected.

$$\begin{eqnarray} &&\fl -{\gamma }_{1}{c}_{1}^{\prime}={\Lambda }_{s}{c}_{1} \end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray} &&\fl -{\gamma }_{1}{a}_{1}^{\prime}+2 \omega {b}_{1}^{\prime}+\mu f{\gamma }_{2}{g}_{1}^{\prime}\nonumber\\ &&=({\Lambda }_{1}-{\omega }^{2}){a}_{1}-{\gamma }_{1}\omega {b}_{1}+\frac{3}{4}\beta (-{a}_{1}^{2}+{b}_{1}^{2}){a}_{3}\nonumber\\ &&-~\frac{3}{2}\beta ({a}_{1}{b}_{1}){b}_{3}-\mu {f}^{2}{g}_{1}+\mu f{\gamma }_{2}\omega {h}_{1} \end{eqnarray} \tag{ 36 }$$

$$\begin{eqnarray} &&\fl -2 \omega {a}_{1}^{\prime}-{\gamma }_{1}{b}_{1}^{\prime}+\mu f{\gamma }_{2}{h}_{1}^{\prime}\nonumber\\ &&={\gamma }_{1}\omega {a}_{1}+({\Lambda }_{1}-{\omega }^{2}){b}_{1}+\frac{3}{2}\beta ({a}_{1}{b}_{1}){a}_{3}\nonumber\\ &&+~\frac{3}{4}\beta (-{a}_{1}^{2}+{b}_{1}^{2}){b}_{3}-\mu f{\gamma }_{2}\omega {g}_{1}-\mu {f}^{2}{h}_{1}-p \end{eqnarray} \tag{ 37 }$$

$$\begin{eqnarray} &&\fl -{\gamma }_{1}{a}_{3}^{\prime}+6 \omega {b}_{3}^{\prime}+\mu f{\gamma }_{2}{g}_{3}^{\prime}\nonumber\\ &&=\frac{1}{4}\beta (-{a}_{1}^{2}+3{b}_{1}^{2}){a}_{1}+({\Lambda }_{3}-9{\omega }^{2}){a}_{3}-3{\gamma }_{1}\omega {b}_{3}\nonumber\\ &&-~\mu {f}^{2}{g}_{3}+3 \mu f{\gamma }_{2}\omega {h}_{2} \end{eqnarray} \tag{ 38 }$$

$$\begin{eqnarray} &&\fl -6 \omega {a}_{3}^{\prime}-{\gamma }_{1}{b}_{3}^{\prime}+\mu f{\gamma }_{2}{h}_{3}^{\prime}\nonumber\\ &&=\frac{1}{4}\beta (-3{a}_{1}^{2}+{b}_{1}^{2}){b}_{1}+3{\gamma }_{1}\omega {a}_{3}+({\Lambda }_{3}-9{\omega }^{2}){b}_{3}\nonumber\\ &&-~3 \mu f{\gamma }_{2}\omega {g}_{3}-\mu {f}^{2}{h}_{3} \end{eqnarray} \tag{ 39 }$$

$$\begin{eqnarray} &&\fl {\gamma }_{1}{a}_{1}^{\prime}-(1+\mu )f{\gamma }_{2}{g}_{1}^{\prime}+2 \omega {h}_{1}^{\prime}\nonumber\\ &&=-{\Lambda }_{1}{a}_{1}+{\gamma }_{1}\omega {b}_{1}-\frac{3}{4}\beta (-{a}_{1}^{2}+{b}_{1}^{2}){a}_{3}\nonumber\\ &&+~\frac{3}{2}\beta ({a}_{1}{b}_{1}){b}_{3}+{\Psi }_{1}{g}_{1}-{\Psi }_{2}{h}_{1} \end{eqnarray} \tag{ 40 }$$

$$\begin{eqnarray} &&\fl {\gamma }_{1}{b}_{1}^{\prime}-2 \omega {g}_{1}^{\prime}+(1+\mu )f{\gamma }_{2}{h}_{1}^{\prime}\nonumber\\ &&=-{\gamma }_{1}\omega {a}_{1}-{\Lambda }_{1}{b}_{1}-\frac{3}{2}\beta ({a}_{1}{b}_{1}){a}_{3}\nonumber\\ &&-~\frac{3}{4}\beta (-{a}_{1}^{2}+{b}_{1}^{2}){b}_{3}+{\Psi }_{2}{g}_{1}+{\Psi }_{1}{h}_{1} \end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray} &&\fl {\gamma }_{1}{a}_{3}^{\prime}+(1+\mu )f{\gamma }_{2}{g}_{3}^{\prime}+6 \omega {h}_{3}^{\prime}\nonumber\\ &&=-\frac{1}{4}\beta (-{a}_{1}^{2}+3{b}_{1}^{2}){a}_{1}-{\Lambda }_{3}{a}_{3}\nonumber\\ &&+~3{\gamma }_{1}\omega {b}_{3}+{\Psi }_{1 s}{g}_{3}-{\Psi }_{2 s}{h}_{3} \end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray} &&\fl {\gamma }_{1}{b}_{3}^{\prime}-6 \omega {g}_{3}^{\prime}+(1+\mu )f{\gamma }_{2}{h}_{3}^{\prime}\nonumber\\ &&=-\frac{1}{4}\beta (-3{a}_{1}^{2}+{b}_{1}^{2}){b}_{1}-3{\gamma }_{1}\omega {a}_{3}\nonumber\\ &&-~{\Lambda }_{3}{b}_{3}+{\Psi }_{2 s}{g}_{3}+{\Psi }_{1 s}{h}_{3} \end{eqnarray} \tag{ 43 }$$

where the following are defined:

$$\begin{eqnarray} {\Lambda }_{s}&=&-1+\beta \left ({c}_{1}^{2}+\frac{3}{2}{r}_{1}^{2}+\frac{3}{2}{r}_{3}^{2}\right ) \end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray} {\Lambda }_{1}&=&-1+\beta \left (3{c}_{1}^{2}+\frac{3}{4}{r}_{1}^{2}+\frac{3}{2}{r}_{3}^{2}\right ) \end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray} {\Lambda }_{3}&=&-1+\beta \left (3{c}_{1}^{2}+\frac{3}{4}{r}_{3}^{2}+\frac{3}{2}{r}_{1}^{2}\right ) \end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray} {r}_{1}^{2}&=&{a}_{1}^{2}+{b}_{1}^{2};\tqs {r}_{3}^{2}={a}_{3}^{2}+{b}_{3}^{2} \end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray} {\Psi }_{1}&=&(1+\mu ){f}^{2}-{\omega }^{2} \end{eqnarray} \tag{ 48 }$$

$$\begin{eqnarray} {\Psi }_{2}&=&(1+\mu )f{\gamma }_{2}\omega \end{eqnarray} \tag{ 49 }$$

$$\begin{eqnarray} {\Psi }_{1 s}&=&(1+\mu ){f}^{2}-9{\omega }^{2} \end{eqnarray} \tag{ 50 }$$

$$\begin{eqnarray} {\Psi }_{2 s}&=&(1+\mu )3 f{\gamma }_{2}\omega . \end{eqnarray} \tag{ 51 }$$

In the same manner as in section 2.2, the nine equations (35)–(43) are solved for steady-state conditions. Equations (40)–(43) are solved for coefficients g₁,h₁,g₃, and h₃ in terms of a₁,b₁,a₃, and b₃. Substitution of these expressions into the preceding equations (36)–(39) is carried out until, as before, equations (36) and (37) are squared and summed. Likewise, the decision regarding interest in studying low- or high-energy oscillations of the bistable device dictates the selection of the parameter c₁ in solving and substituting equation (35) into this last squared–summed equation. However, unlike in section 2.2, there are two terms ${r}_{1}^{2}$ and ${r}_{3}^{2}$ to compute, which necessitates a second equation. Thus, the expressions for a₃ and b₃ are squared and summed to provide a second equation. Equations (52) and (53) represent the two coupled nonlinear response equations.

$$\begin{eqnarray} 2 5 6{\Delta }^{2}{p}^{2}&=&{r}_{1}^{2}K \end{eqnarray} \tag{ 52 }$$

$$\begin{eqnarray} 1 6{\Delta }^{2}{r}_{3}^{2}&=&{\beta }^{2}({\Gamma }_{1}^{2}+{\Gamma }_{2}^{2}){r}_{1}^{6} \end{eqnarray} \tag{ 53 }$$

where K,Δ,Γ₁, and Γ₂ are provided in the appendix.

It is found that equations (52) and (53) are nonlinearly coupled by the terms ${r}_{1}^{2}$ and ${r}_{3}^{2}$. Therefore, a numerical solution procedure is required. We utilize the MATLAB software command fsolve, as has been successfully employed in recent nonlinear vibration studies [41–43]. Having solved the equations for a given set of system parameters, the individual coefficients may be computed from the expressions provided in the appendix, and the stability of the solution may be checked via a perturbation method as in section 2.2.

The average power across the harvesting circuit resistive load R is computed from P = R|i|²/2. By equation (4) which relates the generated current to the harvester response and according to the solution forms equations (9) and (33), the amplitude of the generated power may be expressed as a sum of the respective harmonic displacement amplitudes, P = P₁ + P₃, where ${P}_{1}\propto {r}_{1}^{2}$ and ${P}_{3}\propto {r}_{3}^{2}$. This practically represents the guiding assumptions of the method of harmonic balance that the total system response is the summation of unique harmonic contributions [25]. Consequently, depending on whether one considers the fundamental or superharmonic response magnitude, equation (47), the unique spectral components of the harvested power are computed. We set m₁ and ω₁ to unity to simplify expressions without loss of generality since other dimensionless groups of equations (7) and (8) are normalized to these parameters. The values of the generated power related to the fundamental and order-3 superharmonic are therefore determined, respectively, by equations (54) and (55).

$$\begin{eqnarray} {P}_{1}&=&(\omega \varepsilon {r}_{1})^{2}/2 \end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray} {P}_{3}&=&(3 \omega \varepsilon {r}_{3})^{2}/2. \end{eqnarray} \tag{ 55 }$$

In the following analytical investigations we utilize the parameters γ = 0.01 and ε² = 0.04 to represent a mechanically lightly damped bistable harvester with an electromagnetic coupling strength typical of many inductive platforms [31, 32]. From equations (54) and (55), it is seen that the power is directly proportional to the response magnitude squared of the bistable harvester and thus the trends of the displacement amplitude correspond to the trends of the generated power. Therefore, for a more meaningful comparison against the bistable energy harvester without auxiliary linear oscillator, we assess the power density, which is the power normalized to the net mass of the system. The power density is computed from P_1,3/(1 + μ) and is employed to evaluate the advantage obtained with the proposed configuration as compared to the individual bistable harvester.

To determine the effects of coupling on the fundamental system response, the roots of equation (24) were computed for the system parameters p = 0.2,β = 1,f = 1,γ₂ = 0.001, and a range of mass ratios, μ = [0.003,0.01,0.1,0.3,1,2]. The harvester nonlinearity strength β may realistically be varied over a large range of values by tailoring the severity of the post-buckled state, and we here employ a value similar to that used in recent bistable energy harvesting studies [27, 28, 35]. The excitation level p used in this section represents a mild level with respect to its normalized definition in equation (7). For a more intuitive interpretation, a recent investigation determined that the normalized excitation level p = 0.2 was comparable to an absolute base acceleration of approximately 2.5 m s⁻² for an individual bistable harvester beam of 8 g mass and linear natural frequency around 15 Hz [25]. It is recognized that high-energy interwell responses are most favorable in a steady-state energy harvesting context as compared to low-energy orbits [14] and, consequently, we focus only on the high-energy dynamics of the systems. The response amplitude, phase, and power density of the bistable inertial generator with and without the attached linear oscillator are shown in figure 2. The individual bistable harvester responses may be computed by setting μ = 0 in the analysis. The displacement magnitude and phase of the linear oscillator are plotted in figure 3.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 2 shows that for the smallest mass ratios, the results intuitively converge onto those computed without the appendage oscillator. As the mass ratio of the add-on linear oscillator increases, the consequence is a uniform amplification of the bistable harvester response magnitude for excitation frequencies less than the natural frequency of the linear oscillator (here f = 1 so that ω₂ = 1). Greater mass addition yields more amplification above the displacement amplitude of the individual bistable harvester, figure 2(a). One trade-off is that the frequency at which enhanced high-energy dynamics are no longer sustainable may be less than that for the system without attached oscillator. (We note for completeness the general rule that where stable high-energy orbits are not predicted, this indicates that the system returns to the low-energy orbit [35].) While the generated power, which is directly proportional to the mechanical response amplitude squared as indicated by equation (54), is amplified uniformly below the linear oscillator natural frequency, the power density is enhanced in a bandwidth upper-bounded by this new cut-off frequency, figure 2(c). The lower bound of this bandwidth of performance improvement appears to be an inexact inflection point across the curve for the individual bistable harvester; in this case, the inflection of the lower bound occurs across a narrow band ω ≈ 0.769 for μ = 0.003 to ω ≈ 0.781 for μ = 2. Since the power density is what defines the advantage over the individual bistable energy harvester, enhancing this metric across the observed bandwidth translates into an improvement of the energy harvesting system.

The corresponding high-energy dynamics of the linear oscillator are shown in figures 3(a) and (b). Recall that in this example, the natural frequency of the appended oscillator occurs at ω = 1 since the tuning ratio is f = 1. Unlike a linear 2DOF system that exhibits two resonant behaviors, the response of the linear oscillator in figure 3(a) appears as if it is directly excited to resonance. This challenges the classical interpretation of a two-mass–spring system exhibiting two unique modes. Clearly, the linear oscillator responses in figure 3(a) suggest a single stable resonant peak. What is further interesting is that the linear oscillator is mostly in-phase with the bistable harvester for excitation frequencies less the linear oscillator natural frequency and approaches a classical resonant 90° phase lag as ω → 1, figure 3(b).

An explanation is proposed for why the high-energy dynamics of the bistable harvester with auxiliary oscillator are magnified but may lose stability at a frequency less than that which destabilizes the individual bistable harvester response. As suggested above, from the perspective of the linear oscillator, the bistable harvester appears to be an equivalent base excitation source judging by the auxiliary oscillator's response magnitude and phase characteristics, figure 3. As a result, the linear oscillator is merely being excited to resonance. However, this equivalent base excitation source may be vulnerable to the vibration of the linear oscillator attachment, and the dynamic interplay between the bistable and linear masses influences the response of the bistable harvester in two ways. Firstly, the response amplitude of the harvester becomes magnified as the linear oscillator is excited at frequencies approaching its resonance, figure 2(a) as ω → 1; this is similar to a positive feedback mechanism. Secondly, after it passes a 90° phase lag, the resonant linear oscillator naturally tends to approach a 180° phase lag with respect to its effective base excitation source (i.e. the bistable harvester). However, this out-of-phase response works directly against the bistable harvester which would otherwise be capable of sustaining high-energy dynamics up to greater frequency were the linear oscillator absent, here approximately ω = 1.7 as indicated by the individual bistable harvester response curve in figure 2(a). Therefore, although it may be an effective base excitation source for the linear oscillator, the bistable harvester in high-energy dynamics represents an excitation source that may be worked against and we find that this influence of the linear oscillator may destabilize the snap through response. Therefore, the bandwidth of enhanced energy harvesting power density appears to be governed in magnitude by the mass ratio μ and upper-bounded in frequency by the tuning ratio f, as will be further explored below.

We now consider the impact of the tuning ratio f and employ irrational numbers to avoid misleading results due to potential internal resonance phenomena [20]. Namely, we consider the case of the individual bistable harvester and the harvester with an auxiliary oscillator of μ = 1 when the tuning ratio is f = π/i, where i = 1,2,...,6. All other system parameters remain the same as from section 3.1.

As was postulated in section 3.1, figure 4 demonstrates that the response magnitude of the harvester with the add-on oscillator is uniformly amplified at frequencies less than the appendage linear oscillator natural frequency. In some cases, the tuning ratio may be used to extend the stable high-energy response bandwidth, as in figure 4 for the largest tuning ratio considered, f = π (red plots). Likewise, as noted in section 3.1, power density improvement occurs within a band of frequencies having as its upper bound the appendage linear oscillator resonance (if the excitation level is sufficient to sustain high-energy response to this frequency). These observations propose an initial framework for the designer seeking improvement upon the performance achieved with an individual bistable energy harvester. Enhanced power harvesting density may be obtained with the addition of the linear oscillator by first understanding the present excitation level p, then tailoring the tuning ratio to set the bandwidth of enhanced response, and lastly modifying the mass ratio according to potential mass constraints so as to change the level of the amplification to meet power harvesting goals.

Zoom In Zoom Out Reset image size

A photograph of the test setup is shown in figure 7. A bistable device is constructed using an inertial mass connected to a translational slide bearing and pre-compressed spring (left of the figure). Two guiding rods are included in the design so as to eliminate the potential for bending or twisting of the compressed spring during snap through. A linear oscillator (right of the figure) on a separate slide bearing could be connected to the bistable mass via a removable linear spring. In this study, two linear coupling springs are employed such that the resulting tuning ratios are f₁ ≈ 4 and f₂ ≈ 2 given that the uncoupled bistable oscillator is observed to have a linear natural frequency of ω₁/2π ≈ 4 Hz. The mass of the appended oscillator is such that μ ≈ 0.75. The system is mounted on an electrodynamic shaker platform (shaking horizontally left and right). The motion of the linear oscillator is isolated from the direct shaker motion. Accelerometers are mounted on the inertial masses while a potentiometer measures the shaker displacement; all measurements are in the horizontal left–right direction. In a similar spirit to the model formulation, we assume that the power generation is proportional to the response amplitude squared of the bistable oscillator given the inclusion of appropriate transduction mechanisms, for example a translational induction mechanism in the bistable oscillator slide bearing or rotational induction generators that constrain the pre-compressed spring.

Zoom In Zoom Out Reset image size

In section 3.2 it is observed that the tuning ratio f may be tailored so as to control the bandwidth and rate of the amplified bistable energy harvester response as compared to the device without the auxiliary element. To validate the analytical predictions, tests are run using slow forward frequency sweeps from 2 Hz at a rate of 0.15 Hz s⁻¹ with the bistable devices initially vibrating in continuous high-energy orbit. The frequency response function (frf) of the bistable oscillator inertial mass acceleration to the shaker input acceleration is evaluated with and without the appendage linear oscillator to serve as a normalized response magnitude metric.

Figure 8 plots the response with and without appendage linear oscillators having tuning ratios f₁ and f₂. The individual bistable device frf (black crosses) exhibits the anticipated trend of gradually increased amplitude towards the frequency at which interwell vibration is destabilized. Once the auxiliary linear element is added to the bistable device, as analytically predicted throughout this paper, the response amplitude (here, acceleration frf) is uniformly amplified below the linear oscillator natural frequencies. In these examples, the linear oscillator natural frequencies are approximately 16 Hz for f₁ and 8 Hz for f₂. As compared to the individual bistable oscillator frf, the addition of the appendage linear oscillator of tuning ratio f₁ (red squares) amplifies the frf by approximately 80% at 4 Hz; when using the auxiliary oscillator of linkage with f₂ (blue circles) 360% enhancement at 3.75 Hz is obtained. These findings are qualitatively in good agreement with those predicted in section 3.2 by the fundamental harmonic balance analysis, figure 4, and demonstrate the utility of the simplified one-term solution approach for this response regime.

Zoom In Zoom Out Reset image size

The analytical predictions of section 4 indicated that at certain frequencies near the point of high-energy orbit destabilization, the system with the appendage oscillator may exhibit two interwell vibration responses uniquely characterized by their proportion of fundamental and superharmonic vibration. The two responses are co-existent and the ultimate determination as to which is attained is dependent on the initial conditions. To evaluate this feature experimentally, we excite the bistable oscillator at 4 Hz with and without the auxiliary element of tuning ratio f₁ ≈ 4 and monitor the steady-state responses.

Figures 9 and 10 respectively plot the acceleration autospectra and time series of the two co-existent forms of interwell vibration exhibited by the bistable oscillators with and without the appendage element. Figure 9 shows the response corresponding to continuous high-energy orbit having a predominant response frequency commensurate with the excitation. The system with auxiliary linear oscillator diffuses approximately one order of magnitude less vibrational energy to the order-1/3 subharmonic (at three times the excitation frequency) as compared with the individual bistable device. The spectral concentration of bistable harvester snap through energy is an advantageous result because harvesting circuitry should optimally switch interface connections at the frequency of greatest harvester response to maximize power flow to the storage element [28, 44, 45]. We find that at the driving frequency, the bistable device autospectrum with appendage oscillator is 48.7% greater than the individual bistable device response: 1016 (m s⁻²)² Hz⁻¹ as compared to 685.1 (m s⁻²)² Hz⁻¹, figure 9(a). Due to the linkage to the appendage oscillator, an additional spectral peak is observed at a frequency of approximately four times the excitation, which corresponds to the natural frequency of the add-on element; however, this peak is several orders of magnitude less than that measured at the driving frequency and does not detract from the overall advantage obtained via the system with auxiliary oscillator.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

When the excitation at 4 Hz is begun with different initial conditions, alternative steady-state vibrational responses are discovered, and these are presented in figure 10. In both cases, sub- and superharmonic responses of fraction- and integer-multiples of order-3 are apparent. However, the individual bistable oscillator exhibits a far more diffuse autospectrum typical of chaotic dynamics; the time series in figure 10(b) likewise reflects aperiodic vibration. In contrast, the system with auxiliary oscillator yields a substantially more concentrated multi-harmonic response, in fact concentrating more of the vibrational energy at 1/3 of the excitation frequency than at the drive frequency itself. The time series of the bistable oscillator having an add-on element in figure 10(b) demonstrates the more organized composition of the vibration. At 1/3 of the excitation frequency, the bistable device autospectrum with the auxiliary element is 49.6% greater than the maximum autospectral peak for the individual bistable device that occurs at the drive frequency itself: 1.328 (m s⁻²)² Hz⁻¹ as compared to 0.8338 (m s⁻²)² Hz⁻¹, figure 10(a).

These results suggest that a bistable energy harvester with appendage linear oscillator is less susceptible to non-ideal excitation conditions than the individual harvester because its vibrational response may become more concentrated to target multi-harmonic spectral components, which is more suitable for efficient power harvesting circuits [28, 44, 45]. This feature was suggested in the analytical predictions of figures 5 and 6 in terms of the multiple solution responses for the coupled system in the band near ω ≈ 0.8. In contrast, the individual bistable device may attain either periodic high-energy orbit or adverse aperiodic response less useful for power harvesting. The experimental results of figures 9 and 10 validate the trends of the analytical predictions and exemplify the robustness of the proposed design methodology.