The dynamics of an omnidirectional pendulum harvester

The pendulum applied to the field of mechanical energy harvesting has been studied extensively in the past. However, systems examined to date have largely comprised simple pendulums limited to planar motion and to correspondingly limited degrees of excitational freedom. In order to remove these limitations and thus cover a broader range of use, this paper examines the dynamics of a spherical pendulum with translational support excitation in three directions that operate under generic forcing conditions. This system can be modelled by two generalised coordinates. The main aim of this work is to propose an optimisation procedure to select the ideal parameters of the pendulum for an experimental programme intended to lead to an optimised pre-prototype. In addition, an investigation of the power take-off and its effect on the dynamics of the pendulum is presented with the help of Bifurcation diagrams and Poincaré sections.


Introduction
We consider a lightly damped spherical pendulum which can be forced by support excitation in three orthogonal directions and which is fitted with an active power take-off. The mathematical modelling of the spherical pendulum has been investigated in the past, and the papers by Miles are particularly noteworthy ( [18][19][20]). In the first article, the author examined the stability of planar oscillations for a spherical pendulum which was excited with small amplitudes close to the natural frequency. He found that for a particular set of parameters the planar dynamics become unstable and this results in stable non-planar solutions [19]. The later work shows bifurcations and chaotic dynamics of the system [18]. His modelling and the resulting equations laid the foundation for much of the following work. The works by Olsson [22,23] are informative because they give a clear overview of the pendulum for readers with different levels of expertise. The first experimental validation of the spherical pendulum equations derived by Miles [18,19] was given by Tritton in 1968 [27], and he concluded that the theoretical end experimental results correlate. In 1999, Aston [2] examined the spherical pendulum with a slight variation of the mathematical model to ensure that there would be no singularity at the rest position. However, the work is mainly concerned with the planar movement of the pendulum with small perturbations from the plane. The article contains an analysis of the bifurcatory responses and Lyapunov exponents of the system. Tritton and Groves [28] under-took similar research to obtain the Lyapunov exponents but used the equations derived by Miles [18,19] instead. They concluded the coexistence of alternative attractors inside and outside the range of non-stable fixed points. Furthermore, they showed that the probability of finding chaotic dynamics is strongly dependent on damping [28]. This dependence on damping is important with regards to this work, even though the energy harvester presented is lightly damped sudden switch-on of the power take-off directly corresponds to a strongly increased damping.
The dynamics of a spherical pendulum with a vertical high-frequency vibrating suspension were investigated by Markeyev [15] in 1999. Leung and Kuang [12] investigated the dynamics of the spherical pendulum which was excited close to the natural frequency in all directions. An experimental and numerical analysis of the spherical pendulum was conducted by Pospíšil, Fischer and Náprstek where they examined the influence of damping on the dynamics of the system [24]. These articles provide important insights into the different methods of excitation of a spherical pendulum.
An article on distinguishing the transition to chaos in a spherical pendulum is attributable to Kana et al. [11]. A comparison between the numerical and experimental results of the chaotic dynamics of a spherical pendulum was done by Cartwright and Tritton [6]. Similarly, Bryant quantified the amount of order in the chaotic dynamics of a damped forced spherical pendulum with the help of Lyapunov exponents [4]. These three works are essential to understand the potential chaotic dynamics of the proposed omnidirectional energy harvester. The discussion of stability and quasi-periodic resonance by Náprstek and Fischer [21] is informative in the context of the present work as it shows the softening and quasi-periodic behaviour of the spherical pendulum that is also visible in this system.
It is important to mention the paper by Litak et al. [14], which focused on the dynamic response of a spherical pendulum. The system is excited horizontally, the difference with the work of other researchers is that it is not exited with a simple sine or cosine function, but instead with a combination of them both. As a result, the trajectory of the excitation has the shape of a Lissajous curve. The authors mainly observe the bifurcation diagrams and Lyapunov exponents of their system.
In terms of energy harvesting, there are four main methods for converting energy from mechanical vibra-tion, these are piezoelectric transduction, electromagnetic induction and electrostatic transduction [9]. Electromagnetic induction is assumed for this work. Since the direction of the motion in the proposed system changes constantly, an electrical rectifier is needed to convert the AC into a usable DC, and this procedure is mentioned in several articles [7,25,29] and [8].
The work reported in [30] applied an on/off mechanical load operating on a torsional parametric oscillator, set up so that the oscillator would be loaded every half cycle, with the other half cycle as a recovery stroke. The power take-off application in [17] was devised specifically for a pendulum converter.
Previous experimental work on the general concept of the planar pendulum harvester was done by Borowiec et al. [3], Marszal et al. [16], Linang et al. [13] and Graves et al. [10]. However, Anurakpandit, Townsend and Wilson [1] noted correctly that this configuration potentially limits the ability of the system to harvest all the mechanical energy available, and so they introduced a gimballed pendulum energy harvester with two degrees of freedom (2-DOF). Interestingly in the operational cases that were shown the 2-DOF system was seen to be less effective than single DOF harvester. A justification for why the 1-DOF system has a higher peak power output than the 2-DOF system was not offered [1].
It is evident that little work has been reported on omnidirectional pendulum energy harvesters in the literature. In addition, the work on the spherical pendulum discussed in the literature to date has mostly been theoretical and mathematical and cannot easily be transferred to an experimental approach for an energy harvester. This paper will addresses this lack of information and in order to do this introduces an omnidirectional pendulum energy harvester. This is investigated numerically here, with a particular focus on optimising the operating points in order to obtain an effective harvesting of energy for the subsequent experimental validation.

