Vibro-impact phenomena in engineering systems, considered an adverse effect in some settings, are an intrinsic part of the mechanism in others. In energy harvesting, a vibro-impact component is often intentionally introduced to increase the power output or the system’s bandwidth. The impacts can be treated as “hard” for instantaneous impacts or “soft” for compliant materials. Since both types of models exhibit complex dynamics, a comparison is non-trivial. We develop a soft impact model for a vibro-impact energy harvester, calibrating it with the relevant hard impact model for large stiffness, and systematically compare the different phenomena and dynamics in various compliant regimes. Numerical results are used in two different parametric analyses, considering the bifurcation diagrams in terms of device size and external forcing parameters. Varying the natural frequency of the membranes that form the impact boundaries, we observe shifts in the bifurcation structure that promote period-1 orbits for increased softness parameters, often generating higher power output, but also introducing parameter sensitivities for increased softness. Complementary analytical results reveal unstable periodic orbits and co-existing behaviors, potentially missed by computational methods, that can influence the bifurcation structure and in turn the energy output. A non-dimensional formulation highlights the significance of ratios of external and natural frequencies in delineating soft and hard impact scenarios parametrically. The soft impact model exhibits new symmetry breaking bifurcations related to key quantities that characterize the soft impact dynamics, such as the effective restitution coefficients, the impact phase, and the contact time interval, not captured by hard impact models.

Impacts are present in a variety of physical systems and are known to introduce nonlinearities that can drastically change the system response. Predicting the influence of impacts requires a careful approach that can handle significant complexity since there are different types of physical phenomena depending on the stiffness, geometry, and relative velocity of the impacting bodies. Despite the prevalence of impacts in several applications and mechanisms, few studies compare different models or study their limitations. Even the comparison between simple compliant (soft) and non-compliant (hard) models is not straightforward due to the complex dynamics presented by them. Thus, we develop a soft impact model for a vibro-impact energy harvester, calibrating it with the relevant hard impact model for large stiffness, and systematically compare the different phenomena and dynamics in various compliant regimes. A non-dimensional formulation highlights the significance of ratios of external and natural frequencies in delineating soft and hard impact models parametrically. The soft impact model exhibits new symmetry breaking bifurcations related to key quantities that characterize the soft impact dynamics, such as the effective restitution coefficients, the impact phase, and the contact time interval, which cannot be represented by hard impact models. We observed shifts in the bifurcation structure that promote period-1 orbits for increased softness parameters, often generating higher power output, but also introducing parameter sensitivities for increased softness, which can dictate optimum values for energy harvesting applications. Furthermore, we develop analytical results that complement computational methods and reveal unstable periodic orbits and co-existing behaviors that can influence the bifurcation structure and energy output. Overall, we manage to provide significant advances in the design of vibro-impact energy harvesters and the different phenomena modeled by hard and soft impact formulations.

## I. INTRODUCTION

Impacts are present in a variety of physical systems and are known to introduce nonlinearities that can drastically change the system response. Predicting the influence of impacts requires a careful approach that can handle significant complexity since there are different types of physical phenomena depending on the stiffness, geometry, and relative velocity of the impacting bodies. Despite the prevalence of impacts in some applications and mechanisms, such as gears, hammer drills, jack hammers, and others, few studies compare different models or study their limitations, while most rigorous explorations have focused primarily on the behavior of specific mechanisms such as rigid or elastic impacts, or specific geometries. For example, even though models of impacts have been studied for decades, very few reports compare soft and hard model predictions for the same impacting system, as in Refs. 1–3.

In mechanical systems, the appearance of impacts may be due to tolerances, specific design, or aging, just to name a few. The dynamics driven by impacts is typically undesirable in some applications such as machining where the introduction of nonlinearities can lead to chattering,^{4} or in rotors where impacts increase wear, leading to increased maintenance^{5} and sometimes causing system failure.^{6} In other settings, impacts can be a design feature such as in seismic mitigation where gaps are introduced to uncouple the structure from its surroundings,^{7} or in energy harvesters where impacts are a source of nonlinearities that can expand the operational bandwidth. Optimizing performance in any of these settings requires a thorough understanding of the resulting nonlinear response.^{8,9}

Some of the first studies of impact phenomena in mechanical engineering were performed by Shaw *et al.*^{10–12} where the authors use an impact oscillator system composed of a mass attached to a spring and dash-pot and an impact beam placed at a certain distance from the mass. In these early studies, impacts are considered non-complaint (hard impacts), where impacts are instantaneous and the relationship between the velocities of the impacting objects is captured via a restitution coefficient.^{13} These models can explain several behavioral patterns of an impacting system like chattering,^{15} grazing,^{13} and other nonlinear phenomena.^{16} However, since they do not describe compliant barriers (soft impacts), several subsequent works were focused on modeling the contact forces as a continuum over a time interval as summarized by Gilard and Sharf.^{14}

Following these initial studies, several authors explored analytical solutions of vibroimpact systems, including analysis of grazing phenomena as in Ref. 17. Others considered grazing numerically as in Ref. 18 reporting hard and soft grazing bifurcations and developing a numerical method to evaluate grazing points.^{19} Grazing was also investigated in the soft impact setting, using both models and experiments for impact oscillators.^{20–23} Others have studied the influence of impacts in stochastic systems,^{24} in the effectiveness of control in impacting systems,^{25,26} and within the context of nonlinear bumpers.^{27}

Recently, a number of studies have analyzed the differences between soft and hard impact models and the relevant limitations of these models. Rebouças *et al.*^{28} calibrated a hard impact model with experimental results, using a relatively large stiffness in order to mimic the hard impact model. They verified that even in this setting, the experimental restitution coefficient was dependent on the impact velocity and on the excitation phase at impact. Blazejczyk-Okolewska *et al.*^{29} calibrated the energy dissipation in a soft impact model by utilizing the restitution coefficient of a hard impact model in a single sided impact oscillator. Their results showed that the soft model reproduces the behavior of the hard model for a range of restitution coefficients. Later, Okolewski *et al.*^{1,2} extended the calibration to all values of restitution coefficients, providing their results in bifurcation diagrams in terms of excitation frequency. However, they did not pursue interpreting the different behaviors and the physical phenomena underlying the bifurcation structure. Chehaibi *et al.*^{30} made a brief comparison on the energy dissipation of linear and nonlinear boundary models but did not include hard impacts.

The use of soft or hard impact models can greatly affect the analysis and design in applied settings. For example, impacts have been introduced in the design of several energy harvesters as a source of nonlinearities to improve the device’s performance. Drezet *et al.*^{31} compare different designs for an energy harvester including hard stops in the system. Likewise, Halim and Park^{8} introduce impacts in their energy harvester design, thus broadening the frequency band for efficient operation. Several designs use impacts and contact in conjunction with different harvesting methods, such as in triboelectric devices,^{32} magnet induction harvesters,^{33} and piezoelectric based^{34,35} or electrostatic harvesters.^{9} These designs can be developed for a wide range of energy sources such as ocean waves, structural vibrations, human locomotion, etc. In energy harvesting of structural vibrations, a series of studies by Serdukova *et al.*^{36,37} employs maps to investigate a two-sided vibro-impact energy harvester (VI-EH). In addition to the analysis, they compare analytical and numerical results describing motions composed of impacts on alternating ends of the harvester, as well as grazing incidences that lead to transitions between different types of impact sequences. They also perform a parametric study that indicates design considerations for optimizing the harvester’s power output. These studies are based on the hard impact model with a constant restitution coefficient, leaving many open questions about the corresponding soft impact model for this system that can translate to real applications.^{28}

The goal of this study is to develop a soft impact model for the vibro-impact energy harvesting pair previously described in Ref. 36 and to systematically compare its dynamics with that of the corresponding hard impact model. In the model, a ball moves freely within an externally excited cylindrical capsule, between impacts with dielectric polymer membranes on either end. Energy is generated as the ball impacts and deforms the membranes, and, thus, an accurate model for these impacts is critical. The comparison between hard and soft impact models is achieved primarily through numerical simulation and maps of the dynamics between impacts. In contrast to Refs. 36–38, the hard impact model is used as a starting point for calibration of the soft impact model with a large stiffness parameter. Then, appropriate parameters are varied in order to compare the differences and similarities of the hard and soft impact models when varying the stiffness of the impact boundary. We also discuss the physical phenomena behind the observed differences, demonstrating how this understanding guides the system design for optimal energy output. Our study has several differences as compared with Refs. 1 and 2; we study a model based on the two-sided impact pair, and we capture the variation in dynamics and bifurcations in terms of several different parameters: device length, stiffness of the impact boundary, and excitation frequency. Through appropriate non-dimensionalizations, we obtain an insight into a number of unexpected results; for example, even with stiff boundaries in both the soft and hard impact models, they can disagree if the excitation frequency is sufficiently high. We also identify scenarios where a softer membrane leads to a sharp increase in energy output, with increased parameter sensitivity near the optimal output. Furthermore, the effective restitution coefficients exhibited by the soft model are dependent on the excitation phase at impact and may undergo symmetry breaking bifurcations that produce asymmetric motions with novel dynamical properties, not observed in the hard impact case. We explore these and other differences within the particular context of VI-EH, given the range of applications recently reported,^{39} and with a particular interest in energy harvesting from mechanical and civil structures.

