Qualitative changes in bifurcation structure for soft vs hard impact models of a vibro-impact energy harvester

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,⁴ or in rotors where impacts increase wear, leading to increased maintenance⁵ and sometimes causing system failure.⁶ 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,⁷ 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.¹³ These models can explain several behavioral patterns of an impacting system like chattering,¹⁵ grazing,¹³ and other nonlinear phenomena.¹⁶ 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.¹⁴ 

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.¹⁹ 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,²⁴ in the effectiveness of control in impacting systems,^25,26 and within the context of nonlinear bumpers.²⁷ 

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.²⁸ 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.²⁹ 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.³⁰ 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.³¹ compare different designs for an energy harvester including hard stops in the system. Likewise, Halim and Park⁸ 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,³² magnet induction harvesters,³³ and piezoelectric based^34,35 or electrostatic harvesters.⁹ 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.²⁸ 

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,³⁹ 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.

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 $β$ from the horizontal one. Here, we focus on the case $β=0$ and include some observations about the $β≠0$ in Sec. III.

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

where $τ$ is the time, $ω$ is the excitation frequency of the main structure, $A0$ is the amplitude of excitation in newtons, and $′$ represents the derivative with respect to time. Assuming that $M≫mb$⁠, 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−Z$⁠,

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

In the hard impact model, the membranes are non-compliant and the impact is treated as instantaneous at the time $τj$ of the $jth$ impact. Hence, the relationship between the instantaneous velocities immediately before [$X′(τ−)$] and after [$X′(τ+)$] 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 $τj<τ<τj+1$⁠, the motion is given by

where $ϕj=mod(ωτj,2π)$ is the excitation phase at $τj$ and $+(−)s/2$ corresponds to the last impact $X(τj)$ located on the right (left) side. Equations (5) and (7) are then used to determine $X′(τj+1−)$ and $τj+1$⁠.

