June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

International Journal of Bifurcation and Chaos

⃝c World Scientific Publishing Company

OPTIMIZING THE ELECTRICAL POWER IN AN ENERGY

HARVESTING SYSTEM

MATTIA COCCOLO
Nonlinear Dynamics, Chaos and Complex Systems Group, Departamento de Fı́sica
Universidad Rey Juan Carlos, Tulipán s/n, 28933 Móstoles, Madrid, Spain
GRZEGORZ LITAK
Department of Applied Mechanics, Lublin University of Technology,
Nadbystrzycka 36, PL-20-618 Lublin, Poland

JESÚS M. SEOANE
Nonlinear Dynamics, Chaos and Complex Systems Group, Departamento de Fı́sica
Universidad Rey Juan Carlos, Tulipán s/n, 28933 Móstoles, Madrid, Spain

MIGUEL A.F. SANJUÁN

Nonlinear Dynamics, Chaos and Complex Systems Group, Departamento de Fı́sica
Universidad Rey Juan Carlos, Tulipán s/n, 28933 Móstoles, Madrid, Spain

Received (to be inserted by publisher)

In this paper, we study the vibrational resonance phenomenon (VR) as a useful mechanism for
energy harvesting purposes. A system, driven by a low frequency and a high frequency forc-
ing, can give birth to the vibrational resonance phenomenon, when the two forcing amplitudes
resonate and a maximum in amplitude is reached. We apply this idea to a bistable oscillator
that can convert environmental kinetic energy into electrical energy, that is, an energy har-
vester. Normally, the VR phenomenon is studied in terms of the forcing amplitudes or of the
frequencies, that are not always easy to adjust and change. Here, we study the VR generated
by tuning another parameter that is possible to manipulate when the forcing values depend on
the environmental conditions. We have investigated the dependence of the maximum response
due to the VR for small and large variations in the forcing amplitudes and frequencies. Besides,
we have plotted color coded figures in the space of the two forcing amplitudes, in which it is
possible to appreciate different patterns in the electrical power generated by the system. These
patterns provide a useful information on the forcing amplitudes in order to produce the optimal
electrical power.

1. Introduction
It is possible, all around us, to see small but powerful electrical devices greedy of energy. In fact, the
performance of the electronic devices has increased rapidly along with their energy needs, while the capacity
of the batteries has become suddenly inadequate [Paradiso & Starner, 2005]. Therefore, powering electronics
devices without depending exclusively on the batteries, but transforming the environmental energy into
electrical energy, has grown as an interesting open research field, called energy harvesting.

1
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

2 M. Coccolo, G. Litak, J.M. Seoane, and M.A.F. Sanjuán

