INVESTIGATION OF STOCHASTICALLY EXCITED SUSPENSION SYSTEMS USING
AN INERTIALLY NONLINEAR VIBRATION ABSORBER WITH ENERGY
HARVESTING CAPABILITIES

By

Joel Cosner

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

Mechanical Engineering—Doctor of Philosophy

2025

ABSTRACT

The push towards electric vehicles in modern times is complemented by mobile and renew-

able power technology.

In automotive vehicles, ride comfort and road handling are also

important metrics to consider. Energy harvesting shock absorbers (EHSA) have potential,

but traditional applications have resulted in tradeoffs between ride comfort, road handling,

and power generation. Hence, research involving new technologies for power generation with

multi-objective performance capabilities is very valuable. This work specifically investigates

a novel approach for vibration suppression and electric power harvesting in systems sub-

jected to stochastic and broadband excitations, focusing on automotive suspensions. The

proposed solution harnesses the concept of inertance to provide mass amplification via ro-

tational inertia and also introduces inertial nonlinearity, enhancing bandwidth and overall

performance.

This work specifically explores vibration suppression, energy transfer, and the qualita-

tive change in the probability density function (PDF) with the application of a device with

nonlinearity due to inertia. This device is an inerter-based pendulum vibration absorber

(IPVA). The change in the PDF is called P-bifurcation. The IPVA is first applied to a single-

degree-of-freedom (SDOF) spring-mass-damper system subjected to white noise excitation.

A perturbation method identifies and tracks bifurcation points, revealing that the marginal

PDF of the pendulum’s angular displacement undergoes a P-bifurcation, transitioning from

monomodal to bimodal behavior. A cumulant-neglect technique predicts the system’s mean

squares, demonstrating significant vibration suppression near the P-bifurcation by transfer-

ring kinetic energy from the structure to the pendulum. Results are validated via Monte

Carlo simulations (MCS), approximating the PDF and mean squares.

The study extends to integrating energy harvesting into the nonlinear device, applying it

to a quarter-car model under class C road conditions (average road, ISO 8608). The impact

of pendulum length on power generation, ride comfort, and road handling is assessed. Near

P-bifurcation, simultaneous enhancements are seen in power output (40%), ride comfort

(60%), and road handling (60%) compared to a linear benchmark. A Wiener path integration

(WPI) method predicts the PDF and its second derivative, enabling efficient detection of

monomodal, bimodal, and rotational PDF regions in the noise intensity-electrical damping

plane. MCS confirm performance improvements of up to 43% in energy transfer and 20% in

power harvested compared to optimized linear systems, alongside at least 59% gains in ride

comfort and road handling. A novel bifurcation detection algorithm reduces computational

demands by linking qualitative PDF changes to performance metrics.

Experimental studies verify P-bifurcation and the energy harvesting IPVA’s effectiveness

in vibration reduction and energy harvesting for a SDOF structure under Gaussian broad-

band base excitation. Various experimental scenarios help identify unknown mathematical

model parameters, minimizing discrepancies between experimental and simulated results for

the pendulum’s velocity, with all RMS velocity differences below 3%. The fitted model

predicts power, vibration suppression, and P-bifurcation boundaries in the noise intensity-

electrical damping plane, corroborated experimentally. Power spectral density analyses re-

veal bimodal and rotational systems outperform monomodal configurations, enhancing power

and suppressing resonant peaks by up to factors of four and two, respectively. Near reso-

nance, mean square relative velocity improves by a factor of two.

These findings inform the design and manufacture of a full-scale IPVA-integrated EHSA

for off-road vehicles. The shock design involves balancing durability, weight, ride comfort,

energy harvesting, and road handling. The shock absorber undergoes sinusoidal testing on a

hydraulic machine to fit a developed mathematical model. A nonlinear least-squares based

optimization routine fits model parameters, yielding two viable parameter sets. Numerical

simulations implement the shock in a quarter car model with class F road excitation (off-road;

ISO 8608) and quantify energy harvesting, ride comfort, and road handling. A correlation

between performance and the PDF shape is finally made.

ACKNOWLEDGEMENTS

I would like to thank my advisor Dr. Wei-Che Tai, committee (Dr. Dan Segalman, Dr.

Zhaojian Li, Dr. Xiaobo Tan), lab peers (Aakash Gupta, Majid Baradaran Akbarzadeh,

Guoxin Li, Jamal Ardister), and family/friends for their continued support and help during

my PhD and for my professional development. Majid Baradaran Akbarzadeh also helped

with the transportation and experimentation associated with the full-scale energy harvesting

shock absorber. I would also like to acknowledge various undergraduate students who aided

me in conducting successful experiments: Sean Lishawa, Max Weidemann, and Olivia Reyes.

iv

TABLE OF CONTENTS

CHAPTER 1

INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . .

1

CHAPTER 2

ANALYSIS WITH THE NONLINEAR VIBRATION ABSORBER 13

CHAPTER 3

ENERGY HARVESTING, RIDE COMFORT AND ROAD HAN-
DLING IN AUTOMOBILE SUSPENSION SYSTEMS . . . . . .

CHAPTER 4

AN EXPERIMENTAL STUDY: P-BIFURCATION, ENERGY
HARVESTING, AND VIBRATION SUPPRESSION . . . . . . .

CHAPTER 5

FUTURE AND ONGOING WORK: AN EXPERIMENTAL STUDY
WITH THE FULL-SCALE ENERGY HARVESTING SHOCK
ABSORBER . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

35

62

89

BIBLIOGRAPHY . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

APPENDIX A: DRIFT AND DIFFUSION COMPONENTS OF MOMENT EQUA-

TIONS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

APPENDIX B: CLOSURE TECHNIQUES . . . . . . . . . . . . . . . . . . . . . . . 127

APPENDIX C: LINEAR SYSTEM RESPONSE . . . . . . . . . . . . . . . . . . . . . 129

APPENDIX D: LINEAR SYSTEM RHI ANALYTICAL SOLUTION . . . . . . . . . 130

APPENDIX E: DAMPING AND FRICTION FITTING IN THE TIME DOMAIN

FOR THE SDOF SYSTEM EXPERIMENT . . . . . . . . . . . . 132

v

CHAPTER 1

INTRODUCTION

1.1 Overview of the Work

This work seeks to provide a unique solution to how one might simultaneously achieve

vibration suppression and harvest electric power when considering an underlying excitation

that is naturally stochastic and broadband, i.e. randomly varying with a characteristically

wide power spectral density (PSD) as applied to suspension systems. The solution at hand

specifically involves a device that uses inertance to provide a mass amplification effect and

inertial nonlinearity, which broadens the effective bandwidth of operation.

The development of a device that leverages inerter characteristics—utilizing rotational

inertia instead of mass to minimize acceleration between two terminals [1]—fits well within

automotive suspension systems, where vehicle size constraints are often inflexible. As further

outlined in Section 1.2, the mass amplification accommodated by an inerter can provide a

dramatic enhancement of vibration suppression capabilities when applied to structures with

stochastic excitation, including vehicle suspension systems. This work specifically utilizes a

ball screw inerter type.

As will be further discussed in Section 1.2.3, one of these inerters equipped with energy

harvesting capabilities will generally fall short of providing optimal vibration suppression

while simultaneously generating optimal power output, if a linear design is used. In the case

of a vehicle suspension system, ride comfort and road handling may suffer simultaneously

with electric power generation. Hence, in this work, a nonlinear design is used. The design

chosen specifically accommodates nonlinear effects due to inertia rather than stiffness, uti-

lizing a horizontally attached pendulum vibration absorber to the ball screw inerter device.

The device is called the inerter pendulum vibration absorber (IPVA) and is coupled with

a DC generator when energy harvesting is implemented. It will be shown that this IPVA

accommodates random pendulum motion around more than one distinct position and with

large jumps in between each position, giving a probability density function (PDF) with two

1

peaks (bimodal). This behavior resembles that of popular bistable harvesters with random

broadband excitation found in literature such as in [2, 3, 4, 5] and [6], while those devices

are stiffness based and more bulky and/or too complicated for vehicle suspension systems.

Furthermore, the IPVA energy harvesting device is shown to accommodate superior vibra-

tion suppression and energy harvesting simultaneously, when compared to the comparable

linear system.

Starting in Chapter 2 and thereafter, it is found that performance is linked deeply with

the qualitative shape of the PDF. Performance is measured in terms of energy harvesting

simultaneously with vibration suppression as applied to a single degree of freedom (SDOF)

structure; as applied to a vehicle suspensions system using a quarter car model, performance

is in terms of energy harvesting, ride comfort, and road handling. Hence, unique methods

are developed to predict the phenomenological bifurcation involving a qualitative change in

the number and locations of the local maxima of the PDF, P-bifurcation [7, 8]. P-bifurcation

is further explained in Section 1.2.1.

The following details the organization of the culmination of this work hereinafter. This

research is comprehensively motivated in Section 1.2. Then in Chapter 2 the IPVA is studied

with applications involving a SDOF structure, comprised of a linear spring-mass-damper

system subjected to white noise base excitation to quantify vibration suppression. The

principle of stochastic bifurcation is explored with exciting connections between bifurcation

criteria and optimal kinetic energy transfer between the structure and the pendulum.

In Chapter 3 the device is integrated into an automotive suspension system whereby a

road profile defined by an ISO 8608 [9] standard power spectral density is used as the excita-

tion. Energy harvesting, ride comfort, and road handling are quantitatively and qualitatively

explored. The multi-objective performance is studied in conjunction with P-bifurcation of

the PDF in an attempt to build a meaningful correlation between performance and bifur-

cation. A unique adaptation of the Wiener path integration (WPI) formulation originally

outlined in [10] is augmented with a calculus of variations approach for the prediction of

2

P-bifurcation.

Finally, in Chapter 4, a culmination of the experimental work associated with the energy

harvesting device applied to a SDOF system subjected to random broadband base excita-

tion is presented and discussed in detail. The analytical treatment and simulation results of

Chapter 2 directly motivate an abundance of experimental qualitative validation and quan-

titative verification which are presented in detail in Chapter 4. The experimental studies as

well as the results of Chapter 3 further serve as a precursor to an experimental study of the

full-scale inertially nonlinear Energy Harvesting Shock Absorber (EHSA) which is addressed

in Chapter 5. The creative and unique experimental works coupled with the theoretical

analysis are believed to solidify a significant and scientifically impactful doctoral thesis.

1.2 Motivation and Background

1.2.1 A P-bifurcation primer with literature suggestions

It is important to give further information on P-bifurcation for a more thorough under-

standing of the underlying theory associated with part of this work. Literature associated

with stochastic dynamics such as [11], which provides an overview of the advancements in

stochastic dynamics by the close of the twentieth century, defines P-bifurcation as qualitative

changes in the stationary PDF and the related Markov process. For a more in-depth under-

standing of the mathematical methods used to examine P-bifurcations, one can refer to works

such as [12], which compiles the stochastic dynamics literature from the late 20th century. A

notable article in the text is one on P-bifurcation in the noisy Duffing-van der Pol equation

[13]. The authors of this work noted that the random motion was fast in nature and occurred

as random oscillations or rotations along the trajectories of the deterministic system. They

further used this fact to justify the use of stochastic averaging [14] and then proceeded to

derive the stationary PDF from the Fokker-Planck-Kolmogorov equation, where the FPE is

a partial differential equation governing the PDF [15]. This allowed them to mathematically

determine when the extrema in the PDF changed, indicating P-bifurcation. On the other

hand, the authors in [16] begin by reviewing deterministic center manifold theory, utilizing

3

unfolding parameters to examine behavior around bifurcations such as Hopf, transcritical,

and pitchfork types. They then derive the FPE and its solution near the center manifold

to identify P-bifurcations. Analyses cover both additive noise, where noise is added to the

state space; and multiplicative noise, where noise multiplies the states in the system. In ad-

dition to the aforementioned analysis methods, the authors of [17] mentioned the notion of

bifurcations in the moment equations governing the statistical moments associated with the

system [14, 15]. One analytical method employed in the current work involves detecting the

change in stability of moment equation fixed points by determining when the determinant of

the Jacobian vanishes. This method is detailed in Chapter 2 and used again in Chapter 4.

Another analytical formulation is developed in this work which utilizes the WPI technique

[18, 19, 20] and exploits knowledge of the curvature of the PDF at the bifurcation. The WPI

method is detailed and used in Chapter 3.

1.2.2 Vibration mitigation application

The inerter is a mechanical device with two terminals, each of which exerts an equal

and opposite inertial force proportional to the relative acceleration between the terminals

[1]. The constant of proportionality in this case is known as the inertance. Furthermore,

the transmission mechanism used by inerters commonly includes the ball screw [21] or rack-

pinion [22]. These devices provide a means for conversion of linear translational motion into

rotational motion. In this way, a light-weight object that rotates with the device may add

an effective translational inertia manifested by the moment of inertia of the object instead

of its mass. As a result, a well-chosen transmission ratio can provide an effect of mass

amplification, multiplying the mass of the lightweight object by a factor of 10 or more.

Therefore, the inerter is a favorable device for vibration suppression in many engineering

structures and components.

Researchers have been attracted to inerters in the context of vibration mitigation in-

volving broadband random excitation. For example, Marian and Giaralis [23] integrated

the inerter with the classical tuned mass damper (TMD), known as the tuned mass damper

4

inerter (TMDI), to suppress the response of stochastically motion-excited structures. They

concluded that the “mass amplification effect” of the inerter can enhance the performance

of the classical TMD. Hu and Chen [24] proposed four types of inerter-based dynamic vibra-

tion absorber (IDVA) to suppress the response of linear spring-mass systems force-excited

by white noise and evaluated their H∞ and H2 performance compared with the classical

TMD. Specifically, the proposed IDVAs incorporated a spring connected with the inerter in

series or parallel, such that the inerter provided another degree of freedom to the DVA. They

concluded that the additional degree of freedom is necessary to enhance the H∞ and H2 per-

formance and enlarge the effective frequency bandwidth of response suppression. Moreover,

Joubaneh and Barry [25] studied the performance of four models of electromagnetic reso-

nant shunt TMDI (ERS-TMDI) on both vibration suppression and energy harvesting and

identified the best model. Their parametric studies showed that increasing the inertance en-

hances the performance of the best model in terms of both vibration mitigation and energy

harvesting. On the other hand, Tai [26] integrated the inerter with a torsional mass damper,

known as the tuned inerter-torsional-mass damper (TITMD), to suppress response of linear

spring-mass-damper systems motion-excited by white noise. In comparison with the TMDI,

the TITMD achieved 20 − 70% improvement in the H2 performance when their weights were

identical.

In this work, the IPVA is in fact a nonlinear dynamic absorber. Nonlinear dynamic

vibration absorbers by virtue of inertial nonlinearity have been studied by many researchers

to suppress vibration of structures. One specific class of nonlinear dynamic vibration ab-

sorbers, known as autoparametric vibration absorbers [27], utilize specific inertial nonlinear-

ity to transfer the kinetic energy of a primary structure requiring reduced vibration to the

vibration absorber when subject to harmonic excitation [28, 29, 30] and random excitation

[31, 32]. By virtue of the energy transfer phenomenon, the autoparametric vibration ab-

sorbers have been shown to achieve vibration mitigation and energy harvesting at the same

time [33, 34, 35, 36, 37].

5

Figure 1.1 Electric vehicle quarterly sales for the years 2021-2024 in millions. The y-axis is
the quantity in millions. Orange is China, blue is Europe, green is the United States, and
grey is the rest of the world. The chart was taken from the International Energy Agency
(IEA) [38].

1.2.3 Energy harvesting considerations

The International Energy Agency (IEA) reported quarterly sales for 2021-2024 and the

result is shown in Fig. 1.1. In subsequent years a clear monotonically increasing trend is

seen. The IEA specifically note a 25% increase from Q1 in 2023 to Q1 in 2024. This growth

is substantiated by the fact that 2023 had a record-breaking 14 million electric vehicles sold,

or 18% of all vehicles sold [39]. Furthermore, the International Council on Clean Trans-

portation (ICCT) has reported governments that intend to phase out vehicles that are not

zero-emissions by a designated year [40] and this is shown in the map of Fig. 1.2. This

6

Figure 1.2 Governments with targets to only sell zero-emission vehicle by a certain date [38].

phasing out process can be coined ”electrification” of the automotive industry. This ”elec-

trification” serves as motivation for research in efficient power generation in electric vehicles

in order to accommodate less stringent battery charging and energy storage requirements.

Other needs for efficient energy harvesting in vehicles can be realized in cases when the mil-

itary wishes to remain untethered to power grids and energy independent in expeditionary

environments. Energy harvesting shock absorbers (EHSAs), researched for decades for the

purpose of converting vibration energy into electrical power [41, 42], could be of significant

use to electrical vehicle manufacturers as well as the military as they can effectively act as

a mobilized power resource. A recent review of research in the field suggests a significant

power generation potential for varying degrees of road classes from ISO 8608.[43].

Recent research in the field of small-scale energy harvesting has given promise to the use

of bi-stable systems to enhance energy harvesting; see the extensive reviews of [44, 45]. Other

recent examples exist such as in [4], where the authors designed an electromagnetic energy

harvester device with stiffness nonlinearity due to the geometry of the configuration which

was shown to outperform the corresponding linear device with filtered and band-limited

7

stochastic excitation. Overall, while much research progress has been made regarding bi-

stable harvesters, stiffness nonlinearities are unfortunately not easy to implement or practical

in the case of vehicle suspension systems.

To understand the importance of nonlinearity in EHSAs, it is crucial to look at the

abundance of research in traditional EHSA designs [46, 47, 48]. Traditional linear designs

[49, 50, 51, 52] cannot achieve optimal vibration suppression and power harvesting simulta-

neously. In fact, research shows that linear designs trade suspension performance for energy

recovery [50, 51, 52]. Additionally, peak energy harvesting performance lies in the narrow-

band linear resonance regime [53] whereas road irregularities are stochastic and broadband

in nature [54]; hence, road tests have frequently reported unsatisfactory and unsteady energy

recovery efficiency of linear EHSA designs [43]. Performance degradation in linear suspen-

sion systems when purposed to harvest energy has motivated research efforts in regard to

bistable suspension systems.

One significant feature of bistable systems which contributes to improved energy har-

vesting has been identified as the frequent transitions between potential wells, also known as

interwell oscillation [2, 3, 5]. When interwell oscillations are permitted, the multi-stability in

energy harvesters in the case of random excitation manifests itself as a multimodal (multiple

peaks) PDF associated with the states of the system [4, 55, 56]. As such, analysis often

involves obtaining the PDF by solving the FPE with the finite element method [56, 57] or

by obtaining an approximation with the WPI method [10, 18, 19, 58]. Otherwise, statistical

moment equations governing the system’s statistics can also be derived [59, 60] as another

means of analysis. A Monte Carlo simulation (MCS) is then often used for verification pur-

poses [58, 60, 61, 62], which allows for the determination of the necessary statistics with

direct numerical simulation.

The device introduced in Chapter 2 is a specific device which in lieu of stiffness nonlin-

earities harnesses inertial nonlinearity and accommodates bimodality. This device is studied

with stochastic excitation and it will be found that when the PDF of the system undergoes a

8

bifurcation from monomodal to bimodal, meaning that two equally probable equilibria exist,

an energy transfer occurs between the structure and the pendulum. To this end, this energy

transfer is responsible for mitigating vibrations of the structure, while the energy in the

pendulum is also increased. While energy harvesting is not incorporated yet in Chapter 2,

the energy transfer found serves to motivate transferring the energy to the electrical domain

in an energy harvesting context.

In [62], the IPVA device was further explored when the device was implemented in a

quarter-car model suspension system whereby it was used in conjunction with a planetary

gear system to successfully drive a generator while simultaneously optimizing ride comfort.

The authors in [62] determined that the nonlinear device outperformed the linear EHSA with

an increase in power of 45% and a 45% reduction in the root-mean-square (RMS) acceleration

of the sprung mass. An exhaustive Monte Carlo simulation (MCS) was used to discover and

illustrate these findings. It is the intent in Chapter 3 to return attention to the quarter car

model with increased suspension stiffness to accommodate off-road vehicle simulation, but

this time investigate the system using knowledge of stochastic P-bifurcation [63]. While this

is explored in Chapter 3, it is not a trivial task. Due to the cost of computing the PDF, P-

bifurcation analysis has focused on systems with five or fewer state variables, [64, 65, 66] for

example. Additional research investigating stochastic bifurcation and analytical formulations

for the PDF has been done such as in [55, 67]. However, these approaches unfortunately rely

on the stochastic averaging approach. The method of stochastic averaging can only produce

accurate results in the case of weak excitation and weak nonlinearity, allowing one to assume

a slow varying response [15]. As the suspension system has more states and weak excitation

and nonlinearity cannot be assumed, an effective analysis approach remains open. To that

end, an analytical formulation is developed in this work which utilizes the WPI technique

[18, 19, 20] and exploits knowledge of the curvature of the PDF at the bifurcation.

The research surrounding the multi-objective performance capabilities of the nonlinear

IPVA when subjected to random excitation such as in [58, 62] is very relevant to modern

9

research. As pointed out in a recent review of dual-objective harvesters [68], the research

surrounding the use of nonlinearity to achieve multi-objective performance is still scarce. It is

even more scarce when considering excitation which is stochastic in nature. The research that

does exist is often associated with the longtime running research topic of nonlinear energy

sinks (NES) [69] which inherently require cubic stiffness nonlinearity via the use of elastic

components or magnetic arrangements. For example, the researchers in [70] theoretically

coupled an NES with a giant magnetostrictive energy harvester (GMEH) and applied the

device to a single degree of freedom structure subjected to Gaussian white noise excitation to

achieve energy harvesting and vibration suppression objectives simultaneously. Regardless,

the current work avoids stiffness nonlinearities which could make the system too bulky and

impractical if applied to an automobile suspension system.

1.2.4 Experimental Investigations of Bimodal Energy Harvesting Systems with

Broadband Random Excitation

The experimental exploration of bi-stable energy harvesters subjected to random exci-

tation has accompanied much of the research in this particular field. For example, in [2]

the authors qualitatively verified simulation results for a Duffing-type harvester subjected to

band-limited noise excitation, investigating how the center frequency and the bandwidth of

the excitation affected the output voltage in the mono-stable and bi-stable configurations.

In [3], the researchers experimentally validated simulation results for a bi-stable energy

harvester under band-limited white noise excitation which was transformed into a flexible

bi-stable energy harvester exhibiting a variable potential energy function due to a magnet

on a flexible beam. Similarly, in [4] the authors experimentally validated simulation results

for a bi-stable two-degree-of-freedom harvester system and investigated the effect of varying

excitation intensity, bandwidth, and center frequency. Likewise, in [5] the authors explored

and experimentally validated the theoretical results for a piezoelectric bi-stable harvester

subjected to wide-band random excitation, while varying excitation intensity and the sep-

aration between equilibrium positions. In short, [2, 3, 4, 5] and others like [6] all share a

10

common conclusion: bi-stable energy harvesting systems will outperform mono-stable coun-

terparts under broadband random excitation precisely when the excitation is large enough to

produce interwell oscillations and they all experimentally validate or verify their theoretical

results. Nevertheless, prediction of the parameter values required for interwell oscillations

has remained open thus far.

To understand when multiple peaks in the PDF can occur, or equivalently when interwell

motion capabilities change for multi-stable systems, one can understand P-bifurcations which

was discussed in the end of Sec. 1.2.2. Some approaches in the past to study bifurcations

have included the numerical computation of the PDF with the finite element method [57]

or deriving estimates for the PDF [57, 64, 65, 66], which have been computationally limited

to systems with less than five states. Then [71] employed a technique based on Shannon

Entropy, whereas this study introduces a WPI-based PDF approximation algorithm tailored

for a system with seven states, as detailed in Chapter 3. With that being said, to the

author’s knowledge, no openly available research attempts to predict P-bifurcation pertaining

to experiments with wide-band random excitation. Lastly, while many researchers have

remedied interwell oscillation difficulties by studying tri-stable [72, 73] or even quad-stable

[74] energy harvesters to boost energy harvesting efficiency, it is hypothesized from evidence

in [62] and is shown in Chapters 2 and 3 of this work that by using inertial nonlinearity

rather than stiffness nonlinearity, a bimodal PDF can still be produced without a multistable

potential and is often associated with large amplitudes of oscillation just as in multi-stable

systems with interwell motion. Chapter 4 exploits such a system and seek to predict changes

in the number of peaks in the PDF in experiments, indicating when the large amplitude

oscillations may occur.

Chapter 4 focuses on the bimodal IPVA energy harvesting device applied to an SDOF

structure subjected to Gaussian broadband excitation. The objective is to experimentally in-

vestigate P-bifurcation of the PDF and the correlation with simultaneous energy harvesting

and vibration suppression. Notably, this study will include the experimental parameter char-

11

acterization for a developed prototype to accommodate accurate predictions of energy har-

vesting potential, vibration suppression, and bifurcation boundaries. Adequate experimental

data will be used to verify predictive capabilities, confirm the bifurcation phenomenon, and

gather information regarding power harvested and vibration suppression.

12

CHAPTER 2

ANALYSIS WITH THE NONLINEAR VIBRATION ABSORBER

2.1 Introduction

In this chapter, we conduct a P-bifurcation analysis of a linear spring-mass-damper sys-

tem (primary structure) incorporating the IPVA, which has white noise base excitation,

made to an equivalent force excited system. The qualitative change in the number and loca-

tions of the maxima of the probability density function (PDF) is thereby examined as well

as the correlation between P-bifurcation and kinetic energy transfer of the linear system to

the pendulum.

The chapter is organized as follows. In Sec. 2.1.1, the stochastic differential equations

of the structure-IPVA system are derived. In Sec. 2.2, stochastic bifurcation is investigated

while implementing a bifurcation detection and tracking procedure. Numerical validation

including realizations in the time domain is imposed with the aid of a Monte Carlo simulation,

while an arc-length continuation scheme is developed to track bifurcation points as a function

of parameters. Energy transfer from the primary structure to the pendulum is investigated in

the neighborhood of bifurcation with the use of a cumulant-neglect procedure. Furthermore,

the relationship between the energy transfer and vibration mitigation is investigated.

In

Sec. 2.3, a parametric study is conducted to study the system at bifurcation while varying

pendulum length, mass ratio, noise intensity, and pendulum damping ratio. A discussion

is then given in Sec. 2.4. Finally, this chapter will serve as a prelude to applications in

automobile suspension systems.

2.1.1 The Inerter Pendulum Vibration Absorber

The design of the IPVA is shown in Fig. 2.1. The IPVA is attached between the ground

and a linear oscillator (primary structure) of mass M , stiffness k, and damping c. Linear

motion x of the mass is converted to angular motion θ of the ball screw with effective radius

R = L/2π, where L refers to the lead value associated with the ball screw as a measure of

13

Figure 2.1 Schematic of IPVA attached to a structure of mass (M) with defined degrees of
freedom θ, and ψ. The second pendulum on the left is assumed to be fixed.

the ratio of linear displacement of the nut to full rotation of the screw. More specifically,

x = Rθ. A pendulum, of length r and mass m, is also attached perpendicular to the screw,

at a radius of Rp from the center of the screw and with an angular displacement of ψ with

respect to the attachment point. Note that two pendulums are shown, but one is to be fixed

in place to avoid rotating unbalance. The free pendulum will oscillate and absorb energy

from the primary structure, with varying magnitudes of energy transfer dependent on the

parameter values of the system. F corresponds to the forcing applied to the structure, and