Model of the spherical pendulum with active power take-off
The geometrical configuration of the spherical pendulum is given in Fig. 1  There are two necessary and sufficient generalised coordinates required to define the motion of the pendulum shown in Fig. 1, and these are given here by θ and φ. Together they describe all possible motions of the pendulum as configured. θ describes the planar oscillation of the pendulum and φ describes the rotation of the plane of oscillation of the pendulum. The length l of the pendulum is notionally defined as being from the pivot to the centre of the end mass, m. The two polar coordinates can be related to the local Cartesian coordinates by means of the kinematic relations given in Fig. 1. Two coordinate frames are used in order to provide a definition for prescribed excitations operating on the pivot, and these are given in three dimensions by the functions u(t), v(t) and w(t), as stated in Eq. (1), where U 0 , V 0 and W 0 correspond to the excitation amplitudes and Ω u , Ω v and Ω w refer to the excitation frequencies.
The excitations u(t) and v(t) will excite the pendulum directly and are simply representative of synchronous forced vibration of the system, whereas the vertical motion w(t) will excite the pendulum through principal parametric resonance if the excitation is close to an integer multiple of the linear free damped natural frequency. The most significant throughput of mechanical energy in this direction will be when that integer mul-tiple is two [5,17,26], as required for principal parametric resonance. This multiplier of two is used for the complete numerical evaluation of this work, and therefore, this work is concerned with generic forcing. The linear free damped natural frequency of the pendulum is given in Eq. (2).
After deriving the kinematic relations with respect to time, the potential (3) and kinetic (4) energies of the system can be constructed.
The potential (3) and kinetic (4) energies are substituted into Lagrange's equations of the second kind, and two ordinary nonlinear differential Eqs. (5) and (6) are obtained. In order to make the equations more descriptive of physical reality, classical linear viscous damping terms are introduced into each one. Then, a power take-off term in the angular direction of θ is added to Eq. (5) to allow conversion of energy from the system. This term is in the form of a load torque which switches on and off dependent on the sign of the velocity of the swing of the pendulum. The power take-off term is discussed in more detail in Sect. 3.1.
In order to avoid scaling effects, the following dimensionless parameters are introduced.
This leads to two dimensionless equations of motion (8) and (9) which can then be solved numerically through a suitably controlled process of numerical integration.