Equations (6) and (7) with $τ=τ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∈(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 $τj$ and $τj+1$ are given by

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

where $Δj$ is the time interval between consecutive impacts. Since Eq. (9) is a transcendental equation for $τ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 $τj+1>τj$⁠. After $τj+1$ is obtained, Eq. (7) can be used to calculate $X′(τ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,…$ 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 $ℓ$⁠:$q$ for $ℓ≠n$ and/or $q≠m$⁠. Later, we also consider higher excitation frequencies $ω$ for which it is possible to have $pT$ periodic orbits with a low number of impacts over $pT$⁠.

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−Z$ 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 $τj,c$ (⁠$τ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→∞$⁠. 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 $τj,c≈τj,f$⁠, we find the value $c$ for which the ratio $ρ$ of velocities before and after contact with the membrane is comparable to $r=−X′(τj+)/X′(τj−)$⁠, that is,

We refer to this ratio $ρ$ 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′(τj,c−)$⁠, $Uc(τ)$ is obtained from Eq. (11),

where $ω1=ω02−(c/2mb)2$⁠. Note that for large $k$⁠, the natural frequency associated with the mass-membrane system, $ω0=k/mb$⁠, is also large. Then, the contact time $Δτj,c$ in Eq. (12) is small, that is,

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

To complete the calibration, we evaluate $Uc′(τ)$ in Eq. (17) at $τ=τj,f$ from Eq. (16), and divide by $U′(τj,c−)$ to get the expression for $ρ$⁠,

Then, we set $ρ=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−1/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≠0$⁠, the excitation phase influences the system dynamics and $Δτj,c$ is no longer negligible in these comparisons.

It is useful to obtain analytical expressions for the relative displacement and velocity over the contact interval $Δτ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(τ)$ and $Uc′(τ)$ 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 $ζ=c/2mb$⁠, $ω12=ω02−ζ2$ and the signs $∓$ and $±$ 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 $ϕj,c=mod(ω1τj,c,2π)$⁠. The transition time $τj,f$ is obtained by setting $Uc(τ)=∓s/2$ and $τj+1,c$ can be determined from the equations for $Mi$ for $i=1,…,4$⁠, analogous to the procedure for hard impacts. Then, we have the set of equations,

where $Δτj,f=τj+1,c−τj,f$⁠, and $ϕj,f=mod(ω1τj,f,2π)$⁠.

It is useful to recall that the maps $Mi$ for $i=1,…,4$ are defined in the same way as in the hard impact model away from impact, with relative velocity continuous across $τj,c$ and $τj,f$ as in Eq. (13). After determining $Δτj,c(Δτ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′(τj,f)$ and the velocity of the next impact $U′(τ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′(τj,c)=Uc′(τj+l,c)$⁠, $Uc(τj,c)=Uc(τj+l,c)$ and $mod(τj,c,2π/ω)=mod(τj+l,c,2π/ω)$ 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.³⁶ 

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.

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 $δ=|x|−s/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 $ϵM$ 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(δmax,j)$ is the maximum capacitance of the polymer reached at the maximum penetration of the ball $δ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 $δmax,j$⁠, appearing in the definition of energy $Ej$ via $As$ [Eq. (30)], is given by the displacement $Uc$ when $Uc′=0$ in the soft model and, thus, is calculated directly from Eq. (21). In contrast, the hard model does not include an expression for $δ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 $δ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 $τ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 $τ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 $τ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.

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.

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 $ω0$ of the membrane, and a second focused on the variation of excitation frequency $ω$ vs natural frequency $ω0$⁠.

First, we define a non-dimensional parameter $d=sω2Aπ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ω2Aπ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ω2Aπ2$⁠, which satisfies

Here, $κ=ω0πω$ is proportional to the ratio of natural frequency $ω0$ and the excitation frequency $ω$⁠, and $˙$ indicates the derivative with respect to the non-dimensional time $t=τωπ$⁠. There are also non-dimensionalized parameters $γ=gsin⁡(β)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 $μ=cπmω$⁠.

An analysis of Eqs. (38) and (37) reveals the parameters that influence the different responses of the soft and hard impact models. Since $μ$ 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 $κ$ is the only remaining parameter through which we can track these differences. By definition, $κ∝ω0/ω$ compares the time scale of the membrane dynamics with that of the excitation. If $κ≫1$⁠, then the contact time $Δτj,c$ for the ball with the membrane is very short, and the soft impact model mimics the hard impact case. In contrast, when $κ=O(1)$⁠, then $Δτ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 $κ=O(1)$⁠. Accordingly, we define the system with $κ≫1$ 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$⁠.

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 $κ$⁠, via variation of the natural frequency $ω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˙$ for the hard impact case are shown in Fig. 3(a), while the bifurcation diagram constructed with the Poincar$e´$ section of $x˙$ based on the excitation phase $ωτj=0$ is shown in Fig. 3(b). As we are considering the case $β=0$⁠, a periodic orbit or aperiodic behavior that is not symmetric with respect to the transformation $x→−x$ [shown in black in panels (a) and (b)] co-exists with its reflected counterpart, shifted by a phase difference of $π$ 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˙j$⁠, 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˙j$ 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°M2°M3$⁠, 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.³⁷ 

The difference in the relative impact velocities shown in (a) for hard impacts is closely related to $ϕ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˙j$⁠, particularly when the velocities of the ball $Zj′$ and the capsule $Yj′$ have opposite sign. In contrast, subsequent impacts on the same side occur when the ball and capsule are close in phase, yielding small $x˙j$⁠. 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⁡(Yj′)=−sgn⁡(Zj′)$⁠. 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˙j$⁠, followed by a second impact on the bottom with $sgn⁡(Yj′)=sgn⁡(Zj′)$⁠, yielding small $x˙j$⁠.

We analyze the results from the soft impact model, first comparing the stiff case $κ=ω0π/ω≫1$ to the hard impact model. Then, we consider smaller values of $ω0$ while fixing the excitation frequency at $ω=π$ 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 $κ=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 $β=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 $ρ$ ratios instead of a fixed $r$ as given in the hard model, with the distribution of the $ρ$ values concentrated at the reference value $r=0.5$⁠, as follows from the calibration of $c$ for large $κ$⁠. The spread in $ρ$ is due to small differences in the excitation phase $ϕj,c$ at each impact, due to the (small) time interval $Δτj,c$ during which the ball is in contact with the membrane as shown in Fig. 5(h). The resulting differences in the phase $ϕ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 $κ=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 $κ=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 $ω0$ (and thus $κ$⁠), with $ω$ fixed. (The influence of different values of $ω$ is studied in Sec. III B.) Figure 6 shows different bifurcation diagrams with $κ$ as the bifurcation parameter at different values of $d$⁠, with other parameters as in Table II. Large values of $κ$ indicate large stiffness, for which the soft impact model results agree with the hard impact case indicated by dashed red lines. Smaller values of $κ$ indicate a soft membrane, for which the behavior naturally differs from that of the hard impact model. The diagrams are shown for $κ$ 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 $κ$ until a rapid decrease to near zero impact velocities at $κ=4.04$⁠. For these lower values of $κ$⁠, we identify co-existence of different stable attractors. The bi-stable orbits at $κ=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 $κ$⁠, 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 $κ$ and transitions to 1:2 and then to a 1:1 periodic solution at $κ=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 $κ$ decreases the penetration of the ball into the membrane increases. Similarly, for the case of $d=0.16$⁠, shown in Fig. 6(c), as $κ$ decreases the system undergoes transitions from a 1:2 periodic orbit to 1:1 behavior at $κ=57.14$⁠. For lower values of $κ$ with $d=0.16$⁠, the 1:1 impact velocities increase to a maximum at $κ=3.84$⁠, followed by a rapid decrease of the impact velocities for smaller $κ$⁠. 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 $κ$ 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 $κ=400$⁠, with its attractor seen in Fig. 5(c). As the value of $κ$ decreases from $κ=400$⁠, the chaotic behavior eventually transitions to a 1:2 orbit at $κ=36.36$ and then to a 1:1 solution at $κ=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 $κ$ 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 $κ$⁠. 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 $κ$⁠.

The influence of $κ$ has implications for energy generation. Decreasing $κ$ leads to a period-$T$ behavior with only high velocity impacts. These impacts have the largest contribution at an optimum value of $κ$ for a given $d$⁠. However, the impact velocity decreases substantially for $κ$ 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 $κ$ decreases past its optimum value, it is important to notice that nonlinear effects persist for smaller $κ$⁠, apparent through the other attracting orbits that coexist with the near-linear response as observed for $d=0.2$ and $κ=3$ in Fig. 7(a). From a design perspective, low values of $κ$ can be beneficial for energy extraction, but the system must be designed in such a way to ensure $κ$ 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 $κ$ 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³⁶ (not shown in Fig. 5), with equal magnitude of impact velocities $|X′(τ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 $β$⁠, including $β=0$⁠. As would be expected, for $β>0$⁠, this type of asymmetry is enhanced throughout the dynamics.³⁸ Figures 6 and 7 indicate that this bifurcation from symmetric to asymmetric 1:1 orbits also occurs for increasing $κ$ as well as for different combinations of $r$ and $ω$ 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 $Δτj,c$⁠. Thus, this second type of asymmetric bifurcation cannot be observed in a model with instantaneous impacts.

Figure 8 illustrates the influence of $κ$ on the bifurcations with $d$ as the bifurcation parameter, comparing these for several values of $κ$⁠. Figure 8(a) shows that as $κ$ decreases, the regions with 1:1 periodic orbits with larger impact velocities $u˙$ are pushed to lower values of $d$⁠, until the 1:1 behavior becomes the attractor for all $d$ as can be seen for $κ=12.65$⁠. For smaller values of $κ$⁠, 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 $κ$⁠, experience substantial qualitative changes as $κ$ decreases. The values of $ρ$ shown in Fig. 8(b) for each impact also yield valuable insights. For larger values of $κ$⁠, $ρ$ varies slightly around its reference value $r$ from the corresponding hard impact model. However, as $κ$ decreases, the variability of values in $ρ$ increases, particularly for chaotic ranges. It is important to emphasize that, for softer membranes and lower impact velocities, there are values of $ρ>1$⁠. This does not violate any energy principle as the greater contact time $Δτ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 $ρ$ that diverge from the reference value $r$ for the hard model. This behavior is observed for $κ=12.65$ and $d=0.142$⁠, where a symmetry breaking bifurcation yields two branches in the diagram for $ρ$ vs $d$⁠, highlighted by a black circle in Fig. 8(b). There is also a slight asymmetry in the contact times $Δτ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 $ρ$ follows from the behavior of the excitation phase $ϕj,c$⁠, which has slight variation around $ϕ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 $ϕj,c$ corresponds to the symmetry breaking point of $ρ$⁠, 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 $κ$ 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 $ρ$⁠, the impact velocities can be even higher for small $κ$⁠, for which the system appears to be nearly optimized in terms of both the contact time $Δtj,c$ and the corresponding excitation phase $ϕj,c$⁠, yields the ideal conditions for energy harvesting.