In this sense, mechanical vibrations are a possible and reliable energy source that can be exploited. In
fact, different ideas have been proposed in order to transform environmental vibration kinetic energy into
electrical energy, like using piezoelectric or electrostatic effects [Anton & Sodano , 2007; Arnold, 2007; Litak
et al., 2010], i.e., coupling a mechanical system as a source to a transduction mechanism. Some numerical
results [Gammaitoni et al., 2009] have shown that the bistable energy harvesters are able to provide more
energy in a frequency broadband than their linear counterparts, initially preferred [Roundy et al. , 2003].
It means that the nonlinear oscillators can exploit a wider spectrum of vibration frequencies, making them
more adaptable and the best choice as shown in [Friswell et al., 2011; Borowiec et al., 2014; Haris et al.,
2014; Pellegrini et al., 2012; Harne & Wang , 2013; Twiefel & Westermann, 2013].
In the scientific literature on vibrations is usual to have the impression that all vibrations are detri-
mental because most publicized work discusses vibration reduction in one form or another. But since we
talk about exploiting vibrations as a source of electrical energy, nonlinear environmental vibrations can be
intentionally introduced into designs or enhanced in amplitude. The final objective is to take advantage of
benefits and then studied by using mathematical tools for modeling and predicting the vibrational behav-
iors. Besides, it is important to state that vibrations are defined to be linear when the corresponding force
of the potential is linear (F = −kx) and nonlinear otherwise. Therefore, we have decided to use nonlinear
vibrations, not only because they are more general and realistic, but also because energy harvesters are
supposed to be small enough to be contained into electronic devices. That means that even small vibration
amplitudes result large in consideration with the system. Moreover, sometimes in kinetic energy harvesting
systems, it can be useful to increment the effect of the vibrations in amplitude in order to produce more
energy [Coccolo et al., 2014].
A phenomenon that can help to achieve this enhancement is the vibrational resonance (VR) [Landa
& McClintock, 2000], where the resonance concept is the tendency of a system to oscillate with greater
amplitude at some frequencies than at others. The VR phenomenon appears when a bistable system with
a low (LF) and a high frequency (HF) forcing gives a response amplitude at the LF that grows until a
maximum and then decreases, while we vary the amplitude of the HF forcing or the frequency. Other
theoretical aspects of the VR have been developed in [Gitterman, 2001; Blekhman & Landa, 2004; Zaikin
et al., 2002]. The VR phenomenon can be explained as an amplification of the LF signal, due to a reduction
of the stiffness of the system induced by the HF force. In other words, it happens when the low frequency
signal is able to induce cross-well transitions. So far, this phenomenon has been thoroughly studied in a
large class of dynamical systems [Gandhimathi et al., 2006; Daza et al., 2013; Jayakumari et al., 2009;
Deng et al., 2009; Rajasekar et al., 2011], among others. In [Coccolo et al., 2014], we have studied the
influence of this phenomenon on a Duffing oscillator driven by a bi-harmonic excitation, focusing our work
on the effects of VR in harvesting systems. Following the same path, we have decided to study a bistable
oscillator, studied in [Cottone et al., 2009; Gammaitoni et al., 2009; Vocca et al., 2012] numerically and
experimentally, driven by a HF and LF forcing, in which the double well potential is given by the repulsion
of an external and an internal magnet situated as shown in Fig. 1. In this case, the two forcings simulate
environmental vibrations and we analyze their interaction in order to exploit the VR phenomenon by tuning
a different parameter values of ∆, which is the distance between the magnets (see Fig. 1). Accordingly,
we have focused our attention in the power generated by the system by producing a two-dimensional map
with an appropriate color code, that associates to every LF and HF forcing pair (F, f ) a different color,
that represents a value of the mean electrical power < P >. As we will show later on, different patterns
can be visualized. To summarize, we will analyze the phenomenon of VR for different values of the distance
between magnets ∆. Furthermore, we compute the electrical power generated for different frequencies and
amplitudes of the forcing terms for various values of ∆.
This paper is organized as follows. In Sec. 2, we describe our prototype model, the bistable nonlinear
oscillator which emulates our mechanical system. The phenomenon of the VR in our system and its impli-
cations in energy harvesting is described in Sec. 3. In Sec. 4, we study the response of the average power
of the system in the vibrational resonance regime. A map of the electrical power response of the system is
described in Sec. 5. A discussion and the main conclusions of this manuscript are presented in Sec. 6.
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

Optimizing the electrical power in an energy harvesting system. 3

Fig. 1. This figure shows the harvester considered. Parameter ∆ is the distance between the magnets, PZT the piezoelectric
patch, R and Cp are the resistance and the capacitor of the coupled circuit, and the beam displacement is along x axis.

Fig. 2. Plot of the potential U (x), given by Eq. 3, for different values of the distance between the magnets ∆ = 0.005, ∆ =
0.007, ∆ = 0.010 and ∆ = 0.02. The rest of the parameters are shown in Table 1, in agreement with [Cottone et al., 2009;
Gammaitoni et al., 2009; Vocca et al., 2012]. All magnitudes follow the International System of units.

Table 1. The simulations parameters

Parameters Description Value

m cantilever equivalent mass 0.018 Kg
γ viscous parameter 0.022 kg/s
Kv piezoelectric coupling constant 0.0011 N/V
Kc position to voltage coupling coeff. 4.15 × 103 V /m
RL load resistance 300 × 103 Ω
C piezoelectric capacitance 112 × 10−9 F
Kef f cantilever effective elastic constant 26.6 N/m
M effective magnetic moment 0.051Am2
µ0 magnetic permeability 4π × 10−7 H/m

2. The model
We have chosen for our study the mechanical system depicted in Fig. 1, which has been studied numerically
and experimentally in [Cottone et al., 2009; Gammaitoni et al., 2009; Vocca et al., 2012]. This system can
be modeled as a bistable potential like and the equation of motion reads:

dU (x)
mẍ + + γ ẋ + Kv v(t) = F cos(ωt) + f cos(Ωt), (1)
dx

