Potential well hopping and performance of ocean energy harvesters

Journal of Sound and Vibration 465 (2020) 115008

Contents lists available at ScienceDirect

Journal of Sound and Vibration journal homepage: www.elsevier.com/locate/jsvi

Potential well hopping and performance of ocean energy harvesters Dane Sequeira ∗ , Brian P. Mann Dynamical Systems Research Laboratory, Department of Mechanical Engineering & Materials Science, Duke University, Durham, NC, 27708, USA

article info

abstract

Article history: Received 19 April 2019 Revised 21 September 2019 Accepted 6 October 2019 Available online 10 October 2019 Handling Editor: L. Huang

This paper examines a gimballed horizontal pendulum and studies how it can be used to harvest energy from ocean waves. The system’s behavior is analyzed under both static and dynamic scenarios. Threshold escape behavior for parametrically excited systems with a time dependent term in the potential energy function is discussed and a criterion is proposed for predicting escape events. Performance metrics are identified to quantify and compare different responses. Numerical and experimental studies are conducted showing how the system can be designed for enhanced performance by altering geometric parameters to suit various excitations. The system’s response to both deterministic single frequency and stochastic multi-frequency excitations is investigated. Design implications are discussed. © 2019 Elsevier Ltd. All rights reserved.

Keywords: Nonlinear energy harvester Threshold escape criterion Time dependent potential energy function Stochastic multi-frequency excitation source

Since the invention of the waterwheel in ancient times, scientists and engineers have sought to harvest the vast amounts of energy available in the natural environment for useful purposes. Solar panels, ocean buoys, and wind turbines are just a few examples that have become commonplace in modern society. Among the various different types of harvesters, vibratory energy harvesters specifically seek to extract the kinetic energy available in surrounding environments and convert that into electricity. Historically, the most common way of designing one of these devices has been to tune the harvester so that its natural frequency matches the frequency of the environmental excitation [1]. This idea is based on fundamental concepts in linear vibration theory, however it is not particularly helpful for designing a system in cases where the excitation frequency is either unpredictable or wide ranging as is often the case in environmental excitation sources [2]. Ocean waves present a particularly intriguing environmental excitation source because of the magnitude and density of energy that they possess, but are generally characterized by broad frequency spectrums that render linear designs sub-optimal [3]. Numerous nonlinear designs have been studied in the literature to better understand and ultimately address this difficulty. Dostal et al. investigated the stochastic nature of ocean wave excitations by analyzing a traditional pendulum that was vertically excited by random heave translations [4]. Yurchenko et al. studied a similar heave driven pendulum but found that performance could be improved by orienting the pendulum on an inclined angle to align the system’s natural frequency with the excitation [5]. Sah and Mann analyzed the response of an inclined pendulum under rotational rolling excitation that mimics the roll of ocean buoys as opposed to a translational heave. This work also showed how the horizontal pendulum’s natural frequency approaches zero as the system’s angle of inclination goes to zero [6]. Finally, Ding et al. combined many of these studies by exploring how excitation and natural frequencies could be matched for a roll-driven underwater mooring platform [7]. Collectively, these studies form a cohesive picture of how ocean energy harvesters have traditionally been designed to extract energy from waves.

∗ Corresponding author. E-mail address: [email protected] (D. Sequeira).

https://doi.org/10.1016/j.jsv.2019.115008 0022-460X/© 2019 Elsevier Ltd. All rights reserved.

2

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

Fig. 1. Image showing both the system’s rotational degrees of freedom and its geometric parameters.

In a recent paper, Sequeira et al. proposed inserting an additional rotational degree of freedom (i.e. gimbal) into the classic horizontal pendulum to create an energy harvesting device that triggers dynamic bifurcations as a function of excitation amplitude rather than frequency [8]. This feature allows the system to operate well under conditions where the amplitude of environmental excitation is known even if the frequency is highly unpredictable. Although the findings showed that the system’s response depends on complex relationships between the geometric parameters, electrical circuitry, and environmental conditions, the higher level phenomenological and design implications were left unstudied. This paper addresses those questions to make general statements about energy harvesters that are designed for threshold escape behavior rather than matching frequencies and applies those conclusions to the gimballed horizontal pendulum as a case study. The remainder of this paper is organized as follows. First, an overview of the gimballed horizontal pendulum is given to describe the system’s geometry as well as some general static and dynamic behavioral trends. This largely recaps findings in the previously published paper and sets the groundwork for the remainder of investigations. Next, the system’s potential and kinetic energies are described and a criterion is developed to determine when potential well hopping will occur. This threshold breakthrough criterion applies not only to this particular system, but to parametrically excited systems with time-dependent potential energy terms in general. A metric is also developed to quantify performance and allow for comparisons across different simulations. Next, these insights are used to develop design considerations and show how geometric system parameters can be adjusted to maximize performance. Both numerical and experimental results are presented. Finally, these results are generalized from a single frequency excitation to multi-frequency stochastic ocean waves. A model is developed relating ocean waves to an excitation input, parameter sweeping simulations are conducted, and high level behavioral characteristics are analyzed. 1. System overview This section gives an overview of the gimballed horizontal pendulum system and identifies some of the key phenomena that the remainder of this paper will expand upon. Fig. 1 shows two different views of the gimballed horizontal pendulum which can be most easily described as two individual bodies, a disc of mass M and a point of mass m, that move under three consecutive rotational degrees of freedom as illustrated in Fig. 1a. The first rotation, 𝜓 , represents a parametric excitation which takes the form of a simple sinusoid. In an ocean energy harvester, this excitation would correspond to a rolling motion of the harvester riding the waves and will be referred to throughout the remainder of this paper as the “tilt angle.” The first portion of this paper will consider a single frequency sinusoid of the form 𝜓 (t) = A sin 𝜔t but will be expanded in the later sections to include multiple frequencies and phase shifts. The second rotation, 𝜃 , describes rotation about the gimballed degree of freedom. The third rotation, 𝜙, represents rotation of the point mass about the harvester’s shaft and is the degree of freedom that would be directly linked to an electrical circuit for harvesting energy. Because the electromechanical transducers that can be used to convert this system’s kinetic energy to useful electricity generally rely on shaft rotation rate, this paper will largely focus on the relationship between 𝜙̇ and 𝜓 (t). The rotations can be mathematically described using the following rotation matrices:

⎡1 ⎢

R(𝜓 ) = ⎢0

{ where O =

⎢ ⎣0

0

0 ⎤

cos 𝜓

sin 𝜓 ⎥ ,

− sin 𝜓 }

î O , ˆjO , ̂ kO

⎥

⎥ cos 𝜓 ⎦

⎡cos 𝜃 ⎢

R(𝜃 ) = ⎢ 0

⎢ ⎣ sin 𝜃

0 1 0

− sin 𝜃 ⎤ ⎥ 0 ⎥, ⎥ cos 𝜃 ⎦

⎡ cos 𝜙 ⎢

R(𝜙) = ⎢− sin 𝜙

⎢ ⎣

sin 𝜙

0⎤

cos 𝜙

0⎥ ,

0

0

⎥

⎥ 1⎦

(1)

represents an inertial reference frame as shown in Fig. 1b and rotations in 𝜓 , 𝜃 , and 𝜙 correspond to

{

three subsequent frames of rotation A =

}

î A , ˆjA , ̂ kA , B =

{

î B , ˆjB , ̂ kB

}

{

and C =

}

î C , ˆjC , ̂ kC , respectively.

Position vectors to the disc and point mass, respectively, can be expressed as rM ∕O = Ĥ kA − b̂ kB

rm∕O = Ĥ kA − 𝓁̂ kB + rˆjC

(2)

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

3

Fig. 2. (a,b) Images showing sample static equilibria as the tilt angle is varied. (c) Three dimensional bifurcation diagram highlighting the critical point 𝜓 = 𝜉 s where the system transitions from having two coexisting stable equilibria to a single trivial one. Solid and dashed lines denote stable and unstable equilibria, respectively (d) Quasi-static numerical simulation (i.e. 𝜔 is small). illustrating how the system goes from small oscillations about a single attractor to continuous rotational motion across stable attractors when A ≥ 𝜉 s . A black dashed vertical line has been superimposed at the point A = 𝜉 s for reference.

where H represents the distance from the axis of excitation to the gimbal axis, b describes the distance from the gimbal axis to the disc’s center of mass, 𝓁 is the vertical distance from the gimbal axis to the point mass, and r represents the horizontal distance from the axis of shaft rotation to the point mass. These geometric parameters are shown in Fig. 1b. Previous work by Sequeira et al. used these descriptions to derive this system’s governing equations of motion which can be found in that article [8]. That paper then used those equations to conduct both a static and dynamic analysis of the system with results shown in Fig. 2. In doing this analysis, a few key insights were identified that will serve as the foundation for this paper moving forward. First, a bifurcation point was discovered that occurred for sufficiently large tilt angles where two neighboring stable equilibria would collide at an intermediate saddle node and give birth to a new, lone stable equilibrium as illustrated by Fig. 2a and b. Fig. 2a shows the system at a tilt angle of 𝜓 = 0 in one of its two stable orientations. Although the other position is not pictured, the reader can visualize how a symmetric stable orientation ( should) exist. This trend continues as the tilt angle

̃ ±𝜙̃ → (0, 0) where 𝜃̃ and 𝜙̃ denote equilibria is increased with the coexisting stable equilibria both approaching zero ±𝜃,

positions in 𝜃 and ( 𝜙,)respectively. At some critical tilt angle |𝜓 | = 𝜉 s , these coexisting solutions converge onto a single trivial

̃ 𝜙̃ = (0, 0) as shown in Fig. 2b. The bifurcation point at which this transition occurs was solved analytically equilibrium at 𝜃, by Sequeira et al. as −1

𝜓 = 𝜉s = cos

[(

𝓁 m + bM mr

)

] sin 𝜉s

(3)

Fig. 2c illustrates these findings through a three dimensional bifurcation diagram that was constructed by varying the tilt angle 𝜓 and calculating the corresponding equilibrium angles 𝜃̃ and 𝜙̃. This figure shows how the system alternates between having multiple and a single stable attractor as a function of 𝜓 .

4

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

Numerical simulation reveals the implications of this finding. Fig. 2d shows the results from a simulation where the rotational amplitude A of the parametric excitation 𝜓 (t) = A sin 𝜔t is linearly increased during the course of the simulation. Values for the response 𝜙 are shown on the range [−𝜋, 𝜋 ] by utilizing the modulo function and a vertical dashed line is superimposed at the point A = 𝜉 s . As shown, the system undergoes a qualitative change in dynamic behavior when the excitation amplitude reaches the bifurcation point. In particular, the system transitions from exhibiting small oscillations about a single stable attractor to continuous rotational motion across stable attractors. While the example shown here is for a quasi-static simulation (i.e. 𝜔 is small), further cases were studied and showed that a dynamic bifurcation point, 𝜉 d , can be shifted for a given geometric configuration by varying either the excitation frequency or the damping. These numerical studies were validated experimentally with results showing good agreement. This work laid a useful foundation for studying the gimballed horizontal pendulum and its ability to act as a vibrational energy harvester that can be designed for threshold escape behavior rather than the traditional method of matching frequencies. The sections that follow will build on that analysis to make general statements about threshold escape behavior in energy harvesters and apply those insights to a case study involving the design of ocean energy converters. 2. Energy investigation This section focuses on the dynamic relationship between kinetic, potential, and dissipative energy in the gimballed horizontal pendulum to better understand how it responds to various excitations. In particular, a criterion is developed to predict the onset of threshold escape behavior, and a performance metric is proposed to assess how this type of motion affects the system’s performance as an energy harvester. This is important because it will determine whether threshold escapes are desirable in this energy harvester and also inform how to either encourage or avoid them.