This work is organized as follows. In Sec. II, we describe the vibro-impact pair system and give the equations of motion for the corresponding hard and soft impact models. We derive the (semi)-analytical solutions for both models using an approach that maps the dynamics between consecutive impacts. This leads naturally to a technique for calibrating the soft impact model with the hard impact model. In addition, we give the electric model, describing the extraction of the power generated by each impact and defining an appropriate power output metric for the device. In Sec. III, we calibrate hard and soft impact models for large stiffness and perform two separate parametric analyses of the soft model response as the stiffness varies. Also, we analyze the energy output of the system in different circumstances, illustrating how these results can guide device design for optimal energy output. Finally, we present the conclusions of the work and future directions.

## II. THE VIBRO-IMPACT MODEL

This study addresses the main differences between soft and hard impact models in the context of a two-sided VI-EH. The mechanism is presented in Fig. 1(a), showing a rigid cylindrical capsule with mass $M$ (gray in the figure) and displacement $Z$. We consider an harmonic excitation $Fext$ of the frame as we consider that vibrations are filtered by an external structure where the harvester is attached. A ball (black) with mass $mb$ is free to move within the capsule along its axis without friction, with displacement $Y$. On each end of the capsule, two dielectric polymer membranes (yellow disks) are placed in such a way that the ball can move a distance $s$ from one membrane to another without contact. The apparatus is also considered to be in a constant gravitational field of acceleration $g$ and at an angle $\beta $ from the horizontal one. Here, we focus on the case $\beta =0$ and include some observations about the $\beta \u22600$ in Sec. III.

The equations of motion when the ball is not in contact with the membrane can be expressed as

where $\tau $ is the time, $\omega $ is the excitation frequency of the main structure, $A0$ is the amplitude of excitation in newtons, and $\u2032$ represents the derivative with respect to time. Assuming that $M\u226bmb$, the ball dynamics do not affect the motion of the capsule, and the equations of motion can be simplified by writing them in terms of relative displacement $X=Y\u2212Z$,

The schematics for the hard and soft impact models are shown in Figs. 1(b) and 1(c), respectively.

### A. Modeling the hard impact dynamics

In the hard impact model, the membranes are non-compliant and the impact is treated as instantaneous at the time $\tau j$ of the $jth$ impact. Hence, the relationship between the instantaneous velocities immediately before [$X\u2032(\tau \u2212)$] and after [$X\u2032(\tau +)$] the impact with the membrane wall is defined by a coefficient of restitution $r$, which together with Eq. (3) completes the hard impact model,

Since the system is linear between impacts, the solution can be written explicitly, by integrating Eq. (4) with respect to time, using the condition of Eq. (5), and the continuity of $X$ across impacts. For example, for $\tau j<\tau <\tau j+1$, the motion is given by

where $\varphi j=mod(\omega \tau j,2\pi )$ is the excitation phase at $\tau j$ and $+(\u2212)s/2$ corresponds to the last impact $X(\tau j)$ located on the right (left) side. Equations (5) and (7) are then used to determine $X\u2032(\tau j+1\u2212)$ and $\tau j+1$.

Equations (6) and (7) with $\tau =\tau j+1$ provide the definition of a map of the dynamics between two consecutive impacts. This representation has the advantage of explicitly stating the system dynamics in closed form for each map, allowing analyses beyond purely numerical integration methods as shown in Refs. 37 and 40. Then, the different types of dynamics in the hard impact model can be expressed via combinations of four different maps, $Mi$, $i\u2208(1,2,3,4)$, as shown in Fig. 2, used to construct the dynamics of the hard model. The maps $M1$ and $M2$ correspond to cases where the consecutive impacts are on different ends of the capsule, with $M1$ from left to right and $M2$ from right to left. Analogously, $M3$ and $M4$ correspond to the evolution between two impacts on the same end of the capsule, on the left for $M3$ and on the right for $M4$. The initial and terminal conditions for the map $Mi$ between two impacts at $\tau j$ and $\tau j+1$ are given by

Evaluating Eqs. (6) and (7) at $\tau =\tau j+1$, the maps are in terms of $X\u2032(\tau j+1\u2212)$ $X\u2032(\tau j\u2212)$ only. Then, without loss of generality, we drop the $\u2212$ superscript throughout the remainder of this manuscript. Given the impact time $\tau j$, with corresponding displacement and velocity $X(\tau j)$ and $X\u2032(\tau j)$, respectively, the time of the next impact $\tau j+1$ for the map $Mi$ is determined from Eq. (6) subject to the initial and terminal conditions of Eq. (8),

where $\Delta j$ is the time interval between consecutive impacts. Since Eq. (9) is a transcendental equation for $\tau j+1$, it has several solutions and in general must be solved numerically. To ensure that the solution is physically relevant, that is, corresponding to consecutive impacts, we seek the smallest $\tau j+1>\tau j$. After $\tau j+1$ is obtained, Eq. (7) can be used to calculate $X\u2032(\tau j+1)$,

We highlight that, by using maps, one can calculate the system orbit for a given initial condition, but this method cannot capture sticking conditions as that would require an infinite sequence of maps.

To distinguish between impacting orbits described by different sequences of $Mi$, we use the notation $n$:$m$/$pT$ to categorize periodic orbits for the VI-EH device with $T$-periodic external excitation, as in Ref. 37. Here, $n$ ($m$) is the number of impacts on the left (right) of the capsule per time interval $pT$, for $p$ an integer. In the case where $p=1$, we simplify the notation to $n$:$m$. The notation captures period doubling sequences, e.g., $n$:$m$/$pT$ for $p=2,4,\u2026$ is a period doubling sequence. Grazing bifurcations, where the trajectory has a zero relative velocity at impact, yield other types of instabilities and transitions;^{36,37} for example, in the VI-EH system, a grazing incidence typically precedes a transition from a $n$:$m$ (or $n$:$m$/$pT$) periodic orbit to a ($n+1$):$m$ orbit, while, more generally in other cases, the transition may be from $n$:$m$/$pT$ to $\u2113$:$q$ for $\u2113\u2260n$ and/or $q\u2260m$. Later, we also consider higher excitation frequencies $\omega $ for which it is possible to have $pT$ periodic orbits with a low number of impacts over $pT$.

### B. Modeling and calibrating the soft impact dynamics

In contrast to the hard impact setting, the soft impact model includes dynamics for impacts with the compliant membranes on both ends of the capsule, with elastic coefficient $k$ and a linear dissipation coefficient $c$, shown in Fig. 1(c). Hence, the soft impact model includes additional equations of motion that account for the motion of the ball while it is in contact with one of the membranes. Specifically, maps $M5$ and $M6$, shown in Fig. 2, must be included in order to describe the dynamics of the ball while in contact with the membrane.

The equations of motion [Eq. (1)] are the same for both the hard and soft impact models, since these capture the dynamics when the ball is not in contact with the membrane. The difference is in the way the two models treat the contact with the membrane. We define $U=Y\u2212Z$ as the relative displacement for the soft impact model when the ball is not in contact with the membrane, to distinguish it from $X$ used for the hard impact model. Furthermore, we use $Uc$ to denote the relative displacement while the ball is in contact with either membrane, corresponding to maps $M5$ and $M6$. Combining the equations of motion for $U$ and $Uc$ yields the system,

In general, there is a finite time interval for contact with the membrane, which we denote as

for $\tau j,c$ ($\tau j,f$) the time at which the ball initiates (terminates) contact with the membrane. During this interval, $Uc$ describes the motion of the system. Furthermore, we assume continuity of displacement and velocity at these times,