where U (x) is the potential energy function of the cantilever with equivalent mass m, x is the vertical
displacement around its mean value, the third term on the left γ ẋ, gives the energy dissipation due to the
bending and in the fourth term we have v(t) the voltage and Kv , the coupling constant of the piezoelectric
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

4 M. Coccolo, G. Litak, J.M. Seoane, and M.A.F. Sanjuán

Fig. 3. Figures (a) and (c) plot Qx vs ∆, for f = 0.1 × 10−2 and f = 0.5 × 10−2 . Figures (b) and (d) plot Qv2 vs ∆,
for f = 0.1 × 10−2 and f = 0.5 × 10−2 . All different curves are plotted for different values of F : F = 0.11 × 10−1 , F =
0.12 × 10−1 , F = 0.13 × 10−1 , and F = 0.14 × 10−1 . Note that there is a value of the ∆ parameter that maximizes the Q
factor, showing the occurrence of VR, although it shifts to smaller values when the F value grows up. In all these curves, we
have used (x0 , y0 , v0 ) = (0.001, 0, 0) as the initial condition.

Fig. 4. Figures (a) and (c) plot the Q factor for the position, Qx vs ∆, for f = 0.1 × 10−2 and f = 0.5 × 10−2 . It is possible
to see that in both figures the curves do not change significantly when the Ω value grows up. Figures (b) and (d) plot the Q
factor for the power Qv2 vs ∆, for f = 0.1 × 10−2 and f = 0.5 × 10−2 . Here, the scale changes and also the shape of the
curves. Moreover, a new peak shows up in Fig. 4 (d). All different curves are plotted for different values of the HF frequency
Ω. In all these curves, we have used (x0 , y0 , v0 ) = (0.001, 0, 0) as the initial condition.

sample, so that Kv v(t) gives the energy transferred to the electric load R with the coupling equation:
v̇ + v(t)/τp − κc ẋ = 0, (2)
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

Optimizing the electrical power in an energy harvesting system. 5

Fig. 5. Figures (a) and (b) plot the average electrical power, < P >∝< V 2 > vs ∆, for f = 0.1 × 10−2 and f = 0.5 × 10−2
for different values of F . Figures (c-d) plot the average electrical power, < P >, vs ∆, for f = 0.1 × 10−2 and f = 0.5 × 10−2
for different values of Ω. Figures (a-c) show a peak for the same ∆ value, but figure (d) maximum is displaced with respect
the others and it is in agreement with Fig. 4 (d). In all these curves, we have used (x0 , y0 , v0 ) = (0.001, 0, 0) as the initial
condition.

Fig. 6. Figures (a) and (c) plot the Q factor of the displacement Qx vs ∆, for f = 0.1 × 10−2 and f = 0.5 × 10−2 . Figures
(b) and (d) plot the the Q factor of the power Qv2 vs ∆, for f = 0.1 × 10−2 and f = 0.5 × 10−2 . All different curves are
plotted for different values of F , larger with respect of Fig. 3 (a-d). The figures show that the peaks of both Qx and Qv2 shift
in function of the value of the F . In all these curves, we have used (x0 , y0 , v0 ) = (0.001, 0, 0) as the initial condition.

where κc is the position to voltage coupling coefficient. The time constant of the piezoelectric dynamics,
τp , is related to the coupling capacitance Cp and to the resistive load R by τp = RCp . Finally, F cos(ωt)
and f cos(Ωt) are the low frequency (LF) and the high frequency (HF) excitations, respectively. Notice
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

6 M. Coccolo, G. Litak, J.M. Seoane, and M.A.F. Sanjuán

Fig. 7. Figures (a) and (c) plot the Q factor of the displacement Qx vs ∆, for f = 0.1 × 10−2 and f = 0.5 × 10−2 . Figures (b)
and (d) plot the the Q factor of the power Qv2 vs ∆, for f = 0.1 × 10−2 and f = 0.5 × 10−2 . All different curves are plotted
for different values of HF frequency Ω, larger with respect of Fig 4(a-d). The figures show that the peaks of both Qx and Qv2
show different behaviors. Below, we plot figure (d) for specific ranges of frequency Ω. We have used (x0 , y0 , v0 ) = (0.001, 0, 0)
as the initial condition.