in the present work it is defined to be a Gaussian white noise with zero mean and standard

√

deviation Γ =

2S, where S is the noise intensity.

In the sections that follow, we will

attempt to quantitatively and qualitatively describe the effect of parameter variation of this

IPVA on the bifurcation of the probability density function as well as the relationship with

vibration suppression.

14

𝜃𝜃=𝑥𝑥−𝑦𝑦𝑅𝑅𝜓𝜓𝑚𝑚𝑥𝑥Attachment to Ground𝑅𝑅𝑝𝑝𝑟𝑟To derive the equations of motion using Lagrange’s equations, we first derive the kinetic

energy, potential energy, and virtual work associated with the system. The kinetic energy

for the system is given by

T =

h

1
2

M ( ˙x + ˙y)2 + J ˙θ2 +

(cid:0)

Jp + mr2

(cid:1)(cid:16)

(cid:17)

(cid:16)

(cid:16)

(cid:17)(cid:17)i

˙ψ + ˙θ

2 + m

R2
p

˙θ2 + 2Rpr cos(ψ) ˙θ

˙θ + ˙ψ

(2.1)
where ˙(·) = d (•) /dt represents the first time derivatives, and J and Jp correspond to the

moment of inertia about the axis of rotation of the screw and the pendulum, respectively.

The potential energy is given by

V =

1
2

k (Rθ)2

(2.2)

Furthermore, a torsional viscous damping coefficient cp is introduced to account for energy

loss at the pivot point of the pendulum. The virtual displacement of the mass is given by

δx = Rδθ and the virtual displacement of the pendulum in the direction of the pendulum

damping torque is δψ.

It follows that the summed virtual work due to the mechanical

damping of the primary structure motion as well as pendulum motion is written as

δQ = −cR2 ˙θδθ − cp

˙ψδψ

(2.3)

With (2.1), (2.2) and (2.3), the equations of motion for the system were derived and are

given here.

(cid:0)

M R2 + J + mR2

p + Jp + mr2 + 2mRpr cos (ψ)

(cid:0)

(cid:1)

¨θ +

mr2 + mRpr cos (ψ) + Jp

(cid:1)

˙ψ

+cR2 ˙θ + kR2θ − 2mRpr ˙ψ ˙θ sin (ψ) − mRpr sin (ψ) ˙ψ2 = F (t)R,
(cid:0)

(cid:1)(cid:1)

(cid:1)

(cid:0)

(cid:0)

mr2 + Jp

¨ψ +

Jp + m

r2 + Rpr cos (ψ)

¨θ + cp

˙ψ + mRpr sin (ψ) ˙θ2 = 0.

(2.4)

where the base excitation is converted to direct forcing mathematically via F (t) = −M ¨y.

In (4.6), ¨(•) = d2 (•) /dt2 represent the second time derivatives. Let us now define ˆM as
R2 )R2 = ˆM R2. Thus,
the moment of inertia term J can be absorbed into the effective mass term with the inertial

the effective structural mass. It then follows that M R2 + J = (M + J

addition given by the inertance J/R2. Evidently, for small R, mass amplification is large.

15

Next note that the pendulum shaft moment of inertia Jp is also negligible relative to the

inertia ˆM R2 and so it can be set to zero without any noticeable change in the resulting

dynamics. Next, ψ = ϕ + ϕ0, where ϕ0 corresponds to an offset pendulum angle and ϕ is the

angular displacement with respect to ϕ0.

To facilitate a parametric study, the following dimensionless parameters were appropri-

ately chosen.

r

µr =

mR2
p
ˆM R2

, ω0 =

τ = ω0t,

˙(·) = ω0 (·)

, ξ =

, η =

r
k
ˆM
Rp
′
, ¨(·) = ω2
0 (•)

′′

c
2 ˆM ω0

, ξp =

cp
2 ˆM R2ω0

(2.5)

The equations of motion for the system are then rewritten in terms of the dimensionless

parameters as follows:

M (ψ) x′′ + Cx′ + Kx + g (x, x′) = f

(2.6)

where

2

6
4

2

6
4

M (ψ) =

C =

1 + µr(1 + η2 + 2η cos (ψ))

µrη(η + cos (ψ))

3

7
5 ,

µrη(η + cos (ψ))

3

2

2ξ

0

7
5 , K =

6
4

0

2ξp

0

1 0

0 0

3

0

µrη2
1

7
5 , f =

B
@

W (τ )

C
A ,

g (x, x′) =

B
@

µrη

−(2ψ′θ′ + ψ′2) sin (ψ)

θ′2 sin (ψ)

0

1

C
A

(2.7)

In (3.2) and (2.7), x = (θ, ψ)T , M is the inertia matrix which is inherently nonlinear due to

the Euler acceleration as a result of a rotating reference frame, C is the damping matrix,

K is the stiffness matrix, and the vector f is associated with stochastic forcing. The force

term W (τ ) = F (τ )/ ˆM Rω2

0 is then defined to be the normalized Gaussian white noise with

zero mean and noise intensity D, while ξ and ξp are the damping ratios associated with the

structure and the pendulum respectively. Furthermore, g represents the nonlinear Coriolis

and centrifugal terms.

16

Note that the equations of motion are written in terms of the normalized time variable,

τ . Any temporal results obtained will thus be on a time scale of τ . The current study is

primarily associated with the response characteristics of the system as a function of η and

µr. The value of η is a measure of the pendulum length, while µr is proportional to the

pendulum-structure mass ratio. Both parameters quantify the nonlinear coupling associated

with the system as can be seen upon inspection of g.

To appropriately analyze the system with stochastic forcing, the equations of motion

(3.2) are transformed into the Itô stochastic differential form [31, 75] as follows.

dy = a (y) dτ + Σ (y, τ ) dB

(2.8)

where

1

C
C
C
C
C
C
C
A

0

B
B
B
B
B
B
B
@

dy1

dy2

dy3

dy4

=

dy =

1

C
C
C
C
C
C
C
A

0

B
B
B
B
B
B
B
@

dθ

dψ

d ˙θ

d ˙ψ

and B is a scalar zero-mean and delta-correlated Wiener vector process, with ⟨dB(τ )⟩ = 0

(cid:10)

(cid:11)

and

dB(τ )dBT(τ )

= 2Ddτ , where ⟨·⟩ represents the expectation operator and D is the

noise intensity. Additionally, a and Σ are a 4×1 drift vector and diffusion vector, respectively,

namely

2

6
4

a (y) = −

2

6
4

Σ (y, τ ) =

3

0

7
5 y −

B
@

1

C
A

02×1

M−1g (x, x′)

02×2

I

M−1K M−1C
√

1

0

3

02×2

M−1

7
5

B
@

2D

C
A

0

(2.9)

The explicit forms of a and Σ are given in Appendix A.

2.2 Stochastic Bifurcation

It is known that the pendulum may oscillate with respect to a nonzero equilibrium point

when subject to harmonic excitation [76]. If such asymmetric oscillations exist when subject

17

to random excitation, these asymmetric oscillations correspond to a nonzero mean response

and indicate a change in the locations (and likely the number) of the maxima of the PDF.

It is shown in Sec. 2.2.3 that such P-bifurcation indeed exists in the system and leads to a

transfer of the kinetic energy of the primary structure to the pendulum.

2.2.1 Bifurcation Detection and Tracking

P-bifurcations of highly nonlinear systems are commonly predicted in the literature by

computing the PDF with numerical techniques, such as the cell mapping method [7], the finite

difference method [8], and the path integration method [77]. The computational cost of these

methods are dependent on the number of the state variables, and can be too expensive when

applied to the current system of four state variables. A perturbation method is employed

instead in this work to detect and track the P-bifurcation of the system, which is based

on the following principle in random vibration. It is known that when a linear oscillator is

subject to zero-mean Gaussian white noises, the response is zero-mean Gaussian [78]. This

implies that when subject to zero-mean Gaussian white noises, the response of a linearized

nonlinear oscillator is zero-mean Gaussian with respect to a local maximum of the PDF

because that is where the response spends relatively much time [79].

In other words, if

the equations of motion (2.8) is linearized with respect to a local maximum, the response

of the linearized system should be zero-mean. Linearization of (2.8) will nevertheless lead

to a degenerate stiffness matrix because the pendulum lacks any stiffness terms, which will

incur computational difficulties. A perturbation method is therefore applied to (2.8), instead

of linearization, which is explained as follows. Because only asymmetric oscillation of the

pendulum is of interest, a small perturbation ϕ = ψ − ϕ0 is introduced to (2.8), where

|ϕ| << 1 and ϕ0 represents the location of an unknown local maximum to be determined.

After expanding with respect to ϕ0 and retaining only linear terms of ϕ, one arrives at

dˆy = ˆa (ˆy, ϕ0) dτ + ˆΣ (ˆy, ϕ0, τ ) dB

(2.10)

18

where ˆy = (θ, ϕ, θ′, ϕ′)T and

ˆa (ˆy, ϕ0) =

(cid:18)

a (y) +

(cid:20)

(cid:20)

ˆΣ (ˆy, ϕ0, τ ) =

Σ (y, τ ) +

(cid:19)

(cid:21)

ϕ

∂a (y)
∂ψ
(cid:18)

∂Σ (y, τ )
∂ψ

ψ=ϕ0
(cid:19)

(cid:21)

ϕ

ψ=ϕ0

(2.11)

Because (2.10) governs small oscillation in the neighborhood of a local maximum, the

mean value ⟨ˆy⟩ should be approximately zero. In other words, one can use (2.10) to find a

nontrivial ϕ0 such that ⟨ˆy⟩ = 0. To this end, the stationary moment equations of (2.10) are

derived. Specifically, a moment generating function mpqrs = ⟨yp

1yq

2yr

3ys
4

⟩ up to order n = 2 is

defined, where p, q, r, s ∈ W with p + q + r + s ≤ n. With the use of mpqrs, one can apply the

Itô formula to (2.10) to generate stationary moment equations for the system [80], namely
(cid:28) (cid:16)

(cid:29)

(cid:28)

(cid:29)

(cid:17)

4X

i=1

∂yp

3ys
4

1yq
2yr
∂yi

(ˆa)i

+

1
2

4X

i,j=1

ˆΣ ˆΣT

ij

3ys
4

∂2yp

1yq
2yr
∂yi∂yj

= 0

(2.12)

where (·)i and (·)ij represent the i-th and ij-th component of a vector and square matrix,
respectively, and ˆa and ˆΣ are defined in (2.11). Because only small oscillation is of interest,

the Gaussian closure technique (outlined in Appendix B) is applied to (2.12), resulting in 14

independent nonlinear algebraic equations. The resultant moment equations are very long

and not given explicitly in this paper; the general form of the equations are

ˆh (m1000, m0100, · · · , m0020, m0002, ϕ0) = 0

(2.13)

where ˆh represent 14 nonlinear functions of the four first order moments and ten second order

moments and the unknown ϕ0. Note that (2.13) are underdetermined (14 equations and 15

variables). Assuming that the mean of the small pendulum angular displacement to be zero,

i.e., m0100 = ⟨y2⟩ = ⟨ϕ⟩ = 0, ϕ0 can be determined. The resultant ϕ0 will approximate the

locations of maxima of the PDF.

The stationary moment equations given by (2.13), provide a means to determine the value

of ϕ0 as a function of system parameters. To accomplish this, an arc-length continuation

technique is employed. Arc-length continuation provides a means for tracking the value of

19

ϕ0 while varying a parameter of interest. Denote by z the vector consisting of ϕ0 and the
13 moments excluding m0100, and by h = ˆh(m0100 = 0) the vector consisting of all the 14

nonlinear equations in (2.13) with m0100 = 0. Applying the arc-length continuation from

[81] to (2.13), one arrives at the following:

8

>>>><
>>>>:

h(z; λ) = 0

Jz(z; λ) dz

∂λ

ds + ∂h(z;λ)
ds )2 = 1

dz
ds