In order to systematically identify the differences between the soft and hard impact models under external excitation, we calibrate parameters in the setting where both models are expected to yield similar results, namely, for large values of $k$, approximating $k\u2192\u221e$. For large $k$, the dissipative effects linked to the linear dissipation coefficient $c$ in the soft model and restitution coefficient $r$ in the hard model should result in the same dynamics for both models. To facilitate our analysis, we calibrate the models in the absence of excitation $A=0$, observing that the properties of the membrane such as the dissipation coefficient $c$ and restitution coefficient $r$ do not change with external influences. We consider simply the dynamical behavior of the membrane with a given initial velocity at the beginning of the contact interval. Then, using the dynamics given by Eq. (11) in the large $k$ limit and $A=0$, where $\tau j,c\u2248\tau j,f$, we find the value $c$ for which the ratio $\rho $ of velocities before and after contact with the membrane is comparable to $r=\u2212X\u2032(\tau j+)/X\u2032(\tau j\u2212)$, that is,

We refer to this ratio $\rho $ as the effective restitution coefficient of the soft model, obtained for all $k$, as studied in Sec. III. For $A=0$ and given the initial contact velocity $U\u2032(\tau j,c\u2212)$, $Uc(\tau )$ is obtained from Eq. (11),

where $\omega 1=\omega 02\u2212(c/2mb)2$. Note that for large $k$, the natural frequency associated with the mass-membrane system, $\omega 0=k/mb$, is also large. Then, the contact time $\Delta \tau j,c$ in Eq. (12) is small, that is,

for large $k$. By differentiating $Uc(\tau )$ in Eq. (15) with respect to $\tau $, we obtain the velocity of the ball while in contact with the membrane,

To complete the calibration, we evaluate $Uc\u2032(\tau )$ in Eq. (17) at $\tau =\tau j,f$ from Eq. (16), and divide by $U\u2032(\tau j,c\u2212)$ to get the expression for $\rho $,

Then, we set $\rho =r$ and solve for $c$ in terms of $r$,

By this calibration of $c$, the soft impact model mimics the impact-induced energy losses of the hard impact model when $A=0$ and $k$ is large. In that setting, Eq. (16) indicates a small $O(k\u22121/2)$ difference between the models, observed for the shift in impact times.

Throughout the comparisons in Sec. III, for a given set of parameters $[mb,r]$ in the hard impact model, together with a given large $k$ in the soft model, the value of $c$ is calibrated as in Eq. (19). Then, $c$ varies with $k$ according to Eq. (19) for all values of $k$ that are studied for the given $mb$ and $r$. The calibrated value of $c$ for $A=0$ facilitates a systematic comparison of the soft and hard impact models under excitation, allowing us to track the differences as the ball impacts the membranes with different values of the elastic coefficient $k$. Since we are interested in the forced system for $A\u22600$, the excitation phase influences the system dynamics and $\Delta \tau j,c$ is no longer negligible in these comparisons.

### C. Analytical results for the soft impact model

It is useful to obtain analytical expressions for the relative displacement and velocity over the contact interval $\Delta \tau j,c$, which provide representations of the dynamics that are complementary to numerical simulations. While a complete analytical treatment is not the focus of this study, we provide some examples in Sec. III A that show how these results may provide insight beyond the computations.

Solving for $Uc(\tau )$ and $Uc\u2032(\tau )$ from Eq. (11) and using continuity as in Eq. (13), we obtain the analytical expression for the displacement and velocity of the ball while in contact with the soft membrane,

where $\zeta =c/2mb$, $\omega 12=\omega 02\u2212\zeta 2$ and the signs $\u2213$ and $\xb1$ indicate whether the ball is in contact with the right (top sign) or left (bottom sign) membrane. The constants $CI$ and $CII$ are given by the particular solution of the system, while $CIII$ and $CIV$ are determined by continuity in Eq. (13), to get

where $\varphi j,c=mod(\omega 1\tau j,c,2\pi )$. The transition time $\tau j,f$ is obtained by setting $Uc(\tau )=\u2213s/2$ and $\tau j+1,c$ can be determined from the equations for $Mi$ for $i=1,\u2026,4$, analogous to the procedure for hard impacts. Then, we have the set of equations,

where $\Delta \tau j,f=\tau j+1,c\u2212\tau j,f$, and $\varphi j,f=mod(\omega 1\tau j,f,2\pi )$.

It is useful to recall that the maps $Mi$ for $i=1,\u2026,4$ are defined in the same way as in the hard impact model away from impact, with relative velocity continuous across $\tau j,c$ and $\tau j,f$ as in Eq. (13). After determining $\Delta \tau j,c(\Delta \tau j,f)$, one can then differentiate the maps from Eqs. (24) to (25) with respect to time and find the velocity at which the ball leaves the membrane $U\u2032(\tau j,f)$ and the velocity of the next impact $U\u2032(\tau j+1,c)$, completing the definition of the maps,

These equations can be used to find the periodic solutions semi-analytically using the periodicity conditions $Uc\u2032(\tau j,c)=Uc\u2032(\tau j+l,c)$, $Uc(\tau j,c)=Uc(\tau j+l,c)$ and $mod(\tau j,c,2\pi /\omega )=mod(\tau j+l,c,2\pi /\omega )$ for some integer $l$. This is done by finding a solution, given by a specific sequence of $l$ maps, that minimizes the cost function,

This approach, which we refer to as the optimization method, can detect any given solution of the system. Its stability can be classified by the matrix product of the corresponding sequence of Jacobians,

where $Ji$ is the Jacobian of the $ith$ map.^{36}

We demonstrate the use of this semi-analytical method in Sec. III A, illustrating how it finds unstable solutions that are not captured by simple numerical integration, but nevertheless influence the dynamics. It is also useful in identifying bi-stable behavior, which may or may not be found via simulation, depending on the choice of initial conditions (IC). Based on other analytical studies,^{36–38,40} we expect this semi-analytical method to be useful for future parametric studies and understanding influence of sequences of bifurcations on the energy output.

### D. Energy harvesting

So far, we have provided the results in terms of the relative displacements and velocities. These results, particularly those for the impact velocities, play a central role in the energy harvesting capabilities of the VI-EH device. By applying an electric potential difference $V$ along the membrane thickness, the deformation by the impacting ball yields a change of the membrane capacitance, consequently generating electrical energy. The change in capacitance can be calculated using the surface area $As$ of the polymer and its thickness $w$. Both of these quantities can be written in terms of the penetration of the ball $\delta =|x|\u2212s/2$, the ball radius $rb$, and the membrane initial radius $R0$, resulting in

with

Assuming that the membrane is incompressible, its change in thickness $w$ due to deformation is given by the conservation of volume,

Finally, when a voltage is applied along the thickness of the membrane, its capacitance $CM$ can be approximated by

where $\u03f5M$ is the membrane permittivity.

The change in capacitance due to an impact produces a change in potential energy that can be harvested using the circuit model described in Ref. 9, with the energy $Ej$ extracted from the device at the $jth$ impact given by

Here, $Cmax,j=CM(\delta max,j)$ is the maximum capacitance of the polymer reached at the maximum penetration of the ball $\delta max,j$, $Cmin=CM(0)$ is the minimum capacitance reached when the membrane has no deformation, and $CT$ is the capacitance of a transfer capacitor used to transfer the energy from the circuit to the battery. It is important to highlight that Eq. (34) assumes that the circuit dynamics is much faster than that of the mechanical system, and, hence, the dynamics of both systems are decoupled from each other.

The maximum displacement $\delta max,j$, appearing in the definition of energy $Ej$ via $As$ [Eq. (30)], is given by the displacement $Uc$ when $Uc\u2032=0$ in the soft model and, thus, is calculated directly from Eq. (21). In contrast, the hard model does not include an expression for $\delta max,j$, since it is based on impact with a rigid wall. Thus, the calculation of $Ej$ for the hard impacts uses an expression that is calibrated to the soft model, similar to the calibration of the linear dissipation coefficient $c$. We recall that for large $k$ and without excitation the soft and hard impact models exhibit similar behavior, so we use the maximum value of $Uc$ without excitation from Eq. (15) to provide $\delta max,j$ for the hard impact model. Specifically, we take the maximum of Eq. (15), based on a calibrated $c$ as in Eq. (19), to obtain the maximum displacement for the hard model,

Finally, $Ej$ can be used to calculate the average energy output $P$ per half period for a specific n:m/$pT$ periodic orbit or, analogously, the time-averaged energy output for aperiodic attractors of the system. We define $P$ as the sum of the energy output over impacts, divided by an appropriate time interval length $\tau P$ for that behavior,

For example, if the orbit has period $pT$, then $Nimp$ is the total number of impacts in $pT$, with $Ej$ the energy generated by the $jth$ impact. In that case, we take $\tau P=2pT$ in order to compare $P$ and $Ej$, the energy for a single impact. For example, for 1:1 orbits, $P$ and $Ej$ should have the same value. If the behavior is aperiodic, then $P$ is computed over a larger time interval $\tau P$ to capture its approximate time-averaged output. Throughout this paper, we calculate $P$ using the electrical and dimensional parameters as provided in Table I.