that we have chosen for convenience F ≫ f as the respective forcing amplitudes, while ω = 0.5 and Ω are
the frequencies, that have always to satisfy the VR condition Ω ≫ ω. The HF frequency value is Ω = 5 if
not indicated otherwise.
The potential energy function in Eq. 1 is given by:
1
U (x) = Kef f x2 + (ax2 + b∆2 )−3/2 , (3)
2
with Kef f , a and b representing constants related to the physical parameters of the cantilever. The param-
eter a is given by a = d2 (µ0 M 2 /2πd)−2/3 where µ0 is the permeability constant, M the effective magnetic
moment and d = 2.97 a geometrical parameter related to the distance between the measurement point and
the cantilever length. The variable b is written as b = a/d2 . The tuning parameter ∆ is used to move from
a monostable potential to a bistable potential (see Fig. 2). In other words, it is the distance between the
magnet on the top of the cantilever and the fixed magnet. Later, we discuss the influence of this parameter
on the electrical power harvested in the system of Eqs. (1-3). We start, from now on, every simulation
from the initial condition (x0 , y0 , v0 ) = (0.001, 0, 0) (y = ẋ), using a Runge-Kutta integrator [Butcher,
1987] of fourth order. The parameter values used for the simulations are described in Table 1, according
to [Cottone et al., 2009; Gammaitoni et al., 2009; Vocca et al., 2012], where a comparison between the
numerical simulations and an experimental apparatus has been done. We do not explicit the units other
than in the Table 1, because they do not change along the article. We have decided to study the harvester
subjected to continuous vibrations for values of the amplitudes and frequencies sufficiently sensitive to a
human being, as we can see in [NSW, 2006].

3. The vibrational resonance

The first step is the numerical study of the response of the system defined by Eqs. (1)-(3) for the displace-
ment variable x, and the electrical power P ∝ V 2 , where V is the voltage generated by the system. When
the system has a LF signal of amplitude F and frequency ω, and a HF signal of amplitude f ≪ F and
frequency Ω ≫ ω, the VR phenomenon can take place.
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

Optimizing the electrical power in an energy harvesting system. 7

Fig. 8. Figures (a-e) plot the curves of Fig. 7 (d) for a specific range of the frequency Ω. It is possible to appreciate that
the scale changes abruptly from the figure (c) to (d), but differently from Fig. 6 (a-d) the peaks are settled around the value
∆ = 0.001. To make the legend easier to read, the position of the Ω values match the ones of Fig. 7 (d).

The usual procedure to search for VR is to compute, for different amplitudes f or frequency Ω [Landa
& McClintock, 2000], the Q factor
√
(Cs2 + Cc2 )
Q= , (4)
F

where
∫ nT
2
Cs = Γ(t) sin(ωt)dt (5)
nT 0
∫ nT
2
Cc = Γ(t) cos(ωt)dt, (6)
nT 0

and Γ(t) is, in our case, the displacement x or the electrical power generated P ∝ V 2 , the number of
complete oscillations of the LF signal is n and T = 2π/ω is its period. The VR occurs, if it is possible
to find a value of f or of Ω that maximizes the Q factor. This means that a particular value of the HF
periodic signal has been found that optimizes the response of the system to the weak LF periodic signal.
Here, we have explored a different scenario, where the values of the frequencies and forcing amplitudes are
fixed by the environmental conditions and we modify only the parameter ∆.
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

8 M. Coccolo, G. Litak, J.M. Seoane, and M.A.F. Sanjuán

Fig. 9. Figures (a-b) plot the average electrical power < P >∝< V 2 > vs ∆, for f = 0.1 × 10−2 and f = 0.5 × 10−2 for
different values of the F . Figures (c-d) plot the average electrical power < P > vs ∆, for f = 0.1 × 10−2 and f = 0.5 × 10−2
for different values of the HF frequency Ω. The figures (a-b) show peaks in the generated electrical power < P > for the same
values for which the VR occurs, as shown in Figs. 7. The figures (c-d) show that the behavior of the electrical power generated
is largely affected by the value of Ω and we lose correlation with the Q factor figures, as we show in detail below. In all these
curves, we have used (x0 , y0 , v0 ) = (0.001, 0, 0) as the initial condition.

Fig. 10. Figures (a-c) show in details Fig. 9 (c). We can observe that the dependence of < P > on ∆ takes larger values for
frequencies Ω in the range 60 ≤ Ω/ω ≤ 80.