· dz
ds + ( dλ

dλ
ds = 0

(2.14)

where s represents the arc-length parameter, Jz is the Jacobian matrix of h, while λ is a

parameter of interest. The arc-length parameter is finally defined in the last of (2.14). The

culmination of (2.14) define fifteen first order, nonlinear differential-algebraic equations that

can be numerically integrated to obtain the maximum location ϕ0(λ) as a function of λ,

given with initial conditions that satisfy (2.13).

According to Sharif-Bakhtiar and Shaw [76], asymmetric oscillation of the pendulum is

induced by a pitchfork bifurcation which belongs to limit point bifurcation. To determine

the limit point, the first equation in (2.14) are differentiated with respect to the bifurcation

parameter λ, resulting in

∂h (z; λ)
∂z

dz
dλ

= − ∂h (z; λ)

∂λ

(2.15)

Limit point bifurcation occurs when any component of dz/dλ → ∞ or when the determinant

|∂h (z; λ) /∂z| = 0. The limit point and critical value of the bifurcation parameter can be

determined by solving the first equation in (2.14) together, i.e.,

8
><

>:

h (z; λ) = 0
(cid:12)
(cid:12)
(cid:12) ∂h(z;λ)
∂z

(cid:12)
(cid:12)
(cid:12) = 0

(2.16)

To track the limit point in the parameter space, (2.16) is augmented with the arc-length

20

continuation scheme, namely

8

>>>>>>><
>>>>>>>:

h(z; λ, γ) = 0
(cid:12)
(cid:12)
(cid:12) ∂h(z;λ,γ)
∂z

(cid:12)
(cid:12)
(cid:12) = 0

Jz(z; λ, γ) dz

dz
ds

· dz
ds + ( dλ

∂λ

ds + ∂h(z;λ,γ)
ds )2 + ( dγ

ds )2 = 1

dλ

ds + ∂h(z;λ,γ)

∂γ

dγ
ds = 0

(2.17)

where γ is a second bifurcation parameter of interest. Since the limit point occurs with

ϕ0 = 0, ϕ0 = 0 is substituted in (2.17), which annihilates three first order moments m1000 =

m0010 = m0001 = 0 and two second order moments m1010 = m0101 = 0. Solving (2.17) yields

the limit point curve ϕ0 (λ, γ) = 0 in the λ − γ plane, which describes the values of λ and γ

for which the limit point bifurcation occurs.

2.2.2 Numerical Validation

Using (2.14), ϕ0(µr) and ϕ0(η) were determined as a function of µr and η, respectively, for

specific parameter values. The results are shown in Fig. 2.2. There are several things worth

noting in Fig. 2.2. First, as seen in Fig. 2.2(a), ϕ0(η) bifurcates into two solution branches

around η = 0.65. For η values smaller than this critical value, only one solution branch ϕ0 = 0

exists, which indicates that the pendulum symmetrically oscillates with respect to ϕ0 = 0.

This is expected because when η or µr are small, the nonlinear terms g in (4.10) are weak;

thus, the system behaves similar to a linear system. As η increases and exceeds the critical

value, the nonlinearity becomes sufficiently strong such that asymmetric oscillation occurs.

In the context of random vibration, coexistence of two local maxima means a bimodal PDF;

that is, the pendulum can symmetrically oscillate with respect to either maximum location

and randomly “jump” between each other. The same bifurcation is observed in ϕ0(µr); see

Fig. 2.2(b).

To validate the prediction in Fig. 2.2, a Monte Carlo simulation (MCS) was used to

numerically compute the marginal PDF of ϕ for various η values. Numerically integration

of (2.8) was specifically done by ItoProcess of Mathematica 12.1, which is a stochastic dif-

ferential equation solver. The time step was chosen to be 10−3 units to ensure accuracy.

21

(a) ϕ0(η) with µr = .3, ξ = .02, ξp = .07 and D = .02.

(b) ϕ0(µr) with η = .65, ξ = .02, ξp = .07 and D = .02 .

Figure 2.2 Bifurcation diagram of ϕ0 against (a) η and (b) µr.

22

An integration time of 500 units was found to be sufficiently long so as to ensure stationar-

ity as it allowed for asymmetric pendulum oscillation to occur. Furthermore, a sufficiently

large number of realizations (around 4 × 105) were recorded to ensure accuracy. SmoothK-

ernelDistribution of Mathematica 12.1 was then used to compute the marginal PDF of the

realizations. The result for various values of η and specified parameter values corresponding

to Fig. 2.2(a) is plotted in Fig. 2.3. There are two things worth noting in Fig. 2.3. First, it is

evident that around η = .65 a bifurcation occurs as the distributions become bimodal after-

wards, which matches the prediction in Fig. 2.2(a). Second, it is apparent that an increase

in the pendulum length parameter η leads to an increase in the distance between the max-

imum locations. The maximum locations predicted by Monte Carlo simulation are plotted

as dots for comparison in Fig. 2.2(a). In comparison with the Monte Carlo simulation, the

prediction by (2.14) matches well until η = 0.8. Afterwards the nonlinearity becomes too

strong such that the perturbation method is no longer accurate. The implications of this

bifurcation are further addressed in the sections that follow.

In order to provide additional insight and verification for the determined marginal PDFs,

realizations in the time domain for the pendulum position and velocity are also shown in

Fig. 2.4. The realizations were produced with η values before, near and after bifurcation.

The duration of time was chosen to be from 500 to 1000 units as it is assumed that the system

is stationary after 500 units of time. The first remark is that the mono-modal oscillation

around ϕ = 0 is apparent before bifurcation. Near bifurcation we also start to see a drift

in the oscillation away from ϕ = 0 accompanied by larger amplitude spikes in the angular

displacement and velocity. Lastly, after bifurcation we see the pendulum tends to oscillate

about one equilibrium angle for a more significant amount of time before making a jump to

another equilibrium angle.

2.2.3 Stochastic Bifurcation and Energy Transfer

In this section, the mean square velocities of the primary structure and pendulum oscilla-

tion are investigated. Specific attention is paid to investigating how the bifurcation influences

23

Figure 2.3 Monte Carlo simulation of the marginal PDF for ϕ evaluated with values of η
around bifurcation and with µr = .3, ξ = .02, ξp = .07, and D = .02.

Figure 2.4 Realizations in the time domain for pendulum angular displacement and velocity
before bifurcation, near bifurcation, and after bifurcation. µr = .3, ξ = .02, ξp = .07, and
D = .02.

the mean square velocity. Eqn. (2.13) allows for the determination of the bifurcation point

of the moment equations and is built upon the assumption of a Gaussian distribution in the

neighborhood of ϕ0. A higher order of nonlinearity is required to accurately predict all mean

(cid:16)

(cid:17)
ϵ¯θ, ϵ ¯ϕ, ϵ ˙¯θ, ϵ ˙¯ϕ

T

= ϵ¯y, where ϵ is a small

square values. To this end, it is assumed that y =

q

scaling parameter. In this work, ϵ =

D

2(ξ+ξp) , which corresponds to the root-mean-square

displacement of a linear system that is obtained by removing the IPVA from the system and

24

 (rad)p()=.7=.65=.6=.8=1()()()assuming the viscous damping ratio is ξ + ξp. Substituting y = ϵ¯y into (2.8) and expanding

in a Taylor series up to order r = 7 with respect to ϵ results in

¯a (¯y, ϕ0) = a (ϵ¯y, ϕ0)|

ϵ=0 +

(cid:18)

rX

ϵi ×

¯Σ (¯y, ϕ0, τ ) = Σ (ϵ¯y, ϕ0, τ )|

i=1

rX

ϵ=0 +

i=1

ϵi ×

(cid:19)(cid:12)
(cid:12)
(cid:12)
(cid:12)

∂ia (ϵ¯y, ϕ0)
i! ∂ϵi

(cid:18)

ϵ=0
∂iΣ (ϵ¯y, ϕ0, τ )
i! ∂ϵi

(cid:19)(cid:12)
(cid:12)
(cid:12)
(cid:12)

ϵ=0

(2.18)

where a and Σ are defined in (4.14). Stationary moment equations with moments mpqrs =

⟨yp

1yq

2yr

3ys
4

⟩ up to order n = 4 are then generated from (2.18) by the cumulant-neglect tech-

nique [82, 83]. All cumulants of higher order than n = 4 are neglected by assuming that

they have negligible contribution. To this end, cumulants of orders higher than n = 4 are

set to zero to obtain expressions of higher order moments in terms of lower orders, resulting

in 69 nonlinear algebraic equations governing all moments up to fourth order, while addi-

tionally having a dependence on ϕ0; see the appendix for more information. The system is

thus underdetermined unless ϕ0 is specified. Using the solution for ϕ0 from (2.14), ϕ0 can

be substituted into the resultant 69 equations to allow for a solution for the 69 unknown

moment variables.

Choosing the parameters associated with Fig. 2.2(a), corresponding to bifurcation oc-

curring at η = .65, the ϕ0(η) values predicted thereof was substituted into the 69 moment

equations to calculate the mean square velocities ˙θ and ˙ϕ of the ball screw and pendulum,

respectively. The results are plotted in Fig. 2.5. Because ˙θ is proportional to the velocity of

the primary structure via ˙x = R ˙θ, it is a measure of the kinetic energy of the primary struc-

ture. As such, we can examine how the energy transfers to the pendulum as a function of η.

To validate the results, the mean square velocities computed by the MCS are also plotted

in Fig. 2.5, which show indiscernible discrepancy in comparison with the results obtained by

the moment equations.

There are several things worth noting in Fig. 2.5. First, Fig. 2.5(a) shows a general

decrease in ˙θ as a result of increasing η.

It specifically reaches a minimum value at the

critical point where the bifurcation occurs (η ≈ 0.65). Meanwhile, Fig. 2.5(b) shows that

25

(a) Ball screw mean square velocity as a function of η compared to linear system response.
µr = .3, ξ = .02, ξp = .07, D = .02, ϵ = 1/3.

(b) Pendulum mean square velocity as a function of η. µr = .3, ξ = .02, ξp = .07,
D = .02, ϵ = 1/3.

Figure 2.5 Comparison of mean squares of nonlinear system obtained by MCS and cumulant-
(a) Ball screw mean square
neglect and mean square of linear system while varying η.
velocity. (b) Pendulum mean square velocity.

˙ϕ reaches a maximum at the bifurcation point. The observed trends clearly demonstrate

that the kinetic energy of the primary structure transfers to the pendulum, which reaches

a maximum at the bifurcation point. To demonstrate the efficacy of response suppression

D

E

via this energy transfer, the mean square velocity

˙θ2
l

of the ball screw without the IPVA

is plotted with solid line in Fig. 2.5(a), which is derived in the Appendix C. For a fair

comparison, the viscous damping ratio of the linear system ˆξ = ξ + ξp is set, which equals

the sum of the two damping ratios of the nonlinear system. As shown in Fig. 2.5(a), the

response of the primary structure with IPVA is significantly lower than that without IPVA,

26

especially at the bifurcation point. Second, to investigate the influence of ϕ0 on the mean

squares, the mean squares were also computed with the proposed cumulant-neglect scheme

using the ϕ0(η) obtained from the MCS (dots in Fig. 2.2(a)), plotted with crosses in Fig. 2.5.

Although the MCS and perturbation method (2.14) show discrepancy when predicting ϕ0(η)

at large η values (see Fig. 2.2(a)), it is interesting to note that the mean squares obtained

by these two methods show little discrepancy in Fig. 2.5. This implies that the discrepancy

is not large enough to cause significant errors in the mean square values, while it should be

noted that higher order moments might lose accuracy with the perturbation ϕ0 solution. This

is, however, beyond the scope of this work. Cumulant-neglect used with the perturbation ϕ0

solution is therefore sufficient for the remaining analysis in this chapter.

By choosing η = .65, which corresponds to the maximum energy transfer in Fig. 2.5, µr is

varied to understand its influence on bifurcation. The mean square ball screw and pendulum

velocities as a function of µr were calculated with the cumulant-neglect procedure with a

ϕ0(µr) substitution obtained from the perturbation method, which are shown in Fig. 2.6a and

Fig. 2.6b, respectively. It is clear that the optimal energy transfer occurs at the bifurcation

point around µr = 0.3. Larger values of µr lead to lower pendulum velocities, while smaller

values will also increase the velocity of the ball screw and hence the response of the primary

structure.

2.3 Parametric Study

Sec. 2.2 shows that the limit point bifurcation promotes the energy transfer and response

suppression for a specific noise intensity level. It is worthwhile to investigate the performance

of response suppression with different noise intensity levels. By setting λ = ξp and γ = D to

be the two bifurcation parameters, (2.17) are used to generate the bifurcation curve in the

λ − γ plane, which determines the values of λ and γ for which the limit point bifurcation

occurs. The damping ξp is chosen as a bifurcation parameter because there exist many

engineering applications that allow for variable damping. For example, Kecik and Borowiec

[84] suggested installing an electromagnetic motor at the pivot point of the pendulum, which

27

(a) Ball screw mean square velocity as a function of µr compared to linear system re-
sponse. η = .65, ξ = .02, ξp = .07, D = .02, ϵ = 1/3.

(b) Pendulum mean square velocity as a function of µr. η = .65, ξ = .02, ξp = .07,
D = .02, ϵ = 1/3.

Figure 2.6 Comparison of mean squares of nonlinear system obtained by cumulant-neglect
and mean square of linear system while varying µr. (a) Ball screw mean square velocity. (b)
Pendulum mean square velocity.

is shunted with resistive loads. As such, the pendulum damping can be passively adjusted

by varying the resistance. Using the pulse width modulated (PWM) set-up chopper, Kim

28

rand Okada [85] showed that semi-active damping control of electromagnetic transducers can

be achieved, suggesting the pendulum damping can be controlled in a semi-active manner.

(a) Bifurcation point for D(ξp) for various η values. µr = .3, ξ = .02.

(b) Bifurcation point for D(ξp) for various µr values. η = .65, ξ = .02.
Figure 2.7 Bifurcation diagram of limit points in the D − ξp plane with varying (a) η and
(b) µr.

Shown in Fig. 2.7a are then the bifurcation curves for a few values of the pendulum

length parameter (η). The bifurcation curves for a few values of µr are additionally shown

29

pD=.4=.6=.8pDr=.3r=.5r=.1in Fig. 2.7b. There are several things worth noting regarding Figs. 2.7a and 2.7b. First, the

bifurcation curves show the values of noise intensity, ξp, µr and η that lead to a pitchfork

bifurcation. It should be noted that for a chosen µr, η and noise intensity, one identifies two

possible choices of pendulum damping on each curve except on the vertex. For example,

as shown in Fig. 2.7a, for µr = 0.3, η = 0.6, and D = 0.05, ξp ≈ 0.01 and ξp ≈ 0.04

each lead to a pitchfork bifurcation. The smaller values of pendulum damping on each

curve lead to intermittent pendulum motion that is characterized by oscillations around

one equilibrium angle marked by sporadic pendulum rotations. As analysis involving such

intermittent rotation is not within the scope of this paper, we restrict our attention to larger

values of pendulum damping, i.e., the right portion of each bifurcation curve in Fig. 2.8.

Second, the bifurcation curves also show that one can adjust the pendulum damping so as

to maintain a bifurcation condition for a wide range of noise intensity, η and µr values.

Shown in Figs. 2.8(a) and 2.8(b) are the mean square velocities corresponding to points

along the right portion of the bifurcation curves in Fig. 2.7a. It generally appears as though

an increase in η reduces the pendulum velocity and only slightly influences the angular veloc-

ity of the screw. Furthermore, an increase in noise intensity D along the bifurcation curves

has the effect of monotonically increasing the total energy transfer. To further quantify vibra-

tion suppression, the mean square velocity associated with the linear system corresponding

to the best response suppression is also shown in Fig. 2.8(a). The value of damping cor-

responding to that obtained by following the η = 0.8 bifurcation curve of Fig. 2.7a was

thus chosen for the linear system. It is interesting to note that the IPVA system always

outperforms the linear system for the particular set of parameters considered herein.

After choosing parameters corresponding to the bifurcation points on the right portion

of the bifurcation curves with a noise intensity of D = 0.02, we then show the variation

in the locations of local maximum ϕ0 as a function of η in Fig. 2.9. The combined results

shown in Fig. 2.8 and Fig. 2.9 show that a bifurcation point corresponding to smaller η

is correlated with higher pendulum mean square velocity. Furthermore, the bifurcation

30

diagrams of Fig. 2.9 show that a bifurcation point corresponding to smaller η leads to a

larger increase per η in the distance between the locations of local maxima.

Moving our attention to the limit point bifurcation diagrams shown in Fig. 2.7b, we can

then repeat the analysis with variation in µr. The results for the mean square ball screw

and pendulum velocities, at various bifurcation points, are shown in Figs. 2.10a and 2.10b,

respectively. It is interesting to note that an increase in µr will lead to lower velocity of the

ball screw and primary mass, while the velocity of the pendulum will first increase and then

decrease. Hence, a larger pendulum-structure mass ratio could be warranted if vibration

suppression is of sole interest.

2.4 Discussion on Findings

In this chapter, the vibration suppression, energy transfer and bifurcation characteristics

of a linear oscillator (primary structure) incorporating an inerter-based pendulum vibration

absorber (IPVA) were investigated as a function of normalized pendulum length η and mass

ratio µr. A perturbation method was introduced in conjunction with arc-length continua-

tion to detect and track the pendulum equilibrium angles corresponding to local maxima

of the marginal probability density function (PDF) in the parameter space. The obtained

pendulum equilibrium angles were implemented in a cumulant-neglect based procedure with

a 7th order Taylor expansion about the equilibrium angles to obtain the mean squares of the

system. It was shown that the marginal PDF underwent a P-bifurcation (mono-modality to

bi-modality) at critical values of η and µr. In comparison with mean squares of the linear

structure without the IPVA, it was shown that the IPVA led to effective vibration mitiga-

tion of the structure via transferring the kinetic energy of the structure to the pendulum.

Furthermore, it was shown that the energy transfer was maximum in the neighborhood of

bifurcation. The results were validated by an MCS that was used to numerically approximate

the marginal PDF for the pendulum angle as well as the mean square values.

31

(a) Ball screw mean square velocity vs. noise intensity via cumulant-neglect for various η values.
µr = .3, ξ = .02.

(b) Pendulum mean square velocity vs. noise intensity via cumulant-neglect for various η values.
µr = .3, ξ = .02.

Figure 2.8 Mean squares of nonlinear system with parameters corresponding points along
the bifurcation curves in Fig. 2.7(a): (a) ball screw mean square velocity and (b) pendulum
mean square velocity.

32

Figure 2.9 Pendulum equilibrium angle of oscillation ϕ0(η) for various bifurcation conditions
with µr = .3, ξ = .02.

33

(a) Ball screw mean square velocity vs. noise intensity via cumulant neglect for various µr values.
η = .65, ξ = .02.

(b) Pendulum mean square velocity vs. noise intensity via cumulant neglect for various µr values.
η = .65, ξ = .02.

Figure 2.10 Mean squares of nonlinear system with parameters corresponding points along
the bifurcation curves in Fig. 2.7(b): (a) ball screw mean square velocity and (b) pendulum
mean square velocity.

34

CHAPTER 3

ENERGY HARVESTING, RIDE COMFORT AND ROAD HANDLING IN
AUTOMOBILE SUSPENSION SYSTEMS

3.1 Introduction to the Energy Harvesting Shock Absorber

In this chapter, the nonlinear device introduced in the previous chapter is implemented

in a quarter-car model suspension system in an EHSA. The questions to be answered are

whether the bifurcation phenomena demonstrated in the last chapter will still show up and

whether we can make a strong correlation between optimal energy harvesting, vibration

suppression and P-bifurcation. To this end, performance metrics such as sprung mass accel-

eration, power harvested and road handling are evaluated near bifurcation.

The rest of this chapter is organized as follows.

In Sec. 3.2, the system design with

equations of motion (EOM) are shown followed by the development of an analytical bi-

furcation theory coupled with a Wiener path integration (WPI) formulation. An efficient

perturbation-based algorithm is then defined which covers a large parameter space.

In

Sec. 3.3, a parameter study is conducted to allow the determination of power harvested, ride

comfort, and road handling metrics as a function of a pendulum length parameter which is

responsible for the nonlinearity. Electrical efficiency and energy flow are also considered.

Finally, the relevant findings are summarized in Sec. 3.5.

3.2 Stochastic Bifurcation Detection Algorithm

In the previous chapter it was shown that P-bifurcation, involving the change in the

maxima of the stationary PDF, was apparent and led to substantial transfer of the kinetic

energy of the primary structure to the pendulum. As the same device is considered in

this chapter concerning the EHSA, this phenomenon is further explored in association with

the quarter-car model. The theoretical development of the method used to predict this

P-bifurcation is presented in this section, where an efficient implementation of the Wiener

path integration formulation originally outlined in [10] is modified to fit the current problem

and augmented with a calculus of variations approach for the prediction of change in the

35

curvature of the PDF.

3.2.1 Equations of motion of the EHSA

Figure 3.1 Schematics of an IPVA-integrated quarter-car model.

The IPVA-integrated quarter-car system is shown in Fig. 3.1. As shown, xs denotes the

position of the sprung mass (vehicle body) Ms, xus the position of the unsprung mass (wheel)

Mus, and xr the road elevation. The top and bottom eyelets in Fig. 3.1 then attach to Ms

and Mus, respectively. This connection then converts the translational motion between the

eyelets to the angular motion of the ball screw inside the case. The holonomic constraint

xs − xus = Rθ then relates the suspension deflection xs − xus to the angular displacement

of the ball screw θ with R = L/2π, where L is the lead value for the ball screw. Note that

attached to the ball screw is then a carrier that carries four pendulums. Each of the planet

gears is mounted to the shaft of the corresponding pendulum so that all pendulums have

the same angular position ϕ relative to the carrier. Furthermore, they each revolve around

the sun gear such that the angular displacement of the latter is given by θ − ϕ. As the sun

gear is mounted to the generator shaft, the angular motion θ − ϕ serves as the input of the

generator.

To simplify a parameter variation study, the following dimensionless parameters and

36

normalized time were chosen in the present work.

r

µr =

u =

Mus
Ms

, ξe =

, ξ =
r

ks
Ms

mR2
p
MsR2 , ωo =
Jgn2
g
MsR2 , ωt =
τ = ω0t,

, µg =

cm
2Msωo
kt
Mus
d
dτ

= ()

′

, ψus =

xus
R

, v =

,

ce
2MsR2ωo
xr
R
r
Rp

, η =

, Ψr =

ωt
ωo

,

(3.1)

Note that Rp and r are the distance between the pendulum pivot point and center of rotation

of the carrier, and the length of the pendulum, respectively. The parameters Jg, ng, ks, and kt

represent the principal moment of inertia of the generator rotor, the generator gearbox ratio,

suspension stiffness, and tire stiffness, respectively. ω0 is then a frequency that is dependent

on the sprung mass and suspension stiffness which is chosen as the characteristic frequency

to define dimensionless time τ . Finally, ce denotes the electrical damping constant due to

the generator [86] and cm is the mechanical damping constant which quantifies mechanical

losses associated with the ball screw motion R ˙θ. The equations of motion for the system

were then derived in terms of the dimensionless parameters and dimensionless time. With

z = [θ, ϕ, ψus]T , the dimensionless EOM are as follows:

Mz′′ + Cz′ + Kz = F (z, z′, Ψr)

(3.2)

where

M =

F (z, z′, Ψr) =

and

2

6
6
6
6
4

0

B
B
B
B
@

3

7
7
7
7
5

, C =

M11 M12

M12 M22

1

0

1

0

u + 1

µrη(2ϕ′θ′ + (ϕ′)2) sin (ϕ)

−µrη (θ′)2 sin (ϕ)

uv2Ψr

2

6
6
6
6
4

1

C
C
C
C
A

2 (ξ + ξe) −2ξe 0

−2ξe

0

2ξe

0

0

0

3

7
7
7
7
5

, K =

2

6
6
6
6
4

1 0

0 0

0

0

0 0 uv2

3

7
7
7
7
5

,

(3.3)

M11 = µg + µr

(cid:0)

η2 + 2η cos(ϕ) + 1

(cid:1)

+ 1, M12 = ηµr(η + cos(ϕ)) − µg, M22 = η2µr + µg

37

(3.4)

Ψr corresponds to the dimensionless road profile acting as the excitation. The EOM for the

road excitation is given here

Ψ′

r + vcΨr = n(τ )

(3.5)

where n(τ ) = W (τ )

ωoR is a dimensionless random disturbance where W (τ ) is a white noise

process. n (τ ) has expectation E[n(τ )] = 0 and auto-correlation function E [n(τ )n(ϵ + τ )] =

dδ (ϵ) and d = 2πGrV
ω0R2

is the dimensionless noise intensity. Furthermore, Gr quantifies the road

roughness and V is the driving speed as described in [86]. E[·] represents the expected value

of quantity “·”. Finally, ωc is the cutoff frequency that allows for the power spectral density

(PSD) of the road excitation to remain bounded when the frequency ω equals zero. Then vc

is a dimensionless cutoff frequency defined by the ratio of the physical cutoff frequency ωc

to ω0. The dimensionless power spectral density (PSD) is given by

SΨr =

d
(ω/ω0)2 + v2
c

(3.6)

3.2.2 The Wiener Path Integration Implementation

Following similar notation and framework as in [18, 19], the transition PDF p(yf , τf |yi, τi)

h

i

(cid:1)′

(cid:0)

zT

T

starting at τi as yi and ending

gives the probability density of all states y =

zT , Ψr,

at τf as yf . It is written as

Z

p(yf , τf |yi, τi) =

C{yf ,τf ;yi,τi}

W [y(τ )][dy(τ )]

(3.7)

where W [y(τ )] is the probability density functional on the space of all possible paths

C{yf , τf ; yi, τi} starting at (yi, τi) and terminating at (yf , τf ) with yi = y (τi) and yf =

y (τf ).

For a white noise process n (τ ), the probability density functional is well known and can

be expressed as

W [n (τ )] = C exp

(cid:20)

Z

−

(cid:21)

τf

τi

n2 (τ )
2d

dτ

(3.8)

38

where C is a normalization constant that induces a total probability of one [18]. Substituting

the expression for the white noise as a function of Ψr and Ψ′

r given in (3.5), one obtains

W [y(τ )] = C exp

"

(cid:20)

−

= C exp

−

Z

τf

Z

τi

τf

τi

#

(Ψ′

r + vcΨr)2
2d

(cid:21)

dτ

,

L (Ψr) dτ

(3.9)

where L is then the Lagrangian of the system and is in general a functional of the states

and state derivatives of the system. It is evident from (3.9) that W [y (τ )] can be maximized

by finding a path y∗ that minimizes L, which is also known as the most probable path.

Substituting the maximized W [y∗ (τ )] into (3.7) and isolating the most probable path will

yield the first-order approximation of the PDF [87]

p(yf , τf |yi, τi) ≈ C exp

−

(cid:20)

Z

(cid:21)

τf

τi

L (Ψ∗

r) dτ

As the definition of L in (3.9) does not explicitly depend on z and z′ in (4.9), it is insufficient

to use L to find y∗. This can be remedied by converting the minimization problem into

a constrained minimization problem. To that end, an auxiliary Lagrangian L∗ (y, λ) =

L (Ψr) + λT (τ )Φ (Ψr, z) is introduced, where the dynamic constraints are

Φ (Ψr, z) = M z′′ + Cz′ + Kz − F (z, z′, Ψr) = 0

(3.10)

and λ(τ ) denotes the time dependent Lagrange multiplier vector [10]. From the calculus of

variations, y∗ that minimizes L∗ can be found by solving the following constrained variational

equation

Z

τf

δ

τi

L∗ (y, λ) dτ = 0

which naturally leads to the Euler-Lagrange equations given here.

∂L∗
∂x
∂L∗
∂λ

− d
dτ

∂L∗
∂x′ +

d2
dτ 2

∂L∗
∂x′′ = 0,

= 0

39

(3.11)

(3.12)

with generalized coordinates x =

(cid:3)

(cid:2)

zT , Ψr

T . The associated boundary conditions are

x(0) = x0,

x′(0) = x′

x(τf ) = xf ,
0, x′(τf ) = x′
f ,

(3.13)

λ(0) = λ0, λ(τf ) = λf ,
0, λ′(τf ) = λ′
λ′(0) = λ′

f

Note that there are a total of seven second-order differential equations in (3.12) and therefore

only fourteen of the boundary conditions in (3.13) can be fully specified. The rest of the

boundary values are determined via the solution process automatically. Finally, the solution

will yield the most probable path of all seven states y from (3.2) through (3.5) in addition

to new pseudo states

Ψ′

r, λT ,

λT

h

(cid:0)

i

(cid:1)′

T

which are solely associated with (3.12). Finally, the

approximate transition PDF for y is given by

p(yf , τf |yi, τi) ≈ C exp

−

(cid:20)

Z

(cid:21)

τf

τi

L (Ψ∗

r) dτ

(3.14)

where Ψ∗

r denotes the most probable path of Ψr. It should be mentioned that the solution of

(3.14) is only the first-order approximation of (3.7) and improvements made by higher order

approximations are left for future work. The second order approximation would involve a

second order functional expansion of the functional in (3.14). Ref.

[87] outlines this well

for the interested reader. Finally, it is worth mentioning that the solution to (3.12) and

(3.13) with traditional boundary value problem solvers such as MATLAB bvp4c will require

additional differentiation if each produced equation does not have the same order of highest

derivative. In such a case, one must turn to other methods to approximate the solution to

(3.11) such as with Ritz method; see [88] for more detailed information. However, for the

case of the current quarter-car model, the highest derivative is uniformly equal to two.

3.2.3 Variation-Based Approach for Curvature Estimation

This section is focused on predicting the P-bifurcation of the PDF. P-bifurcation involves

the change in location of PDF extrema and thus in the present study, the P-bifurcation is

explored near the stationary point ϕf = 0 due to symmetry in the probability distribution.

40

The formulation for the proposed algorithm is outlined in this section. Given the knowledge

of the IPVA system dynamics demonstrated in the previous chapter while taking into account

that only the pendulum will oscillate locally around nonzero values, it is apparent that the

present system will have a monomodal PDF with respect to all states except ϕ if bifurcation

is achieved. As such, one can look at the change in sign of the second derivative ∂2p
∂ϕ2
f

, given the

first derivative ∂p
∂ϕf

= 0, where the approximate value of p is given in (3.14) and ϕf = ϕ (τf ).

In a system with symmetry, one can look solely at the PDF around ϕf = 0. Otherwise, one

may choose a general ϕf of interest in which the second derivative should vanish.

To compute the derivatives, however, the WPI method would involve computing a sur-

rogate model of p [89] where one would solve (3.12) for numerous paths, which can result

in considerable computational costs. In this work, a more efficient approach based on the

calculus of variations for determining the derivatives is developed, which of the formulation

is derived as follows.

First, let us define all pseudo and true states as a vector q =

rewrite (3.12) and (3.13) as

h

yT , Ψ′

r, λT ,

i

T

(cid:1)′

(cid:0)

λT

. Then

q′ = g (q) ,

q (τf ) = qf ,

q (0) = qi

(3.15)

where g is a nonlinear function of q and only 14 of the boundary conditions can be specified.

Then taking variation δ of (3.15) with respect to the boundary qf gives

∂q′
∂qf

δqf =

∂g (q)
∂qf

δqf

After interchanging derivatives ∂/∂qf and ∂/∂τ it becomes

u′ =

∂g (q)
∂q

u

(3.16)

(3.17)

where u = [u1, u2, · · · , un] = ∂q
∂qf

, with ui ∈ Rn×1 for i = 1, 2, · · · , n. Then u ∈ Rn×n given

q ∈ Rn×1, where n = 14 for the current problem. The necessary boundary conditions are

41

then simply

u (τf ) = I

(3.18)

where I denotes the identity matrix of size n. Similarly, taking variation of (3.17) one obtains
(cid:19)

(cid:18)

w′ =

i

h

T

∂
∂qf

⊗ ∂g (q)
∂q

(cid:20)

(I ⊗ u) +

w

∂g
∂q
(cid:21)

(3.19)

with qf =

f , q(2)
q(1)
∂
∂qf
value of the i-th final state, w = ∂
∂qf

f , · · · , q(n)

,

f

=

∂
∂q(1)
f

,

∂
∂q(2)
f

, · · · ,

∂
∂q(n)
f

where q(i)
f

corresponds to the

⊗ u, w ∈ Rn×n2 and ⊗ is the Kronecker tensor product.

The right boundary condition is then

w (τf ) = 0

(3.20)

where 0 denotes the zero matrix of size n×n2. Solving (3.15)–(3.20) together, one can obtain

the most probable path and its first and second variations.

Note that in the present case, the interest is only in ∂2p
∂ϕ2
f

. Suppose the i-th component

qi of q corresponds to the pendulum angle ϕ, i.e., qi = ϕ, then set ui = v and only solve

v′ = ∂g

∂q v in (3.17) and set ˆw = ∂v
∂ϕf

. Then in indicial notation, the variational equations are

given by

v′
l =

ˆw′

l =

vk

∂gl (q)
∂qk
∂2gl (q)
∂qj∂qk

vkvj +

∂gl (q)
∂qj

ˆwj, j, k, l = 1, 2, · · · , n

(3.21)

where the subscript l denotes the l-th component. One then has the following boundary

conditions

vl (τf ) = δli

ˆwl (τf ) = 0, l = 1, 2, · · · , n

(3.22)

where Einstein summation is employed, i is the index corresponding to ϕf and δli is the

Kronecker delta function. Finally, taking the derivative of (3.14) with respect to ϕf , one

42

obtains

∂p
∂ϕf

= −p

Z

τf

0

∂L
∂ϕf

dτ

, while the second derivative gives

∂2p
∂ϕ2
f

= p

"(cid:18)Z

τf

∂L
∂ϕf

(cid:19)

2

dτ

−

#

Z

τf

0

∂2L
∂ϕ2
f

dτ

+ Ψ′

r
d

∂2Ψ′
r
∂ϕ2
f

. Note that Ψ′

r, ∂Ψ′

r
∂ϕf

0

(cid:17)

2

where ∂L
∂ϕf

= Ψ′

r
d

∂Ψ′
r
∂ϕf

and ∂2L
∂ϕ2
f

= 1
d

(cid:16)

∂Ψ′
r
∂ϕf

(3.23)

(3.24)

, and ∂2Ψ′
r
∂ϕ2
f

are solutions

obtained from solving (3.15), (3.21) and (3.22) together. Also, we have the cutoff frequency

vc = 0 as it is unnecessary in the context of computing the PDF approximation and deriva-

tives. When the PDF is stationary one can further note that ∂p
∂ϕf

ϕf =0 ≈ 0 as the current

system has a symmetric PDF about the origin and a maximum or minimum is bound to

(cid:12)
(cid:12)

occur at the center of symmetry. Furthermore, if ∂2p
∂ϕ2
f

β =

(cid:12)
(cid:12)
(cid:12)
(cid:12)

ϕf =0

ϕf =0

∂p
∂ϕf
∂2p
∂ϕ2
f

(cid:12)
(cid:12)

ϕf =0 ̸= 0, the ratio

(3.25)

can be used as a measure of convergence where smaller |β| implies better convergence.

Now that the equations are defined. Fig. 3.2 is the pseudocode for the proposed WPI

algorithm. There are a few things to note in the proposed algorithm. First, (3.12) will

naturally give a system in which negative damping proportional to ξe and ξ exists due to

the term − d
dτ

∂L∗
∂x′ . This makes the system unstable with forward marching schemes such as

one used with the shooting method [90]. Additionally, it had proven to be difficult to find

a suitable initial solution guess that allowed for sufficient convergence using the collocation

method with solvers such as MATLAB bvp4c.

The authors have overcome this difficulty with the following procedure. Perform a back-

ward integration initially with ξe = 0. For the initial conditions, choose all final states

y (τf ) = 0 because the interest is in examining the PDF when all final states are trivial as the

current system has a symmetric PDF. For other cases, the final states should be chosen to re-
r (τf ) , λT (τf ) , λ′T (τf )],

flect the PDF location of interest. Out of seven final pseudo-states, [Ψ′

43

choose one nonzero but small value such that the solution remains bounded on the time in-

terval [0, τf ]. The remaining final pseudo-states are set to zero. Ultimately, for sufficiently

long integration, it is assumed that the response process is well mixed such that the initial

states become all immaterial and do not affect the stationary PDF. They are consequentially

left to be determined from the solution process. Finally, this solution is used as an initial

guess for the solution corresponding to small, but nonzero ξe using the collocation method

with an appropriate boundary value problem (bvp) solver.

Next, the algorithm takes a perturbation-based approach which accelerates the solution

convergence. As such, following a small perturbation of parameters, the preceding solution

is used as the initial guess for the current solution. This is done in an iterative fashion.

Additionally, it was determined that some dependent variables may grow unbounded and

result in failure or inefficiency in computation. This was remedied by setting a maximum

dependent variable value of σ and then multiplying the solution time vector by some factor

(cid:12)
(cid:12)

N < 1 before passing the solution as the initial guess for the next iteration wherein the

integration time is shortened by the same factor N . Lastly, the algorithm only saves ∂2p
∂ϕ2
f

ϕf =0
if |β| is smaller than some specified tolerance βub, implying sufficient convergence. In general,
βub ≪ 1. In practice, it is recommended to initially store all values of ∂2p
∂ϕ2
f

ϕf =0 to allow

(cid:12)
(cid:12)

for the flexibility in changing βub based on the qualitative trends in the solutions. If many

outliers exist, the tolerance is decreased. If not enough solution points meet the tolerance

constraint, the tolerance is increased. A compromise is thereby determined.

3.3 Parametric Study and Performance Evaluation

3.3.1 One-Parameter P-bifurcation

With the curvature estimated method formulated in the previous section, the focus now

moves on to the parametric study wherein the pendulum length parameter η is varied and the

performance metrics given by power harvested, sprung mass acceleration and road handling

are quantified.

For the variation in η study, start by assigning the dimensionless parameters given in

44

v,0]T

v = [ξe, αT
v ]