2.1. Threshold escape behavior In order for interaction to occur across potential wells (commonly called “potential well hopping”), the system must have sufficient energy to overcome the barrier which exists between neighboring stable solutions and the unstable saddle point separating them. Consequently, there are two key factors that must be considered in determining whether a potential well escape can occur: (1) the height of the potential barrier Ub , and (2) the amount of the energy, E, in the system at any given time as a sum of both kinetic and potential energy. Focusing first on the height of the potential energy barrier, Fig. 3 utilizes contour [ ] lines to show how the structure of the system’s potential energy evolves as the tilt angle is varied on the range 𝜓 = 0, 𝜉s . When 𝜓 is small, the system has two coexisting potential energy wells that are separated by a saddle point. As ( 𝜓 increases, )

̃ 𝜙̃ = (0, 0) the difference in energy levels between these points decreases until finally the stable positions converge onto 𝜃,

for 𝜓 ≥ 𝜉 s , where there is no longer any energy required to transition between stable points because they are coincident. Remembering that the system oscillates in time as 𝜓 (t) = A sin 𝜔t, there are two options available to the system as the tilt angle decreases back towards 𝜓 < 𝜉 s and the potential wells separate: either the system will go back into the well that it started in or it will transition into the neighboring well. Which of these two outcomes will occur depends on the phase of the response relative to that of the excitation, but in either scenario, the height of the potential barrier will increase until 𝜓 = 0 and the cycle starts again. Therefore, the potential energy barrier varies periodically as a function of 𝜓 and reduces to zero when 𝜓 ≥ 𝜉s. Turning to the amount of energy in the system at any given time, it is important to understand how damping and excitation frequency affect the way systems absorb energy from parametric excitations. For the quasi-static case, the system will respond by slowly transitioning along stable solutions and the only energy exchange that will occur is in the change in potential energy associated with changing equilibrium positions, which may be zero. If the frequency of oscillation is increased, the system may begin to oscillate about its equilibria as it absorbs energy from the excitation in the form of kinetic energy. If damping is introduced, then the system will lose some of this energy in the response as it is dissipated to the surrounding environment. Therefore, the system will be able to extract the largest amount of energy from a parametric excitation when the frequency is large and the damping is small. Combining these insights on potential barriers and energy absorption with the results in Fig. 2 illustrates why there are two distinct types of motion that the system undergoes (i.e. small oscillations about a single stable attractor and continuous rotational motion across stable attractors) and how it is able to transition between these states under different parametric excitations. The system bifurcates in a quasi-static simulation when the amplitude is allowed to reach the critical point even though there is no kinetic energy in the response because the energy required to do so becomes zero. As the excitation frequency is increased, this bifurcation occurs at lower excitation amplitudes because the energy in the response increases and at some point becomes large enough to overcome the nonzero potential barrier. Contrastingly, when the damping is increased this bifurcation occurs at larger excitation amplitudes because some of the energy in the response is being lost through dissipation. This can be expressed by introducing a dynamic bifurcation point 𝜉 d , which shifts relative to the static bifurcation point based on the particular set of frequency and damping values being examined. This qualitative understanding of why potential well escapes occur at different amplitudes for different combinations of frequency and damping is interesting, however it would be more useful to be able to predict when they will occur. Numerous methods have been used in the literature to study potential well escapes for a broad class of systems. Virgin et al. considered har-

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

5

Fig. 3. Plot of potential energy as a function of 𝜓 showing how a saddle node has collided with neighboring stable equilibria to give rise to a single stable equilibrium when 𝜓 ≥ 𝜉 s . The red cross denotes the unstable saddle point and black dots mark the stable equilibria. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

monically forced systems and determined the approximate conditions for potential well escape by applying harmonic balance to obtain approximate analytical solutions and comparing the energy in the response to the potential barrier [9]. Mann built upon this method by studying a parametrically excited system and showing that this escape criterion needed to be modified to accurately predict potential well escapes under parametric excitation [10]. Specifically, he found that the parametric excitation term should be included in the calculation of the kinetic energy to determine the maximum energy in the system throughout excitation. The gimballed horizontal pendulum system differs from both of these examples due to the fact that the potential energy function includes a time varying term. Specifically, the potential energy U(𝜃, 𝜙, t), takes the form

[

]

U (𝜃 (t ), 𝜙(t ), t ) = g mr cos 𝜙(t ) sin[A sin 𝜔t ] + cos[A sin 𝜔t ] (H(m + M) − (𝓁 m + bM) cos 𝜃 (t ) + mr sin 𝜙(t ) sin 𝜃 (t )) .

(4)

The potential barrier can be calculated from this by evaluating the potential energy function at the stable and unstable equilibrium positions and subtracting. Sequeira et al. calculated the stable equilibrium positions as functions of the tilt angle and found that unstable saddle points exist at (𝜃 = 0, 𝜙 = n𝜋 ) [8]. Using this information, the potential barrier can be expressed as