3.1. Vibrational resonance for small variations in the forcing amplitudes and
frequencies ratio
We have started our analysis with a value of the HF forcing amplitude f = 0.1 × 10−2 (as shown in
Figs. 3 (a-b)). In these figures, the different curves related to different values of the LF forcing amplitude
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

Optimizing the electrical power in an energy harvesting system. 9

Fig. 11. Figures (a-d) show in details Fig. 9 (d). In this case, the dependence of < P > on ∆ takes larger values for frequencies
Ω in the range 60 ≤ Ω/ω ≤ 90.

0.1 × 10−1 < F < 0.14 × 10−1 show how the peaks detach around the value of ∆ = 0.01. The same behavior
happens when we change the amount of the HF forcing amplitude f = 0.5 × 10−2 (see Figs. 3 (c-d)), only
a new peak detach in the curve of F = 0.011 for a slightly larger ∆ value. Also, we have computed the
same curves for different amounts of the frequency Ω, as shown in Figs. 4 (a-d). In Figs. 4 (a) and (c) we
can see that the peaks detach for the same ∆ value. However, for a HF forcing value f = 0.5 × 10−2 , as
shown in Fig 4 (d), the Q factor for the electrical power Qv2 shows a different behavior with respect to the
Qx and another peak pop up for a slightly larger value of ∆. The analysis of the previous figures shows
us that a variation of the experimental parameter ∆ can generate the VR phenomenon for which the role
of the parameter f is crucial. In the situation in which VR exists, we have calculated the value of ∆ for
which we have a maximum in the average electrical power generated by the mechanical system, namely
< P >∝< V 2 >. This maximum can be observed in Fig. 5 (a) where we show the values of < P > for the
case with f = 0.1 × 10−2 . On the other hand, the average electrical power generated by the mechanical
system, for the case with f = 0.5 × 10−2 , is shown in Fig. 5 (b), where qualitatively its behavior is the same
that in the case of the previous figure. We can see, as shown in Fig. 3, that the value of the parameter ∆
for which the peaks are maximal coincides with the value for which the Q factor is also maximum. The
last figure, Fig. 5 (d), shows us that the optimal distance ∆ between the magnets for the VR to occur can
depend on the LF forcing, in fact a new maximum shows up while Ω is increasing.

3.2. Vibrational resonance for larger variations in the forcing amplitudes and
frequencies ratio
In the previous subsection, we have shown that for a small difference of the amplitude F and the frequency
Ω of the forcing, the peak of the VR occurs for almost the same ∆ value. We have also shown that the
maximum in electrical power and the peaks in the Q factor are related. In this subsection, we investigate
the occurrence of the VR peak when the variation of the ratio between the two amplitudes and frequencies
of the forcing is larger, although the vibration values remain in a humanly acceptable scale. We study
the possibility to establish that there are values of the ∆ parameter that maximize the system response
in electrical power for given environmental vibrations, even when the vibrational conditions vary on a
larger scale. The Q factor for the displacement x and the electrical power < P > in the cases of different
amplitudes F can be seen in Figs. 6 (a-d). In these figures, we have a similar behavior between them,
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

10 M. Coccolo, G. Litak, J.M. Seoane, and M.A.F. Sanjuán

Fig. 12. Figures (a-e) show the average electrical power < P >∝< V 2 > vs the ratio Ω/ω and F/f . It is possible to see that
in figures (a)-(b) we have a maximum for values of the HF frequency Ω in the range 80 ≤ Ω/ω ≤ 100. In figures (c)-(d) we
observe that a HF frequency value Ω = 80ω detach for higher values of the amplitudes F . Figure (e) is a zoom of figure (d),
in which the curve Ω = 90ω gives the maximum electrical power for smaller values of the amplitudes.