Description . | Symbol . | Value . | Unit . |
---|---|---|---|

Initial thickness of the membrane | w_{0} | $0.00025$ | m |

Initial radius of the membrane | R_{0} | 0.006 3 | m |

Membrane permittivity | ε_{M} | 4.15 × 10^{−11} | F/m |

Circuit transfer capacitance | C_{T} | 99.35 | nF |

Applied voltage | $V$ | 2000 | V |

Radius of the ball | r_{b} | 0.005 0 | m |

Description . | Symbol . | Value . | Unit . |
---|---|---|---|

Initial thickness of the membrane | w_{0} | $0.00025$ | m |

Initial radius of the membrane | R_{0} | 0.006 3 | m |

Membrane permittivity | ε_{M} | 4.15 × 10^{−11} | F/m |

Circuit transfer capacitance | C_{T} | 99.35 | nF |

Applied voltage | $V$ | 2000 | V |

Radius of the ball | r_{b} | 0.005 0 | m |

Another important measure in energy harvesting is the efficiency of the device. In this study, we focus on output energy since the comparisons in this work are restricted to cases that are subjected to the same input energy. Thus, an improvement in output energy will consequently mean an improvement in efficiency.

## III. PARAMETRIC ANALYSIS OF HARD VS SOFT IMPACTS

We analyze the effects of key parameter combinations, whose variation yields different dynamics in the soft and hard impact models. Two different non-dimensionalizations of the equations facilitate these comparisons: variation of the system parameters of capsule size $s$ or excitation amplitude $A$ vs natural frequency $\omega 0$ of the membrane, and a second focused on the variation of excitation frequency $\omega $ vs natural frequency $\omega 0$.

### A. Influence of effective size vs natural frequency

First, we define a non-dimensional parameter $d=s\omega 2A\pi 2$ based on the ratio between the capsule length and excitation amplitude, thus representing the effective capsule size. Following the form of $d$, the non-dimensional displacement is defined as $x=X\omega 2A\pi 2$ for the hard impact model, which satisfies the non-dimensionalized equations of motion,

For the soft impact model, we define the non-dimensional displacement as $u=U\omega 2A\pi 2$, which satisfies

Here, $\kappa =\omega 0\pi \omega $ is proportional to the ratio of natural frequency $\omega 0$ and the excitation frequency $\omega $, and $\u02d9$ indicates the derivative with respect to the non-dimensional time $t=\tau \omega \pi $. There are also non-dimensionalized parameters $\gamma =gsin\u2061(\beta )A$, capturing the influence of gravity via the ratio $g/A$ of gravitational acceleration to excitation amplitude, and the non-dimensionalized dissipation coefficient of the membrane $\mu =c\pi m\omega $.

An analysis of Eqs. (38) and (37) reveals the parameters that influence the different responses of the soft and hard impact models. Since $\mu $ is proportional to $c$, which is calibrated in the soft impact model to reproduce the dynamics of the hard impact model for large stiffness, then $\kappa $ is the only remaining parameter through which we can track these differences. By definition, $\kappa \u221d\omega 0/\omega $ compares the time scale of the membrane dynamics with that of the excitation. If $\kappa \u226b1$, then the contact time $\Delta \tau j,c$ for the ball with the membrane is very short, and the soft impact model mimics the hard impact case. In contrast, when $\kappa =O(1)$, then $\Delta \tau j,c$ is longer, and the dynamics of the soft impact model exhibits substantial differences from the hard impact model. Here, the non-dimensionalization allows us to track these differences, even in the case of a membrane with large stiffness (often assumed as large $k$), where intuition would suggest similar behavior for both models. Equation (38) indicates that for large stiffness, the soft model dynamics can diverge from that of the hard one when the excitation frequency is high enough so that $\kappa =O(1)$. Accordingly, we define the system with $\kappa \u226b1$ as stiff, corresponding to the scenario where the soft impact model closely approximates the hard impact results.

Based on this preliminary analysis of the non-dimensionalized models, we numerically explore the system dynamics by integrating Eq. (38) using a fourth order Runge–Kutta method. The parameters of the model are taken from Refs. 36 and 9, with the values listed in Table II unless stated otherwise. Bifurcation diagrams in terms of the capsule length $d$ illustrate the system response. They are constructed using a continuation-like approach: for each value of $d$, we discard the first $375$ excitation periods to avoid transient behavior, then graph the impact velocities vs $d$ for the following $25$ excitation periods. This provides the impact velocities for the attracting dynamical behavior over a sufficiently long time. The result also gives the initial condition for the subsequent value of $d$, decreased by $0.002$, and the procedure is repeated. In this case, decreasing $d$ is achieved by decreasing $s$ but a similar result can be obtained by increasing $A$.

Symbol . | Value . | Unit . | Symbol . | Value . | Unit . | Symbol . | Value . |
---|---|---|---|---|---|---|---|

m_{b} | 0.0035 | kg | k | 560.0 | N/m | κ | 400.0 |

A | 1.0 | m/s^{2} | β | 0.0 | rad | γ | 0 |

g | 9.8 | m/s^{2} | r | 0.5 | … | μ | 172.363 |

s | 0.3 | m | ω | π | rad/s | d | 0.3 |

Symbol . | Value . | Unit . | Symbol . | Value . | Unit . | Symbol . | Value . |
---|---|---|---|---|---|---|---|

m_{b} | 0.0035 | kg | k | 560.0 | N/m | κ | 400.0 |

A | 1.0 | m/s^{2} | β | 0.0 | rad | γ | 0 |

g | 9.8 | m/s^{2} | r | 0.5 | … | μ | 172.363 |

s | 0.3 | m | ω | π | rad/s | d | 0.3 |

#### 1. Influence of effective size *d* in the hard impact model

Before comparing the soft and hard impact models, we first analyze the computational results for the bifurcations in the hard impact model in terms of the non-dimensionalized capsule size $d$. These results illustrate the different behaviors appearing for large (infinite) stiffness of the membranes, which are of interest when considering variable stiffness in the soft impact model below. Subsequently, we consider the variation of $d$ vs $\kappa $, via variation of the natural frequency $\omega 0$ of the membrane.

Several different dynamical behaviors are identified in the bifurcation diagrams for the hard impact model in Figs. 3(a) and 3(b), as illustrated by the corresponding phase planes in panels 3(c)–3(f). The relative impact velocities $x\u02d9$ for the hard impact case are shown in Fig. 3(a), while the bifurcation diagram constructed with the Poincar$e\xb4$ section of $x\u02d9$ based on the excitation phase $\omega \tau j=0$ is shown in Fig. 3(b). As we are considering the case $\beta =0$, a periodic orbit or aperiodic behavior that is not symmetric with respect to the transformation $x\u2192\u2212x$ [shown in black in panels (a) and (b)] co-exists with its reflected counterpart, shifted by a phase difference of $\pi $ rad (shown in purple). For longer capsule length with a larger distance separating the membranes, $0.3>d>0.186$, the attracting system behavior is a 1:1 periodic orbit (alternating impacts on left and right ends of the capsule) depicted in the phase plane of Fig. 3(f). For $d=0.186$, there is a grazing incidence of the inner loop of the 1:1 periodic orbit creating a 2:4/2$T$ stable periodic orbit [Fig. 3(e)]. Co-existing with the stable orbit is a 1:2 unstable periodic orbit, depicted in red on the diagram and on the phase planes of Fig. 3(e). This orbit and its instability are obtained semi-analytically via the optimization method (28) and (29) based on the maps $Mj$. The instability of this orbit underlies the influence of grazing, which then pushes the system to a stable orbit with both a longer period and an increased number of impacts. For smaller values of $d$, there is a series of period doubling bifurcations that is quenched by the grazing incidence of one of the unstable periodic orbit inner loops, generating a transition to a 2:3/2$T$ stable orbit, shown in Fig. 3(d). The same behavior continues down to $d=0.124$ where further decreases of $d$ yield aperiodic, apparently chaotic behavior, with some periodic windows containing high period orbits. Figure 3(c) shows the values in the phase plane corresponding to the Poincaré section of zero phase difference from panel (b) for $d$ in one of the chaotic windows.