Numerical results
In this section, the numerical results with the nondimensional ordinary differential Eqs. (8) and (9) are examined. For the numerical calculations, if not indicated otherwise, the parameters are set to experimentally pragmatic values, as follows, In the following plots, we consider ranges for a which correspond to a = a u = a v = a w and β which correspond to β = β u = β v = 0.5β w . Likewise, the initial conditions (ICs) are set to the following dimensionless practically achievable values: These were selected after an extensive study of the initial conditions. The aim was to avoid two critical points in this calculation. Firstly, the calculation cannot be started at the equilibrium point because of the singularity at θ = 0. The other problem point would be when the initial conditions for θ andθ are too high, and in this case θ would displace through 2π and φ through π when the excitation amplitudes are also high. This already happens in the transient response and is therefore directly influenced by the initial conditions. Large angular displacements of this scale would not be physically achievable in practice and would be restricted by the necessary supporting structure. Therefore, a reasonable value is in the middle of these extreme values, and θ(0) = 0.42 was therefore arbitrarily chosen and for simplicity also applied for the other initial conditions. In the following, these values are referred to as default parameters. The system is excited in three orthogonal directions in order to model realistic conditions. Furthermore, the sampling of the points for the following plots start for a value of τ higher than 800 in order to ensure that the transient part of the response has fully decayed.

Functionality of the power take-off
The power take-off has been designed to load the pendulum every half cycle, similar to the approach taken in the concepts of [30] and [17]. In order to model this behaviour, a power take-off term is introduced into the equation of motion (5) and in the dimensionless equation of motion (8). It is important to mention that for simplicity this term does not assume any of the Fig. 2 Velocity of θ and the power take-off torque over time energy losses through the process of conversion that would occur in a real experimental setting. A suitable mathematical representation for the nominally square wave on/off power take-off uses the sign function or an arctangent switching function, such as that used by McRobb [17], for which an additional degree of tuning is introduced by means of a slight radiusing on the discontinuity in the function. The tuning parameter can generally be chosen freely between 0 and 1. However, since it was necessary to implement a power take-off that was of a form close to that was of a square wave, to maintain consistency of approach with previous work done, it would not be reasonable to select a tuning parameter close to 1 since the power take-off function would then be in the form of an arctan function rather than a square wave. The tuning parameter was therefore set to 0.01 as a compromise. Fig. 2 shows how the rounded square wave power take-off works. It can be seen that the velocity of θ and the rounded square wave are anti-phased and therefore always have opposite signs, whereby a persistent mechanical load is presented to the pendulum, with the potential for an efficient conversion of its kinetic energy. The converted mechanical energy is a function ofθ(τ ) and in order to maximise the efficiency of this conversion, high angular velocityθ(τ ) and a high take-off torque P θ are required. Therefore, in this study only excitation frequencies close to the natural frequency have been applied in order to be able to obtain high angular velocities so that the energy conversion is optimised.