insofar the scale of the F increases. On the other hand, in Figs. 7 (a-d), we have changed the frequency
ratio, showing a strong change in the shape of the curves shown in Figs. 7 (c) and 7 (d) with respect to
the ones shown in Figs. 7 (a) and 7 (b). In Figs. 8 (a-e) we attempt to clarify Fig. 7 (d) by plotting in
the different panels the shape of the curves for a shorter range of Ω. In these last figures, we can see that
the peaks are localized in a narrow area of the figure, typically around the value ∆ = 0.001. When we
study the average electrical power < P > shown in Figs. 9 (a-b), we can see that the maximum electrical
power follows the peak of the Qv2 as F increases, while otherwise in the Figs. 9 (c-d) it does not happen.
In particular we can appreciate besides the difference in scale, that the maxima in the electrical power are
displaced with respect to the Qv2 shown in Fig. 7 (d). To better visualize the curves of Figs. 9 (c-d), we
have split the two figures in Figs. 10 (a-c) and Figs. 11 (a-d) respectively. Analyzing these figures, we can
appreciate some similarity between them, in particular that for different values of the ∆ we have different
maximal responses in electrical power depending on the HF frequencies Ω. An analysis of the figures from
smaller to bigger ∆ values can be fruitful. It is possible to see that just before the value ∆ = 0.001 the
maximum is reached by the curve Ω = 100ω, and for ∆ = 0.001 for lower frequencies, ω ≤ Ω ≤ 50ω.
Then, for 0.001 < ∆ < 0.013 the maximum is reached by the curve Ω = 60ω, for 0.013 ≤ ∆ ≤ 0.016 by
the curve Ω = 70ω, and finally for ∆ ≥ 0.016 by the curve Ω = 80ω. A discrepancy can be detected for
Ω = 50ω, in fact in Fig. 11 (b) a higher peak detaches for a different value of ∆ than in Fig. 10 (b). Indeed,
another difference is the scale of the average electrical power generated. We want to underline how the f
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

Optimizing the electrical power in an energy harvesting system. 11

Fig. 13. The figure represent the average electrical power < P >∝< V 2 > generated by the harvester for different values
of the forcing amplitudes f, F for ∆ = 0.007. The color gradient indicates the amount of power generated in function of the
forcing amplitudes, from blue to red, respectively from the smallest amount of electrical power to the highest. It is possible to
see that the bigger the two amplitudes the bigger is the electrical power generated. We have used (x0 , y0 , v0 ) = (0.001, 0, 0) as
the initial condition.

value can be decisive in order to harvest more energy. We have generated all these figures with an initial
condition (x0 , x˙0 , v0 ) = (0.001, 0, 0). We have also studied the average electrical power generated by 300
different initial conditions, then we have calculated the mean value and plotted it in function of ∆. The
figures generated were not only similar, but equal to the figures above. This tells us that the system has a
strong robustness towards variations of the initial conditions.

4. Numerical analysis of the average electrical power of the system

After showing the response of the system in function on the parameter ∆, we study the corresponding
response in function of both the forcing amplitude F/f and the frequency Ω/ω ratios since these parameters
are crucial in the energy harvested by the mechanical system.
In Figs. 12 (a-b) it is possible to observe a maximum for 80 < Ω/ω < 100, also in this case the
magnitude of the peaks of the electrical power response is very different, due to the value of the HF forcing
amplitude f . Figures 12 (c-e) show the electrical power response of the system for larger amplitudes where
every curve is plotted for different Ω values. In Fig.12 (c) the response grows in an oscillatory way and
only for F/f ≥ 300 the curve related with Ω = 80ω detaches and indicates more generated electrical power
than the others. On the other hand, in Fig. 12 (d) we can see that for large values of F/f , the curve of
Ω = 80ω gives a higher response in the electrical power with respect to the others. But, if we take a look
at Fig. 12 (e), that is a zoom of the previous figure in the region F/f < 100, for values of the forcing
F/f < 50, the maximum generated electrical power is given by the curve of Ω = 90ω. To summarize,
the largest value of the response of the system takes place for f = 0.5 × 10−2 as we can easily observe in
Figs. 12(b) and 12(d) where the maximum occurs for Ω = 80ω.
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

12 M. Coccolo, G. Litak, J.M. Seoane, and M.A.F. Sanjuán

Fig. 14. The figure represent the average electrical power < P >∝< V 2 > generated by the harvester for different values of
the forcing amplitudes f, F for ∆ = 0.01. The color gradient indicates the amount of electrical power generated in function of
the forcing amplitudes, from blue to red, from the smallest amount of electrical power to the highest respectively. It is possible
to see that some structures show up and we have regions of F, f values for which the electrical power is maximum. We have
used (x0 , y0 , v0 ) = (0.001, 0, 0) as the initial condition.

5. The map of the electrical power response of the system