Figure 3(g) shows the energy $Ej$ for each impact and the average energy generated $P$ for orbit (36), thus quantifying the energy harvesting performance of the hard impact system at different $d$. In general, there are two different types of sequences of impact velocities corresponding to different ranges of $P$. Periodic 1:1 orbits with large relative impact velocities $x\u02d9j$, corresponding to alternating maps $M1$ and $M2$ in Fig. 2, generate energy in the range of 3 mJ to 4 mJ of energy per impact. In contrast, those periodic orbits or attractors with a combination of high and low values of $x\u02d9j$ generate lower $P$. This second type of behavior corresponds to 1:$n$ periodic orbits with $n>1$, which include subsequent impacts on the same side of the capsule, corresponding to maps $M3$ and $M4$ in Fig. 2. For example, 1:2 periodic motion is defined by the composition of maps $M1\xb0M2\xb0M3$, with $M3$ yielding a low energy output below 0.3 mJ. Then, the 1:1 periodic solution in Fig. 3(e) yields a higher power output $P$ than for other 1:$n$ orbits, noting also that $P$ increases with $d$. It may seem surprising that there are no drastic jumps in $P$ that might be expected in transitions from 1:1 to 1:$n$ behavior for $n>1$, with only the slope of $P$ changing with these types of transitions. This is partly due to the choice of parameters and reflects the choice of $P$ as an averaged energy output over time, in contrast to averaging over impacts. For $n>1$, the additional impacts per period in 1:$n$ orbits tend to have low impact velocity, so that transitions to 1:$n$ orbits typically yield smaller increases in energy per time interval, as compared with the decreases in energy per impact.^{37}

The difference in the relative impact velocities shown in (a) for hard impacts is closely related to $\varphi j$ as shown in Fig. 3(g), which is the excitation phase at impact that captures the phase shift of the impact relative to the external oscillation of the apparatus. Larger phase shifts result in larger $x\u02d9j$, particularly when the velocities of the ball $Zj\u2032$ and the capsule $Yj\u2032$ have opposite sign. In contrast, subsequent impacts on the same side occur when the ball and capsule are close in phase, yielding small $x\u02d9j$. The time series depicted in Fig. 4 illustrate this contrast, where panel (a) shows that period 1:1 orbits are usually phase locked, having a fixed phase shift with $sgn\u2061(Yj\u2032)=\u2212sgn\u2061(Zj\u2032)$. In contrast, Fig. 4(b) illustrates 2:1 periodic orbit, where a reduced phase shift of the first impact on the bottom yields lower $x\u02d9j$, followed by a second impact on the bottom with $sgn\u2061(Yj\u2032)=sgn\u2061(Zj\u2032)$, yielding small $x\u02d9j$.

#### 2. Influence of effective size *d* vs stiffness *κ* in the soft impact model

*κ*

We analyze the results from the soft impact model, first comparing the stiff case $\kappa =\omega 0\pi /\omega \u226b1$ to the hard impact model. Then, we consider smaller values of $\omega 0$ while fixing the excitation frequency at $\omega =\pi $ as used in the calibration of the soft model (see Table II).

Figure 5 compares the results for the hard and soft impact models for $\kappa =400$ ($k=560$ N/m). By comparing both bifurcations diagrams and phase planes, it is clear that the models agree well with each other. The differences between the bifurcation diagrams of Fig. 5(b) are related to co-existing reflected orbits for $\beta =0$, with the soft impact model in close agreement to the left/right reflection of the hard impact model results. Given the sensitivity near bifurcation points (in this case a bifurcation from 1:1 to either 2:1 or 1:2 behavior), small changes in the parameters can result in transfers to reflected orbits. This can be verified by comparing Figs. 5(c)–5(f) with Figs. 3(c)–3(f), showing the matching periodic orbits for the two models. We emphasize that the matching orbits were obtained by simulating the system with different ICs so that the results correspond to the attracting system behavior. Another difference in the soft impact model is the shift in the range of $d$ values where periodic windows between chaotic behavior are observed. This is expected as the narrow periodic windows of hard impacts are sensitive to parameter or model perturbations. As shown in Fig. 5(g), the soft impact model yields a narrow range of $\rho $ ratios instead of a fixed $r$ as given in the hard model, with the distribution of the $\rho $ values concentrated at the reference value $r=0.5$, as follows from the calibration of $c$ for large $\kappa $. The spread in $\rho $ is due to small differences in the excitation phase $\varphi j,c$ at each impact, due to the (small) time interval $\Delta \tau j,c$ during which the ball is in contact with the membrane as shown in Fig. 5(h). The resulting differences in the phase $\varphi j,c$ at impact can be observed by comparing Fig. 5(j) with Fig. 3(h). The energy output at each impact $Ej$ and average energy output $P$ are also the same for both models, as shown by comparing Figs. 5(i) and 3(g). The close comparisons for the soft and hard impact models across all performance variables in the case of $\kappa =400$ demonstrate that the calibration based on Eqs. (14)–(19) successfully captures the dynamics of the hard impacts for large membrane stiffness.

Once the calibration is established at large $\kappa =400$, it is possible to explore other parameters ranges for the soft impact model, to identify differences with hard impacts. We consider different values of $\omega 0$ (and thus $\kappa $), with $\omega $ fixed. (The influence of different values of $\omega $ is studied in Sec. III B.) Figure 6 shows different bifurcation diagrams with $\kappa $ as the bifurcation parameter at different values of $d$, with other parameters as in Table II. Large values of $\kappa $ indicate large stiffness, for which the soft impact model results agree with the hard impact case indicated by dashed red lines. Smaller values of $\kappa $ indicate a soft membrane, for which the behavior naturally differs from that of the hard impact model. The diagrams are shown for $\kappa $ on a $log$ scale in order to cover a sufficiently large parameter range. Figure 6(a) depicts the results for $d=0.2$, where a 1:1 orbit is maintained for most values of $\kappa $ until a rapid decrease to near zero impact velocities at $\kappa =4.04$. For these lower values of $\kappa $, we identify co-existence of different stable attractors. The bi-stable orbits at $\kappa =3$, shown in the phase plane in Fig. 7(a), are a 1:1 orbit with low impact velocities, shown in blue, and a 2:2/3$T$ orbit with higher impact velocities, shown in black. The bi-stability indicates that for lower values of $\kappa $, the nonlinear effects due to impacts can have little influence on the system behavior via the low energy orbits. Yet, for the same parameters, the nonlinear effects still play a significant role in the overall system dynamics, since co-existing attractors can be accessed at higher energy levels. We note that the optimization method (28) and (29) is helpful in identification of the bi-stable orbits since they may not be observed numerically in the absence of the necessary ICs.

Figure 6(b) shows the results for $d=0.18$, with a 2:4/2$T$ periodic solution for large $\kappa $ and transitions to 1:2 and then to a 1:1 periodic solution at $\kappa =28.57$. Complementary phase planes are shown in Fig. 7(b) for a sequence of periodic orbits corresponding to Fig. 6(b). From these phase planes, we see that as $\kappa $ decreases the penetration of the ball into the membrane increases. Similarly, for the case of $d=0.16$, shown in Fig. 6(c), as $\kappa $ decreases the system undergoes transitions from a 1:2 periodic orbit to 1:1 behavior at $\kappa =57.14$. For lower values of $\kappa $ with $d=0.16$, the 1:1 impact velocities increase to a maximum at $\kappa =3.84$, followed by a rapid decrease of the impact velocities for smaller $\kappa $. Figure 7(c) shows that at higher impact velocities, the periodic orbits are nonlinear and stay in contact with the membrane for longer times. However, as $\kappa $ decreases further, so does the influence of the impacts, as the periodic orbits are then ellipsoidal with small amplitudes corresponding to a weaker influence of the nonlinearities on the dynamics. The nearly linearly behavior coincides with increased softness of the membrane, which in turn yields low impact velocities in which the influence of the impact nonlinearity is reduced. This behavior extends across a range of parameters, as verified in the diagrams studied throughout this section. Finally, Fig. 6(d) shows the case where chaos is present for large $\kappa =400$, with its attractor seen in Fig. 5(c). As the value of $\kappa $ decreases from $\kappa =400$, the chaotic behavior eventually transitions to a 1:2 orbit at $\kappa =36.36$ and then to a 1:1 solution at $\kappa =12.12$.

Considering these examples together, we observe the following trends: in parameter regimes where the stiff system exhibits 1:$n$ periodic behavior for larger $n$, reducing the stiffness via decreasing $\kappa $ has a noticeable impact, leading to 1:$n$ behavior for reduced $n$ and eventually 1:1 behavior. Since bifurcations to 1:$n$ behavior for larger $n$ occur for smaller values of $d$, it is in this range of parameters where we see the greater influence of reduced $\kappa $. This trend follows intuitively from the fact that the differences between the soft and hard impact models are in the modeling of the impacts so that behavior with a larger number of impacts per period are more sensitive to variations in $\kappa $.