Bifurcation diagrams
In order to optimise the power output, an ideal torque has to be selected which, depending on the deflection, gives an ideal operating point for the power output. Therefore, in the following, the bifurcation diagrams for the angular position θ are considered for different excitation frequencies and amplitudes to give a broad overview of the dynamical system. This is performed without a power take-off and then with a switched on power take-off. In Fig. 3, the bifurcations of the angular position θ for three different excitation frequencies β are plotted as a function of the excitation amplitude a without the power take-off. Within the bifurcation diagram for β = 1.0, a stable region exists starting at around a = 0.162. In this region, the value for θ drops drastically and the dynamics of the system become ordered. Due to the general softening property of this sort of pendulum, the oscillations at an excitation frequency 10 percent lower than the natural frequency at β = 0.9 are periodic until the excitation amplitude exceeds the value 0.18. With the highest excitation frequency at β = 1.1, it is noticeable that the stable region, which is present for β = 1.0, is no longer existing. Fig. 4 zooms in on the 3D bifurcation diagram of Fig. 3 for an area around a = 0.16 for the excitation frequency β = 1.0. This area is of particular interest since the values of θ decrease here, which results in a less effective energy conversion. For a range of a = 0.161 to 0.164, there is a region of order within the bifurcation diagram. This periodic regime comes with a drop of θ max from 2.15 rad to 1.21 rad. But even before the ordered area appears, the beginning of a Fig. 4 Bifurcation diagram for the angular position θ with the excitation amplitude (a) as control parameter. The power takeoff torque P θ is set to zero, and the excitation frequency β is set to 1.0 bifurcation forms within the chaotic looking behaviour starting at a value of a = 0.158. After the region of order disappears at around a = 0.164, an area emerges that suggest chaotic behaviour. However, within there are areas in which the points are denser (indicated by the rectangles), these areas indicate unstable periodic orbits which remain up to a value of around a = 0.168.
In Fig. 5, the bifurcation diagram for the three different excitation frequencies β and the excitation amplitude a with a switched on power take-off torque is shown. It is to be expected that the angular displacement of θ decreases with an increase in the power take-off since this directly corresponds to a greater damping. As expected, there is a large decrease in the motion of the bob. It can be seen that although the general behaviour of the three frequency ranges explored has become more periodic or quasi-periodic, see also the Poincaré sections in Sect. 3.3, bifurcations do still take place. For the excitation frequency β = 0.9, this increase in the power take-off does not show a big effect on the dynamics of the system. However, when comparing the 3D bifurcation diagram without and with a power take-off the start point where θ is initially deflected decreases from a = 0.05 to a = 0.10. The dynamics for the frequency β = 1.0 become generally more periodic up to a value of a = 0.18 when the power take-off is switched on. Also, the previously observed drop of θ for the excitation frequency β = 1.0 is visible and is starting at around a = 0.147. For the frequency β = 1.1, the dynamics drastically change when the power take-off is switched on. For a value of a from 0.1 to 0.135, the system behaves periodic. After the  Bifurcation diagram for the angular position θ with the excitation amplitude (a) as control parameter. The power takeoff torque P θ is set to 0.10, and the excitation frequency β is set to 1.0 fixed point, bifurcations occur, but they are generally more bounded compared to the bifurcations in Fig. 3. Fig. 6 zooms in on the 3D bifurcation diagram of Fig. 5 for a value of β = 1.0 and starting at an area around a = 0.147. It is noticeable that the value of the fixed point drops from θ max = 1.36 to θ max = 1.01 at a = 0.147. After the drop, 15 limit cycles occur and they merge into 5 limit cycles at a = 0.152. After the excitation amplitude exceeds the value of a = 0.155, the behaviour of the harvester changes back to a fixed point which indicates periodic behaviour.

Poincaré sections
The Poincaré section for different excitation frequencies and no power take-off is shown in Fig. 7, for the default parameters that are defined at the beginning of Sect. 3. The points were sampled at the local maxima of the excitation parameter u(τ ) [2]. In order to confirm the dynamic behaviour at β = 0.9, it is zoomed in on the area and the points are slowly converging with each excitation cycle to the right top corner. Therefore, it can be assumed that despite the long time period ( τ = 800) that is truncated at the start, the transient part of the response has not fully decayed for this set of parameters. However, a large proportion of the points can be found in the upper right corner so it can be assumed that it is a fixed point and the other points indicate the transient response that is close to completion. At β = 1.1, the system shows potential chaotic dynamics and at an excitation frequency of β = 1.0 the system behaves quasi-periodically. This shows that even small changes in the excitation frequency can change the qualitative behaviour of the system.
After introducing the power take-off, the Poincaré section in Fig. 8 changes drastically, as it is expected from Fig. 5. This means that β = 0.9 is still a fixed point, but the transient part of the response has fully decayed compared to the plot with switched off power take-off in Fig. 7. For β = 1.0 the dynamics change from quasi-periodic dynamics in Fig. 7 to a fixed point. The possible chaotic behaviour of β = 1.1 in Fig. 7 becomes quasi-periodic with a switching on of the power take-off in Fig 8. This confirms that the overall