We have seen in a previous section that an enhancement of the average electrical power generated by the
system has taken place when the VR phenomenon occurs. Now, we explore this by using a (F, f ) plot, which
we call forcing space. For that purpose, we have carried out different simulations for three different values
of ∆, while the other parameters do not change so that ω = 0.5, Ω = 10ω and (x0 , x˙0 , v0 ) = (0.001, 0, 0).
The figures that we have generated show a color gradient, from blue to red in function of the electrical
power, for different values of both HF and LF forcing amplitudes. In Fig.12 and in Fig. 14, even if the
scale of the electrical power is different, a common pattern is observed. We can see that the higher the
two forcing amplitudes, the higher the electrical power generated. However, we have a completely different
pattern in Fig. 13, where more complicated structures are shown. It is possible to distinguish several red
regions, that show us the values of the forcing amplitudes that allow us to harvest the maximum electrical
power. These regions are surrounded by yellow ones, in which the electrical power that we can harvest is
smaller. Finally, we can observe some blue regions for very small forcing amplitudes f , where the electrical
power generated is even smaller or almost zero. Note that the color gradient scale of Figs. 13, 14 and 15 are
very different. In the Fig. 13, the electrical power generated by the system is much higher than in the other
cases. With the analysis of these last figures we have shown that for the case of ∆ = 0.001, an enhancement
of the electrical power generated by the system appears. Besides, we can visualize as well the regions in
the forcing space for which this occurs when the pattern is more clear. We think that this kind of graphical
tools might help experimentally to obtain an optimal enhancement of the average electrical power. We
think that, this is crucial for a better understanding of the electrical power associated to different values
of the forcing amplitudes. In fact, it gives us the opportunity to adjust the harvester to the environmental
conditions in order to harvest the exact amount of energy that we need from the system.
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

Optimizing the electrical power in an energy harvesting system. 13

Fig. 15. The figure represent the average electrical power < P >∝< V 2 > generated by the harvester for different values of
the forcings amplitudes f, F for ∆ = 0.015. The color gradient indicates the amount of electrical power generated in function of
the forcing amplitudes, from blue to red, respectively from the smallest amount of electrical power to the highest. It is possible
to see that the bigger the two amplitudes the bigger is the electrical power generated. We have used (x0 , y0 , v0 ) = (0.001, 0, 0)
as the initial condition.

6. Conclusions
In this work, we have studied a mechanical system which has been modeled as a bistable oscillator designed
to harvest energy driven by a HF and a LF harmonic external forcing. We have shown how the distance
between the magnets, ∆, may be crucial for the occurrence of the VR phenomenon, and as a matter of
fact we have computed the values for which VR occurs. Then, we have related the VR with the optimal
electrical power responses. We have also computed the average electrical power for small and large variations
in the amplitudes and forcing frequencies. Furthermore, we have calculated the average electrical power
< P >∝ V 2 generated for different values of F and f where further information related to the different
regions of the forcing space can be derived as the optimal situation for which the electrical power generated
reaches its maximum and minimum. Some interesting patterns have been observed mainly for the case
∆ = 0.001. We think that these graphical study constitutes a useful tool to show for which values of the
forcing amplitudes the electrical power generated is enhanced, as well as to gain a better knowledge about
its prediction. To summarize, we believe that a complete study of this kind should give the possibility
to set the harvester so that it can respond, maybe automatically, in the best way to the environmental
conditions.

7. Acknowledgements
We acknowledge financial support from the Spanish Ministry of Economy and Competitivity under Project
No. FIS2013-40653-P. MC acknowledges financial support for a research stay from the Universidad Rey
Juan Carlos, under grant number PPR-2012-28, and the warm hospitality received at Lublin University
of Technology, where part of this work was carried out. GL was supported by the Polish National Science
Center under the grant agreement No. 2012/05/B/ST8/00080.
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

14 REFERENCES