The influence of $\kappa $ has implications for energy generation. Decreasing $\kappa $ leads to a period-$T$ behavior with only high velocity impacts. These impacts have the largest contribution at an optimum value of $\kappa $ for a given $d$. However, the impact velocity decreases substantially for $\kappa $ below this optimum, with the possibility for abrupt changes in the system behavior as shown in Fig. 6(a). As discussed above, this behavior follows from the softness of the membrane and low impact velocities that reduce the influence of the impact nonlinearity. While the orbit becomes more and more linear as $\kappa $ decreases past its optimum value, it is important to notice that nonlinear effects persist for smaller $\kappa $, apparent through the other attracting orbits that coexist with the near-linear response as observed for $d=0.2$ and $\kappa =3$ in Fig. 7(a). From a design perspective, low values of $\kappa $ can be beneficial for energy extraction, but the system must be designed in such a way to ensure $\kappa $ does not decrease past the value corresponding to optimal impact velocity or be perturbed to any co-existing low impact velocity behavior.

Before further exploring the influence of $\kappa $ on the bifurcation structure, we note that Figs. 3–7 all show the prevalence of symmetry breaking bifurcations. Specifically, for larger $d$, e.g., $d=0.275$ in Fig. 5(a), the stable 1:1 orbits are symmetric in the phase plane^{36} (not shown in Fig. 5), with equal magnitude of impact velocities $|X\u2032(\tau j|$ for left and right impacts. As $d$ decreases, there are symmetry breaking transitions from symmetric 1:1 orbits to asymmetric 1:1 orbits. The asymmetry in the 1:1 behavior is clear from phase planes such as Fig. 3(f), 5(f), and 7)(b), with increased asymmetry in the trajectories between impacts leading to grazing and transitions to 1:$n$ or $n$:1 behavior. This type of symmetry breaking bifurcation has been systematically explored in the hard impact model,^{36–38,40} occurring for decreasing $d$ as also observed here for both hard and soft impacts. It can appear for any value of inclination angle $\beta $, including $\beta =0$. As would be expected, for $\beta >0$, this type of asymmetry is enhanced throughout the dynamics.^{38} Figures 6 and 7 indicate that this bifurcation from symmetric to asymmetric 1:1 orbits also occurs for increasing $\kappa $ as well as for different combinations of $r$ and $\omega $ as apparent from Figs. 11 and 16. As we see below, the soft impacts allow a second type of symmetry breaking bifurcation, associated with the dynamics over the finite contact time $\Delta \tau j,c$. Thus, this second type of asymmetric bifurcation cannot be observed in a model with instantaneous impacts.

Figure 8 illustrates the influence of $\kappa $ on the bifurcations with $d$ as the bifurcation parameter, comparing these for several values of $\kappa $. Figure 8(a) shows that as $\kappa $ decreases, the regions with 1:1 periodic orbits with larger impact velocities $u\u02d9$ are pushed to lower values of $d$, until the 1:1 behavior becomes the attractor for all $d$ as can be seen for $\kappa =12.65$. For smaller values of $\kappa $, there is a range of apparently chaotic behavior for larger values of $d$, with its phase plane dynamics shown in Fig. 11(a) for $d=0.2$. In Fig. 8(a), we observe the same patterns as in Fig. 6, noting that orbits with larger numbers of impacts per period, e.g., for smaller $d$ and larger $\kappa $, experience substantial qualitative changes as $\kappa $ decreases. The values of $\rho $ shown in Fig. 8(b) for each impact also yield valuable insights. For larger values of $\kappa $, $\rho $ varies slightly around its reference value $r$ from the corresponding hard impact model. However, as $\kappa $ decreases, the variability of values in $\rho $ increases, particularly for chaotic ranges. It is important to emphasize that, for softer membranes and lower impact velocities, there are values of $\rho >1$. This does not violate any energy principle as the greater contact time $\Delta \tau j,c$ in these cases enables energy to be transferred from the capsule to the ball, which is not possible for hard impact models. Also, another type of symmetry breaking orbit can appear in the soft impact model, characterized by asymmetric values of $\rho $ that diverge from the reference value $r$ for the hard model. This behavior is observed for $\kappa =12.65$ and $d=0.142$, where a symmetry breaking bifurcation yields two branches in the diagram for $\rho $ vs $d$, highlighted by a black circle in Fig. 8(b). There is also a slight asymmetry in the contact times $\Delta \tau j,c$ on the left and right ends. Note that this symmetry breaking bifurcation follows from the dynamics during the finite contact, rather than from the dynamics in between impacts and, thus, is not observed in the hard impact model. The asymmetry in $\rho $ follows from the behavior of the excitation phase $\varphi j,c$, which has slight variation around $\varphi j$ from the hard impact model, diverging from it only when there is qualitatively different behavior for the soft impacts. Figure 8(c) indicates that this change in $\varphi j,c$ corresponds to the symmetry breaking point of $\rho $, shown in panel (b). The asymmetric 1:1 periodic orbit is shown in Fig. 11(b) for $d=0.12$ and has a signature in the bifurcation structure vs $\kappa $ shown in Fig. 6(d). As $d$ decreases further, these branches show additional bifurcations, with jumps in impact phase, velocity, and energy. For larger values of $\rho $, the impact velocities can be even higher for small $\kappa $, for which the system appears to be nearly optimized in terms of both the contact time $\Delta tj,c$ and the corresponding excitation phase $\varphi j,c$, yields the ideal conditions for energy harvesting.

Finally, Fig. 8(d) shows $\Delta tj,c$ (non-dimensionalized $\Delta tj,c$), the time that the ball stays in contact with the membrane at the $jth$ impact. This value has substantial variation when the system dynamics is aperiodic or chaotic for the softest membrane shown with $\kappa =4$. In contrast, for periodic behavior, $\Delta \tau j,c$ is small for larger $\kappa $ steadily increasing as $\kappa $ decreases. This can be seen also from the phase planes in Fig. 11(d) for $d=0.2$, showing the sequence of 1:1 orbits corresponding to the $\kappa $ values from the bifurcation diagrams in Fig. 8(a). This sequence illustrates that the penetration into the impact region of the orbit, and consequently $\Delta \tau j,c$, increases as $\kappa $ decreases, while the system is in the nonlinear regime.

After exploring the influence of $\kappa $ on the bifurcation diagram vs $d$, we investigate also the influence of the reference restitution coefficient $r$. Figures 9 and 10 show bifurcations vs the parameter $d$ for small $r=0.1$ and large $r=0.95$, respectively. In each case, we recall that the calibration of the soft impact model is applied for $mb=0.0035$ kg and $\kappa =400$ so that Eq. (19) gives $c$ using the respective values of $r$. Then, the soft impact model closely reproduces the behavior of the hard impact model for the largest value of $\kappa $ shown, which is the basis for comparison as we vary $\kappa $ ($\omega 0$) and thus also $c$ via Eq. (19).

For small and large values of reference value $r$, we observe similar influences of $\kappa $ as previously described for $r=0.5$. As $\kappa $ decreases, we observe the same shift of the bifurcations of 1:1 periodic behavior to a range of smaller values of $d$. The diagrams show that the 1:1 behavior is dominant over much of the parameter range shown, except for small $d$ where chaotic behavior dominates for $r=0.95$ (see Fig. 10), and for small $d$, large $\kappa $, and $r=0.1$, where 1:$n$ orbits survive. One prominent difference for $r=0.1$ is shown in Fig. 9(b), where we observe symmetry breaking in the values of $\rho $ for larger values of $\kappa $ as compared with $r=0.5$. This observation follows intuitively, since stronger energy dissipation in such systems facilitates stable asymmetric orbits, consequently breaking the symmetry of $\varphi j,c$ which in term yields different values of $\rho $ at each impact. We also observe bifurcations to asymmetric 1:1 behavior and asymmetric values of $\rho $, for reference $r=0.95$ and $\kappa =12.65$. For smaller $d$ on this branch for larger $r$, the symmetry breaking leads to chaotic behavior (see Fig. 10). A second difference is the chaotic behavior that is observed for $r=0.95$ over a larger range of $\kappa $ for small $d$, as shown, e.g., in Fig. 6(d). This behavior is consistent with the chaotic behavior observed for large $r$ and small $d$ over wide parameter ranges in previous studies of the hard impact model,^{36} suggesting that smaller values of $\kappa $ are required to induce periodic behavior for larger $\rho $ and reference $r$.

For $r=0.1$, small values of $\kappa $ and larger values of $d$ result in nearly non-impacting orbits, where very soft impacts yield relatively slow motion of the ball. Then, the ball rarely impacts the larger capsule ends, as indicated by the sparse relative impact velocities in Fig. 9(a) captured only sporadically within the time window used computationally. In contrast, for smaller $d$ and small $\kappa $, the ball regularly impacts the membrane in its 1:1 orbit, collecting energy from the capsule motion via a shift in $\varphi j,c$, the effective restitution coefficient $\rho $, contact time $\Delta \tau j,c$, and the energy output increase.