Average power output
The average dimensionless power output P avg is a useful quantity for determining the efficiency of an energy harvester. Thus, in this section it is analysed for a variation of other parameters of the spherical pendulum. The actual power output P act in Watts can be calculated using Eq. (10). The part m l 2 ω 2 0 comes from the nondimensionalisation of the power take-off torque, and the additional ω 0 derives from the dimensionless time.
A comment on the average power output is made here; firstly, the absolute value of the power is taken; in an experiment, this would be achieved either mechanically or electrically with a rectifier, and secondly, the root mean square dimensionless power is used to determine the power output for a time period from τ = 800− 900.
In this paper, we refer to this as the average power output. The late reading of the power guarantees that the transient part of the response has fully decayed. Since the system is highly nonlinear with a singularity at θ = 0, numerical errors are inevitable. One numerical error that appeared is a stop of the calculation between τ = 800 and 900, and thus, the displayed result of P avg is lower than it actually should be. Another numerical error is where the results for the average power are disproportionately too high, but these can be filtered out. The described errors seem to be features of the numer- Fig. 9 Average power output with the dimensionless power takeoff torque P θ as control parameter for different excitation frequencies ical calculation rather than features of the dynamics of the system, but they have been minimised as far as practically possible by means of suitable filtering. The average dimensionless power output (P avg ) is plotted against the dimensionless power take-off torque P θ for different values of β and shown in Fig. 9. The other parameters are set to the default values as defined in sect. 3.1. For β = 0.9, the most power can be extracted until P θ exceeds a value of 0.16. After this, the converted power drops abruptly to a value close to zero. When the power take-off, torque gets higher than this the excitation at the natural frequency converts more power and drops to zero at a value of P θ = 0.22. The lowest power output is extracted for β = 1.1 and stops converting power at P θ = 0.24. Interestingly, the graph is somewhat different compared to those for the other two frequencies. The graph reaches a local maximum for a value of P θ = 0.22 and then drops slightly with increasing power take-off torque before it comes to the drastic drop to zero.
The reason why no more power can be converted after a certain point is that the pendulum bob movement stops because the power take-off torque is too high. It can be seen that this usually results in narrow operating areas where the energy conversion is effective. These two observations can be seen in the following plots as well.
In Fig. 10, the average dimensionless power output P avg over the dimensionless power take-off torque P θ is shown for a variation of the excitation amplitude a. The average power output hardly changes when changing the excitation amplitude until P θ reaches a value of 0.03. However, as the torque increases the graphs diverge from each other. Generally, it can be said that the bigger the excitation amplitude, the more power can be converted. However, as in the previous Fig. 9 it can be seen that the power take-off drops abruptly to zero when the power take-off torque gets to big. For an excitation amplitude of a = 0.2, power can be converted until P θ = 0.215, for a = 0.16 the power conversion ends at a value of P θ = 0.22. The least power and for the shortest period can be drawn from an excitation amplitude of a = 0.1, and here, it reaches its maximum at P θ = 0.136.
In Fig. 11, the average power output is plotted against the excitation amplitude a for different excitation frequencies β. In the beginning, the plot-lines start at a distance to each other but converge with an increase in the excitation amplitude. For the range shown, the most power can be extracted with β at 0.9 followed by 1.0 and 1.1. This corresponds with the softening properties of the pendulum. The average power P avg is plotted over the excitation frequency β for different excitation amplitudes a in Fig. 12. Generally, it can be said that the larger the excitation amplitude a is, the more power can be extracted and this is in correlation with observations of Fig. 10. Because of the softening character of the pendulum, the maximum value for P avg is not at the linear natural frequency but instead below at approximately β = 0.9. In addition, it can be seen that only an operating range of β = 0.9 to 1.1 is feasible for the chosen parameters. This confirms the previous limitation of the excitation frequencies close to the natural frequency of the system.

Characteristic plots
A overview of the entire system is given by the characteristic plots in Fig. 13. Here, P avg is plotted over the excitation amplitude a and the excitation frequency β for a variation of the power take-off torque P θ . For the lowest power take-off torque of P θ = 0.05 given in Fig. 13a, the least amount of energy can be converted. However, the working range is given over the entire area shown. If the power take-off torque is increased further, the power that can be extracted also increases for most of the areas as can be seen in Fig. 13b. The first sign of a decrease in the operating area is visible in the right bottom corner. Furthermore, it is noticeable that for a range around β = 1.0 and a = 0.15 less power can be extracted. This is already expected in the bifurcation diagrams Figs. 3, 4 and 5, as there the angle θ is smaller, and thus also less velocity is to be expected which consequently leads to less energy that can be converted. This phenomenon can also be seen in the following two plots in varying degrees of intensity.
In Fig. 13c, the power take-off torque is increased to a value of 0.15 and it shows that the converted power continues to increase as well, but it should be noted that the working range of the pendulum energy harvester is greatly reduced. This means that in regions where no power can be extracted, the excitation amplitude is too low and/or the excitation frequency is unfavourable to overcome the load applied by the power take-off torque, and as a result the pendulum bob does not move and therefore no power can be extracted out of the system. Furthermore, the area with the lower power output in the middle of the effective operating area at around a = 0.16 and β = 1.0 increased further in size.
In the last Fig. 13d, the power take-off torque is further increased to a value of 0.2. This results in a further restriction of the operation range of the energy harvester. It is no longer possible to convert energy when the excitation amplitude is too low. Likewise is it not possible any more to convert energy for the excitation frequencies below β = 0.95. This shows that the softening effect of the pendulum slowly disappears with an increase in the power take-off torque. However, within the boundaries where energy can be converted, more energy is converted compared to all the following plots. Like before, the area with the lower power output at around a = 0.17 and β = 1.025 increased further in size, but at the same time drifted slightly in the direction of the top right hand corner.