References
Anton, S.R. & Sodano, H.A. [2007] “A review of power harvesting using piezoelectric materials (20032006),”
Smart Mater. Struct. 16, 1.
Arnold, D.P. [2007] “Review of microscale magnetic power generation,” IEEE Trans. Magn. 43, 3940.
Blekhman, I.I., & Landa, P.S. [2004] Conjugate resonances and bifurcations in nonlinear systems under
biharmonical excitation. Int. J. Non linear Mech. 39, 421.
Borowiec, M., Rysak, A., Betts D.H., Bowen, C.R, Kim, H.A., & Litak, G. [2014] “Complex response of
the bistable laminated plate: Multiscale entropy analysis,” Eur. Phys. J. Plus 129, 211.
Butcher, J. C. [1987] “The numerical analysis of ordinary differential equations, Runge-Kutta and general
linear methods,” Wiley, Chichester and New York.
Coccolo, M., Litak, G., Seoane, J.M., & Sanjuán, M.A.F. [2014] “Energy harvesting enhancement by
vibrational resonance” Int. J. of Bifurcat. Chaos 24, 1430019 .
Cottone, F., Vocca, H., & Gammaitoni, L. [2009] “Nonlinear energy harvesting,” Phys Rev Lett 102,
080601.
Daza, A., Wagemakers, A., Rajasekar, S. & Sanjuán, M.A.F. [2013] “Vibrational resonance in a time-
delayed genetic toggle switch,” Commun. Nonlinear Sci. Numer. Simulat. 18, 411–416.
Deng, B., Wang, J., & Wei, X. [2009] “Effect of chemical synapse on vibrational resonance in coupled
neurons,” Chaos 19, 013117 .
Friswell, M.I., Ali, S.F., Adhikari, S., Lees, A.W., Bilgen, O. & Litak, G. [2011] “Nonlinear piezoelectric
vibration energy harvesting from an inverted cantilever beam with tip mass,” J. Int. Mat. Syst. Struc.
23, 1505.
Gammaitoni, L., Neri, I., & Vocca, H. [2009] “Nonlinear oscillation for vibrational energy harvesting,”
Appl. Phys. Lett. 94, 164102.
Gandhimathi, V.M., Rajasekar, S. & Kurths, J. [2006] “Vibrational and stochastic resonances in two
coupled overdamped anharmonic oscillators,” Phys. Lett. A 360, 279–86.
Gitterman, M. [2001] “Bistable oscillator driven by two periodic fields,” J. Phys. A Math. Gen. 34, 355.
Haris, P., Bowen, C.R., & Kim, H.A. [2014] “Manufacture and characterisation of piezoelectric broadband
energy harvesters based on asymmetric bistable laminates,” J of Multifunc. Compos. 3, 113.
Harne, R.L., & Wang, K.W. [2013] “A review of the recent research on vibration energy harvesting via
bistable systems,” Smart Mater. Struct. 22, 023001.
Jayakumari, S., Chinnathambi, V., Rajasekar, S. & Sanjuán, M.A.F. [2009] “Vibrational resonance in
groundwater-dependent plant ecosystems,” Chaos 19, 043128.
Landa, P.S. & McClintock, P.V.E. [2000] “Vibrational resonance,” J. Phys. A Math. Gen. 33, 433–38.
Litak, G., Friswell, M.I. & Adhikari, S. [2010] “Magnetopiezoelastic energy harvesting driven by random
excitations,” Appl. Phys. Lett. 96, 214103.
Paradiso, J.A. & Starner, T. [2005] “Energy scavenging for mobile and wireless electronics,” IEEE Pervasive
Comput. 4, 18–27.
Pellegrini S.P., Tolou N., Schenk, M., & Herder, J.L. [2012] “Bistable vibration energy harvesters: A review
journal of intelligent material systems and structures,” J. Intell. Mater. Syst. Struct. 24, 1303.
Rajasekar, S., Abirami, K., & Sanjuán, M.A.F. [2011] “Novel vibrational resonance in multistable systems,”
Chaos 21, 033106.
Roundy, S., Wright, P.K., & Rabaey, J. [2003] “A study of low level vibrations as a power source for wireless
sensor nodes,” Comput. Commun. 26, 1131–1144.
The data and tables on the forcing amplitudes and frequencies can be found in:
Department of Environment and Conservation [2006] “Assessing vibration: a technical guideline,”
http://www.environment.nsw.gov.au/resources/noise/vibrationguide0643.pdf
Twiefel, J., & Westermann H. [2013] “Survey on broadband techniques for vibration energy harvesting,”
J. Intell. Mater. Syst. Struct. 24, 1291.
Vocca, H., Travasso, F., Neri, I., & Gammaitoni, L. [2012] “Kinetic energy harvesting with bistable oscil-
lators,” Appl. Energ. 97, 771.
Zaikin, A.A., López, L., Baltanás, J.P., Kurths, J. & Sanjuán, M.A.F. [2002] “Vibrational resonance in
June 22, 2015 11:37 ”Colisesa˙IJBC˙2015 (1)”

REFERENCES 15

noise induced structures,” Phys. Rev. E 66, 011106.