### B. Excitation frequency vs natural frequency of the membrane

In this section, we consider the influence of the membrane stiffness $k=mb\omega 02$ on bifurcations in terms of excitation frequency $\omega $ so that the non-dimensionalization does not vary with the bifurcation parameter $\omega $, and we choose the capsule length $s$ and the natural frequency of the polymer $\omega 0$ as the reference parameters. This choice yields displacements $\upsilon =U/s$ for the soft model and $\chi =x/s$ for the hard model and time variable $t\u2217=\tau \omega 0$. Then, the non-dimensionalized equations of motion are

where $g\u2217=gsin\u2061(\beta )s\omega 02$, $A\u2217=As\omega 02$, $\omega \u2217=\omega \omega 0$, and $c\u2217=cmb\omega 0$. Note that the parameter $\omega \u2217=\pi /\kappa $, recalling $\kappa $ as the non-dimensional parameter from Sec. III A that characterizes the soft membrane. Then, an increase of $\omega \u2217$ has the same influence as reducing $\kappa $, namely, that the soft impact dynamics diverge from the hard impact dynamics. Since we are interested in the influence of $\omega 0$ on the variability of the bifurcations in terms of $\omega $, in this section, we show the results with $\omega $ and $k=m\omega 02$ on separate axes, rather than in terms of $\omega \u2217$.

As in Sec. III A, we calibrate the soft impact model with the hard impact model for $k$ large (and $\omega 0$ large) so that the corresponding bifurcation diagram serves as a reference for other values of $\omega 0$. Figure 12(a) shows the bifurcation diagram for the dimensional relative impact velocities $U\u2032(\tau j,c)$ (magenta) and $X\u2032(\tau j)$ (black) vs $\omega $, while Fig. 12(b) shows the bifurcation diagram based on the Poincar$e\xb4$ section of $U$ and $X$ when the excitation phase is zero. The bifurcation results show several windows of periodic behavior and chaos for different values of $\omega $, illustrated in the corresponding phase planes. For low values of excitation frequency, the system presents a period-$T$ solution with multiple impacts, including both higher and lower relative impact velocities [panel (c)]. As $\omega $ increases, this orbit undergoes several grazing bifurcations until a chaotic region is reached. Beyond the chaotic response, a 1:1 periodic orbit emerges at $\omega =3.12$ rad/s [panel (d)]. This initiates a window of 1:1 behavior, which is the largest periodic window for this parameter combination. At $\omega =6.89$ rad/s, another chaotic behavior appears sharply, continuing up to $\omega =9.97$ rad/s [panel (e)]. This complex behavior gives way to a period-3$T$ solution [panel (f)], following a narrow sequence of periodic orbits with a successively shorter period. Finally, at $\omega =11.93$rad/s, the system is again chaotic, which continues for larger values of $\omega $.

The curious sequence of a number of impacts per excitation period $T$, identified from the periodic windows in panels (a) and (b), has implications for the energy output as shown in Fig. 12(g). For small values of $\omega $, the system impacts each side of the capsule multiple times before the period of excitation $T$ is reached. In that case, the time scale associated with $\omega $ is long relative to that of the system dynamics, associated with $\omega 0$. Since both large and small impact velocities $Uj,c\u2032$ with large and small $Ej$ are included in the orbit, the average power $P$ remains relatively high but fluctuates with $\omega $. For slightly larger $\omega $, there is a mismatch between the period of the excitation and the system dynamics, yielding a chaotic response, dominated by low velocity impacts and thus low energy output. Increasing $\omega $ further, the excitation frequency is large enough to match the system dynamics, yielding period-$T$ 1:1 behavior over a large window of $\omega $ values. The 1:1 orbits provide the most consistent large energy output for $Ej$ and $P$. For $\omega $ increased further, the system experiences several periods of excitation for each impact, leading to orbits whose periods are composed of multiple excitation periods. Not only is the energy per impact $Ej$ lower than for the 1:1 orbits, but the time-averaged power $P$ is even lower, given the less frequent impacts. The periodic windows for these orbits are interspersed with windows of chaotic behavior, again with low energy output.

Comparing the results for both models, we observe complete agreement between soft and hard impacts for $k=560$ N/m, thus establishing a good reference scenario for exploring other values of $k$ and, thus, $\omega 0$. Likewise, the success of the calibration is observed from the values of $\rho $ for each impact, highly clustered around the reference value of $r=0.5$ as shown Fig. 12(i). The slight variation in $\rho $ around the reference value $r=0.5$ appears to decrease with increasing $\omega $ for the periodic responses, while in chaotic regions, it spreads out due to the different impact phases for these dynamics, as shown in Fig. 12(h). The same type of variation can be seen for $\Delta \tau j,c$ in Fig. 12(j), where in the chaotic regions it is larger for some values than the reference values.

Following the calibration of the soft model with the hard model at $k=560$ N/m, we consider the bifurcation diagrams as $k$ varies. Figures 13–15 show the $Uj\u2032$ vs$\omega $ bifurcation diagrams for several different values of $k$ and three different reference values $r=0.5,0.1,0.95$. For the mid-range reference value $r=0.5$, we see a similar sequence of bifurcations as in Fig. 12 for all but the smallest value of $k$ (softest membrane). Multi-impact T-period orbits are observed for small $\omega $, and larger windows of 1:1 orbits are followed by a transition to chaos for intermediate $\omega $. Periodic windows with 1:1/$nT$ behavior for $n>1$, occur for larger $\omega $, also shown in Fig. 16(b), with chaos dominating as $\omega $ increases further. For the smallest value of $\omega 0=4$ rad/s, i.e., the softest membrane considered, 1:1 orbits are observed for a limited range of smaller $\omega $, which exhibit larger average energy output $P$. For this value of $\omega 0$ and larger values of $\omega $, we again observe behavior with sporadic or no impacts over multiple excitation periods, as also observed in Fig. 9 for larger $d$. Figure 16(a) shows the 1:1 behavior in the phase plane for the values of $k$ as in Fig. 13 and $\omega =\pi $, in agreement with Sec. III A.

Looking to the other measures, panel (b) shows a spread of $\rho $ in chaotic regions or high periodic asymmetric behavior, otherwise concentrated near the reference value $r=0.5$. As seen previously in Fig. 10, for low impact velocities in the chaotic regime, there are cases where $\rho >1$, which indicates that the membrane gives energy to the ball as the interaction between the capsule and the ball depends on the impact phase and the time spent in contact. Furthermore, there is a slight decrease in $\rho $ for smaller values of $k=0.56$ and $0.056$ N/m and low $\omega $. In panel (d), larger $\Delta \tau j,c$ is observed over a larger range of $\omega $ as $k$ decreases, as observed in Sec. III A for softer membranes. In panels (e) and (f), increases in (average) impact energy $Ej$ ($P$) are consistent with increased $Uj,c\u2032$ and increased $\Delta \tau j,c$. In general, energy produced by the membrane increases as $k$ decreases up to $k=0.56$ N/m, due to the change in excitation phase $\varphi j,c$. Increased $Ej$ and $P$ for the softest membrane $k=0.056$ N/m is limited to low $\omega $, corresponding to 1:1 behavior. Chaotic behavior observed for larger $k$ tends to reduce mean energy output, as it includes impacts with small velocities.

For the small reference value $r=0.1$, decreasing $k$ limits the chaotic regions to higher $\omega $, with periodic behavior dominating over much of the range of $\omega $ shown. For larger and medium values of $k$, this includes $T$-periodic orbits with multiple impacts per period including low velocity impacts with $Uj,c\u2032$ near zero. The range of this behavior increases with decreasing $k$, extending over the entire range of $\omega $ for $k=5.6$ N/m. This asymmetric, multiple impact behavior is captured by the multiple branches of small $\rho $ in panel (b), and a branch of $\Delta \tau j,c$ in panel (d) that increases with decreasing $\omega $. For smaller $k=0.56$ N/m, this behavior is replaced by 1:1 orbits, reflected in $\rho $ near the reference value $r=0.1$, and $\Delta \tau j,c$ increasing with decreasing $k$ for softer membranes. For the softest membrane $k=0.056$ N/m, as seen for $r=0.5$, the 1:1 behavior is limited to low values of $\omega $, once again with larger values of $Ej$ and $P$. As $\omega $ increases in this case, $\Delta \tau c,j$ decreases and $\rho $ increases, corresponding to nearly linear behavior illustrated in Fig. 16(c). This behavior is analogous to observations from Fig. 6, where an optimal orbit is reached for a small value of $k$, below which the system exhibits nearly linear behavior. Here, the membrane has less time to dissipate the energy of the ball, yielding larger $\rho $, but limited $\Delta \tau c,j$ results in near-linear dynamics and reduced energy output. Higher $\omega $ in this softest case again yields sporadic or no impact behavior over multiple excitation periods.