Optimisation process
As can be seen from the previous sections, the operating range and the effectiveness of the energy harvesting are highly dependent on the excitation amplitude, excitation frequency and the power take-off torque. Thus, a flow chart is presented here in order to find the optimal operation point, as shown in Fig. 14. After the start, the parameters such as initial conditions, damping, mass and length are sensibly selected based on the given system. In order to get a first overview, we will now take a closer look at the 3D bifurcations as shown in Fig. 3 and 5. In the following, a decision is made as to which parameter should be preselected for further calculation. In an experimental setting, these parameters would be dictated by the design of the construction and by the level of excitation at which the system operates. After Characteristic plots for β over a for the average power output P avg for different power take-off torques P θ the decision on the parameters is made, two sub-loops are formed. On the left sub-loop, the power take-off torque has been chosen. This, for example, is the case when the load generator has a pre-defined torque that it will apply. In order to find the ideal operation point for the selected torque, we need to plot the P avg over a given in Fig. 11 and P avg over β as can be seen in Fig. 12. From these plots, optimal excitation amplitudes and frequencies are chosen.
The right sub-loop can be chosen if the excitation amplitude and excitation frequency of the systems are known. Then, P avg is plotted over P θ referring to Fig. 9 and 10 for different excitation values close to the value previously selected, and this results in the selection of an ideal torque. In the following, it is checked whether it is iterated enough in order to get optimal operating parameters. If another iteration cycle is needed, this decision gives the opportunity to jump back to the selection of the starting conditions. When the process is complete, the ideal points for the system can be selected. These are then compared with the characteristic plots, see Fig. 13 for details, in order to double-check the results and avoid errors, so it is possible to return to the beginning. If everything is correct, the process stops and so the experiments can be initiated with the selected values.

Conclusions
The aim of this paper consists of two parts: firstly, it has been the aim to provide a theoretical analysis of a omnidirectional pendulum energy harvester in order to remove the limitations that previous simple pendulum harvesters have operated under. Additionally, the second aim of this work has been to propose an optimisation procedure to select the ideal parameters of the pendulum for an experimental programme intended to lead to an optimised pre-prototype.
The energy harvester presented here has then been investigated for its dynamic behaviour using bifurcation diagrams and Poincaré sections. In general, the system behaves more periodically when power is extracted from it. The amount of power that can be realised has been investigated by plotting the average power for various parameters of the pendulum. These plots show the influence of the power take-off torque, excitation amplitude and excitation frequency on the average power output. Additionally, the results also show that the operating ranges strongly depend on these variables and that the energy conversion can stop abruptly if the parameters are not chosen with caution. In order to show the operating ranges of the omnidirectional pendulum, energy harvester characteristic diagrams have been created. In general, it can be said that the average power output increases with an increasing torque, but the operation range decreases likewise. Finally, a flow chart has been created that illustrates the selection of an ideal operating point depending on the predefined conditions.
Further research will focus on the experimental side of the problem. An analytical evaluation of the problem using the asymptotic approach of the perturbation method of multiple scales is underway.

Conflict of interest
The authors declare that they have no conflict of interest.

Ethical standards
The authors ensure the compliance with ethical standards for this work.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/ by/4.0/.