r(τf ), λT (τf ), λ′T (τf )]T = [Ψ′

Define variable parameters: ΓT
Set varying parameter initial values: Γv,0 = [0, αT
Set domain: Γv ∈ [Γv,L, Γv,R], Γv,L ≤ Γv,0 ≤ Γv,R
Set fixed parameter values: Γf = Γf,0
Set y(τf ) = 0
[Ψ′
Solve (3.15), (3.21), and (3.22) together using backwards integration.
sol = solution to boundary value problems
sol.z = state solution vector
sol.time = time vector
BC = True States and Variations at τ = τf and Pseudo-States at τ = 0
while Γv,L ≤ Γv ≤ Γv,R do

r(τf )f , λT (τf )f , λ′T (τf )f ]T ,

1: function WPIAlgorithm
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:

Γv = Γv + ϵ
N ∈ (0, 1)
σ = maximum desired value for any dependent variable
while solution not obtained do

▷ ϵ is a small perturbation in parameters
▷ scaling factor

Solve (3.15) and (3.21) with BC and sol as guess for collocation method with

any bvp solver, e.g., MATLAB bvp4c

if bvp solver fails (singularity in Jacobian of collocation equations) or does

not finish in a reasonable time then
sol.time = N · sol.time
solution not obtained

end if
end while
sol = solution obtained
if σ ≤ max(sol.z) then

sol.time = N · sol.time

19:

20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:

31:

end if
save sol
set βub ≪ 1
if |β|< βub then
save ∂2p
∂ϕ2
f

(cid:12)
(cid:12)
(cid:12)

ϕf =0

end if
32:
end while
33:
34: end function

▷ Save the solution
▷ Upper bound of |β|

▷ Save the second derivative value,

Figure 3.2 Pseudocode for WPI Algorithm.

(3.1) numerical values which can be reasonably correlated with a passenger car and a random

excitation intensity correlated with a class C road as defined in [86]. As such, all dimension-

less parameter values with the exception of µg given in Table 3.1 have been derived from

recent work in [62]. Additionally, in this chapter, V = 113 km/hr which serves as a reason-

45

able upper bound for a passenger car driving on an average road. Also, note that the value

of dimensionless µg was determined by considering the Maxon 353300 DC motor which has

a rated power of 250W with a gearbox giving a 3.8:1 ratio, so that the speed would be in the

nominal range. The motor also has an internal resistance Rint = 1.06Ω which will be needed

when accounting for electrical efficiency.

Some important physical parameters were then

Parameters, Dimensionless Value

µr
µg
ξ
ξe
u
v
d

0.1521
0.0194511
.02
0.0337
0.14
4.4137
.0806

Table 3.1 Dimensionless Parameters of IPVA system.

determined by assuming the ball screw lead value of L = 100 millimeters. From companies

such as THK Co., Ltd., this appears to be the largest value immediately offered. A large

lead value will help keep the dimensionless intensity at a minimum when rougher roads are

considered later on. Some physical parameters are given in Table 3.2.

Starting with the

Parameters Value

Meaning

Ms
Mus
kt
ks
Rp
m
cm
V
Gr
R
κ
Rint

Sprung Mass
Unsprung Mass
Tire stiffness
Suspension stiffness
Carrier radius
Mass of pendulum
Mechanical damping
Driving speed

374 kg
52.36 kg
248,085 N/m
90,963 N/m
.05 m
3.74 kg
233.74 Ns/m
31.3 m/s
64/(2π)2 × 10−6 m· cycle Road roughness coefficient
Effective ball screw radius
.1/(2π) m
Motor torque constant
0.207 Nm/A
Motor internal resistance
1.06 Ω

Table 3.2 Phyical Parameters of IPVA system.

WPI method introduced in the previous section, sgn

(cid:16)

(cid:17)

(cid:16)

∂2p
∂ϕ2
f

log10

1 +

(cid:17)

(cid:12)
(cid:12)
(cid:12)

(cid:12)
(cid:12)
(cid:12) 1
β

is plotted against

46

dimensionless pendulum length η in Fig. 3.3a. This function was chosen to clearly show the

qualitative characteristics of the second derivative of the marginal PDF as the magnitudes

of the second derivative vary drastically. Additionally, β−1 is used in place of the second

derivative in order to eliminate the unknown integration constant C.

In Fig. 3.3a there is a clear indication that the sign turns positive as η is increased past

0.84. This indicates an increase in curvature and more importantly a bimodal probability

density function. This is verified with the marginal PDF plot for ϕ in Fig. 3.3b, obtained

via an MCS with ∼ 7 · 104 realizations. Note that a marginal PDF is generally obtained by

integrating the PDF with respect to all variables besides the variable of interest. The PDF

is then normalized with respect to the entire integral such that the total probability is unity.

The marginal PDF of ϕ is chosen in this work as this is where the bimodal distribution is

clearly seen. The marginal PDF then simply gives the probability distribution for the pen-

dulum angular position relative to the carrier. In this case, it is clear that the PDF becomes

bimodal after η = 0.7. For η = 0.9, the PDF has two distinct peaks.

It should also be

noted that a 50-second integration time was chosen to ensure that the approximate solution

for the second derivative was close to stationary and reasonably accurate in a qualitative

sense. While lower integration times of only 10 seconds were tried initially (not shown in

this work), they were qualitatively accurate only for subzero curvature solutions. However,

this fact has been exploited by the authors in the bifurcation analysis later on in Sec. 3.4.

i.e. the initial integration time is decreased in the case of low electrical damping as this

most likely will result in rotation of the pendulum and thus near-zero and possibly negative

curvature. This helps immensely with the efficiency of the algorithm. Next, we will examine

how the P-bifurcation influences the suspension performance.

3.3.2 Multi-objective Suspension Performance

The power from the IPVA into the electrical domain (Pnl) is directly proportional to

generator velocity squared, where the proportionality constant is the electrical damping

coefficient ce = κ2n2

g/Re. Here κ is the torque constant of the motor, ng is the gear ratio

47

(a)

(b)

Figure 3.3 (a) Logarithm of the scaled and shifted magnitude of the second derivative of the
marginal PDF vs. η. The integration time is 50 seconds. (b) Marginal PDF for pendulum
angle ϕ for various η values.

and Re is the total electrical resistance [91]. The expected value of the electrical power is

then just

and in dimensionless form it is

(cid:20)

(cid:16)

Pnl = E

ce

˙θ − ˙ϕ

(cid:21)

(cid:17)

2

h

i

ˆPnl = E

ξe (θ′ − ϕ′)2

(3.26)

(3.27)

Upon consideration of losses due to internal resistance, one can write the expected value

of the harvested dimensionless power ˆPnl,H as

ˆPnl,H = E

(cid:20)

Rl
Rl + Rint

(cid:21)

ξe (θ′ − ϕ′)2

(3.28)

where ξe now represents the electrical damping coefficient with internal electrical resistance

Rint as well as load resistance Rl. Rl can further be written in terms of ξe so that (3.28)

becomes

(cid:20)(cid:18)

ˆPnl,H = E

ξe − 2R2ω0RintMs

κ2n2
g

(cid:19)

(cid:21)

ξ2
e

(θ′ − ϕ′)2

(3.29)

It is apparent from (3.29) that losses can be reduced with a decrease in internal resistance

Rint which should have a large influence on the choice of generator.

48

The formula for the linear system mean power corresponding to the case of locked pen-

dulum such that ϕ = ϕ′ = 0 is identical to that derived in [62] and is given here.

Pl =

ceπGrV kt
ce + cmR2

(3.30)

In a similar fashion to the nonlinear case, the linear system dimensionless mean power can

be written with electrical efficiency considerations as

ˆPl,H = ηe

ξeduv2
4 (ξe + ξm)

where electrical efficiency ηe is given by

ηe =

Rl
Rl + Rint

(3.31)

(3.32)

For the linear system, there exists an optimal load resistance giving optimal harvested power
l,H . For a given set of parameters, this can be determined via ∂ ˆPl
ˆP opt

= 0. Substituting the

∂Rl

optimal load resistance in (3.31), one arrives at

dκ2uv2n2
g

ˆP opt

l,H =

p

κ2n2

g + 2R2ω0RintξmMs

(A) (B)

where

√

(cid:16)q

p

(cid:17)

A = 2

2

B = κ2n2

g +

g + 4ξmR2ω0RintMs + 2R
2κ2n2
q
√

(cid:0)

ξmω0RintMs
(cid:1)

2R

ξmω0RintMs

κ2n2

g + 2ξmR2ω0RintMs

+ 2ξmR2ω0RintMs

(3.33)

(3.34)

Now the mean electrical power is plotted against η in Fig. 3.4a, normalized with respect

to the optimal linear system. Given an optimized linear system with a mean power of

Pl = 21.0932 Watts, a clear increase of 40% is shown in Fig. 3.4a for η = 0.9, in the bimodal

regime.

In order to quantify ride comfort, Fig. 3.4a also shows the root mean square (RMS) sprung

mass acceleration normalized with respect to the linear system. The detailed derivation of

the root mean square sprung mass acceleration for the linear system with no pendulum

49

(a)

(b)

Figure 3.4 (a) Harvested power and RMS sprung mass acceleration vs. η for the proposed
IPVA-integrated system, normalized with the optimal linear benchmark system. Note that
P opt
l,H = 21.0932 W, and the electrical efficiency of the nonlinear system is 83%. (b) Road
handling index vs. η for the proposed IPVA-integrated system normalized with respect to
linear system γRHI,l. Solid vertical line indicates P-bifurcation.

is given in [62]. Note that as the present work considers a locked pendulum rather than

no pendulum at all, the solution is slightly different than in [62]. However, the solution

procedure is the same and the solution is omitted here for brevity. It is clear from Fig. 3.4a

that the ride comfort improves with an increase in pendulum length parameter η. There

is an approximate 60% reduction for η = 0.9.

It is finally noted that the simultaneous

optimization of power and ride comfort is a major hallmark of the current design; the non-

linearity allows for simultaneous objectives which are otherwise competing when the linear

system is used [92].

Another metric worth exploring is the road handling index (RHI)[86] which will be given

the symbol γRHI in this paper. This is the ratio of root mean square dynamic tire force to

static weight from the car. Lower numbers indicate safer driving conditions as the car makes

better contact with the road. This mathematically is stated as

γRHI = RMS

(cid:18)

ktR (ψus − Ψr)
(Mus + Ms) g

(cid:19)

(3.35)

where RMS(·) denotes the root mean square value of “ · ” and g is the acceleration due to

gravity. This quantity normalized with respect to the linear system value γRHI,l is plotted vs

50

η in Fig. 3.4b. Note that the detailed formulation of the linear system road handling index is

given in Appendix D. With the IPVA system, there is apparently simultaneous improvement

in γRHI for larger η when the system is bimodal.

It should finally be noted that the P-

bifurcation condition with bifurcation parameter η does not guarantee optimal harvested

power, sprung mass acceleration, or road handling. However, it serves as a suitable starting

point for further optimization with respect to η. Note that the power harvested seems to

approach a horizontal asymptote near this P-bifurcation in Fig. 3.4a. In this case, it is quickly

discovered that larger choices of η are desired as the other metrics continue to improve.

Returning to the subject of power, let us investigate power generation further with power

spectral density (PSD) analysis. The PSD associated with the generator velocity is plotted

vs. normalized frequency ω/ω0 for the monomodal system in Fig. 3.5a and the bimodal

system in Fig. 3.5b. Note that the two natural frequencies ωn1 and ωn2 of the linear system

with locked pendulum η = 0.9 have been determined to be 0.74ω0 and 2.54ω0, respectively. It

is then noted, that for both η cases the introduction of the IPVA has the effect of annihilating

ωn2 which simultaneously has the effect of broadening the effective bandwidth of operation.

Also, as the power harvested is directly proportional to the area under the curve, both cases

show improvement over the linear system while power is notably increased in the bimodal

case. Specifically, for ω/ω0 < 1 and ω/ω0 > 3 the total area is noticeably larger. Therefore,

these high-frequency and low-frequency areas are responsible for significant improvements in

energy harvesting. In fact, the total area increase for the frequency range shown is 23% for

the monomodal case and 34% for the bimodal case.

3.3.3 Energy Flow

The power input to the IPVA, each component in the IPVA, and the electrical domain

are analyzed in this section to better understand the energy flow in the system as a function

of η.

Starting with the power input to the IPVA (PIP V A), the result derived in [62] is used.

51

(a)

(b)

Figure 3.5 (a) Power spectral density of the generator velocity. η = 0.5. (b) Power spectral
density of the generator velocity. η = 0.9. ω0 =

ks/Ms.

p

Namely,

(cid:16)

(cid:17)

PIP V A = −

Ms ¨xs + cmR ˙θ + ksRθ

R ˙θ

(3.36)

To determine the power input to the pendulum (Tc) first write the energy in physical

units as

(cid:18)

Tc =

1
2

m

(cid:16)

R2
p

˙θ2 + r2

˙θ + ˙ϕ

(cid:17)

2

(cid:16)

(cid:17)(cid:19)

+ 2Rpr cos (ϕ) ˙θ

˙θ + ˙ϕ

(3.37)

Then the power input to the pendulum, Pp = dTc

dt , is simply

h

(cid:16)

(cid:17) (cid:16)

Pp = m

r2

+rRp

˙θ cos(ϕ)

˙θ + ˙ϕ
(cid:16)

(cid:17)
¨θ + ¨ϕ

(cid:17)
¨θ + ¨ϕ

− rRp

˙θ ˙ϕ sin(ϕ)
(cid:17)
(cid:16)

(cid:16)

(cid:17)

˙θ + ˙ϕ

+ rRp

¨θ cos(ϕ)

˙θ + ˙ϕ

+ R2
p

i
˙θ ¨θ

The mechanical energy of the generator rotor is then

and the power is

Tg =

(cid:16)

Jgn2
g
2

(cid:17)

2

˙θ − ˙ϕ

(cid:16)

Pg = Jgn2
g

˙θ − ˙ϕ

(cid:17) (cid:16)

(cid:17)
¨θ − ¨ϕ

(3.38)

(3.39)

(3.40)

Finally, the expected value of power into the electrical domain is given in (3.26). Note that

the second derivatives in (3.36), (3.38) and (3.40) are written in terms of the system states as

determined by (3.2-3.4) before they are made dimensionless. The substitutions are omitted

52

Figure 3.6 Mean input power into IPVA, generator rotor, pendulum, and electrical domain
vs. η.

here as they are too lengthy. For a range of η, the mean powers are then estimated with

the statistical expectation derived from an ensemble of ∼ 7 · 104 realizations. The results

are shown in Fig. 3.6. It is apparent that for all η, most of the energy that is input to the

IPVA flows into the electrical domain. The total power also increases with an increase in the

dimensionless pendulum length η. Additionally, there is practically zero net energy flow in

the pendulum and generator rotor. These components are acting to only transfer the power

into the electrical domain.

3.4 Two-parameter P-bifurcation and Performance Evaluation Around the Bi-

furcation Boundary

As demonstrated in Sec. 3.3, the P-bifurcation is critical to the suspension performance.

Therefore, we study it in-depth in this section. Now that the WPI algorithm efficiently
(cid:17)

(cid:16)

estimates ∂2p
∂ϕ2
f

, a simple bifurcation detection algorithm which checks for change in sgn

∂2p
∂ϕ2
f

is implemented. Note that for this study, the dimensionless pendulum length is set to

η = 0.97 as it was shown in the preceding section that large pendulum lengths might be

53

correlated with energy transfer. Also, to allow for rotation one must have η < 1. For η > 1,

the pendulum will collide with the ball screw. Also, the sign change is determined as a

Parameters Value

Meaning

Ms
Mus
kt
ks
Rp
m
cm
R
κ
Rint

Sprung Mass
Unsprung Mass

374 kg
52.36 kg
300263 N/m Tire stiffness
110095 N/m Suspension stiffness
.062 m
3.74 kg
256.67 Ns/m Mechanical damping
.1/(2π) m
0.207 Nm/A Motor torque constant
1.06 Ω

Carrier radius
Mass of pendulum

Motor internal resistance

Effective ball screw radius

Table 3.3 Physical Parameters of the IPVA System on Rough Class F Road.

function of dimensionless noise intensity (d) and electrical damping ratio ξe. For this study,

some physical parameters are also given in Table 5.1. Note that the suspension stiffness and

tire stiffness has been increased in this case to allow for d < 0.9 when traveling on a rough

class F road with a velocity of 5 m/s while keeping all dimensionless parameters (excluding

d and ξe) the same as in Table 3.1.

Using the developed WPI-based curvature estimation algorithm, a bifurcation curve has

been produced and is shown in Fig. 3.7a. Three distinct regions have been determined

based on the WPI approximation and verified with an MCS of ∼ 3.6 · 104 realizations and a

simulation time of τ = 1000 for a few different car velocities and damping values. Given a

parameter set in region I, the system is expected to primarily exhibit rotation and thus the

PDF will have near-zero curvature at ϕ = 0. Region II and III then correspond to monomodal

and bimodal PDFs, respectively, while it should be noted that rotations in region II are still

possible. In the case of rotations in region II, oscillations are likely to occur near ϕ = 2πk

following rotation, where k is an integer.

It should finally be noted that the dashed line

was estimated with the use of the MCS. The dashed line prediction was necessary as no

detectable change in the sign was determined with the WPI approximation and convergence

54

criteria used in the algorithm. Additionally, the noise intensity stops at around d = 0.9 as it

was proven difficult to observe any further qualitative distinctions past this point given the

fact that rotation is highly probable.

To demonstrate these assertions, let us define the value of ξe lying on the bifurcation

curve to be ξe,bif and ξe ∈ [0.5ξe,bif , 1.5ξe,bif ]. Then the marginal PDFs for pendulum angle

ϕ for three distinct velocities corresponding to bifurcation points in the bifurcation diagram

are shown in Figs. 3.7b-3.7d. Note that the PDFs extend to very large angles in the

positive and negative directions. The domain has been cropped only for visualization near

the origin and the area under all PDFs for the entire domain is in fact equal to unity. Next,

note that Figs. 3.7b and 3.7c correspond to the case of the PDF transitioning between

bimodal and monomodal (region III to region II) when the driving speed is V = 0.1973 m/s

and V = 0.5371 m/s respectively. In both cases, a PDF with near zero curvature at the

approximate bifurcation point is apparent. A major difference is then in the bimodal PDFs

where the larger velocity leads to a larger distance between peaks. Moving to Fig. 3.7d one

sees the case of the PDF transitioning between monomodal and pure rotation (region II to

region I) when the driving speed is V = 2.6432 m/s. In this case, it is interesting to note

that the PDF at the approximate bifurcation point actually takes the form of what appears

to be a monomodal type whose curvature is nearly zero at the origin. While it has proven

to be too difficult to numerically verify this curvature with the use of the MCS for sufficient

accuracy, one can easily see that the MCS results in a PDF at the predicted bifurcation

point which indeed is qualitatively different than the others, indicating a transition. Hence,

WPI has given a reasonable bifurcation prediction.

3.4.1 Performance Around P-bifurcation

In this section, we look at performance as a function of driving speed on a class F road

d = 2πGrV

in the neighborhood of P-bifurcation. As a reminder, the dimensionless noise intensity is
ω0R2 where for a class F road, one has Gr = 4096/(2π)2 × 10−6 m· cycle [86]. As a
result, one may vary d by varying the driving speed V . It should be noted that d is also

55

(a)

(b) V = 0.19728 m/s

(c) V = 0.53718 m/s

(d) V = 2.6432 m/s

≈ 0 with certain
Figure 3.7 (a) Bifurcation boundaries in d vs. ξe plane. Region I: ∂2p/∂ϕ2
f
rotation. Region II: ∂2p/∂ϕ2
f < 0, monomodal oscillations and possible rotation. Region III:
∂2p/∂ϕ2
f > 0, bimodal and rotation is less likely. The dashed line is an MCS prediction. (b)
PDF transition from region III to II with V = 0.19728 m/s. (c) PDF transition from region
III to II with V = 0.53718 m/s. (d) PDF transition from region I to II with V = 2.6432 m/s.

linearly proportional to Gr and so one might interpret the results in this section as a function

of road class with constant velocity if they wish. Furthermore, note that an upper limit of

driving speed V = 6 m/s is chosen because a class F road means a very rough off-road style

terrain. Driving speeds are generally low on off-road terrain for safety and to prevent damage

to the vehicle suspension.

Using (3.33) one can obtain the physical optimal power P opt

l,H with

P opt
l,H = 2Msω3

0R2 ˆP opt

l,H

(3.41)

56

(a)

(b)

(c)

(d)

Figure 3.8 Performance metrics around bifurcation with 50% variation in ξe. (a) Optimal linear
system power into the electrical domain.
(b) Nonlinear system mean power into the electrical
domain normalized with the optimal linear system. (c) Nonlinear system RMS sprung mass accel-
eration normalized with the optimal linear system. (d) Nonlinear system RHI normalized with the
optimal linear system.

The total power into the electrical domain for the optimized linear system is then

P opt

l =

P opt
l,H
ηe

(3.42)

The optimized linear system power that flows into the electrical domain is then shown in

Fig. 3.8a as a function of driving speed. This optimal system will now be used in normaliza-

tion to evaluate the performance of the nonlinear system around the bifurcation curve shown

in Fig. 3.7a. Let us now restrict our electrical damping domain to be ξe ∈ [0.5ξe,bif , 1.5ξe,bif ].

Moreover, by taking the ratio of of (3.26) to (3.42) at these ξe values, one can quantify

57

(a)

(b)

(c)

Figure 3.9 (a) Electrical efficiency for the nonlinear system around bifurcation.
system mean power harvested, normalized with the optimal linear system.
system harvested power.

(b) Nonlinear
(c) Optimal linear

the normalized expected value of power that flows into the electrical domain relative to the

optimized linear benchmark system when the system is near bifurcation. This is plotted vs.

ξe/ξe,bif on a class F road for six different speeds in Fig. 3.8b. It is finally apparent that for

all driving speeds shown, the energy transferred into the electrical domain is between 130%

and 143% of the optimal linear system when electrical damping is chosen to accommodate

P-bifurcation. Note that it is also apparent that the power into the electrical domain appears

to vary less in magnitude as electrical damping is increased past the predicted bifurcation

point. This might signify that an optimal solution is near.

58

In Fig. 3.8c is the RMS acceleration of the sprung mass plotted using the same parameter

values, normalized with respect to the linear system that was optimized for power. A key

result here is that the nonlinear system outperforms the linear system and allows for at least

59% better ride comfort for all driving speeds near bifurcation.

It is also noted that for

most driving speeds, the best ride comfort occurs when ξe is near the predicted bifurcation

curve. However, the case of V = 0.19728 m/s and V = 2.6432 m/s seem to break this

trend. The former and latter show that the best performance in terms of ride comfort

occurs when the system is bimodal and when the system is monomodal with intermittent

rotation, respectively. While this means that simultaneous performance objectives cannot

be completely met in the case of the V = 0.19728 m/s case, it is actually apparent that the

power into the electrical domain decreases by at most 2%. In the case of V = 2.6432 m/s,

simultaneous objectives are met if electrical damping is larger than the bifurcation value.

In Fig. 3.8d, is then the road handling index normalized with respect to the same op-

timal linear system and for the given parameter set. It is finally interesting to note that

the performance in this metric somehow qualitatively mirrors that of the ride comfort which

means the road handling performance objective is automatically met if the ride comfort per-

formance objective is met. Quantitatively, as much as 65% improvement in this metric over

the linear system is seen near the bifurcation point. Excellent road handling can therefore

be expected with the proposed nonlinear EHSA.

To further understand the power potential of the nonlinear EHSA, the electrical efficiency

determined via (3.32) is plotted vs. driving speed on a class F road in Fig. 3.9a. The result

from (3.29) is then divided by (3.33) and one can show that

Pnl,H
P opt
l,H

=

ˆPnl,H
ˆP opt
l,H

(3.43)

The corresponding harvested power is then plotted in Fig. 3.9b. However, first note that the

zero values of power and electrical efficiency correspond to when the electrical damping ξe

is greater than the upper limit ξe,max that is achieved with only internal resistance; ξe,max =

κ2n2
g
2Msω0R2Rint

. Simply put, no power can be harvested when there is zero load resistance. Using

59

parameters from Table 5.1 one finds that ξe,max = 0.1791. While energy transfer might be

large near bifurcation, the small electrical efficiency corresponding to large electrical damping

results in significant decreases in power harvested. This implies that simply choosing a point

near P-bifurcation curve will not guarantee the good power harvested; it is also important

to choose an extremely efficient generator. In the best case of ξe = 0.5ξe,bif , one can expect

around a 20% to 40% improvement over the linear system. Given the result in Fig. 3.9c, this

correlates with around 12 Watts per shock at a driving speed of 0.19728 m/s and around

412 Watts per shock at a driving speed of 5.6757 m/s.

3.5 Discussion on Findings

In this chapter, analysis was done which pertained to the performance of an IPVA-based

EHSA connected to a quarter-car suspension system while varying the dimensionless pendu-

lum length η. The performance of the device was measured relative to a linear benchmark

with locked pendulum whereby electrical efficiency was considered and the linear system

was thereby optimized for power. Excellent power harvesting, sprung mass acceleration, and

road handling of the device was determined to occur for large η when the system PDF was

bimodal, prompting the development of a Wiener path integration formulation to determine

approximately when the PDF is bimodal. Some PSD analysis showed that for normalized

frequencies of ω/ω0 < 1 and ω/ω0 > 3 the total power was much larger than that of the

linear benchmark system. Energy flow analysis also showed that more energy flowed into

the electrical domain as the pendulum length was increased.

A variation-based formulation was then adopted which, in addition to the value of the

PDF, allowed for a prediction for the estimated sign of the second derivative of the PDF at

any point of interest. An efficient bifurcation detection algorithm was thereby developed,

which resulted in the prediction of three distinct regions in the noise intensity and electri-

cal damping plane. These regions corresponded to monomodal, bimodal, and finally flat

(rotation) probability distributions. The algorithm results were verified with an MCS.

Next, the performance of the EHSA was measured against the optimized linear system

60

driving on a class F road and varying driving speeds with electrical damping which put

the system at or near the bifurcation curve. It was determined that energy transfer into the

electrical domain may be as much as 43% larger than that of the optimized linear benchmark

system when the nonlinear system is in the neighborhood of P-bifurcation. Furthermore,

harvested power can be around 20% improved, while ride comfort and road handling may be

improved by at least 59% for all investigated driving speeds. Electrical efficiency was shown

to diminish the power harvested from the nonlinear EHSA for high electrical damping while

ride comfort and road handling were much less sensitive to electrical damping variation.

Therefore, it is first important to select the most efficient generator.

In the case that a

generator is not very efficient, it is recommended to select electrical damping values less than

required for bifurcation when power harvesting is the most important metric considered.

While the bifurcation detection algorithm did not allow for complete optimization in

the performance, it serves as a means for choosing initial parameter values to be used in

a less exhaustive MCS with a full parameter sweep across the reduced realistic domain of

values that will ultimately optimize the system parameters.

It should be noted that the

proposed algorithm could be used in any case where the qualitative nature of the probability

distribution is of interest and there is uniformity in the highest derivatives of the produced

Euler-Lagrange equations. Otherwise, other approximate methods such as the Ritz method

might be substituted in the algorithm to solve the Lagrangian optimization problem.

61

CHAPTER 4

AN EXPERIMENTAL STUDY: P-BIFURCATION, ENERGY
HARVESTING, AND VIBRATION SUPPRESSION

This chapter focuses on the bimodal IPVA energy harvesting device with an application to

a SDOF structure subjected to Gaussian broadband excitation. The objective is to conduct

an experimental investigation into P-bifurcation of the device and the correlation with si-

multaneous energy harvesting and vibration suppression. Notably, this study includes the

experimental parameter characterization for a developed prototype to accommodate accurate

predictions of energy harvesting potential, vibration suppression, and bifurcation boundaries.

Adequate experimental data is used to verify predictive capabilities, confirm the bifurcation

phenomenon, and gather information regarding power harvested and vibration suppression.

4.1 The IPVA Device and System Model

In this section, the IPVA device, which was built and designed for experiment, and the

corresponding system model are introduced. The model for the system is shown in Fig. 4.1.

The system has a suspended mass Ms with a degree of freedom xs, stiffness ks, mechanical

damping cm, and an attached IPVA device with planetary gear coupling to a DC generator.

The system also has a base motion xr.

A holonomic constraint is imposed by the ball screw giving xs − xr = Rθ, where R is the

lead value of the screw divided by 2π. Figures 4.1b and 4.1c then show the bottom and top

(a)

(b)

(c)

Figure 4.1 (a) Model of the suspended mass system with IPVA device. (b) Bottom view of
carrier. (c) Top view of pendulum and gear arrangement with connection to ball screw.

62

of a carrier that holds two pendulums, each of mass mp and principal inertia Jp about their

centers of mass. Additionally, the center of mass of each pendulum is displaced at a distance

r from the respective pivot points which are at a distance Rp from the center of the carrier.

These pendulums are free to rotate with an angle ψ relative to the carrier, while the carrier

with principal moment of inertia Jc about its vertical axis is free to rotate with an angular

degree of freedom θ about the ball screw’s vertical axis. Also attached to each pendulum is

then a planetary gear of mass mgear and principal moment of inertia Jpg. These planetary

gears are then meshed with a sun gear connected to the input of a gearbox whose output

is connected to the input of a DC generator. The gearbox shaft has a degree of freedom

ψg = θ − ψ due to the holonomic constraint imposed by the planetary gear system. The

motor rotor then has a degree of freedom ngψg, where ng is the gear ratio. The principal

moment of inertia of the sun gear and generator rotor are Jsg and Jr, respectively.

4.1.1 Equations of Motion

To derive the equations of motion using Lagrange’s equations, we first derive the kinetic

energy, potential energy, and virtual work associated with the system. The kinetic energy

for the system can be broken down into energy associated with each component. The energy

associated with the suspended mass and carrier is

T1 =

1
2

Ms(R ˙θ + ˙xr)2 +

1
2

˙θ2

Jc

Then the energy associated with the pendulums and planetary gears is
#

" 

!

T2 =

1
2

2X

2X

(mpr2 + Jp) +

Jpg

˙ψ + ˙θ

2X

(cid:16)

(cid:17)

2

i=1
(cid:16)

i=1

(cid:16)

(cid:17)(cid:17)

mp

R2
p

˙θ2 + 2Rpr cos(ψ) ˙θ

˙θ + ˙ψ

+

+

i=1
2X

(cid:0)

i=1

mgearR2
p

(cid:1)

˙θ2

The energy associated with sun gear and generator rotor is given by

T3 = 1/2

(cid:0)

Jrn2

g + Jsg

(cid:1) (cid:16)

(cid:17)

2

˙θ − ˙ψ

(4.1)

63

Finally, the total kinetic energy is given by

T = T1 + T2 + T3
" 

2X

(cid:0)

!

(cid:1)

+ Jc

mgearR2
p

1
2

2X

i=1
(cid:0)

=

+

+

+

2X

i=1

i=1

(mpr2 + Jp) +

(cid:16)

2X

i=1

mp

R2
p

˙θ2 + 2Rpr cos(ψ) ˙θ
(cid:17)
(cid:1) (cid:16)

2

1
2

˙θ2
!

(cid:16)

(cid:17)

2

Jpg

˙ψ + ˙θ

(cid:16)

(cid:17)(cid:17)

˙θ + ˙ψ

(cid:21)

Jrn2

g + Jsg

˙θ − ˙ψ

+

Ms(R ˙θ + ˙xr)2

(4.2)

where ˙(•) = d (•) /dt represents the first time derivative. The potential energy is given by

V =

1
2

ks (Rθ)2

(4.3)

Furthermore, a torsional viscous damping coefficient cp is introduced to account for en-

ergy loss at the pivot point of the pendulum due to bearings. Assuming negligible internal

inductance such as in [51, 62, 86], one can also assume a generator electrical damping torque

opposing the flow of current due to Faraday’s law and proportional to the generator ro-

tor velocity with coefficient ce. Specifically, define the electrical damping imposed by the

generator as

ce =

gκ2
n2
t
Ri + Rl

(4.4)

where Ri and Rl are the internal and load electrical resistance, respectively. κt is the torque

constant of the generator. Additionally, mechanical losses in the direction opposing the sun

gear and generator rotor motion ψg can be approximated with a viscous damping torque

with coefficient cgm and coulomb friction torque in the same direction. This accounts for

mechanical energy dissipation within the generator housing as well as dissipation in the

gears. Finally, note that the virtual displacement of the mass Ms is given by δxs = Rδθ, the

virtual displacement of the pendulum in the direction of the pendulum damping torque is δψ

and the virtual displacement of the sun gear is δψg = δ (θ − ψ). It follows that the summed

64

 
virtual work due to the mechanical damping, electrical damping, and frictional torque is

written as

δQ = −cmR2 ˙θδθ − (ce + cgm)

(cid:16)

(cid:17)

(cid:16)

(cid:17)

˙θ − ˙ψ

δ (θ − ψ) − cp

˙ψδψ

− Tf

˙θ − ˙ψ

(cid:16)

(cid:17)

|

˙θ − ˙ψ

|

δ (θ − ψ)

(4.5)

where −Tf

( ˙θ− ˙ψ)
|( ˙θ− ˙ψ)| δ (θ − ψ) accounts for the friction - based losses in generator rotor and gear
motion. With (4.2), (4.3) and (4.5), the equations of motion for the system were derived

using Lagrangian methodology and are given here.

M R2 + Jc + Jrn2

g + Jsg +

2X

mpR2

p +

!

mpr2

¨θ

2X

2X

Jp +

i=1

i=1

!

¨θ

!

Jpg − Jrn2
g

− Jsg

¨ψ

2X

i=1

2X

Jpg + 2

mpRpr cos (ψ)

i=1

i=1

2X

2X

mpRpr cos (ψ) +

Jp +

(cid:1)

i=1

i=1

˙θ − (ce + cgm) ˙ψ + kR2θ

2X

i=1
2X

i=1

mgearR2

p +

mpr2 +

2X

i=1

ce + cgm + cmR2
2X

+

+

+

(cid:0)

− 2

mpRpr ˙ψ ˙θ sin (ψ) −
(cid:17)
(cid:16)

i=1

= −M ¨xrR − Tf

˙θ − ˙ψ

(cid:16)

|

˙θ − ˙ψ

(cid:17)

,

|

2X

i=1

mpRpr sin (ψ) ˙ψ2

!

2X

mpr2 +

2X

2X

Jp +

Jpg + Jrn2

g + Jsg

¨ψ

i=1

2X

i=1
2X

i=1
2X

(cid:0)

Jp +

Jpg +

mp

i=1

i=1

i=1

r2 + Rpr cos (ψ)

(cid:1)

!

− Jrn2
g

− Jsg

¨θ

+

+

2X

i=1

mpRpr sin (ψ) ˙θ2

+ (ce + cp + cgm) ˙ψ − (ce + cgm) ˙θ = Tf

(cid:16)

(cid:16)

˙θ − ˙ψ

(cid:17)

(cid:17)

|

˙θ − ˙ψ

|

(4.6)

In (4.6), ¨(•) = d2 (•) /dt2 represent the second time derivative. Let us next define ˆM as

65

 
 
 
 
 
the effective structural mass constituted by components only multiplying ¨θ. Then,

2X

i=1
(cid:18)

mgearR2

p + M R2 + Jc

=

2mgearR2

p/R2 + M +

(cid:19)

Jc
R2

R2 = ˆM R2

(4.7)

where it should be noted that here 2mgear constitutes the mass of planet gears, vertical

shafts and bearings at the pendulum pivot point combined. Thus, these moment of inertia

terms can be absorbed into the effective mass term without loss of generality. In order to

simplify parameter identification, the effects of friction torque and viscous damping torque

are combined and approximated by a single viscous damping torque in the direction opposing

the generator motion and with coefficient cf . Additional terms can be combined and the

following set of dimensionless parameters are defined.

µr =

µg =

ξe =

2mpR2
p
MsR2 , µp =
n2
gJr + Jsg
MsR2
ce
2MsR2 ω0

, α =

, ξf =

,

, η =
r

r
Rp
ks
Ms
, ξp =

2 (Jp + Jpg)
MsR2
ˆM
Ms

, ω0 =
cf
2MsR2 ω0
′′
= ω2
0

, ξ =

cm
2Ms ω0

,

cp
2MsR2 ω0

τ = ω0t, (•)

′

= ω0

˙(•), (•)

¨(•)

(4.8)

The equations of motion for the system are now rewritten in terms of the dimensionless

parameters as follows:

M (ψ) x′′ + Cx′ + Kx + g (x, x′) = f

(4.9)

where

M(ψ) =

2

6
4

M11 M12

M12 M22

3

7
5 ,

M11 = 1 + α + µp + µg + µr

(cid:0)

1 + η2 + 2η cos (ψ)

(cid:1)

M12 = µrη(η + cos (ψ)) − µg + µp

M22 = µrη2 + µg + µp

66

2

6
4

2

6
4

C =

K =

2(ξ + ξe + ξf ) −2(ξe + ξf )

−2(ξe + ξf )
0

3

2(ξe + ξf + ξp)

1

3

7
5 ,

1 0

7
5 , f =

B
@

W (τ )

C
A ,

0 0
0

0

g(x, x′) = µrη

−(2ψ′θ′ + ψ′2) sin (ψ)

B
@

θ′2 sin (ψ)

1

C
A

(4.10)

In (4.9) and (4.10), x = (θ, ψ)T , M is the dimensionless inertia matrix, C is the di-

mensionless damping matrix, K is the dimensionless stiffness matrix, and the vector f is

associated with stochastic excitation. The excitation term W (τ ) = −¨xrM/M R ω2

0 is then

defined to be the normalized Gaussian white noise with zero mean and constant two-sided

power spectral density D, given that D = d/( ω3

0R2) where d is the physical noise intensity in

units of squared acceleration per frequency. ξ and ξp are the damping ratios associated with

the structure and the pendulum, respectively. α quantifies added mass due to inertance,

while µr and µp quantify pendulum mass and moment of inertia relative to structural mass.

η quantifies the length of the pendulum. Furthermore, g represents the nonlinear Coriolis,

centrifugal terms, and friction. It is finally apparent that µr and η are the quantifiers of

non-linearity in the system. Additionally, note that the equations of motion are written in

terms of the normalized time variable, τ = ω0t, where ω0 is the characteristic frequency

associated with the suspended mass and spring stiffness.

4.2 The Experimental Setup and Characterization

The experimental setup correlating with Fig. 4.1 is now shown in Fig. 4.2. The setup

consists of the IPVA, a generator, a mass suspended with eight coil springs and clamps,

a signal conditioner, an electrodynamic shaker, an amplifier, a vibration controller, and

a pendulum tracking system. In this section, the measured or experimentally determined

values for all parameters outlined in the previous section are given.

67

Figure 4.2 The Experimental System. (1) Pendulum tracking system. (2) Electrodynamic
shaker. (3) Amplifier, signal conditioner, and vibration controller. (4) DC generator. (5)
Pendulums on the carrier. (6) Ball screw and nut with the coupler to suspended mass. (7)
Camera and light for pendulum tracking. (8) Suspended mass.

4.2.1 Fitting the Frequency Response Function

Mechanical damping ξ, characteristic frequency ω0, and rotational inertia of ball screw

and carrier were fit in the frequency domain by fitting the norm squared frequency response

function (FRF) of the relative acceleration (¨xs − ¨xr) for the linear system with pendulums

and gears removed. For a constant acceleration base excitation experiment steady-state

response, ¨xr = A · exp j ωt, xs = Xs (j ω) · exp j ωt and xs − xr = RΘ (j ω) · exp j ωt. Then

the norm squared FRF takes the form given here.

(R|Θ|)2 ω4
A2

=

(cid:0)

r4

1 − r2 (1 + αl)2

(cid:1)2

+ 4ξ2r2

(4.11)

where r = ω
ω0

and αl = α − 2mgearR2

p/MsR2. A sine sweep test with 0.15g base acceleration

(g is the acceleration due to gravity) was done and the values of ω0, ξ, and α were fit using a

nonlinear least squares algorithm in MATLAB. The close agreement between the theoretical

68

Figure 4.3 Norm squared relative acceleration frequency response function experimental fit.

and experimental FRF is shown in Fig. 4.3.

4.2.2 The Experimental Methods

With the experimental setup shown in Fig. 4.2, a three-hour random excitation experi-

ment was done for the case of four different load resistances attached to the generator, giving

four different electric damping (ξe) values with the use of (4.4) and (4.8). This is outlined

in Table. 4.1. The resistor values were chosen to incorporate a sufficiently large electrical

damping range while using resistors that were imediately available. Note that open circuit

implies infinite load resistance. Also, note that a total of thirty minutes of data was re-

moved in order to eliminate transient beginning and ending dynamics which might have a

misleading effect when interpreting some statistical information. The total duration of the

analyzed data in each case was then about two and a half hours. The excitation used in the

69

Load Resistance (Rl) Electric Damping (ξe)
Open Circuit
10 Ω
3.9 Ω
1 Ω

0
.02
0.037
0.0637

Table 4.1 Electric load resistances and associated electrical damping value.

experiment was such that the base acceleration had a root mean square (RMS) value of 0.8g

where g is the acceleration due to gravity and the corresponding PSD is shown in Fig. 4.4.

Note that this PSD was constructed to be broadband on the frequency band between two

and thirteen cycles per second while also including the resonant frequency determined in the

FRF fit and also within the voltage, current, and stroke limits of the electrodynamic shaker.

Shaker current, voltage, and minimum frequency limitations are directly proportional to the

acceleration, velocity, and displacement output limitations, respectively.

For the purpose of verifying bifurcation predictions, the same experiments were also done

for the case of 0.6g RMS base acceleration as well as a 0.8g RMS base acceleration case with

a shorted circuit (zero load resistance). The corresponding PDFs will be analyzed and it shall

be determined whether the bifurcation boundaries can be used to predict their qualitative

nature.

4.2.3 Characterization of Experiment

Various fitted, measured, and generator parameters used in this study are presented in

Table 4.2 and Table 4.3. It should be noted that µr was directly measured, while µp and η

were computed using SolidWorks. Additionally, ce varies with load resistance as described by

(4.4), but once a load resistance is selected, it is calculated using the reported values for the

Maxon 110207 4.5 Watt DC motor with a gear ratio of ng = 5.2. Finally, although mechanical

damping ξ was determined through the FRF fit, it was observed that it significantly differs

from the value obtained in the actual experiment described in Sec. 4.2.2; hence, it is discarded

hereinafter. The value is included in Table E.1 of Appendix E. It should finally be noted that

stiffness ks and mass Ms in ω0 were chosen to keep the linear resonant frequency low and

70

Figure 4.4 PSD for 0.8g RMS base acceleration excitation experiment.

accommodate a lower excitation frequency upper bound while also keeping the mechanical

damping ratio ξm as low as possible with the springs and mass which were available on

hand. Pendulum parameters µr and η were chosen after numerical simulation predicted that

they were sufficiently large to accommodate the bimodal phenomenon of interest while also

remaining relatively compact. The choice of generator torque constant κt, internal resistance

Ri, and gear ratio ng was governed by a combination of what was commercially available

and the intent to limit generator inertia µg, limit internal resistance, and use a gear ratio

which would put the generator speed near its nominal rated value in the documentation.

71

Dimensionless Parameters Value

µr
η
µp
µg
α

2.43
0.45
0.38
0.28
5.93

Table 4.2 Dimensionless Parameters of Experimental IPVA system.

Physical Parameters Value

κt
Ri
ω0
Ms
R

0.0168 Nm/A
3.17 Ω
12.26 Hz
9.39 kg
.03
2π m

Table 4.3 Physical Parameters of Experimental IPVA system.

4.2.3.1 Parameter Fitting with Frequency Domain Optimization

To further accommodate complicated dynamics not modeled in (4.9) and (4.10) which

may only be noticeable in the experiment with random excitation, Ξ = (ξf , ξp, ξ)T has been

left as a vector of variables to be tuned in an optimization routine. The new approximate

model is then a function of unknown linear damping coefficients Ξ. Also, in this case, W (τ )

in (4.10) is the dimensionless form of the exact excitation input to the shaker which was also

used to generate Fig. 4.4. The optimization routine is now defined for the determination

of the most suitable Ξ that forces the RMS velocities for the ball screw, pendulum, and

generator obtained from the simulation of the EOM and from the experiment outlined in

Sec. 4.2.2 to be within three percent agreement. To accurately capture pendulum dynamics,

especially when there is rotation, the objective function is chosen to minimize the error in

the mean square pendulum velocity on the frequency band below two cycles per second. The

72

optimization problem is formally defined as

(cid:12)
(cid:12)
(cid:12)
(cid:12)1 −
(cid:12)

X

0≤fk<2

(cid:12)
(cid:12)2

(cid:12)
(cid:12)Vψexp(fk)
|Vψsim(fk)|2

(cid:12)
(cid:12)
(cid:12)
(cid:12)
(cid:12)

min
Ξ

subject to

max
i=1,2,3

|Fi| < 0.03

(0, 0, 0)T ≤ Ξ ≤ (0.01, 0.03, 0.15)T ,

M(ψ)x′′ + C(ξe, Ξ)x′ + Kx + g(x, x′) − f = 0,

(4.12)

where

x(0) = x0,exp, x′(0) = v0,exp,
(cid:18)

F =

1 −

RMS(θ′
RMS(θ′

exp)
sim)

, 1 −

RMS(ψ′
RMS(ψ′

exp)
sim)

, 1 −

RMS(ψ′
RMS(ψ′

g,exp)
g,sim)

(cid:19)

.

where RMS(·) is defined as the RMS value of (·). The symbols (·)exp and (·)sim represent

the experimental and simulated values of (·), respectively. Additionally, x0,exp and v0,exp

are the initial position and velocity in the experiment, respectively.

It is also important

to note that the time period used for both the simulation and experimental data in the
(cid:17)

(cid:16)

optimization process was the same. Furthermore, Vψ (fk) =

P

N −1
n=0

˙ψ (n) exp

−i2π fk
fs

n

denotes the discrete Fourier transform of the pendulum velocity, where n is the sample

number, fk = kfs/N is the chosen frequency in cycles per second, fs is the sampling rate

and N is the total number of samples. The objective function is then the absolute value of

the error in the mean spectral power of the pendulum velocity below two cycles per second.

The upper limit of two cycles per second was chosen to account for the pendulum dynamics

associated with frequencies below the lowest excitation frequency, including the rotational

component at zero cycles per second. This should help fit more closely potential rotational

dynamics as well as possible subharmonic features in the response.

To execute the optimization described in (4.12), an initial three-dimensional grid of finite

parameter values was established. Initial fitting in the time domain detailed in Appendix

E helped narrow down the parameter space to the range specified in (4.12). Next, parallel

computing leveraging high-performance computing resources with 4 × 128 CPU cores and

73

921.6 gigabytes of memory per 128 cores to solve the EOM in (4.12) with ode45 in MAT-

LAB. Consequentially, the objective function value was determined. The minimal objective

function that satisfied the constraint was then determined and the optimized damping pa-

rameters for all resistant loads along with the corresponding maximum RMS velocity percent

errors are presented in Table 4.4. It should be noted that each resistance load case was as-

sociated with a unique optimal Ξ. Throughout the rest of this study, Ξ is considered as a

function of ξe (dependent on load resistance). This approach is necessary because the model

described in (4.12) does not fully capture the complex dynamics involved. For instance,

intricate dynamics associated with ball screw damping, discussed in [93, 94], are omitted.

Following the optimization routine of (4.12), a comparison between experimental and

simulated PSDs is displayed in Fig. 4.5. Note that the maximum RMS velocity error is

below 3% across all PSDs and the low frequency contributions are accurate qualitatively,

while the peak PSD values do not align as closely. It was proven difficult to obtain complete

agreement. While this is deemed acceptable given the primary focus on the RMS error in

this research, it highlights some flaws in the mathematical model.

In concluding this section, it is noted that all RMS velocities and power, with power

being directly related to the square of the generator velocity as shown in previous chapters

and in [62], were within a 3% error margin. Consequently, it is anticipated that predictions

regarding vibration suppression (reduction in relative velocity of Ms) and energy harvesting

will have a maximum of 3% error, at least when electrical damping and noise intensity are

relatively close to the experimental values in magnitude. Additionally, it is important to

highlight that the optimization that initially only minimized RMS velocity error resulted

in RMS error as low as 0.62%. However, predictions related to the rotational dynamics

were poor. Hence, this work uses the objective function in (4.12) for enhanced accuracy in

predicting complicated dynamics and improving the likelihood of forecasting P-bifurcation

in the following section.

74

Load Resistance
(Rl)

ΞT = (ξf , ξp, ξ)

Open Circuit
10 Ω
3.9 Ω
1 Ω

(0.005, 0.004, 0.0821)
(0, 0, 0.1093)
(0, 0, 0.1093)
(0.0025, 0, 0.1093)

Maximum
RMS Velocity
Error
2.84%
2.10%
1.37%
2.74%

Table 4.4 Fitted damping values for all load resistance cases and corresponding maximum
% RMS velocity error.

(a) Rl = 1 Ω

(b) Rl = 3.9 Ω

(c) Rl = 10 Ω

(d) Open circuit.

Figure 4.5 Fitted simulation PSD (solid line). Experiment PSD (dashed line). Each figure
column corresponds to a different load resistor case.

4.3 Verification of P-bifurcation with Power Harvested and Vibration Suppres-

sion Exploration

The bifurcation tracking formulated in Chapter 2 will now be used to predict the critical

values in a two-parameter plane for which a P-bifurcation, or change between the total

number of PDF maxima, may occur. In this section, the stochastic bifurcation analysis is

briefly reintroduced and the prediction is compared with experimental results.

75

4.3.1 Moment Equations And P-bifurcation

Just as in Chapter 2, but choosing this time to use symbolic matrix notation rather than

summation notation [14], the equations governing the steady state solution to moments of

all orders are

(cid:28)

(∇gm)T ˆa +

ˆΣ ˆΣ

T ∇ (∇gm)T

(cid:17)(cid:29)

= 0

(cid:16)

1
2

tr

(4.13)

where ⟨·⟩ is the expectation operator, gm = ˆyi

variables in ˆy = (θ, ϕ, θ′, ϕ′)T up to a desired order n = i + j + k + l, ∇ = ∂

4 is a polynomial consisting of the state
∂ ˆy , and tr(·)
represents the trace of the matrix (·). Note that ϕ is the pendulum angle of oscillation

2 ˆyk

3 ˆyl

1 ˆyj

relative to some fixed angle ϕ0. ˆa and ˆΣ are the deterministic and stochastic contributions,

respectively. Specifically, they are defined by first writing the stochastic differential form

[15] of (4.12) as follows.

dy = a (y) dτ + Σ (y, τ ) dB

(4.14)

where y = (θ, ψ, θ′, ψ′)T and

2

6
4

a (y) = −

3

0

7
5 y −

B
@

1

C
A

02×1

M−1g (y)

02×2

I2×2

M−1K M−1C
√

0

1

3

2

6
4

Σ (y, τ ) =

02×2

M−1

7
5

B
@

D

C
A

0

(4.15)

(cid:10)

dB(τ )dBT(τ )

(cid:11)

= Ddτ .

and in this case dB is a scalar differential Wiener process, with

Just as in Chapter 2, the drift and diffusion vectors are now first-order Taylor expanded

about some pendulum angle ϕ0 to obtain

ˆa (ˆy, ϕ0) =

(cid:18)

a (y) +

(cid:20)

(cid:20)

ˆΣ (ˆy, ϕ0, τ ) =

Σ (y, τ ) +

(cid:19)

(cid:21)

ϕ

∂a (y)
∂ψ
(cid:18)

∂Σ (y, τ )
∂ψ

ψ=ϕ0
(cid:19)

(cid:21)

ϕ

ψ=ϕ0

(4.16)

To alleviate the infinite hierarchy of dependence of moments on higher order moments, as

detailed in Chapter 2 with reference to [14, 15], n is set equal to two, and Gaussian closure

76

is used. Next, the Jacobian of the moment equations in (4.13) with respect to moment

variables and ϕ0, is evaluated at ϕ0 = 0 and the parameter values required for a vanishing

Jacobian determinant are determined. Then with the algorithm developed in Chapter 2 all

P-bifurcation solution candidates are readily found.

4.3.1.1 P-bifurcation Verification with Experiment

With the constructed excitation PSD shown in Fig. 4.4, the experiment excitation is only

an approximation to the theoretically assumed Gaussian white noise. Given the theoretical

dimensionless noise D = d/( ω3

0R2), d was chosen so that the integral of the PSD of the

white noise on the frequency band from f1 = 2 Hz to f2 = 13 Hz was equal to the mean

square experiment acceleration value of (0.8g)2/2. The division by two is necessary to keep

total power conserved since one must integrate over the negative frequencies as well in the

case of a symmetric two-sided spectral density such as with the theoretical Gaussian white

noise. Finally, to obtain d, with Parseval’s theorem one finds

Z

f2

f1

d · df = (f2 − f1)d =

(0.8g)2
2

d =

(0.8g)2
2(f2 − f1)

(4.17)

The dimensionless value is then D = 0.2687.

Using the noise intensity, ξe, Ξ(ξe) from Table 4.4, and the bifurcation tracking algorithm

in Chapter 2, a bifurcation diagram was generated for each resistor case in Table 4.1 with

the RMS acceleration excitation of 0.8g. The maximum and minimum boundary for a

given electrical damping was then determined and the resultant boundary region is shown

in Fig. 4.6. In this figure are also markers indicating the experimental damping and noise

intensity values, including the values corresponding to the experiments only used to test the

boundary accuracy. It is assumed that the region between the curves in Fig. 4.6 serves as

an estimate for where p-bifurcation can occur. In other words, an electrical damping-noise

intensity pair outside of this region is assumed to be not correlated with a bifurcation point.

77

Figure 4.6 Lower and upper bifurcation estimation in the D − ξe plane. Point markers and
square markers indicate the 0.8g RMS excitation and 0.6g RMS excitation experiment values,
respectively. The right-most 0.8g marker corresponds to Rl = 0.

The marginal PDFs for the pendulum angle corresponding to each experimental param-

eter set were also obtained via a Monte Carlo simulation with ∼ 4 · 104 realizations, a

dimensionless simulation time of 104 and a time step of 10−3. Note that this simulation

time corresponds to only ∼ 130 seconds in physical time which is drastically smaller than

the experimental time of two and a half hours. As a result, any dynamics that develop on

a timescale larger than 130 seconds might not be accurately captured. However, given the

number of realizations needed to reach stationarity, running the MCS with a simulation time

of two and a half hours is impractical within a reasonable timeframe. With this disclaimer

aside, the resultant marginal PDFs are shown in Fig. 4.7 and the corresponding Marginal

PDFs from the experiment are shown in Fig. 4.8. Note that the pendulum will undergo

oscillations centered around ψ = 2πk radians for any integer k with possible intermittent

rotations in between. As such, the pendulum angle ψ has been replaced in the PDFs with

ψ mod 2π and further restricted to the domain of [−π, π] which is all accomplished using

78

Figure 4.7 Marginal probability density functions for the pendulum angle from simulation
using experiment parameters and 0.8g RMS excitation. Each PDF correlates with a point
marker in Fig. 4.6.

the ’wrapToPi’ function in MATLAB. This does have the side effect of the separation of the

PDFs from the ψ axis near ψ = ±π when rotation exists.

Initially, it is important to recognize that the bifurcation diagram in Fig. 4.6 delineates

both a right and left boundary region. Considering the relationship between damping and

resistance in Table 4.1, transitioning from a 3.9 Ω to a 1 Ω load resistor necessitates entering

the right bifurcation boundary region from the left and ending near the outer edge. With

that being said, the original bifurcation boundary generated with the parameters from the

1 Ω load resistor case actually had an outer edge that was to the left of the experimental

data point. In this case, it is assumed that bifurcation has in fact occurred. Notably, the

PDF in Fig. 4.7 associated with the 1 Ω case has negative curvature at its origin and is

distinctly monomodal, while the PDF for the 3.9 Ω case on the left side of the boundary

region shows positive curvature at the origin and is distinctly bimodal. The experimental

PDFs in Fig. 4.8 exhibit similar characteristics in terms of curvature at the origin, with more

79

Figure 4.8 Marginal probability density functions for the pendulum angle from the experi-
ments with 0.8g RMS excitation. Each PDF correlates with a marker in Fig. 4.6.

Figure 4.9 Marginal probability density functions for the pendulum angle from the experi-
ments with 0.6g RMS excitation. Each PDF correlates with a square marker in Fig. 4.6.

80

pronounced bimodality. Furthermore, the experimental short circuit case is shown in Fig. 4.8,

corresponding to the marker with the highest damping value in Fig. 4.6. The PDF is indeed

more distinctly monomodal, indicating that the transition has fully occurred from bimodal to

monomodal. The right bifurcation boundary has been predicted reasonably well. Moreover,

if the left boundary region is accurate, the 10 Ω case is likely to display a bimodal PDF as

well. Both simulation and experimental PDFs in Figs. 4.7 and 4.8 reveal this bimodality, also

with rotation indicated by the separation from the ψ axis at ψ = ±π, attributed to reduced

electrical damping. Crossing the left boundary region to examine the open circuit scenario,

one might anticipate another qualitative shift in the PDFs. In this instance, Figs. 4.7 and

4.8 suggest a possible flattening of the PDFs, while they show an apparent change from

positive to negative curvature at the origin, indicating bifurcation. The precise nature of

this change remains uncertain, yet a comparable qualitative change was noted in the Chapter

3, where the left boundary was presumed to lead to pure rotation. Lastly, Fig. 4.9 shows

the experimental PDFs for 0.6g RMS acceleration. The boundary would predict that all of

these should be bimodal, with the exception of the open circuit case which lies in the region

of uncertainty. After closer examination, the three bimodal predictions are in fact verified,

albeit more asymmetric due to frictional effects or not enough time to become stationary.

The PDF with a prediction in the bifurcation region of uncertainty then appears be possibly

locally bimodal with a somewhat flattening of the PDF. This indicates a potential transition

towards pure rotation.

Some concluding observations pertaining to this subsection are now provided. The PDFs

derived from the MCS were produced using numerous realizations in the time domain, in

contrast to the PDFs associated with the experiment, which were derived from a single two-

and-a-half-hour experiment. Consequently, it is highly probable that the PDFs depicted in

Fig. 4.7 are closer to being stationary compared to those in Fig. 4.8, as suggested by the

greater symmetry in Fig. 4.7. However, there remains a possibility that the duration of

the simulation was insufficient to achieve fully developed PDFs, even if they are relatively

81

symmetric compared to the experimental PDFs. Furthermore, it is important to highlight

that generating these PDFs took around four hours, leveraging high-performance computing

resources with 8 × 128 CPU cores and 921.6 gigabytes of memory per 128 cores. Addition-

ally, it should be noted that the Gaussian white noise excitation utilized in the analysis

is merely an approximation of the band-limited noise employed in the experiment. Lastly,

the precision of the pendulum tracking system was constrained by a sampling rate of five

hundred samples per second, yet its measurement accuracy was confirmed to be within a

couple of degrees. Nevertheless, the qualitative change in the PDF was predicted with the

use of the approximate bifurcation boundary region.

4.3.2 Experimental Exploration of Power Harvested and Vibration Suppression

In this section, the PSDs are presented for the relative velocity of the suspended mass and

harvested power computed with the data obtained from the experiments discussed in Section

4.2.3.1. Additionally, the same experiments were done for a linear benchmark system where

the pendulums are removed, but equivalent inertia is added so that the natural frequency

remains about the same for a fair comparison. However, note that the pendulums act as

vibration absorbers and so their removal results in very high displacement amplitudes of the

suspended mass. For safety and to preserve the integrity of the experiment, the linear system

was excited at only half of the root mean square acceleration excitation and the resultant

amplitude was scaled by two. This is simply utilizing the principles of linearity; i.e. a force

twice as large will give a response twice as large. Ultimately, after this scaling, the PSDs

from the linear system experiment will be compared to those from the experiment with the

addition of the IPVA.

First, Fig. 4.10 shows the PSD for the electric power harvested across the resistant load.

Note that the power in this work was considered as the power across the resistive load which

was connected in between the generator terminals. The power harvested was calculated

directly from measured voltage (Vm) as P = V 2

m/Rl. In this scenario, it becomes evident that

the 1 Ω resistance setup correlates with the least electric power near the resonant frequency,

82

Figure 4.10 PSD for harvested power for the three finite resistant load configurations.

Figure 4.11 PSD for the relative suspended mass velocity for four configurations.

83

while it still outperforms all linear system cases in terms of mean power, calculated as the

area under the PSD. A closer look at the PDF in Fig. 4.8 reveals that this 1 Ω resistance

setup is also a monomodal configuration. The correlation between a monomodal PDF and

suboptimal energy harvesting exists. For Rl = 10 Ω, the mean power is more than a factor

of four times that of the linear system, and the power near the resonant frequency is more

than double compared to the best linear system (Rl = 3.9 Ω). A significant portion of the

mean power can be attributed to the substantial low-frequency response below the excitation

frequency, including the DC component in the PSD due to the pendulums’ rotation. Note

that the actual amount of power is not discussed in this work as the concern is primarily

with how the power correlates with the characteristics of the PDF and how the nonlinear

system performance compares to that of the linear benchmark system.

Moreover, Fig. 4.11 shows the PSDs for the relative suspended mass velocity. In this

case, the nonlinear system results are only compared with the linear system which uses the

3.9 Ω resistor. This is because this is the configuration in which the linear system performs

the best in terms of power harvesting. In fact, one can derive an optimal load resistance for

the linear system. Noting that the linear system would be the same as the linear system

in Chapter 2, where the linear closed-form solution for the dimensionless mean square ball

screw velocity was determined and is given by

⟨θ′2⟩ =

D
4(ξ + ξe)

(4.18)

The formula for dimensionless mean power harvested ⟨ ˆP ⟩ then just as in Chapter 3 is

Finally, taking ∂⟨ ˆP ⟩
∂Rl

⟨ ˆP ⟩ =

Rl
Ri + Rl

ξe⟨θ′2⟩ =

Rl
Ri + Rl

ξeD
4(ξ + ξe)

= 0, one can find that the optimal load resistance

√

Ri

Ropt

l =

q

κ2
t n2
√

g + 2Msξ ω0Ri
2Msξ ω0

(4.19)

(4.20)

Using the parameters from Table 4.3, one finds when rounded to the hundredths place,

Ropt

l = 3.17 Ω, which is very close to the experimental value of 3.9 Ω.

84

It is clear in Fig. 4.11 that the nonlinear system once again outperforms in terms of

vibration suppression near the resonant frequency except when there is an open circuit

condition and electrical damping is zero. A reduction by a factor of four in the relative

velocity PSD value is seen near the resonant frequency and an overall halving in the mean

squared value of the relative velocity has been computed for the case of Rl = 1 Ω in the

system with the IPVA. As long as power harvesting is a priority and the circuit is not

open, it appears as though the IPVA will allow for superior vibration suppression and power

harvesting over the linear system. This possibility for simultaneous objectives somewhat

echoes the result of utilizing the device in a two-degree-of-freedom quarter car model such

as in [62] and Chapter 3.

A final note is that the nonlinear system with 10 Ω resistance actually had 6% higher

mean square velocity and hence worse vibration suppression. However, a linear system fully

optimized for power using (4.20), given Ropt

l = 3.17 Ω, will have a much larger mean squared

value for the relative velocity of Ms. Using (4.18) one can then find an expression for the

ratio of optimized mean square relative velocity R2⟨θ′2⟩opt to mean square relative velocity

with another load resistance R2⟨θ′2⟩. One finds

⟨θ′2⟩opt
⟨θ′2⟩ =

√

t n2
κ2

g + 2MsξR2ω0 (Ri + Rl)
Msξω0(κ2

r

t n2

2R (Ri + Rl)

g+2MsξR2ω0Ri)
Ri

(4.21)

and with Tables 4.2, 4.3, and 4.4, a 26% increase in the linear system mean squared relative

velocity has been determined which would put it 20% above that of the IPVA system with

10 Ω resistance. In this case, one can similarly use (4.19) and (4.20) to determine a corre-

sponding 8% increase in power for the optimized linear system. This still gives a power that

is less than 30% of the power harvested with the IPVA system and 10 Ω resistance.

In concluding this section, it is important to solidify the hypothesis for why energy har-

vesting appears to be optimal when the nonlinear system PDF is bimodal. More specifically,

the most superior power was found in the 10 Ω case. A closer look at Figs. 4.7 and 4.8 reveals

that this case was also associated with a smaller decrease in PDF magnitude in between the

85

peaks when compared to the bimodal PDF associated with the 3.9 Ω resistance case.

It

is hypothesized that this smaller decrease in PDF magnitude between peaks is responsi-

ble for more frequent large amplitude jumps between the pendulum angles specified by the

PDF extrema. The work in Chapter 3, which was associated with the same energy har-

vesting IPVA device implemented in a quarter car suspension system, also concluded that

energy harvesting performance was generally better when near bifurcation, where the PDF

magnitude reduction is minimal. Furthermore, this hypothesis would be analogous to how

a shallow well in double potential-based systems allows for enhanced interwell oscillations

which is known and exploited by many researchers such as in [6].

4.4 Discussion on Findings

In this chapter, an inertially nonlinear device was equipped with a DC generator to con-

duct experiments on and verify predictions regarding the modality of the marginal PDF

of the pendulum angle of oscillation as well as the device’s ability to suppress vibrations

and harvest energy when applied to a single degree of freedom structure under Gaussian

broadband base excitation, treated analytically as white noise. The experimental excita-

tion, closely approximating Gaussian white noise within experimental constraints, involved

a broadband excitation with a root mean square acceleration of 0.8g and a bandwidth of

eleven cycles per second, incorporating the system’s resonant frequency and a minimum

value of two cycles per second.

To facilitate experimental verification, the system’s unknown parameters were formally

characterized with frequency domain optimization, where the optimization routine was de-

signed to minimize the mean squared error in the pendulum velocity on the frequency band

below two cycles per second while constraining the RMS velocity discrepancy between the

simulations and actual experiments to be below 3%. All parameters whether measured,

calculated, or characterized via the optimization routine were documented. Experiments

utilized four different load resistance settings: open-circuit, 1 Ω, 3.9 Ω, and 10 Ω. Each

resistance setting was associated with a unique set of damping parameters determined from

86

optimization. The RMS velocities and the corresponding mean power (proportional to the

squared RMS velocity) were shown to be predicted within a margin of error of 3%. This

signifies predictive capabilities for both power generation and vibration reduction. Addition-

ally, a visual comparison between the simulation and the power spectral densities (PSDs)

for all velocities demonstrated sufficient agreement.

A numerically determined P-bifurcation boundary in the plane of noise intensity versus

electrical damping was derived for each experiment at 0.8g RMS acceleration excitation,

excluding a extra test scenario with a short circuit, and a bifurcation region was established

based on all four boundaries. Two clear qualitative changes were predicted and experimen-

tally verified when noise intensity remained constant and electrical damping varied. The

right boundary region indicated a transition from a bimodal to monomodal PDF, while the

left boundary region suggested a potential flattening of the PDF. The validity of the bifur-

cation boundary was enhanced with the addition of the short-circuited experiment at 0.8g

and four additional 0.6g experiments. There was significant qualitative agreement between

experimental and simulation PDFs.

The power spectral density for electric power in the IPVA system was compared to a

linear benchmark system. Results showed that power harvested by the IPVA was double

that of the linear benchmark system near the resonant frequency, and mean power was

a factor of four times that of the best linear system tested. Additionally, an analytically

optimized linear system was shown to still produce less than 30% of the power harvested by

the IPVA system. The IPVA configuration with the lowest performance in terms of power

harvesting had the lowest load resistance of 1 Ω, associated with a monomodal PDF.

Finally, the power spectral density for the relative velocity of the suspended mass was

assessed and compared to the most effective experimentally tested linear system for power

harvesting, which had a load resistance of 3.9 Ω. All configurations except the open circuit

configuration outperformed the linear system, by as much as a factor of four at the peak PSD

value. Then all except the open circuit and 10 Ω resistance configuration had as little as half

87

of the mean square relative velocity associated with the linear system. It was determined

that the 10 Ω resistance configuration had 6% higher mean square velocity than the tested

linear system with 3.9 Ω resistive load. However, further analysis showed that a linear system

that was fully optimized for power had 20% higher mean square velocity and still less than

30% of the power when compared to the 10 Ω resistance IPVA configuration.

88

CHAPTER 5

FUTURE AND ONGOING WORK: AN EXPERIMENTAL STUDY WITH
THE FULL-SCALE ENERGY HARVESTING SHOCK ABSORBER

The full-size prototype has been built and is shown in Fig. 5.1 and the corresponding model in

the quarter car is shown in Fig. 5.2. Measured and known parameters are given in Table 5.1.

The objective of this chapter is to introduce the prototype, including design constraints and

considerations, and fully characterize mechanical damping associated with the ball screw,

generator, and pendulum motion. The MTS 810 hydraulic testing machine, equipped with

the available software, lacks the capability to perform random testing. Therefore, sinusoidal

testing is used and damping is determined with a parameter fitting in the time domain. Force-

velocity and force-displacement curves, as are used ubiquitously in EHSA literature such as

in [86, 95], are used to qualitatively observe the effective damping ratio and the total energy

dissipated. The characterized model is then used in a quarter-car model simulation with a

road profile defined in the same manner as in Chapter 3 with a noise intensity correlated with

a class F road (rough road) according to ISO 8608 standards [9]. Finally, power harvested,

road handling, and ride comfort performance metrics are thoroughly investigated with the

fitted model simulation results.

Figure 5.1 Manufactured IPVA integrated EHSA.

89

Figure 5.2 Quarter car model with IPVA integrated EHSA.

Parameters Value

Meaning

Ms
Mus
kt
ks
Rp
m
R
κ
Jg
Rint
Mgenerator
Mtotal

Sprung Mass
Unsprung Mass
Tire stiffness
Suspension stiffness
Carrier radius
Mass of pendulums
Effective ball screw radius
Motor torque constant

324 kg
50 kg
300263 N/m
110095 N/m
.072 m
3.48 kg
0.06/(2π) m
0.245 Nm/A
1340(10−7) kg m2 Generator Inertia
1.41 Ω
4.4 kg
15.71 kg

Motor internal resistance
Generator Mass
Prototype Mass

Table 5.1 Physical Prototype Parameters of the IPVA EHSA.

5.1 Physical Realization of the Prototype, Important Constraints, and Neces-

sary Design Considerations

As demonstrated in Sec. 3.3, the P-bifurcation is critical to the suspension performance.

In Sec. 3.4 the WPI algorithm efficiently estimated ∂2p
∂ϕ2
f

and a simple bifurcation detection

90

algorithm was developed to check for change in sgn

(cid:16)

(cid:17)

∂2p
∂ϕ2
f

.

Ideally, the design of a full-

scale prototype EHSA is designed as close to the model in Chapter 3 as physically possible,

associated parameter values are fitted, and then the bifurcation detection algorithm and

MCS are used to predict bifurcation and performance metrics near bifurcation for a few

driving speeds and as a function of electrical damping on a class F-road. However. some

parameters defined for the shock absorber in Chapter 3 have been deemed either impractical

or not physically realizable. As one example, η = 0.97 is not physically realizable unless the

carrier that holds the pendulum is unrealistically wide and most mass is concentrated near

the end.

Moreover, to keep the width and height down to a size that can convincingly fit inside a

passenger car wheel well, a constraint on the physical design is imposed to keep the shock

diameter (ds) below 300 mm.

The choice of the ball screw was not a trivial one to make as well. As pointed out

in section 3.3, a larger lead value allows for a lower dimensionless noise intensity. This is

apparent given the characteristic radius R = Lead/2π and dimensionless intensity formula

d = 2πGrV
ω0R2 .

In numerical simulation, a larger lead has also been correlated with better

performance. However, it has been determined that the ball screw lead value of 100 mm

leads to a very massive ball screw and nut configuration, which significantly increases the

total mass of the shock absorber and decreases the specific power (power per unit mass) of

the device. Another important factor to consider is the strength of the ball screw. Care was

taken to choose a diameter of the ball screw which was large enough to not break in operation

with potentially imperfectly aligned torsion and axial loads. Overall, a compromise between

large lead, sufficient strength, and low mass has led to the choice of a 60 mm lead ball screw.

The decrease in lead value unfortunately makes the analysis from the Chapter 3 not

directly usable even in the qualitative sense since µg is thereby amplified by 1002/602 times.

A change in DC motor (used as generator) from the Maxon 353300 motor with a nominal

speed of 2960 RPM to the to Maxon 353301 motor with a nominal speed of 2480 RPM has

91

also been done. This is because the system at hand has shown mean generator velocities

( ˙ψg = ˙θ − ˙ϕ) values on the order of only hundreds of RPM and a nominal operating speed

which is closer to the generator input speed results in superior energy harvesting. The gear

ratio has also changed from ng = 3.8 to 3.7 since this is what was immediately available;

note that future work will likely involve a more fine tuned custom generator. Finally, these

changes affect the value of µg and it is rightfully noted that the system dynamics have a

non-trivial dependence on µg. This is apparent by looking at the dependence in the mass

matrix from Chapter 3. Specifically, we have

M11 = µg + µr

(cid:0)

η2 + 2η cos(ϕ) + 1

(cid:1)

+ 1, M12 = ηµr(η + cos(ϕ)) − µg, M22 = η2µr + µg(5.1)

For example, investigation of M12 alludes to the fact that if µg > max (ηµr(η + cos(ϕ))) there

are always negative off-diagonal inertia terms. Otherwise, for η < 1, as is the requirement

for no collisions of the pendulum, the off-diagonal inertia terms may be positive for certain

pendulum angles of oscillation. These off-diagonal terms couple the pendulum and carrier

motion. Therefore, the change in sign can play a significant role in the overall dynamics.

In order to keep the dynamics of the designed prototype similar to that of the system

analyzed in the Chapter 3, an optimization routine has been developed to keep µg small

relative to η2µr. The optimization takes into account the following constraints: ds = 254 mm,

η = r/Rp < 1, the mass of the pendulums m < 4 kg giving a mass ratio m/Ms ≤ 4/324 <

1.5%, and the height of the pendulums (h) is constrained by |h|< 50.8 mm in order to keep

the total length of the shock absorber down so as to fit it in a passenger vehicle and be

around the 70 cm determined from an existing shock absorber that has been measured in

the lab. To achieve the pendulum weight with minimal volume, tungsten with a density of

17.99g/cm3 is used.

92

5.2 Characterization of Energy Harvesting Shock Absorber

5.2.1 Energy method with the linear system

The equations of motion associated with the shock detached from the car for experimen-

tation are then derived using the Euler-Lagrange formulation. Specifically, they are derived

in terms of dimensionless parameters but with physical time t = τ /ω0, pendulum angle rel-

ative to the carrier ϕ (as in Chapter 3), and ball screw angular rotation about its vertical

axis θ.

(cid:18)

α + µp + µgear + µg + µr

MtranR2
MsR2
˙θ
+ (ηµr(η + cos(ϕ)) − µg + µp + µgear) ¨ϕ + 2(ξ + ξe + ξf )ω0

η2 + 2η cos(ϕ) + 1

+

(cid:0)

(cid:1)

(cid:19)

¨θ

− 2 (ξe + ξf ) ω0

˙ϕ − µrη(2 ˙ϕ ˙θ + ( ˙ϕ)2) sin(ϕ) =
(cid:0)

(ηµr(η + cos(ϕ)) − µg + µp + µgear) ¨θ +

Tnl (t)
MsR2 ,
η2µr + µg + µp + µgear

(cid:1)

¨ϕ

− 2ξeω0

˙θ + (2ξe + 2ξp + 2ξf ) ω0

˙ϕ + µrη( ˙θ)2 sin(ϕ) = 0

(5.2)

Mtran is quantified as the mass which translates with the MTS hydraulic machine actuator.

µgear quantifies the summation of all four planetary gear principal moment of inertia about

their rotational axis. α quantifies the moment of inertia about the ball screw vertical axis

for all components which rotate with the ball screw angular displacement, not including the

pendulums. µp then quantifies the principal moment of inertia of the pendulums about the

vertical axis through their center of mass. Tnl (t) /R is the imposed force on the nonlinear

shock along the ball screw rotational axis in order to prescribe the motion xs − xr = Rθ. All

other parameters have been defined previously in chapter 3.

For the locked pendulum, this time with force Tl (t) /R, we have
(cid:18)

(cid:19)

α + µp + µgear + µg + µr(η2 + 2η + 1) +

MadapterR2
MsR2

¨θ + 2(ξ + ξe)ω0

˙θ =

Tl (t)
MsR2

(5.3)

Prescribing ˙θ = −ωΘ cos ωt in the case of harmonic excitation and with µtran = Mtran
Ms

we

have

Tl (t) = ΘMsR2ω(2ω0 cos(ωt)(ξ + ξe + ξf ) −

(cid:0)

α + µp + µtran + µg + µr (η + 1)2

(cid:1)

ω sin(ωt))

93

We may also integrate (5.3) with respect to angular displacement θ over one cycle T = 2π/ω

to obtain
Z

h(cid:0)

T

µg + µr(η2 + 2η + 1) + µtran

(cid:1)

i

¨θ + 2(ξ + ξe)ω0

˙θ

dθ =

0

Then using

R

T
0

¨θdθ =

R

T
0

¨θ ˙θdt =

R

(cid:16)

(cid:17)

T
0

1
2

d
dt

˙θ2

Z

T

0

Tl(t)
MsR2 dθ

(5.4)

dt = ˙θ2|T

0 = 0 a relationship between work

input to the system and the energy dissipated by mechanical and electrical damping is

derived. Specifically, with

R

T
0

˙θ2dt = πωΘ2 we have the following relationship between

damping parameters and work input Win into the system

ξ + ξe

=

≈

1
2πωω0Θ2MsR2

1
2πωω0Θ2MsR2

Z

T

0
N −1X

i=1

Tl (t) ˙θdt

RFl(ti)(θi+1 − θi) =

1
2πωω0Θ2MsR2 Win

(5.5)

Theoretically speaking, experiments can be defined and used to determine a range of values
for ξ assuming known ξe = ce/(2ω0R2Ms) based on the the relationship ce = κ2n2
g
Ri+Rl

and

Table 5.1. The determined parameters can be used for further fitting of ξp and ξf with

the use of 5.2. Note the linear system characterization has not been done for the sake of

preserving limited time and funding; only the characterization of the nonlinear system has

been done. The equations have been derived merely for completeness and to guide the

proposed characterization in future work.

5.2.2 Experimental configuration, constraints, and challenges

Load Resistance, Rl Frequency, ω (Hz) Amplitude, RΘ (mm)
Open circuit
39Ω
22Ω
10Ω
4.6Ω
1Ω

10, 5, 2
10, 5, 2
10, 5, 2
10, 5, 2
10, 5, 2
10, 5, 2

2,6,8
2,6,8
2,6,8
2,6,8
2,6,8
2,6,8

Table 5.2 Experiment values for the linear system limited by the maximum current and
generator speed limitations as well as a maximum input force which remains mostly periodic.

94

Note that the values in Table 5.2 defining the experiments obey the generator speed

and maximum continuous current limitation set by the manufacturer. The experiments are

designed so that the constraints of (5.6) are satisfied; the constraint equation assumed is

based on the linear system velocity, but is used to at least determine a rough estimate for

the upper limit of input displacement. The constraints are derived by taking the circuit of

Fig. 5.3 and neglecting inductance (Li = 0) due to low frequency operation of the generator,

giving the circuit of Fig. 5.4.

il =

κngθω
Ri + Rl

< 3.53 A

˙ψg = ngθω < 5500 RPM

Li

(5.6)

Rl

Ri

il

+ κng
−

(cid:16)

(cid:17)

˙θ − ˙ϕ

Figure 5.3 Circuit diagram with inductance.

Ri

il

Rl

+ κng
−

(cid:16)

(cid:17)

˙θ − ˙ϕ

Figure 5.4 Circuit diagram neglecting inductance.

Other limitations incorporated are the MTS 810 hydraulic testing machine force limit

of 25 kN and maximum peak to peak stroke of around 80 mm. Additionally, a plethora

95

of experimental testing showed that the force imposed (or at least the sensor reading) was

in fact irregular for imposed amplitudes beyond those of Table 5.2. An example data set

which demonstrates the erratic force readings is shown in Fig. 5.5; this example data set

specifically corresponds to experiments with 2 Hz input displacement frequency. While the

velocity reading is regular across all five experiment data sets, the force reading breaks down

into irregularity for the cases of larger input velocity/displacement.

Figure 5.5 Force and velocity corresponding to some experiments with erratic force readings.

In fact, it was proven impossible to reproduce the irregular force using measured voltage

and velocities and an extensive grid of damping parameters with the mathematical model

defined next in the Sec. 5.2.3; this model incorporates a twist in the actuator of the machine

which is introduced in the end of this section. The mathematical model is simply not

capable of capturing some unknown dynamics either associated with the prototype or the

MTS machine itself. The values in Table 5.2 were chosen based on trial and error and visually

inspecting the force data for expected periodicity. The experiments which were kept had force

sensor readings which were at least close to periodic. Experiments with aperiodic/chaotic

96

force readings were stored, but not used for further analysis. It is the opinion of the author

of this work that irregularity of the force measurement could have been due to air in the

system or possible sensor issues, but will remain unclear until further investigated in future

work. Future work will attempt to identify any unmodeled dynamics in the machine itself or

the prototype. If the machine itself or the sensors are responsible for unexpected dynamics,

a different machine and/or sensors may be sought, allowing for higher excitation profiles.

Since the use of the MTS hydraulic testing equipment available incurred hourly charges, the

current research scope was not able to accommodate further investigation.

Figure 5.6 Open Circuit Linear System Displacemnt PSD.

Moreover, Fig. 5.6 shows a representative linear system PSD for the car body velocity

when under random road excitation. It was generated with the assumption of two percent

mechanical damping (ξ) and the parameters of Table 5.1.

It is clear that the prominent

harmonic contributions are around two and six Hertz. The harmonic fitting experiments

therefore will include these frequencies in order to better represent the system dynamics,

albeit with an excitation which is still not very close to the real-world excitation.

Equation (5.2) gives the idealistic view of the equations of motion associated with the

system connected to the hydraulic testing machine shown in Fig. 5.7. However, the testing

machine shown unfortunately is free to rotate at the bottom where the actuator connects to

the shock, further complicating the experiment dynamics. As a result of this, the rotation

97

Figure 5.7 View of experimental setup including the laser, power resistor and data acquisition
software.

also rotates the case which is connected to the generator housing. The mathematical model

associated with this respective system including the twist is derived next in the Sec. 5.2.3

5.2.3 The experimental equations of motion including the twist degree of free-

dom

Figure 5.8 Bottom (generator side) view CAD model of the prototype with angles given.

98

𝝍𝐠=𝐱𝐑	−𝝓𝝓𝑹𝒑𝒓*̂,̂𝒙𝑹Figure 5.9 CAD model side view for the prototype with voltage divider circuit connected to
DAQ.

The equations of motion now have a new angular variable γ introduced. The equations

of motion are derived as follows. The position vector for the pendulum center of mass rp is

rp =

(cid:16)

(cid:16)

x
R
(cid:16)

Rp cos
(cid:16)

+

Rp sin

(cid:17)

x
R

(cid:17)(cid:17)

(cid:16)

+ r cos
(cid:17)

ϕ +
(cid:16)

x
R

+ r sin

ϕ +

x
R

i
(cid:17)(cid:17)

j

(5.7)

where i and j are orthonormal unit vectors in the horizontal plane. This is illustrated in

Fig. 5.8 which shows a cut of the prototype viewed from the bottom and looking at the

generator. x = Rθ is the displacement coordinate for the shock relative to the fixed end

where R = (Ball Screw Lead) /(2π). Note that the pendulum angle ϕ is a measured relative

to the carrier arm which rotates with angle x/R. Due to the planetary gear assembly, the

generator rotor angular rotation with respect to the generator housing is then ψg = x/R − ϕ.

The global experimental view is then shown in Fig. 5.9 which shows the model for shock

99

connected to a voltage divider circuit for safe voltage measurement as well as a disk for

velocity measurement with a laser doppler vibrometer (LDV).

The four pendulums then have a combined mass m and each one has a moment of inertia

Jp about their vertical principal axis. The gears have mass mgear and moment of inertia Jgear

about their vertical principal axis. The case that rotates with the machine has a moment of

inertia of Jcase about the vertical axis going from the bottom connection to the top. Then

all components that translate with the case have mass which sum to Mtran. With this in

mind, the kinetic energy is written as

T =

(cid:16)

1
2

mvp · vp + (4Jp + 4Jgear)
(cid:18)

(cid:19)

˙x
R

− ˙ϕ

2

(cid:0)

+

+ (Jg + Jgear)

(cid:19)

2

(cid:18)

˙ϕ +

˙x
R

Jc + R2
(cid:17)

pmgear

(cid:18)

(cid:1)

(cid:19)

2

˙x
R

+Mtran ˙x2 + (Jcase + Jgen,h) ˙γ2

(5.8)

where vp = d

dt rp is the velocity vector for the pendulum center of mass. As a reminder,

as shown in Fig. 5.8, Rp is the distance from the center of the carrier to the pendulum

attachment point. Jg represents the moment of inertia for the generator rotor around its

axis of rotation, Jgen,h denotes the moment of inertia for the generator housing around its

axis of rotation, and Jc represents the carrier’s moment of inertia with respect to its axis of

rotation.

Since there is no mechanical potential energy in the system, the next step is to write

the virtual work. For simplicity in this work, only linear damping is considered; the virtual

work is specifically defined to be the virtual work done by damping forces in the direction

of generator motion with coefficient cg, in the direction of pendulum motion with coefficient

cp, and in the direction of vertical translation with losses due to the ball screw and nut

interface with damping coefficient cm. Note that this cm also includes the damping in the

experiment which is inherent in the hydraulic machine actuator; the damping determined

will be an overestimate. Unfortunately, determining the machine damping was not deemed

possible due to limited time and is left for future work. Note there is additionally the work

100

due to the Lorentz force in the direction opposing generator rotor motion with coefficient

ce. Lastly, there is an assumed damping which opposes the twisting motion of the case γ

with coefficient ct, and there is a virtual work contribution from the force F (t) input to the

system from the hydraulic testing machine. The virtual work is written as

δW = − (cg + ce)

˙x/R − ˙γ − ˙ϕ

· δ (x/R − γ − ϕ)

(cid:16)

(cid:17)

+ (F (t) − cm ˙x) · δx − cp

˙ϕ · δϕ − ct ˙γ · δγ

(5.9)

With the Euler-Lagrange equations, the equations of motion may be written as

0

1

(cid:18)

(cid:19)

=

−M−1 (C ˙x + g − F)

¨x

C
B
C
B
C
B
¨ϕ
C
B
@
A

¨γ

(5.10)

where x = (x, ϕ, γ)T and

M11 =

µp + µgear + α + µexp + η2µr + µg + µr + 2ηµr cos(ϕ(t))
R

M12 = µp + µgear + η2µr − µg + ηµr cos(ϕ(t))
(cid:1)

(cid:0)

M22 = R

η2µr + µg
Jcase + Jgen,housing
MsR

M33 =

M (ϕ) =

2

6
6
6
6
4

M11 M12

M12 M22

3

7
7
7
7
5

0

0

0

0 M33

2

6
6
6
6
4

C =

2ω0(ξe+ξf +ξm)
R
−2ω0 (ξe + ξf )

−2ω0 (ξe + ξf )

−2ω0 (ξe + ξf )

2Rω0 (ξe + ξf + ξp)

2Rω0 (ξe + ξf )

3

7
7
7
7
5

−2Rω0 (ξe + ξf )
0

2Rω0 (ξe + ξf )
1

2Rω0 (ξe + ξf + ξt)

0

1

(cid:16)

(cid:17)

g

ϕ, ˙ϕ, ˙x

=

B
B
B
B
@

C
C
C
C
A

, F (t) =

C
C
C
C
A

B
B
B
B
@

F (t)
MsR

0

0

− ηµr ˙ϕ sin(ϕ)(R ˙ϕ+2 ˙x)
R

ηµr ˙x2 sin(ϕ)
R

0

101

Unfortunately, the angle γ is not easily measured and the angle ϕ can be measured using

the camera tracking system, but it was not feasible to modify the experimental setup within

a reasonable timeline and budget; the camera may be installed in future work associated

with the prototype. With that being said, the voltage V induced in the DC generator is

related to the velocities with Faraday’s law as given here

(cid:16)

˙γ =

˙x/R − ˙ϕ

(cid:17)

− V
κng

(5.11)

where κ and ng are the torque constant and gear ratio, respectively, as given in Table 5.1.

The relationship of (5.11) substituted into (5.10) gives a set of differential equations for x and

ϕ as a function of the measured voltage (V) and input force (F). Sparing algebraic details,

it is written generically in terms of some function f as

0

1

(cid:16)

B
@

¨x

C
A = f
¨ϕ

(cid:17)

ϕ, ˙x, ˙ϕ, V, F

In dimensionless time τ = ω0t, it is written as

0

1

B
@

x′′

ϕ′′

C
A =

1
ω2
0

f (ϕ, ω0x′, ω0ϕ′, V, F )

(5.12)

(5.13)

5.2.3.1 Fitting the model

In order to finally identify the damping parameters ξ, ξf , and ξp, an optimization routine

whereby the collected data over sixty cycles (from τ = τi to τ = τf = τi + 60 × 2π ω0

ω ) of input

displacement x = RΘ sin ωt is compared to the values obtained in simulation for different

electrical damping ratio values ξe; an additional twenty cycles are then used to test the fitted

model. The optimization routine is defined as

0

B
@

x′′

ϕ′′

s.t.

||x′

min
β

measured
||x′

− x′
sim
||2

||2

+ w1ξ2 + w2ξ2
f

measured
β = [ξ, ξf , ξp, ϕ′
1

0, ϕ0, µtran]

C
A = 1
ω2
0

f (ϕ, ω0x′, ω0ϕ′, V (τ ) , F (τ ))

102

(5.14)

0

B
B
B
B
@

x′ (0)

ϕ (0)

ϕ′ (0)

1

C
C
C
C
A

=

0

B
B
B
B
@

x′
0,measured

ϕ0

1

C
C
C
C
A

ϕ′
0
τi ≤ τ < τf

0 ≤ ξ ≤ 1; 0 ≤ ξf ≤ 1; 0 ≤ ξp ≤ 2

−54 ≤ ϕ′
0

≤ 54

−π ≤ ϕ0 ≤ π

0.08 ≤ µtran ≤ 0.11

Parameters Value

µg
µp
µgear
µr
α
η

0.065
0.05
.015
0.61
0.7
0.47

Table 5.3 Dimensionless Prototype Parameters.

In this work, the optimization was done with the nonlinear least squares algorithm in

MATLAB; within the optimization routine the MATLAB ode45 RK4(5) integrator was used

to numerically solve the coupled nonlinear odes with a relative and absolute error tolerance

set to 10−12. The low tolerance was deemed necessary after trial and error; it gave the most

stable and consistent solution convergence.

There are few details worth mentioning regarding the optimization routine in (5.14).

First, all parameters except µtran and the damping parameters were measured with a scale

and measuring tools or computed with the help of SolidWorks CAD software; in the case

of the generator parameters, the documentation was used. See Table 5.3 for the parameter

values. µtran, the dimensionless quantification of translational mass, includes the mass of

components that are not easily measured such as the actuator of the hydraulic machine; it

was therefore left as a fitted parameter with an educated estimate of the bounds based on

103

the known translational masses. Next, the pendulum initial conditions are fitted as there

was no collection of pendulum rotation data. As such, care must be taken when choosing

the optimal parameter values from the gradient based approach in the case of multiple

local minima. This was somewhat remedied in this work by adding the weighted squared

damping ratio values w1ξ2 and w2ξ2

f to the objective function; the weights were selected

to keep the weighted squared values on a slightly lower order of magnitude than the term

involving the error in velocity so as to not to take priority in the optimization process.

This serves the purpose of forcing the optimization routine to favor solutions with lower

damping. Lastly, while pendulum damping was fitted, it was determined that the the fitted

values were unreasonably high; the low excitation led to very small pendulum motion, which

accommodates a high damping solution. Therefore, for simulation purposes in Section 5.3,

pendulum damping ξp is set to zero as the value determined in the small-scale experiment

in Chapter 4 was negligible or zero depending on load resistance.

5.2.4 The fitted model and simulation

After the fitting process governed by the optimization routine in (5.14) a total of eighty

cycles of data were used to test the associated fit; the total data set consisted of twenty
cycles added to the previous sixty cycles. ϵ = RMS(x′

was then recorded and the

measured)

sim
RMS(x′

−x′
measured)

fitted parameters were recorded as candidates only for the lowest two ϵ cases. These cases

are given in Table 5.4 and with corresponding fitted ball screw velocity ˙θ = ˙x/R in Fig. 5.10.

The measured value matches the fitted value reasonably well; error was found to be around

10% − 17% times the RMS measured velocity amplitude. At this point, it is also reiterated

that the ξ damping ratio values are likely overestimates, which include the hydraulic machine

damping; the values of ξ = 0.25 and ξ = 0.56 also serve as some sort of bounds on the actual

damping, albeit overestimated.

104

Resistance (Ω) Frequency (Hz)

ξ

1Ω
1Ω

2 Hz
4 Hz

0.25
0.56

ξf
0
0

ξp
1.2
2.0

ϵ

0.1730
0.1062

Table 5.4 Damping solutions and corresponding ϵ for the two lowest error cases.

(a) 1 Ω, 2 Hz

(b) 1 Ω, 4 Hz

Figure 5.10 Ball screw velocity; Fitted: dashed black line; Measured: solid red line.

Since it was determined that only two of the experimental data sets were fitted with

the current mathematical model reasonably well, the seemingly erratic force was further

explored. As a start, the measured velocity and force corresponding to the two and four

Hz cases are shown in Figs. 5.11 and 5.12, respectively.

The 1 Ω case clearly exhibits

the most periodic forcing; it was actually determined that it was not possible to reproduce

the force in the other cases via simulation with the current mathematical model, within

some reasonable error bounds. Additional unmodeled dynamics could be prominent. To

further demonstrate that the measured force F(t) is not reproducible, F(t) and ϕ are herein

simultaneously obtained via numerical simulation while using the optimized parameters in

β and experimentally determined quantities ˙x, ¨x, and V (t). Specifically, (5.10) and (5.11)

are combined and rearranged to obtain

0

1

F (t)

B
@

¨ϕ

C
A = −A−1 (ϕ) (B ˙y + g2 − F2)

(5.15)

105

Figure 5.11 Example experiment measured velocity and force - 2 Hz displacement frequency.

Figure 5.12 Example experiment measured velocity and force - 4 Hz displacement frequency.

106

with

2

6
4

B =

0

B
B
B
B
@

˙y =

˙x

˙ϕ

1

C
C
C
C
A

˙x
R
2

6
4

− ˙ϕ − V
κng

1
MsR

−M12

0 M22

3

7
5

A (ϕ) =

− 2ω0(ξe+ξf +ξm)

R

2ω0 (ξe + ξf )

2ω0 (ξe + ξf )

3

7
5

−2ω0 (ξe + ξf ) 2Rω0 (ξe + ξf + ξp) 2Rω0 (ξe + ξf )
1

0

1

0

(cid:16)

g2

ϕ, ˙ϕ, ˙x

(cid:17)

B
@

=

ηµr ˙ϕ sin(ϕ)(R ˙ϕ+2 ˙x)
R

ηµr ˙x2 sin(ϕ)
R

C
A , F2 =

B
@

M11 ¨x

C
A

−M21 ¨x

Finally, the resulting simulated vs. experimentally measured force for the two cases from

Table 5.4 as well as two additional cases which highlight a large discrepancy between the

measured and simulated force is shown in Fig. 5.13; the cases with large discrepancies corre-

spond to the 39 Ω cases with 2 Hz and 4 Hz displacement frequencies. The large discrepancy

between the experimentally measured force and the modeled force is evident in the 39 Ω

cases, while the 1 Ω cases exhibit reasonable agreement.

It is also important to quantify the twist rate due to the imperfect machine. As such,

the fitted pendulum angle ϕ and twist rate ˙γ along with the measured force and voltage

used as inputs in the EOM are shown in Fig. 5.14. It is evident that due to low excitation

conditions, the pendulum velocity relative to the carrier is relatively low, while the simulated

twist variable due to the unwanted twisting in the machine actuator was on the same order

of magnitude as the ball screw velocity shown in Fig. 5.10. This is undesirable and is not

representative of the real application in a vehicle suspension as the twist would surely be

zero in that case.

It is therefore difficult to make any remarks regarding performance of

the shock absorber; only the determined ξ and ξf values from Table 5.4 are assumed to be

somewhat accurate and are used in the next section for simulation purposes.

107

(a) 1 Ω, 2 Hz

(b) 39 Ω, 2 Hz

(c) 1 Ω, 2 Hz

(d) 39 Ω, 2 Hz

Figure 5.13 Examples of measured vs.
Experiment: solid red line.

simulated force. Simulation: dashed black line;

Lastly, shown in Fig. 5.15 are the associated force-velocity and force-displacement curves

for the nonlinear system. The fitted force-velocity curve is relatively accurate, but this is

expected since force was measured and the velocity from simulation was already shown to

match the velocity from experiment with ϵ < 0.18. More importantly, the force-velocity slope

can give an indication of the effective linear damping constant while the area enclosed in

the force-displacement curve gives an indication of total energy dissipated. The slope of the

force-displacement curve is also related to the effective stiffness (positive slope) and inertia

(negative slope). It is apparent, as would be expected, that the system is inertia dominant.

With that being said, the final remark is that the force-velocity and force-displacement curves

cannot truly give any meaningful effective damping/stiffness/inertia information regarding

108

(a) 1 Ω, 2 Hz

(b) 1 Ω, 4 Hz

Figure 5.14 Fitted pendulum angle ϕ and twist rate ˙γ. Measured force and voltage were
used as inputs in the EOM during fitting.

(a) 1 Ω, 2 Hz

(b) 1 Ω, 4 Hz

Figure 5.15 Force vs. velocity (top) and force vs. velocity (bottom)

the actual system of interest as the twist rate ˙γ was too large. The plots are therefore shown

mostly just for completeness. A thorough study involving the curves will be left for future

work.

5.3 Simulation of the Quarter Car Model with Identified Damping

With z = [θ, ϕ, ψus]T , the dimensionless EOM for the shock in the quarter car model

(assuming no twist) shown in Fig. 5.2 are essentially the same as in Chapter 3 with the

addition of the α, µgear, and µp terms since these are not negligible inertia terms. Pendulum

damping is also assumed to be negligible and is not included as mentioned in the previous

109

section. The EOM are:

Mz′′ + Cz′ + Kz = F (z, z′, Ψr)

(5.16)

where

M =

F (z, z′, Ψr) =

and

2

6
6
6
6
4

0

B
B
B
B
@

3

7
7
7
7
5

, C =

M11 M12

M12 M22

1

0

1

0

u + 1

µrη(2ϕ′θ′ + (ϕ′)2) sin (ϕ)

−µrη (θ′)2 sin (ϕ)

uv2Ψr

2

6
6
6
6
4

1

C
C
C
C
A

2 (ξ + ξe) −2ξe 0

−2ξe

0

2ξe

0

0

0

3

7
7
7
7
5

, K =

2

6
6
6
6
4

1 0

0 0

0

0

0 0 uv2

3

7
7
7
7
5

,

(5.17)

M11 = α + µgear + µp + µg + µr

(cid:0)

η2 + 2η cos(ϕ) + 1

(cid:1)

+ 1

M12 = µgear + µp + ηµr(η + cos(ϕ)) − µg

M22 = µgear + µp + η2µr + µg

Ψr corresponds to the dimensionless road profile acting as the excitation. The EOM for the

road excitation is given here

Ψ′

r + vcΨr = n(τ )

(5.18)

where the dimensionless white noise is defined thoroughly in Chapter 3. Note that in this

work, cutoff frequency vc is actually set to zero as it was shown to speed up the simula-

tion. The role of the cutoff frequency is to keep the zero-frequency component for the road

displacement bounded; the RMS road displacement will otherwise grow unbounded. How-

ever, the author, Cosner, has determined via direct simulation that the cutoff frequency has

negligible effects on the overall dynamics of the system and so may be safely set to zero.

Finally, the equations of motion were solved using the parameters of Table 5.3 and Table

5.4 with ξp = 0. u = Mus/Ms and v =

were calculated with parameters in Table

ks/Ms
5.1. For a driving speed of 5 m/s on an extremely rough road (off-road), the dimensionless

√
√

kt/Mus

110

(a) ξ = 0.25

(b) ξ = 0.56.

Figure 5.16 Power into electrical domain versus electrical damping ξe.

noise intensity is also taken as d = 1.94. This would correspond to a reasonable driving

speed on off-road style terrain. The numerical simulation was done with a Euler-Maruyama

based stochastic differential equation (SDE) solver found in the MATLAB Financial Toolbox

with a dimensionless time step set to 10−3 and total dimensionless time of 1000 units. The

MCS consisted of 35(103) simulations for a range of ξe values on the interval [0, 0.5]; the

highest value that the generator can accommodate is only very slightly higher than ξe = 0.5.

It should lastly be noted that the results to be introduced in this section should give an

indication of the possible bounds of performance, but they cannot be used to draw any

precise conclusions since the ξ = 0.25 fitted model and the ξ = 0.56 model may lead to very

different system dynamics.

5.3.1 Quarter car model simulation results

Power into the electrical domain (i.e. with no losses due to internal resistance) is shown

for the ξ = 0.25 case in Fig. 5.16a and for the ξ = 0.56 case in Fig. 5.16b. The peak power

occurs near ξe = 0.2 in both cases at around 250 W for ξ = 0.25 and slightly under 200 W

for ξ = 0.56.

Furthermore, the optimal linear system power may be calculated and remains the same as

(3.33) in Chapter 3. With electrical efficiency considerations (resistive losses) in the nonlinear

system, nonlinear system harvested power is normalized with respect to the optimized linear

111

(a) ξ = 0.25; Optimal linear power is 136 W.

(b) ξ = 0.56. Optimal linear power is 81 W.

Figure 5.17 Power harvested normalized with the optimized linear system versus electrical
damping ξe.

system and is shown for ξ = 0.25 in Fig. 5.17a and for ξ = 0.56 in Fig. 5.17b; a value greater

than unity means that the power exceeds that of the optimized linear benchmark system.

In fact, peak power harvested seems to exceed the linear benchmark system by around 40%

for the ξ = 0.25 case and around 60% for the ξ = 0.56 case. This amounts to a power

of around 190 Watts and 130 Watts respectively; the optimal damping has also shifted to

around ξe = 0.1. It is apparent that the choice of the generator led to subpar efficiency.

A better implementation in the future will most likely involve a customized generator that

may cater to high efficiency needs.

The sprung mass acceleration normalized with respect to the linear benchmark system

optimized for power is shown for the ξ = 0.25 case in Fig. 5.18a and for the ξ = 0.56 case in

Fig. 5.18b. Notably, peak performance occurs again near ξe = 0.1 with around a 20% sprung

mass acceleration reduction in comparison with the linear benchmark system.

The RHI, the ratio of the RMS dynamic load at the road to the vehicle weight (lower

is better for road handling), is normalized with respect to the linear benchmark system

optimized for power is shown for the ξ = 0.25 case in Fig. 5.19a and for the ξ = 0.56 case

in Fig. 5.19b. Interestingly, the normalized RHI is very similar (to the naked eye) to the

normalized sprung mass acceleration.

To gain more insight into the performance in terms of the relationship with the pendulum

112

(a) ξ = 0.25

(b) ξ = 0.56

Figure 5.18 Sprung mass acceleration normalized with the optimized linear system versus
electrical damping ξe.

(a) ξ = 0.25

(b) ξ = 0.56

Figure 5.19 RHI normalized with the optimized linear system versus electrical damping ξe.

dynamics, the marginal PDF (all other states are integrated out) for the pendulum angle

ϕ is shown for a few ξe values in Fig. 5.20. It should be noted that the wrapToPi function

in MATLAB 2023b was used to keep the pendulum angle on the domain of [−π, π]; this

leads to the notable separation at ϕ = ±π when large rotations exist. Also, for each ξ case,

the marginal PDFs are separated into two plots in order to more clearly observe qualitative

changes with varying ξe. The first thing to look at is the qualitative nature of the PDFs

with ξe ≤ 0.24 in Figs. 5.20a and Figs. 5.20b. Once again, the optimal performance in terms

of power, ride comfort, and road handling was found to occur at around ξe = 0.1. This

113

(a) ξ = 0.25. ξe ≤ 0.24.

(b) ξ = 0.56. ξe ≤ 0.24.

(c) ξ = 0.25. ξe ≥ 0.28.

(d) ξ = 0.56. ξe ≥ 0.28.

Figure 5.20 Marginal PDF for the pendulum angle ϕ for various ξe.

corresponds to between the triangle and diamond marked PDFs in the figures.

In the case of ξ = 0.25, the probability density functions around ξe = 0.1 are clearly

bimodal. Additionally, the probability densities at the origin seems to remain the same,

while the curvature at origin decreases with the ξe increase; this is accompanied by a slightly

broadened set of peaks. This potentially indicates improvements are based on enhanced

large jumps in the pendulum angle of oscillation.

In the case of ξ = 0.56, the probability density functions around ξe = 0.1 are perhaps

bimodal or even have more peaks than two (multimodal). In this case, it is a little less clear.

It is possible that the dynamics have become significantly more complicated at high ξ, or

there were not enough simulations for convergence to stationary PDFs. Given the fact that

114

the total simulation was well over 2.5 hours using over 7 × 900GB of memory and 7 × 128

cpu cores on the Michigan State University ICER HPCC, more simulations would likely not

be justifiable. With that aside, the most objectively determined change is that of increased

probability density in the region between inner peaks for an increase in electrical damping.

Finally, in the case of ξ = 0.25 and ξ = 0.56 with ξe ≥ 0.28, the main notable qualitative

change is that of increasing probability density at the origin, indicating a potentially more

monomodal (single peak) system.

5.4 Discussion on Findings

In concluding this chapter, it is first noted that the simulation results seemed to agree

qualitatively with the overlying theme associated with the conclusions in Chapters 2, 3,

and 4. That is, performance in terms of energy transfer from the structure in question to

the pendulums and/or electrical domain seems to be enhanced when the system is near

P-bifurcation or bimodal with small variation between peaks. Road handling specifically,

quantified by the RMS tire dynamic force to gravity force ratio (lower is better) was shown

to also be enhanced simultaneously in the quarter car model. This bimodality/P-bifurcation

that is associated with multi-objective performance was specifically made possible with the

introduction of inertial nonlinearity; in fact, the energy harvesting absorber device is free of

mechanical potential energy.

Next, it is noted that future work is needed for any absolute conclusions with high

objectivity regarding the quarter car model simulation results in this chapter. There are a

few reasons for this. One is that the inability of the hydraulic testing machine to impose

vertical displacement motion with high accuracy without inducing the added twist motion

makes the results questionable. The twisting motion has a notable effect on the voltage from

the generator; the voltage is in turn coupled with the rest of the system. A better/more

stable machine which is more suitable for cyclic or random displacement testing will ensure

higher confidence in experimental results. Another reason is that the linear damping ratios

fitted to the mathematical model may have been greatly overestimated due to the hydraulic

115

machine damping. In the future, it may be necessary to characterize the damping from the

machine independently. Also, the excitation was rather low; the dynamics might have been

friction dominated and no explicit friction models were used in the mathematical model.

Furthermore, either a random excitation experiment with an appropriate broad frequency

bandwidth or a set of experiments with a more comprehensive set of frequencies which

cover a broad bandwidth may lead to different fitted mathematical models. As an example,

the work in Chapter 4 included the fitting of a FRF model using data from a sine sweep

experiment. In this case, the damping ratio value that was fitted to the FRF model, albeit

an almost perfect fit, was very different than fitted in the optimization routine involving

the random excitation over a long runtime; as a reminder, the latter was based on fitting

mean-square values between specified frequency bands. It is the opinion of the author, that

this inconsistency in fitted damping ratio values may be due to the fact that complicated

nonlinearities such as ball screw damping were not present in the assumed modeled. As

such, the best fit to the approximate mathematical model may be more accurate when

the excitation is close to that of real operation or when it is fit with the use of mean

square quantities such as was done in Chapter 4. In future work, it might serve useful to

add nonlinear dynamics associated with the ball screw to the mathematical model and/or

implement random excitation and characterize the system in a manner similar to how it was

done in Chapter 4. Lastly, the ability to perform an experiment with random excitation

which is consistent with the simulation will add incredible value to the simulation results in

the future.

116

BIBLIOGRAPHY

[1] Smith, M. C., 2020. “The inerter: A retrospective”. Annual Review of Control, Robotics,

and Autonomous Systems, 3, pp. 361–391.

[2] Masana, R., and Daqaq, M. F., 2013. “Response of duffing-type harvesters to band-
limited noise”. Journal of Sound and Vibration, 332(25), Dec., pp. 6755–6767.

[3] Zhou, Z., Qin, W., Du, W., Zhu, P., and Liu, Q., 2019. “Improving energy harvesting
from random excitation by nonlinear flexible bi-stable energy harvester with a variable
potential energy function”. Mechanical Systems and Signal Processing, 115, pp. 162–
172.

[4] Li, X., Zhang, J., Li, R., Dai, L., Wang, W., and Yang, K., 2021. “Dynamic responses of
a two-degree-of-freedom bistable electromagnetic energy harvester under filtered band-
limited stochastic excitation”. Journal of Sound and Vibration, 511, Oct., p. 116334.

[5] De Paula, A. S., Inman, D. J., and Savi, M. A., 2015. “Energy harvesting in a nonlinear
piezomagnetoelastic beam subjected to random excitation”. Mechanical Systems and
Signal Processing, 54, pp. 405–416.

[6] Lan, C., and Qin, W., 2017. “Enhancing ability of harvesting energy from random
vibration by decreasing the potential barrier of bistable harvester”. Mechanical Systems
and Signal Processing, 85, Feb., pp. 71–81.

[7] Han, Q., Xu, W., and Sun, J.-Q., 2016. “Stochastic response and bifurcation of period-
ically driven nonlinear oscillators by the generalized cell mapping method”. Physica A:
Statistical Mechanics and Its Applications, 458, pp. 115–125.

[8] Wei, W., Xu, W., and Liu, J., 2021. “Stochastic p-bifurcation analysis of a class of
nonlinear markov jump systems under combined harmonic and random excitations”.
Physica A: Statistical Mechanics and its Applications, 582, p. 126246.

[9]

International Organization for Standardization, 1995. ISO 8608:1995 - Mechanical Vi-
bration - Road Surface Profiles - Reporting of Measured Data. Standard.

[10] Petromichelakis, I., Psaros, A. F., and Kougioumtzoglou, I. A., 2018. “Stochastic re-
sponse determination and optimization of a class of nonlinear electromechanical energy
harvesters: A Wiener path integral approach”. Probabilistic Engineering Mechanics,
53, June, pp. 116–125.

[11] Arnold, L., 2001. “Recent progress in stochastic bifurcation theory”. In IUTAM Sym-
posium on Nonlinearity and Stochastic Structural Dynamics, S. Narayanan and R. N.
Iyengar, eds., Springer Netherlands, pp. 15–27.

[12] Crauel, H. a., 1999. Stochastic Dynamics. New York, NY : Springer New York.

[13] Liang, Y., and Namachchivaya, N. S., 1999. P-Bifurcations in the Noisy Duffing-van

der Pol Equation. Springer New York, New York, NY, pp. 49–70.

117

[14] Roberts, J., and Spanos, P., 2003. Random Vibration and Statistical Linearization.

Dover Civil and Mechanical Engineering Series. Dover Publications.

[15] Soong, T., and Grigoriu, M., 1993. Random Vibration of Mechanical and Structural

Systems. PTR Prentice Hall.

[16] Knobloch, E., and Wiesenfeld, K. A., 1983. “Bifurcations in fluctuating systems: The

center-manifold approach”. Journal of Statistical Physics, 33(3), Dec., pp. 611–637.

[17] Sri Namachchivaya, N., 1990. “Stochastic bifurcation”. Applied Mathematics and Com-

putation, 38(2), pp. 101–159.

[18] Kougioumtzoglou, I. A., and Spanos, P. D., 2014. “Nonstationary stochastic response
determination of nonlinear systems: A wiener path integral formalism”. Journal of
Engineering Mechanics, 140(9), p. 04014064.

[19] Kougioumtzoglou, I., and Spanos, P., 2012. “An analytical Wiener path integral tech-
nique for non-stationary response determination of nonlinear oscillators”. Probabilistic
Engineering Mechanics, 28, Apr., pp. 125–131.

[20] Kougioumtzoglou, I. A., and Spanos, P. D., 2009. “An approximate approach for non-
linear system response determination under evolutionary stochastic excitation”. Current
Science, 97(8), pp. 1203–1211.

[21] Ikago, K., Saito, K., and Inoue, N., 2012. “Seismic control of single-degree-of-freedom
structure using tuned viscous mass damper”. Earthquake Engineering & Structural
Dynamics, 41(3), pp. 453–474.

[22] Wang, F.-C., Liao, M.-K., Liao, B.-H., Su, W.-J., and Chan, H.-A., 2009. “The perfor-
mance improvements of train suspension systems with mechanical networks employing
inerters”. Vehicle System Dynamics, 47(7), pp. 805–830.

[23] Marian, L., and Giaralis, A., 2014. “Optimal design of a novel tuned mass-damper–
inerter (tmdi) passive vibration control configuration for stochastically support-excited
structural systems”. Probabilistic Engineering Mechanics, 38, pp. 156–164.

[24] Hu, Y., and Chen, M. Z., 2015. “Performance evaluation for inerter-based dynamic
vibration absorbers”. International Journal of Mechanical Sciences, 99, pp. 297–307.

[25] Joubaneh, E. F., and Barry, O. R., 2019. “On the improvement of vibration mitiga-
tion and energy harvesting using electromagnetic vibration absorber-inerter: Exact h2
optimization”. Journal of Vibration and Acoustics, 141(6).

[26] Tai, W.-C., 2020. “Optimum design of a new tuned inerter-torsional-mass-damper pas-
sive vibration control for stochastically motion-excited structures”. Journal of Vibration
and Acoustics, 142(1).

[27] Haxton, R. S., and Barr, A. D. S., 1972. “The autoparametric vibration absorber”.

Journal of Engineering for Industry, 94(1), 02, pp. 119–125.

118

[28] Hatwal, H., Mallik, A., and Ghosh, A., 1982. “Non-linear vibrations of a harmonically
excited autoparametric system”. Journal of Sound and Vibration, 81(2), pp. 153–164.

[29] Song, Y., Sato, H., Iwata, Y., and Komatsuzaki, T., 2003. “The response of a dynamic
vibration absorber system with a parametrically excited pendulum”. Journal of Sound
and Vibration, 259(4), pp. 747–759.

[30] Warminski, J., and Kecik, K., 2009. “Instabilities in the main parametric resonance area
of a mechanical system with a pendulum”. Journal of Sound and Vibration, 322(3),
pp. 612–628.

[31] Ibrahim, R. A., and Roberts, J. W., 1976. “Broad band random excitation of a two-
degree-of-freedom system with autoparametric coupling”. Journal of Sound and Vibra-
tion, 44(3), Feb., pp. 335–348.

[32] Ibrahim, R. A., and Heo, H., 1986. “Autoparametric vibration of coupled beams under
random support motion”. Journal of Vibration, Acoustics, Stress, and Reliability in
Design, 108(4), 10, pp. 421–426.

[33] Yan, Z., and Hajj, M. R., 2015. “Energy harvesting from an autoparametric vibration

absorber”. Smart Materials and Structures, 24(11), p. 115012.

[34] Yan, Z., and Hajj, M. R., 2017. “Nonlinear performances of an autoparametric vibration-
based piezoelastic energy harvester”. Journal of Intelligent Material Systems and Struc-
tures, 28(2), pp. 254–271.

[35] Kecik, K., 2018. “Assessment of energy harvesting and vibration mitigation of a pendu-
lum dynamic absorber”. Mechanical Systems and Signal Processing, 106, pp. 198–209.

[36] Felix, J. L. P., Balthazar, J. M., Rocha, R. T., Tusset, A. M., and Janzen, F. C., 2018.
“On vibration mitigation and energy harvesting of a non-ideal system with autopara-
metric vibration absorber system”. Meccanica, 53(13), pp. 3177–3188.

[37] Tan, T., Yan, Z., Zou, Y., and Zhang, W., 2019. “Optimal dual-functional design
for a piezoelectric autoparametric vibration absorber”. Mechanical Systems and Signal
Processing, 123, pp. 513–532.

[38] IEA, 2024. Quarterly electric car sales by region, 2021-2024. https://www.iea.org/data-

and-statistics/charts/quarterly-electric-car-sales-by-region-2021-2024.

[39] IEA, 2024. Global ev outlook. https://www.iea.org/reports/global-ev-outlook-2024.

[40] International Council on Clean Transportation, 2024. Zero-emission vehicle phase-ins:
Passenger cars and vans/light trucks (July 2024). https://theicct.org/zero-emission-
vehicle-phase-ins-passenger-cars-and-vans-light-trucks-july-2024/.

[41] Goldner, R., Zerigian, P., and Hull, J., 2001. “A preliminary study of energy recovery in
vehicles by using regenerative magnetic shock absorbers”. SAE Technical Paper(2001-
01-2071).

119

[42] Karnopp, D., 1989. “Permanent Magnet Linear Motors Used as Variable Mechanical
Dampers for Vehicle Suspensions”. Vehicle System Dynamics, 18(4), Jan., pp. 187–200.

[43] Abdelkareem, M. A., Xu, L., Ali, M. K. A., Elagouz, A., Mi, J., Guo, S., Liu, Y.,
and Zuo, L., 2018. “Vibration energy harvesting in automotive suspension system: A
detailed review”. Applied Energy, 229, Nov., pp. 672–699.

[44] Daqaq, M. F., Masana, R., Erturk, A., and Dane Quinn, D., 2014. “On the Role
of Nonlinearities in Vibratory Energy Harvesting: A Critical Review and Discussion”.
Applied Mechanics Reviews, 66(4), July, p. 040801.

[45] Harne, R. L., and Wang, K. W., 2013. “A review of the recent research on vibration
energy harvesting via bistable systems”. Smart Materials and Structures, 22(2), Feb.,
p. 023001.

[46] Chen, C., Chan, Y. S., Zou, L., and Liao, W.-H., 2018. “Self-powered magnetorheolog-
ical dampers for motorcycle suspensions”. Proceedings of the Institution of Mechanical
Engineers, Part D: Journal of Automobile Engineering, 232(7), June, pp. 921–935.

[47] Singh, S., and Satpute, N. V., 2015. “Design and analysis of energy-harvesting shock
absorber with electromagnetic and fluid damping”. Journal of Mechanical Science and
Technology, 29(4), Apr., pp. 1591–1605.

[48] Zuo, L., Scully, B., Shestani, J., and Zhou, Y., 2010. “Design and characterization
of an electromagnetic energy harvester for vehicle suspensions”. Smart Materials and
Structures, 19(4), Apr., p. 045003.

[49] Jin-qiu, Z., Zhi-zhao, P., Lei, Z., and Yu, Z., 2013. “A review on energy-regenerative
suspension systems for vehicles”. In Proceedings of the world congress on engineering,
Vol. 3, pp. 3–5.

[50] Li, P., and Zuo, L., 2017. “Influences of the electromagnetic regenerative dampers
on the vehicle suspension performance”. Proceedings of the Institution of Mechanical
Engineers, Part D: Journal of Automobile Engineering, 231(3), Feb., pp. 383–394.

[51] Abdelkareem, M. A., Xu, L., Guo, X., Ali, M. K. A., Elagouz, A., Hassan, M. A., Essa,
F., and Zou, J., 2018. “Energy harvesting sensitivity analysis and assessment of the
potential power and full car dynamics for different road modes”. Mechanical Systems
and Signal Processing, 110, Sept., pp. 307–332.

[52] Huang, B., Hsieh, C.-Y., Golnaraghi, F., and Moallem, M., 2015. “Development and
optimization of an energy-regenerative suspension system under stochastic road excita-
tion”. Journal of Sound and Vibration, 357, pp. 16–34.

[53] Casavola, A., Di Iorio, F., and Tedesco, F., 2018. “A multiobjective H∞ control strategy
for energy harvesting in regenerative vehicle suspension systems”. International Journal
of Control, 91(4), Apr., pp. 741–754.

120

[54] Andren, P., 2006. “Power spectral density approximations of longitudinal road profiles”.

International Journal of Vehicle Design, 40(1-3), pp. 2–14.

[55] Li, M., Liu, D., and Li, J., 2023. “Stochastic analysis of vibro-impact bistable energy
harvester system under colored noise”. Nonlinear Dynamics, 111(18), July, pp. 17007–
17020.

[56] He, Q., and Daqaq, M. F., 2015. “New Insights Into Utilizing Bistability for Energy
Harvesting Under White Noise”. Journal of Vibration and Acoustics, 137(2), Apr.,
p. 021009.

[57] Mirzakhalili, E., and Epureanu, B. I., 2019. “Probabilistic Analysis of Bifurcations in
Stochastic Nonlinear Dynamical Systems”. Journal of Computational and Nonlinear
Dynamics, 14(8), 06, p. 081009.

[58] Cosner, J. A., and Tai, W.-C., 2023. “P-Bifurcation Analysis of a Quarter-Car Model
With Inerter-Based Pendulum Vibration Absorber: A Wiener Path Integration Ap-
proach”. Journal of Computational and Nonlinear Dynamics, 19(2), 12, p. 021004.

[59] Daqaq, M. F., 2012. “On intentional introduction of stiffness nonlinearities for energy
69(3), Aug.,

harvesting under white Gaussian excitations”. Nonlinear Dynamics,
pp. 1063–1079.

[60] Cosner, J. A., and Tai, W.-C., 2022. “Stochastic bifurcation and energy transfer in an
inerter-based pendulum vibration absorber”. Journal of Computational and Nonlinear
Dynamics, 17(8), p. 081003.

[61] Psaros, A. F., Kougioumtzoglou, I. A., and Petromichelakis, I., 2018. “Sparse repre-
sentations and compressive sampling for enhancing the computational efficiency of the
Wiener path integral technique”. Mechanical Systems and Signal Processing, 111,
Oct., pp. 87–101.

[62] Hajidavalloo, M. R., Cosner, J., Li, Z., Tai, W.-C., and Song, Z., 2022. “Simultaneous
suspension control and energy harvesting through novel design and control of a new
nonlinear energy harvesting shock absorber”. IEEE transactions on vehicular technology,
71(6), pp. 6073–6087.

[63] Li, Y., Wu, Z., Zhang, G., Wang, F., and Wang, Y., 2019. “Stochastic P-bifurcation
in a bistable Van der Pol oscillator with fractional time-delay feedback under Gaussian
white noise excitation”. Advances in Difference Equations, 2019(1), Dec., p. 448.

[64] Zhang, Y., Jin, Y., Xu, P., and Xiao, S., 2020. “Stochastic bifurcations in a nonlinear
tri-stable energy harvester under colored noise”. Nonlinear Dynamics, 99, pp. 879–897.

[65] Xu, Y., Gu, R., Zhang, H., Xu, W., and Duan, J., 2011. “Stochastic bifurcations in a
bistable duffing–van der pol oscillator with colored noise”. Physical Review E, 83(5),
p. 056215.

121

[66] Kumar, P., Narayanan, S., and Gupta, S., 2017. “Bifurcation analysis of a stochasti-
cally excited vibro-impact duffing-van der pol oscillator with bilateral rigid barriers”.
International Journal of Mechanical Sciences, 127, pp. 103–117.

[67] Sun, Y.-H., Sun, Y., Yang, Y.-G., and Xu, W., 2023. “Bifurcation and stability anal-
ysis of a hybrid energy harvester with fractional-order proportional–integral–derivative
controller and Gaussian white noise excitations”. Probabilistic Engineering Mechanics,
73, July, p. 103464.

[68] Cai, Q., and Zhu, S., 2022. “The nexus between vibration-based energy harvesting
and structural vibration control: A comprehensive review”. Renewable and Sustainable
Energy Reviews, 155, Mar., p. 111920.

[69] Vakakis, A. F., 2001. “Inducing Passive Nonlinear Energy Sinks in Vibrating Systems”.

Journal of Vibration and Acoustics, 123(3), July, pp. 324–332.

[70] Xue, J.-R., Zhang, Y.-W., Niu, M.-Q., and Chen, L.-Q., 2023. “Harvesting electricity
from random vibrations via a nonlinear energy sink”. Journal of Vibration and Control,
29(23-24), Dec., pp. 5398–5412.

[71] Kumar, P., Narayanan, S., and Gupta, S., 2016. “Investigations on the bifurcation
of a noisy duffing–van der pol oscillator”. Probabilistic Engineering Mechanics, 45,
pp. 70–86.

[72] Zhou, S., Cao, J., Inman, D. J., Lin, J., Liu, S., and Wang, Z., 2014. “Broadband
tristable energy harvester: Modeling and experiment verification”. Applied Energy,
133, Nov., pp. 33–39.

[73] Zhang, T., and Jin, Y., 2024. “An improved coupled tri-stable energy harvesting sys-
tem with spring stops for passive control”. Communications in Nonlinear Science and
Numerical Simulation, 135, Aug., p. 108050.

[74] Zhou, Z., Qin, W., and Zhu, P., 2017. “A broadband quad-stable energy harvester
and its advantages over bi-stable harvester: Simulation and experiment verification”.
Mechanical Systems and Signal Processing, 84, Feb., pp. 158–168.

[75] Daqaq, M. F., 2011. “Transduction of a bistable inductive generator driven by white and
exponentially correlated Gaussian noise”. Journal of Sound and Vibration, 330(11),
May, pp. 2554–2564.

[76] Sharif-Bakhtiar, M., and Shaw, S. W., 1992. “Effects of nonlinearities and damping
on the dynamic response of a centrifugal pendulum vibration absorber”. Journal of
Vibration and Acoustics, 114(3), 07, pp. 305–311.

[77] Huang, Z., Zhu, W., and Suzuki, Y., 2000. “Stochastic averaging of strongly non-linear
oscillators under combined harmonic and white-noise excitations”. Journal of Sound
and Vibration, 238(2), pp. 233–256.

122

[78] Soong, T. T., and Grigoriu, M., 1993. “Random vibration of mechanical and structural

systems”. NASA STI/Recon Technical Report A, 93, Jan., p. 14690.

[79] Namachchivaya, N. S., 1990. “Stochastic bifurcation”. Applied Mathematics and Com-

putation, 38(2), pp. 101–159.

[80] Bover, D., 1978. “Moment equation methods for nonlinear stochastic systems”. Journal

of Mathematical Analysis and Applications, 65(2), Sept., pp. 306–320.

[81] Higgins, B. G., and Binous, H., 2010. “A Simple Method for Tracking Turning Points in
Parameter Space”. JOURNAL OF CHEMICAL ENGINEERING OF JAPAN, 43(12),
pp. 1035–1042.

[82] Sun, J.-Q., and Hsu, C. S. “Cumulant-neglect closure method for nonlinear systems

under random excitations”. pp. 649–655.

[83] Wojtkiewicz, S., Spencer Jr, B., and Bergman, L., 1996. “On the cumulant-neglect
closure method in stochastic dynamics”. International journal of non-linear mechanics,
31(5), pp. 657–684.

[84] Kecik, K., and Borowiec, M., 2013. “An autoparametric energy harvester”. The Euro-

pean Physical Journal Special Topics, 222(7), pp. 1597–1605.

[85] Okada, Y., 2002. “Variable resistance type energy regenerative damper using pulse

width modulated step-up chopper”. J. Vib. Acoust, 14, pp. 110–115.

[86] Zuo, L., and Zhang, P.-S., 2013. “Energy Harvesting, Ride Comfort, and Road Handling
of Regenerative Vehicle Suspensions”. Journal of Vibration and Acoustics, 135(1), Feb.,
p. 011002.

[87] Psaros, A. F., and Kougioumtzoglou, I. A., 2020. “Functional series expansions and
quadratic approximations for enhancing the accuracy of the wiener path integral tech-
nique”. Journal of Engineering Mechanics, 146(7), p. 04020065.

[88] Petromichelakis, I., Psaros, A. F., and Kougioumtzoglou, I. A., 2020. “Stochastic re-
sponse determination of nonlinear structural systems with singular diffusion matrices:
A Wiener path integral variational formulation with constraints”. Probabilistic Engi-
neering Mechanics, 60, Apr., p. 103044.

[89] Kougioumtzoglou, I. A., Di Matteo, A., Spanos, P. D., Pirrotta, A., and Di Paola, M.,
2015. “An efficient wiener path integral technique formulation for stochastic response
determination of nonlinear mdof systems”. Journal of Applied Mechanics, 82(10),
p. 101005.

[90] Sundararajan, P., and Noah, S. T., 1997. “Dynamics of Forced Nonlinear Systems Using
Shooting/Arc-Length Continuation Method—Application to Rotor Systems”. Journal
of Vibration and Acoustics, 119(1), 01, pp. 9–20.

123

[91] Li, Z., Zuo, L., Luhrs, G., Lin, L., and Qin, Y.-x., 2013. “Electromagnetic Energy-
Harvesting Shock Absorbers: Design, Modeling, and Road Tests”. IEEE Transactions
on Vehicular Technology, 62(3), Mar., pp. 1065–1074.

[92] Huang, B., Hsieh, C.-Y., Golnaraghi, F., and Moallem, M., 2015. “Development and
optimization of an energy-regenerative suspension system under stochastic road excita-
tion”. Journal of Sound and Vibration, 357, pp. 16–34.

[93] Liu, X., Li, Y., Cheng, Y., and Cai, Y., 2023. “Sparse identification for ball-screw
drives considering position-dependent dynamics and nonlinear friction”. Robotics and
Computer-Integrated Manufacturing, 81, June, p. 102486.

[94] Bowen, L., Vinolas, J., and Olazagoitia, J. L., 2019. “The influence of friction parameters
in a ball-screw energy-harvesting shock absorber”. Nonlinear Dynamics, 96(4), June,
pp. 2241–2256.

[95] González, A., Olazagoitia, J. L., Viñolas, J., Ulacia, I., and Izquierdo, M., 2022. “An in-
novative energy harvesting shock absorber system for motorbikes”. IEEE/ASME Trans-
actions on Mechatronics, 27(5), pp. 3110–3120.

124

APPENDIX A: DRIFT AND DIFFUSION COMPONENTS OF MOMENT
EQUATIONS

The third and fourth components of the drift vector in (4.15) were derived and are given

here.

a3 =

−

(cid:16)

(cid:17)

(η + cos (ψ))

2 ˙ϕξp + η ˙θ2µr sin (ψ)

(cid:16)

η

θ + 2 ˙θ

(cid:16)

A
ξ − η ˙ϕµr sin (ψ)

(cid:17)

(cid:17)

− η ˙ϕ2µr sin (ψ)

A

a4 = (−µrηA)

−1 [2 ˙ϕξp + 2η2 ˙ϕξpµr + 4η ˙ϕξpµr cos (ψ)

+ 2 ˙ϕξpµr − 2η2 ˙θξµr + η3 ˙θ2µ2

r sin (ψ)

+ 2η3 ˙θ ˙ϕµ2

r sin (ψ) + η2 ˙θ2µ2

r sin (2 (ψ))

+ η2 ˙θ ˙ϕµ2

r sin (2 (ψ)) + η3 ˙ϕ2µ2

r sin (ψ) − 2η ˙θξµr cos (ψ)

+ η2 ˙ϕ2µ2

r sin (ψ) cos (ψ) + η ˙θ2µ2

r sin (ψ) + η ˙θ2µr sin (ψ) − ηθµr (η + cos (ψ))] (A.1)

where A = −η (µr cos2 (ψ) − µr − 1). The drift vector, a, is then given by

1

C
C
C
C
C
C
C
A

0

B
B
B
B
B
B
B
@

˙θ

˙ϕ

a3

a4

a =

(A.2)

Likewise, Σ in (4.15) is a 4 × 1 diffusion matrix consisting of the stochastic components

of the system. It is found by collecting coefficients of W and multiplying by the standard

√

deviation, σ, of the white noise. Namely, σ =

2D, where D is the intensity of the white

noise; that is,

0

B
B
B
B
B
B
B
@

Σ =

1

C
C
C
C
C
C
C
A

0

0

σ
µr sin2(ψ)+1
− σ(η+cos(ψ))
η+ηµr sin2(ψ)

125

(A.3)

Eqn. (4.14) can then be integrated numerically with appropriate solvers to generate realiza-

tions in the time domain.

126

APPENDIX B: CLOSURE TECHNIQUES

B.1 Gaussian Closure

The generated moment equations (2.12) are dependent on moments of higher order than

n = 2. One course of action is to use a technique called Gaussian Closure in order to relate

the higher order moments to moments of first and second order in order to “close” the system.

We more precisely assume

⟨yiyjykyl⟩ ≃ −2 ⟨yi⟩ ⟨yj⟩ ⟨yk⟩ ⟨yl⟩ + ⟨yiyl⟩ ⟨yjyk⟩

+ ⟨yiyk⟩ ⟨yjyl⟩ + ⟨yiyj⟩ ⟨ykyl⟩

⟨yiyjyk⟩ ≃ ⟨yi⟩ ⟨yl⟩ ⟨yjyk⟩ − 2 ⟨yi⟩ ⟨yj⟩ ⟨yk⟩ + ⟨yk⟩ ⟨yiyj⟩

+ ⟨yj⟩ ⟨yiyk⟩

(B.1)

Substituting (B.1) into (2.12) leads to a set of fourteen nonlinear equations coupled through

first and second-order moments.

B.2 Cumulant neglect

For the computation of mean-squares, a cumulant neglect procedure was used to derive

stationary moment equations with moments mpqrs = ⟨yp

1yq

2yr

3ys
4

⟩ up to order n = 4. We

accomplish this by first writing the characteristic function Fy(u) in terms of cumulants

λk(p, q, r, s) of order k. Specifically,

Fy(u) = exp

∞X

k=1

ik
k!

nX

p,q,r,s=1

λk(p, q, r, s)up

1uq

2ur

3us
4

!

(B.2)

√

where i =

function by

−1 and p + q + r + s = k. Moments mpqrs are then related to the characteristic

(cid:20)

mpqrs = (−i)ζ

∂ζFy(u)
1∂uq
2∂ur

3∂us
4

∂up

(cid:21)

u=0

(B.3)

where ζ = p + q + r + s. We can then set cumulant orders λk = 0 for k > 4 allowing us to

establish a set of equations relating moments of fifth order and higher to lower orders. The

127

 
resultant relationships and (2.18) are substituted back into stationary moment equations

(2.12), giving 69 algebraic equations.

128

APPENDIX C: LINEAR SYSTEM RESPONSE

In order to quantify the vibration suppression imposed by the IPVA, the closed form

solution for the system without the IPVA is derived. The moment differential equations are

used to determine the solution for the mean square velocity. The equation of motion for the

linear system is given by

θ′
l =

Ω

Ω′ = W − 2 ˆξΩ − θl

(C.1)

Additionally, the white noise input follows a Gaussian distribution with covariance matrix
0

1

ΣΣT =

0

B
@

0

C
A

0 σ2

Lastly, from (C.1), the drift components are

α1 =

Ω

α2 =

−2 ˆξΩ − θl

(C.2)

(C.3)

One can finally input (C.2) and (C.3) into (2.12) to obtain the closed-form solution for

the dimensionless mean square velocity given here.

(cid:11)

(cid:10)

θ′2
l

=

D
2 ˆξ

(C.4)

129

APPENDIX D: LINEAR SYSTEM RHI ANALYTICAL SOLUTION

In this section, the closed-form linear system solution for the road handling index is de-

rived using power spectral density integration. To this end, the tire displacement R (ψus − ψr)

is transformed into the frequency domain, and the ratio of tire displacement to dimensionless

input road velocity is determined. Note that for the remainder of this section, we will set

ˆω = ω/ω0. Then,

with

R (ψus (ˆω) − Ψr (ˆω))
r(ˆω)

Ψ′

= R

N1 + N2 + N3
D1 + D2 + D3 + D4 + D5

,

(D.1)

N1 = (−iu − i) ˆω,

N2 = ˆω2 (2uξ + 2uξe + 2ξ + 2ξe) ,

(cid:0)

N3 = ˆω3

iη2µr + 2iηµr + iµg + iµr + iη2µru + 2iηµru + iµgu + iµru + iu

(cid:1)

,

D1 = uv2,
(cid:0)

D2 = ˆω

D3 = ˆω2

2iuv2ξ + 2iuv2ξe
(cid:0)

(cid:1)

,

−η2µruv2 − 2ηµruv2 − µguv2 − µruv2 − uv2 − u − 1

(cid:1)

,

D4 = ˆω3 (−2iuξ − 2iuξe − 2iξ − 2iξe) ,

(cid:0)

D5 = ˆω4

η2µr + 2ηµr + µg + µr + η2µru + 2ηµru + µgu + µru + u

(cid:1)

√

with i =

−1. Then using the dimensionless power spectral density of the road velocity as

SΨ′

r = d with the principle of H2 norm and integration formula from [14], the road handling

index becomes

Finally,

γRHI =

R
g (Ms + Mus)

s

Z ∞

−∞

(cid:12)
(cid:12)
(cid:12)
(cid:12)

d

1
2π

(ψus (ˆω) − Ψr (ˆω))
r(ˆω)

Ψ′

2

(cid:12)
(cid:12)
(cid:12)
(cid:12)

dˆω

(D.2)

γRHI =

h

(cid:16)

R
g (Ms + Mus)

1
4u2v4(ξ + ξe)
+ u2(3 + 2v2( − 1 + 4(ξ + ξe)2 − 2µg − 2(1 + η)2µr)

1 + u(3 + v2(4(ξ + ξe)2 − 2µg − 2(1 + η)2µr))

130

+ v4(µg + (1 + η)2µr)(1 + µg + (1 + η)2µr))

+ u3(1 + v4(1 + µg + (1 + η)2µr)2

+ v2(4(ξ + ξe)2 − 2(1 + µg + (1 + η)2µr)))

(cid:17)i 1

2

(D.3)

131

APPENDIX E: DAMPING AND FRICTION FITTING IN THE TIME
DOMAIN FOR THE SDOF SYSTEM EXPERIMENT

Pendulum bearing damping cp was also determined using a fitting procedure, but this

time in the time domain. To do this, the gears were removed and the pendulums were simply

given an initial velocity while keeping the carrier fixed and the time response was measured.

The equation of motion for the pendulums in this case is

(Jpendulum) ¨ψ + cp

˙ψ = 0

(E.1)

where Jpendulum is the pendulum’s principal moment of inertia about the pivot point of the

pendulums. For simplicity, only the pendulum vertical shafts were used in the experiment

and then Jpendulum was equal to the moment of inertia of the vertical shaft about its vertical

axis. The equation used for fitting with MATLAB was then

(cid:20)

˙ψ (t) = ˙ψ (0) exp

−

(cid:21)

cpt
(Jpendulum)

(E.2)

A similar experiment was done for the generator mechanical damping cgm. This time

the generator and sun gear were disconnected from the prototype and a time response to an

initial condition was measured. The solution in this case is

˙ψg (t) = ˙ψg (0) exp

(cid:21)

(cid:20)

− cgmt
n2
gJr + Jpg

(E.3)

Finally, with pendulums and planetary gear reattached and the carrier held fixed, the

same experiment was done to determine the friction term Tf . The EOM in this case is given

by

(cid:0)

2Jpg + Jsg + n2

gJr + 2Jpendulum

(cid:1)

¨ψ + (cp + cgm) ˙ψ = −Tf

˙ψ
| ˙ψ|

For a monotonically decreasing ˙ψ, the solution to (E.4) is then

˙ψ (t) =

(cid:18)

Tf
(cp + cgm)

(cid:19)

(cid:20)

+ ˙ψ (0)

exp

− (cp + cgm) t
(Jtotal)

(cid:21)

−

Tf
(cp + cgm)

(E.4)

(E.5)

132

Figure E.1 Raw data for twelve sequential experiments (solid line). Fit based on (E.5)
(dashed line).

with

Jtotal = 2Jpg + Jsg + n2

gJr + 2Jpendulum

Twelve sequential experiments of this kind were done to accommodate potential variability

in the data and the resultant mean fit using (E.5) is shown in Fig. E.1.

Finally, fitted parameters and measured parameters used in this work are given in Table

E.1. Note that µr was measured directly and SolidWorks was used to compute µp as well as

η. Note that ξp was determined to be negligible and was set to zero.

Dimensionless Parameters Value
µr
η
µp
µg
ξ
ξgm
α
ˆTf

2.43
0.45
0.38
0.28
0.063
0.004
5.93
0.004

Table E.1 Dimensionless Parameters of Experimental IPVA system.

133