Considering larger $r=0.95$, for larger values of $\omega 0>4$, the bifurcation sequence is relatively simple, with chaotic behavior for lower $\omega $ and 1:1 behavior or 1:1/$nT$, $n>1$, for higher $\omega $. For periodic orbits, $\rho $ is clustered near the reference value of $r$ and $\Delta \tau j,c$ is small. As $\omega 0$ decreases, we again see the influence of high $\omega $, for which the periodic behavior shifts to 1:1/$nT$ orbits with smaller impact velocities over multiple excitation periods. Figure 16(d) highlights this transition via phase plane comparisons. For the smallest $\omega 0=4$, there is an increased range of chaotic behavior with low impact velocities for high $\omega $. The behavior of $\rho $ and $\Delta \tau j,c$ track these transitions as observed above in Fig. 12. Likewise the largest energy output is observed for 1:1 behavior, increasing with the softer membranes. However, the range of this behavior is limited to lower $\omega $ for softer membranes, which exhibit low velocity impacts and reduced energy output for higher $\omega $.

## IV. CONCLUSIONS

We develop a soft impact model for the vibro-impact energy harvesting pair, previously studied in the context of hard (instantaneous) impacts. In the model, a ball moves freely along the axis of an externally excited capsule, with dielectric polymer membranes at each end. Energy is generated as the ball impacts and deforms the membranes, converting mechanical to electrical energy^{9} using the variable capacitance principle. While the hard impact model is based on instantaneous impacts, the soft impact model includes additional equations for the motion when the ball is in contact with the membranes, modeled as a linear spring.^{41,42}

In order to systematically compare its dynamics with that of the corresponding hard impact model, we establish a consistent calibration procedure for the soft impact model in the setting of large elastic coefficient $k$ for the membrane. This procedure matches the reference restitution coefficient $r$ of the hard impact model to the ratio $\rho $ of relative velocities before and after a soft impact with the membrane, yielding a relationship for the linear dissipation coefficient $c$ in terms of $r$, $k$, and the mass of the ball. This relationship is then used to explore the dynamics of soft impacts over a range of values of $k$, including small $k$ corresponding to soft impacts on compliant membranes with low natural frequencies $\omega 0$.

While the bulk of the comparisons of soft and hard impact dynamics is carried out numerically, we also provide analytical descriptions, based on maps of the dynamics in-between impacts for the hard model, supplemented by maps in the soft impact model that describe the dynamics when the ball is in contact with the membrane. The analytical results provide complementary information about influential unstable and bi-stable behaviors that can be missed within numerical simulations. We capture the variation in dynamics and bifurcations in terms of several different parameters: device length, stiffness of the impact boundary, excitation frequency, and reference restitution coefficient.

We explore two main bifurcation scenarios motivated by two different non-dimensionalizations. The first considers effective capsule size $d$, proportional to capsule length $s$ and inversely proportional to the excitation amplitude $A$, with fixed excitation frequency $\omega $. The second considers variable excitation frequency $\omega $ and fixed capsule length. We studied these two scenarios for different reference restitution coefficients, $r=0.1,0.5,0.95$, comparing changes in the bifurcation sequences for different natural frequencies $\omega 0$ of the soft membrane.

For variable $d$ and fixed $\omega $, softer membranes tend to promote regular 1:1 orbits (alternating left/right impacts) and other periodic behavior, shifting grazing and period doubling bifurcations with more complex behavior to smaller values of $d$. This is particularly apparent for smaller $r$. Within this scenario, we found that typically the optimal energy outputs are achieved for soft membranes, with the output above that predicted by the hard impact model. However, for stiffness parameters slightly below this optimum, impact velocity and energy outputs drop sharply, yielding low output, sometimes with co-existing complex behavior. This comparison suggests that soft impacts may be desirable, but they are also sensitive to parameter fluctuations near optimal energy output. In particular, for small $\kappa $, we observed substantial ranges of $d$ for which the orbits exhibited sporadic or weak impacts or no impacts at all. This absence of energy generation occurred for a range of smaller reference values $r$ (corresponding $\rho $ in the soft impact model), which increased with excitation frequency $\omega $.

The bifurcation sequence for variable excitation frequency $\omega $ shows several types of periodic behavior in between windows of chaotic behavior. Multi-impact $T$-periodic behavior appears for lower $\omega $, 1:1 orbits appear for mid-range $\omega $, and $pT$ $p>1$ periodic orbits with few low impact velocity impacts occur for larger $\omega $. The largest energy output generally occurs for 1:1 orbits, since the multi-impact orbits and chaotic behavior include low velocity impacts, while the $pT$ periodic orbits occur less frequently. Decreasing the elastic coefficient $k$ and, thus, $\omega 0$, tends to expand the range of 1:1 periodic behavior for all $r$ (to which $\rho $ is calibrated) and, thus, increase the energy output. However, for higher $\omega $ and softer membranes, other less desirable behaviors appear: low impact velocity orbits occur for larger $r$, and orbits with sporadic/infrequent impacts occur for mid-range and smaller $r$. The parameter range of multi-impact $T$-periodic orbits also expands for reduced $r$, and chaotic behavior is sustained for large $r$ and low $\omega $ in all but the softest membranes.

In comparing the dynamics of the soft impact model with that of the hard impacts, we found a number of unexpected results. The soft and hard impact models can disagree in the setting of stiff membranes with large $k$, where one might expect them to agree when the excitation frequency is sufficiently high. Then, for a given excitation frequency $\omega $, the calibration of the models requires a sufficiently large $k$. In addition, Figs. 6–10 and 12–15 demonstrate that the divergence of soft and hard impact model solutions also depends on the number of impacts of the behavior, the excitation frequency, and the device length. Therefore, it is not possible to pin-point a value of $\kappa $ where the hard and soft models diverge for all scenarios. The analysis using the non-dimensional parameter $\kappa \u221d\omega 0/\omega $ provides a useful guide, as in Sec. III A, indicating that the soft and hard models agree for sufficiently large $\kappa $, where “sufficiently large” depends on other parameters. Furthermore, in addition to symmetry breaking bifurcations related to the motion between impacts, also observed for hard impact models, we found additional symmetry breaking bifurcations driven by the dynamics while the ball is in contact with the membrane. These occur in the variable effective restitution coefficient values $\rho $ from the soft model, which are dependent on the excitation phase at impact, shifted relative to the excitation. The bifurcations produce asymmetric motions with different energy generation and stability properties that are not captured in the hard impact model with a constant $r$.

We focus here on the horizontal case $\beta =0$. While we have not explored the setting of the inclined energy harvester $\beta \u22600$, nevertheless preliminary results (not shown) show qualitatively similar bifurcation sequences and effects of soft impacts as observed here for $\beta =0$. As has been observed for hard impact models, $\beta \u22600$ removes symmetric phase plane behavior and tends to encourage symmetry breaking bifurcations as well as grazing bifurcations from the basic 1:1 alternating left-right periodic behavior, which generates higher energy outputs. However, given the complex influence of $\beta $ on the sequence of different continuous and non-smooth bifurcations,^{38} we leave a thorough investigation of the influence of $\beta $ within the soft impact model to future work.

From a practical point of view, the polymer-based membranes are rather soft and highly stretchable. This implies that the hard impact model, being easier to study analytically and implement numerically, may not correctly capture the device dynamics when the excitation frequency is not large relative to the natural frequency of the membranes. In this case, a soft impact model should be used to better predict the energy harvesting device power output. Moreover, in some cases (e.g., small $k$ and relatively high $\omega $), the hard model may provide unrealistic results, where it demonstrates consistent power output, in contrast to low or no power output within the soft impact model. Furthermore, the use of the hard impact model may lead to underestimated energy output values, coming from lower estimated impact velocity compared to the corresponding soft model value, as presented in Fig. 6.

## ACKNOWLEDGMENTS

The authors gratefully acknowledge partial funding for this work from NSF-CMMI (No. 2009270) and EPSRC (No. EP/V034391/1).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Dimitri Costa:** Formal analysis (equal); Investigation (equal); Methodology (equal); Validation (equal); Writing – original draft (lead). **Rachel Kuske:** Conceptualization (equal); Funding acquisition (lead); Supervision (lead); Writing – original draft (equal); Writing – review & editing (lead). **Daniil Yurchenko:** Conceptualization (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.