̃ 𝜙, ̃ t) Ub (t ) = U (0, n𝜋, t ) − U (𝜃,

[ ( ) ( )] ̃ sin[A sin 𝜔t] + cos[A sin 𝜔t] −𝓁 m − bM + (𝓁 m + bM) cos 𝜃̃ − mr sin 𝜙̃ sin 𝜃̃ . = g −mr ±1 + cos 𝜙

(5)

This potential barrier function will oscillate periodically while the trough decreases for increasing amplitude until eventually it reaches zero at A = 𝜉 s . Fig. 4a illustrates this result by showing how the instantaneous potential barrier oscillates between a constant maximum value at 𝜓 (t) = 0 and a varying minimum value at 𝜓 (t) = A. This represents an added complexity over the previously presented potential well escape examples which is not only true of the system being studied here but can be generalized to any parametrically excited system that has time dependent terms in the potential energy function. ̇ 𝜃, 𝜃̇ ), during amplitude sweeps for varying excitation Fig. 4b expands on this analysis by tracking the system’s energy E(𝜙, 𝜙, frequencies. Using a similar metric to Mann, a potential well escape should be predicted when condition E ≥ Ub is satisfied [10]. As expected, the system experiences a potential well escape at smaller values of A as 𝜔 increases. This figure also shows that the

6

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

Fig. 4. (a) Plot showing how the potential barrier oscillates in time while the trough decreases with increasing amplitude where the thick solid line denotes the minimum potential energy barrier, the thin solid line indicates the instantaneous barrier, and the thick dashed line marks the maximum potential energy barrier. (b) Plot illustrating how the system is allowed to escape from its potential well when condition E ≥ Ub is satisfied. Red lines denote the system’s energy throughout simulation, the black line indicates the minimum potential energy barrier, and the grey area marks space where escapes have occurred. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

intra-well response becomes increasingly more energetic as it approaches the potential energy barrier, which is consistent with Mann’s study. This result expands on previous work by offering insight on the nature of potential well escapes in parametrically excited systems and provides a useful metric that can aid in design considerations for the gimballed horizontal pendulum as an energy harvester. 2.2. Energy harvester performance In order to compare how the harvester performs under various conditions and use that analysis to inform design, it is important to define what constitutes a desirable outcome and link that to a quantifiable metric. Knowing that 𝜙 is the degree of freedom that would be directly linked to an electrical circuit and that shaft rotation rate squared is proportional to the power transferred to an electrical circuit (regardless of whether electromagnetic inductance or electrostatic capacitance is being used), it makes sense that a performance metric would take the form

Γ=

1 ̇2 c 𝜙 2 e

(6)

where Γ denotes performance and ce represents the effective electrical damping coefficient that couples the mechanical harvester with an electrical circuit. Because this paper is focused on the mechanical system rather than the specifics of the electromechanical transducer, the leading coefficients in this expression can be neglected so that performance is measured as Γ = 𝜙̇ 2 . One additional complexity that needs to be considered is how to account for the fact that 𝜙̇ is constantly changing in response to the parametric excitation. While peak rotation rate and other metrics can offer useful insight, this paper will focus on the mathematical average over a period of time as

Γ = 𝜙̇ 2 =

1

t2

t 2 − t 1 ∫t1

𝜙̇ 2 dt

(7)

where 𝜙̇ 2 represents the time averaged shaft rotation rate squared from an initial time t1 to a sufficiently later time t2 which was determined through experimenting with various incremental values. While there are numerous ways to justifiably measure harvester performance, this is the performance metric that the remainder of the paper will use. Until this point, it has been assumed that potential well hopping should be encouraged for the energy harvester because these escapes are commonly associated with large velocities and amplitude displacements that consequently lead to an increased amount of energy being scavenged from the environment [9,10]. This assumption can now be tested by applying the performance metric to the numerical results discussed earlier and comparing whether the harvester performs better in a state of small oscillations about a potential well or one in which potential well escapes are occurring. These results are shown in Fig. 5. Fig. 5a shows the results for an amplitude sweep similar to the one discussed previously where a vertical black dashed line has been superimposed over the point A = 𝜉 s and the system clearly undergoes dynamic bifurcation at approximately this point. Fig. 5b and c, respectively, show the instantaneous and time averaged shaft rotation rate squared. These figures clearly indicate that the performance of the harvester improves drastically when the system is allowed to transition across stable attractors rather than

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

7

Fig. 5. Plot showing the performance of the harvester as measured by 𝜙̇ 2 during an amplitude sweep for the slow wave, low powered circuit case showing that there is a large increase in kinetic energy for continuous rotational motion. A vertical dashed line has been superimposed at A = 𝜉 s for reference.

(

Fig. 6. Diagrams showing how the critical tilt angle changes as a function of various geometric parameters. (a,b,c) Results showing how 𝜉 s varies on the interval 0, a function of r, m, and 𝓁 , respectively.

𝜋 2

) as

oscillating about a single attractor. Furthermore, these plots highlight the importance of relying on the time averaged results to gauge performance as the large fluctuations in the instantaneous results can lead to misleading interpretations. 3. Design study (single frequency) This section uses the analysis conducted previously to design an energy harvester that operates well under a set of single frequency excitation conditions. Because the relationship between the excitation amplitude and the static bifurcation point has been shown to have a significant impact on performance, design considerations will center around how the system’s geometry can be altered to achieve a situation where the static bifurcation angle is smaller than the excitation amplitude. This can be accomplished by manipulating the geometric parameters in Eq. (3) to suit the excitation source and thereby dictate the type of response that occurs. Fig. 6 shows how the critical tilt angle varies as a function of three different geometric parameters: r, m, and 𝓁 . Dashed 𝜋 horizontal lines have been superimposed at 𝜉 s = 0 and to illustrate the two extreme scenarios in which the system behaves 2 like a simple pendulum rotating in only one of its two degrees of freedom. Concentrating on Fig. 6a for example, if r is small, then the system acts like a simple pendulum rotating in 𝜃 about the point 𝜃 = 0 with little regard for the coupling effects from 𝜙. Alternatively, if r is large, then the system acts like a pendulum rotating in 𝜙 while being held at the fixed orientation 𝜃 = − 𝜋2 . Between these two extreme cases, the system transitions through points where 𝜉 s takes on non-trivial values and the system’s coupling becomes prevalent, leading to more complex dynamic behavior. These results are illustrated through Fig. 6a, b, and c and highlight how the system can be tuned using r, m, and 𝓁 , respectively, to alter 𝜉 s . Leveraging these insights, simulations were conducted to illustrate how the performance of the harvester can be varied by adjusting the geometric parameters in accordance with the excitation source. Fig. 7 shows the results from simulation for two cases with identical excitation sources, A = 0.62 rad and 𝜔 = 0.63 rad/s, but using different geometric parameters, specifically different values of r. The top row of each of these cases shows the instantaneous value of 𝜙̇ 2 throughout the simulation while the bottom two rows show the phase space responses. The phase portraits for Fig. 7a clearly indicate that the system is being confined within a single basin of attraction as neither 𝜃 nor 𝜙 is allowed to cross over the zero point in the phase diagram and instead trace out classic forced trajectories about a stable attractor. This notion was further confirmed through an amplitude

8

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

Fig. 7. Computational results using the excitation parameters A = 0.62 rad and 𝜔 = 0.80 rad/s for cases where the system (a) oscillates about a stable attractor and (b) transitions across stable attractors demonstrating how the latter performs approximately 330% better.

sweep indicating that 𝜉 d ≈ 0.74 rad for this combination of excitation frequency and damping so that A < 𝜉 d . The results from Fig. 7b, however, show more interesting behavior. By simply reducing the radius from r = 0.04 m to r = 0.03 m and consequently reducing the dynamic bifurcation point to 𝜉 d ≈ 0.57 rad so that A > 𝜉 d , the system is allowed to jump freely across basins of attraction as it interacts simultaneously with both stable attractors. Results confirm that this case performs significantly better than the previous one as indicated by the approximately 330% increase in 𝜙̇ 2 from 2.13 to 9.16 (rad/s)2 . To confirm that these findings hold true beyond simulation, studies were conducted to validate them experimentally. Fig. 8 shows the setup that was constructed and used for experimental testing. This system was driven in 𝜓 by an external motor at a specified amplitude and frequency and the response in both 𝜃 and 𝜓 were measured using potentiometers. Because of the practical difficulty in reducing the radius without affecting the mass, both parameters were varied slightly to achieve the overall goal of shifting the system’s dynamic bifurcation point. Fig. 9 shows the results obtained using the excitation parameters A = 0.19 rad and 𝜔 = 3.14 rad/s. Looking at the phase portraits clearly illustrates that the system is confined to small oscillations about the nearest equilibrium point for the case where 𝜉 d > A whereas it is allowed to achieve continuous rotational motion when the dynamic bifurcation point is reduced so that 𝜉 d < A with an approximately 670% associated increase in performance. This result proves that the motion, and consequently the energy extracted, from this system can be significantly altered by adjusting simple geometric parameters to best suit the environmental excitation source. At this point it is useful to compare the performance of the gimballed horizontal pendulum to similar energy harvesters to evaluate both advantages and disadvantages. The two most obvious comparisons are between the roll driven vertical and horizontal pendulum system’s shown in Fig. 10. These are natural comparisons because they appear in the literature and also because of the way they mirror the behavior of the gimballed horizontal pendulum, which behaves like a vertical pendulum when 𝓁 → ∞ and a horizontal pendulum when r → ∞ as evidenced by Fig. 6. Numerical simulations were conducted using the same corresponding geometric and excitation parameters tested previously (i.e. r, 𝓁 , m, H, A, and 𝜔) and the resulting

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

9

Fig. 8. Image of the gimballed horizontal pendulum used for experimental testing.

Fig. 9. Experimental results using the excitation parameters A = 0.19 rad and 𝜔 = 3.14 rad/s illustrating that behavior can be manipulated by changing geometric parameters while maintaining the same wave input to achieve an approximately 670% increase in performance.

performance was measured to be 𝜙̇ 2 = 0.12 and 14.52 (rad/s)2 for the vertical and horizontal pendulums, respectively. The first thing to note is that the vertical pendulum performs substantially worse, which occurs largely because of the misalignment of excitation and natural frequencies. Due to the relationship between a vertical pendulum and its linear natural frequency, √ 𝜔n = g∕𝓁 , these harvesters must be designed unrealistically long to elicit any significant response for characteristically low frequency ocean waves [11]. Contrastingly, the horizontal pendulum actually performs noticeably better for the same conditions. This is likely due to the horizontal pendulum’s interesting characteristic where both the linear natural frequency and potential energy barrier approach zero as the tilt angle goes to zero [6]. Despite this theoretically improved performance, the horizontal

10

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

Fig. 10. Images of roll driven (a) vertical and (b) horizontal pendulums as studied in the literature.

Fig. 11. Plot showing the performance of the horizontal pendulum as a function of an offset tilt angle to illustrate the disadvantage of a buoy which floats at a skewed equilibrium orientation. A dashed line has been included to indicate the gimballed horizontal pendulum’s performance under a similar excitation for comparison.

pendulum will never actually achieve this performance in an ocean wave setting due to the fact that the asymmetrical mass distribution causes any floating body to equilibrate about an offset angle which dampens out regular buoy oscillations [8]. This was observed experimentally by Sequiera et al. and largely motivated their first investigations into the gimballed horizontal pendulum. Fig. 11 shows how the performance of the horizontal pendulum worsens as the offset angle increases and eventually dips below the gimballed horizontal pendulum at an offset angle of approximately 0.26 rad. While this result applies specifically to the parameters used in this study, it highlights why the horizontal pendulum is generally unsuitable for use in ocean energy harvesters. In summation, the gimballed horizontal pendulum offers unique advantages over both the vertical and horizontal pendulums for harvesting energy from low frequency ocean waves. 4. Design study (ocean waves) This section expands the analysis conducted thus far for single frequency deterministic excitations to consider ones that are both multi-frequency and stochastic in order to better emulate a real ocean environment. As in the single frequency case, the goal is to analyze the harvester under various conditions to determine how it should be designed to maximize performance. There are three steps that need to be taken in order to make this extension: (1) a model must be constructed to describe sea surfaces under various conditions, (2) a method must be determined for converting sea surface to the parametric excitation ̈ (t), and (3) the parameter space must be studied under different ocean conditions to identify desirable terms 𝜓 (t), 𝜓̇ (t ), and 𝜓 design criteria.

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

11

Fig. 12. (a) Pierson-Moskowitz spectrum for various wind speeds where black dots denote the peak value for each wind speed and dashed lines indicate the maximum and minimum frequencies being considered. From top to bottom the wind speeds used correspond to U19.5 = 20, 17.5, 15, 12.5, and 10 m/s. (b) Example of a sea surface across space at an instant in time for U19.5 = 10 m/s.

4.1. Generating sea surface states The Pierson-Moskowitz spectrum is a well known method for determining sea surfaces that assumes a fully developed sea (i.e. winds have been blowing at a constant rate for a long period of time) where the wind is the only external force generating waves [12]. It approximates the wave spectral density S, as a function of wave frequency and the wind speed at 19.5 m above the surface of the water U19.5 , as

[

(

𝛼 g2 g S(𝜔, U19.5 ) = 5 exp −𝛽 𝜔 𝜔U19.5

)4 ]

(8)

[

]

where 𝛼 = 8.1 × 10−3 , 𝛽 = 0.74, g = 9.81 m/s2 is the gravitational constant, and 𝜔 ∈ Ω1 , Ω2 represents the range of excitation frequencies being considered in the spectrum. Fig. 12 shows the Pierson-Moscowitz spectrum for five different wind speeds where black dots have been included to denote the peak values for each individual spectrum. From this spectrum, the sea surface W(x, t), can be calculated as a spatially and temporally varying function by dividing the frequency spectrum into N discrete frequencies with corresponding amplitudes and then summing them in a Fourier Series as W ( x, t ) =

N ∑

(

ai sin ki x + 𝜔i t + 𝛾i

i=1

√

where ai (𝜔i , U19.5 ) =

2S(𝜔i , U19.5 )

[

)

Ω2 −Ω1 N −1

(9)

] , ki =

𝜔i 2 g

for waves in deep water (as will be assumed throughout the remainder

of this paper), and 𝛾 i is a randomly chosen phase shift associated with each individual frequency to give the model its stochastic nature. Fig. 12b shows one realization of this process for U19.5 = 10 m/s, however infinitely many possible realizations exist. 4.2. Multi-frequency stochastic excitation source To relate information about the ocean waves to the current study, a model must be established for converting sea surface states into parametric excitation of the harvester. One reasonable method that has been used in the literature is to assume that the device is mounted to a platform or buoy that floats in the water and remains tangent to the surface of the waves as they pass underneath it (i.e. tangent plane assumption) as shown in Fig. 13a [13]. This assumption loses validity near the troughs of waves in circumstances where the length of the platform, D, is large relative to the wavelength, but can be improved by enforcing a constraint on the minimum wavelength (or similarly the maximum frequency) to be included in the Pierson-Moskowitz sea surface formulation for a given value √ of D. Fig. 13b illustrates how this issue is mitigated throughout the remainder of the paper by applying the constraint Ω2 = 𝜋 g∕D to the maximum frequency, which is the equivalent of stating that the minimum allowable wavelength is 8 times the length of the floating platform. Combining Eq. (9) with the tangent plane assumption, ocean waves can be converted into parametric excitation by first computing the spatial derivative of the sea surface

( ) 𝜕 W ( x, t ) ∑ = ai ki cos ki x + 𝜔i t + 𝛾i . 𝜕x i=1 N

(10)

12

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

Fig. 13. (a) Illustration showing how the wave characteristics gathered from the Pierson-Moskowitz spectrum were converted into a rotational amplitude for prescribed motion. (b) Zoomed in image showing the importance of including a maximum wave frequency for a given buoy platform diameter.

Using the deep water assumption described previously to get ki in terms of 𝜔i and the point source approximation which assumes the platform remains at a fixed location x = 0, this equation can be simplified to.

( ) 𝜕 W (0, t) ∑ 𝜔i 2 = ai cos 𝜔i t + 𝛾i 𝜕x g i=1 N

(11)

which represents the slope of the line tangent to the sea surface as a function of time. Finally, the resulting tilt angle can be calculated from this line using trigonometry as −1

𝜓 (t) = tan

[N ∑ i=1

ai

𝜔i 2 g

(

cos 𝜔i t + 𝛾i

)

]

.

(12)

̈ (t) (see Appendix). Fig. 14 Time derivatives can be taken from this equation to get similar expressions for both 𝜓̇ (t ) and 𝜓 shows an example of a time series generated from sea surface states that takes U19.5 as an input to produce one randomized ̈ (t ). representation of the resulting parametric excitations 𝜓 (t), 𝜓̇ (t ), and 𝜓 4.3. Results This section studies the parameter space in order to determine how the system should be designed to maximize performance under various ocean conditions. Similarly to the single frequency case, simulations are conducted in which the excitation source is swept in time and the corresponding performance is measured. To further explore the design space, this process is repeated for varying geometric parameters to identify how performance varies both as a function of environmental conditions (i.e. wind speed) and design parameters (i.e. 𝓁 , r, and m.) These simulations are run numerous times and then averaged to account for the stochastic nature of the ocean waves injected into the simulation through randomized 𝛾 terms. Fig. 15 illustrates how the excitation source is swept during simulation. This can be achieved by combining Eqs. (12), (14) and (17) with the expression for ai (𝜔i , U19.5 ) and allowing the wind speed to vary linearly throughout the simulation as U19.5 = Ur t where Ur is the desired sweep rate. This rate should be kept small to reduce unwanted transients in the simulation results. For the following results, a sweep rate of Ur = 5 × 10−4 m/s was used with a time period of 20, 000 s. Random values of 𝛾 are assigned at the beginning of the simulation and maintained throughout the sweep. Fig. 15a helps visualize this process by showing how the intensity function is reevaluated during simulation as the wind speed is increased and then broken into N discrete frequencies with corresponding intensity values. Fig. 15b shows one example of how the excitation parameters change during a sweep. This sweep, and all of the other ones conducted in this paper, were done using N = 50. Performance results vary for individual frequency sweeps depending on the randomly chosen values for 𝛾 as illustrated through the 10 different sets of results shown in Fig. 16a. To account for these fluctuations, swept simulations were conducted numerous times for a single set of geometric parameters and the results were averaged to generate a new stochastic performance metric 𝜙̇ 2 , as shown in Fig. 16b. A statistical metric known as the relative standard error (RSE) was used to determine when enough trials had been conducted to accurately capture the underlying phenomena across different sweeps as

𝜎

RSE = √x x n

(13)

where x is the mean, 𝜎x represents the standard deviation from the mean, and n is the number of trials being considered. The maximum RSE value is plotted as a function of n in Fig. 16c to show how the average solution converges with an increasing number of trials. The remainder of this paper averages results over 10 trials for each parameter combination because it was observed that the maximum RSE value changed very little beyond this point and was deemed to be sufficiently small (i.e. ≪ 1%). Fig. 17 shows the results for average performance as a function of both wind speed and system geometry. To generate these results, a base set of values is first assigned to the geometric parameters as shown in Table 1 in the Appendix. Then, a design range is given to three of these geometric parameters, r, m, and 𝓁 , that will be the range over which the following studies will

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

13

Fig. 14. Example illustrating how the sea surface state at a given point in space is converted into 𝜓 (t), 𝜓 (̇ t) , and 𝜓 ̈(t).

Fig. 15. (a) Plot illustrating how frequency spectrum varies with different wind speeds U19.5 . (b) Example showing how 𝜓 (t), 𝜓 (̇ t), and 𝜓 ̈(t) change as the wind speed is swept.

measure performance. Starting with r, a simulation sweeping wind speed is conducted using the first value in the design range for r and the base values for the other parameters while the corresponding performance is measured. This is carried out 10 times for each parameter value and the results are averaged to give a stochastic performance. This same process is then repeated across the entire range of r values to generate Fig. 17a and the same study is conducted for both m and 𝓁 to obtain the results shown in Fig. 17b and c, respectively. Finally, these results are mapped from their respective x-axes to the corresponding 𝜉 s values using

14

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

Fig. 16. Performance results during simulations sweeping wind speed where (a) 10 independent trials with random 𝛾 values are conducted, (b) these results are averaged, and (c) the relative standard error is tracked as a function of the number of trials considered.

Fig. 17. Illustration showing how the harvester’s performance changes as a function of wind speed and various geometric parameters. Results reveal desirable design regions and suggest that both 𝜉 s and individual geometric values are important in optimizing design.

Eq. (3) to allow for direct comparison and generate Fig. 17d–f. Hence, Fig. 17a, b, and c represent the same simulation results as Fig. 17d, e, and f, respectively, after being mapped nonlinearly to a different x-coordinate system. Collectively, these results make up a small portion of the complex, multi-dimensional design space which governs the harvester’s performance as a function of wind speed and geometric parameters. Because of the fact that so much of the actual design space is left unstudied, it is impossible to offer any definitive statements about optimal design choices for this system; however, four key trends can be confidently identified from these results that place useful limitations on the design space. These are: 1. The harvester does not perform well when U19.5 is too slow, as the waves provide insufficient excitation. This occurs at approximately U19.5 = 2 m/s. 2. The harvester does not perform well when 𝜉 s is too large because the system loses the opportunities for enhanced performance associated with potential well escapes. This occurs at approximately 𝜉 s = 0.4 rad.

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

15

3. The harvester does not perform well when 𝜉 s is too small. The reason for this drop in performance varies depending on the particular combination of parameters (e.g. when m is too small the rotating point mass is unable to provide the torque necessary to overcome damping), but appears to be consistent in all cases. This occurs at approximately 𝜉 s = 0.02 rad. 4. Beyond a certain threshold, increases in wind speed no longer provide additional increases in harvester performance. This is because the energy distribution described by the Pierson-Moskowitz spectrum shifts towards large amplitude waves at a relatively low frequency as U19.5 → ∞ and the rolling excitation is less sensitive to waves in this frequency range. This occurs at approximately U19.5 = 6 m/s. These findings agree with qualitative intuition about how the harvester should operate and offer useful quantitative limitations on the system’s complex design space. 5. Conclusion This paper performed an in depth analysis of a vibratory energy harvester that can be engineered for threshold escape and investigated its use in a simplified ocean setting to provide both phenomenological and design insights. In doing so, a fundamental explanation of the system’s static and dynamic behavior was given. This included a discussion of bifurcation diagrams, time dependent potential energy surfaces, and threshold breakthrough behavior. A criterion was developed to predict the onset of a potential threshold escape for parametrically excited systems with a time varying term in the potential energy function and a performance metric was described for comparing the desirability of different harvester responses. These tools were leveraged to conduct a design study for the deterministic single frequency excitation with both computational and experimental results showing that a harvester could be designed for superior performance by altering geometric parameters to suit the environmental excitation source. This analysis was then expanded to a multi-frequency stochastic case to emulate an environment in which the excitation is provided by ocean waves. Sea surface states were modeled using the Pierson-Moskowitz spectrum and simulation results provided limitations on the complex multi-dimensional design space. These insights are useful for numerous real world applications. First, they expand on previously existing methods for predicting potential well hopping to include systems with a time varying term in the potential energy function. As demonstrated throughout this paper, these types of systems offer various advantages over traditional systems in energy harvesting because of their ability to bifurcate as a function of excitation amplitude. This is particularly useful when attempting to harvest energy from environmental sources where the frequency of excitation is either unpredictable or wideband as is the case in ocean waves. Furthermore, the simple performance metric described here can be used to compare results across the wide breadth of harvesters already existing in the literature. This type of comparison is rarely included in studies on ocean energy harvesters but could represent a meaningful step towards advancement and collaboration in the field. Finally, the concrete limitations placed on the system’s parameter space enable future studies to conduct in depth analysis on design and performance through optimization. Without some practical constraints, optimization studies for a complex design space such as the one governing this system would be extremely computationally expensive. By restricting some of the key parameters, this computational burden can be reduced significantly. Future work could expand on the studies presented here by including heave and sway dynamics as well as the coupling between the buoy platform and the gimballed horizontal pendulum. These effects could have a significant impact on the harvester’s performance and accounting for them would help to inform design.

̈ (t ): Appendix Equations for converting sea surface conditions into 𝜓̇ (t ) and 𝜓 𝜓̇ (t) =

T B

(14)

where T =−

N ∑

ai

𝜔i 3 g

i=1

(

sin 𝜔i t + 𝛾i

)

(15)

and

B=1+

[N ∑ i=1

(

𝜓̈ (t) =

𝜕T 𝜕t

ai

𝜔i 2 g

)

( B−

(

cos 𝜔i t + 𝛾i 𝜕B 𝜕t

)

]2

.

(16)

) T

B2

(17)

where

∑ 𝜔4 ( ) 𝜕T =− ai i cos 𝜔i t + 𝛾i 𝜕t g i=1 N

(18)

16

D. Sequeira and B.P. Mann / Journal of Sound and Vibration 465 (2020) 115008

and

[

N ∑ ( ) 𝜔2 𝜕B = −2 ai i cos 𝜔i t + 𝛾i 𝜕t g i=1

][ N ∑ i=1

ai

𝜔i 3 g

(

sin 𝜔i t + 𝛾i

] )

.

(19)

Table 1 Base Geometric Parameters. Parameter

m (kg)

𝓁 (m)

r (m)

M (kg)

b (m)

Value Parameter Value

0.02 H (m) 0 .1

0.04 R (m) 0 .1

0.04 h (m) 0.005

0.02 c𝜃 (kg ·m2 /s) 5 × 10−5

0.03 c𝜙 (kg ·m2 /s) 5 × 10−6

References [1] M.F. Daqaq, R. Masana, A. Erturk, D.D. Quinn, On the role of nonlinearities in vibratory energy harvesting: a critical review and discussion, Appl. Mech. Rev. 66 (2014) 1–23. [2] S.C. Stanton, C.C. McGehee, B.P. Mann, Nonlinear dynamics for broadband energy harvesting: investigation of a bistable piezoelectric inertial generator, Physica D 239 (2010) 640–653. [3] C. Zheng, L. Shao, W. Shi, Q. Su, G. Lin, X. Li, X. Chen, An assessment of global ocean wave energy resources over the last 45 a, Acta Oceanol. Sin. 33 (2014) 92–101. [4] L. Dostal, K. Korner, E. Kreuzer, D. Yurchenko, Pendulum energy converter excited by random loads, J. Appl. Math. Mech. (2017) 349–366. [5] D. Yurchenko, P. Alevras, Parametric pendulum based wave energy converter, Mech. Syst. Signal Process. 99 (2018) 504–515. [6] S.M. Sah, C.C. McGehee, B.P. Mann, Dynamics of a horizontal pendulum driven by high-frequency rocking, J. Sound Vib. 332 (2013) 6505–6518. [7] W. Ding, B. Song, Z. Mao, K. Wang, Experimental investigations on a low frequency horizontal pendulum ocean kinetic energy harvester for underwater mooring platforms, J. Mar. Sci. Technol. 21 (2016) 359–367. [8] D. Sequeira, J. Little, B.P. Mann, Investigating threshold escape behavior for the gimballed horizontal pendulum system, J. Sound Vib. 450 (2019) 47–60, https://doi.org/10.1016/j.jsv.2019.03.008. [9] L.N. Virgin, R.H. Plaut, C. Cheng, Prediction of escape from a potential well under harmonic excitation, Int. J. Non-Linear Mech. 27 (3) (1992) 357–365. [10] B. Mann, Energy criterion for potential well escapes in a bistable magnetic pendulum, J. Sound Vib. 323 (3) (2009) 864–876. [11] D. Yurchenko, P. Alevras, Dynamics of the N-pendulum and its application to a wave energy converter concept, International Journal of Dynamics and Control 1 (2013) 290–299. [12] L.H. Holthuijsen, Waves in Oceanic and Coastal Waters, Cambridge University Press, 2007, https://doi.org/10.1017/CBO9780511618536. [13] C.C. McGehee, Dynamics of an Ocean Energy Harvester, Ph.D. thesis, Duke University, 2013.