## Abstract

The spatial variation of the coefficient of restitution for frictionless impacts along the length of a circular beam is investigated using a continuous impact model. The equations of motion are obtained using the finite element method, and direct time integration is used to simulate the collision on a fast time scale. For collision of a pinned beam with a fixed cylinder, the spatial variation of the coefficient of restitution, impulse magnitude, duration of collision, energetics, and the role of damping are investigated. In the absence of significant external damping, the kinematic and kinetic definitions of the coefficient of restitution provide identical results. Experiments validate the results from simulation which indicate that the coefficient of restitution is sensitive to the location of impact.

## 1 Introduction

Multibody systems frequently have non-smooth dynamics due to impacts or collisions between the bodies [1]. Impacts last for a very short duration in comparison to the overall motion of the system and may involve multiple sub-impacts. In many systems, it is important to model the collision accurately for predicting the behavior of the bodies post-collision; examples are juggling systems [2–4] and hopping locomotion [5,6].

In modeling systems with impact, it is often the rigid body dynamics which is of interest. Impacts result in energy being removed from the rigid body motion of the bodies. Depending on the geometric and material properties of the colliding bodies, the configuration of the collision, and the relative velocity of the bodies before impact, the loss of energy may be associated with one or more of the following:

A single dimensionless parameter, the *coefficient of restitution*, is used to quantify the energy losses during an impact; this parameter also relates to the velocities associated with the rigid body motion before and after collision.

In the literature [9,15], impact is usually treated in two ways, depending on the problem being considered. They are as follows:

Instantaneous impact models, or

*algebraic models*[16,17], in which an impact is considered to be of infinitesimal duration; the positions of the colliding bodies do not change, and there is an instantaneous jump in their velocities. The magnitude of the jump in velocity is quantified by the coefficient of restitution. These models are simpler and computationally much cheaper than the alternative, discussed next.Continuous impact models, or

*incremental models*[11], in which an impact has short but finite duration, and the dynamics of the system during collision is explicitly considered. These models make use of contact force laws, the simplest of which is the purely elastic Hertzian contact law [18]. Viscous and hysteretic dissipation may additionally be incorporated into the contact force laws—see Refs. [10,11,19–21], for example.

Previous works have considered impact on flexible bodies assuming prior knowledge of the coefficient of restitution, and a non-smooth velocity of the flexible body post-impact [22–30]. The focus here, however, is on computation of the coefficient of restitution for such impacts. Stoianovici and Hurmuzlu [31] used kinematic data to obtain the coefficient of restitution for collisions of steel bars, observing that internal vibrations and multiple sub-impacts occur during collision. The multi-scale simulation approach proposed by Seifried and collaborators [15,32–37] is capable of incorporating both elastic and plastic material behavior. The simulation approach was used to compute the coefficient of restitution for the impacts of spheres on rods and beams for different relative velocities. The problem of multiple impacts at the same location, which leads to a gradual evolution in the value of the coefficient of restitution due to permanent deformation at the contact point, was also considered. The simulation results were validated with experiments employing high-bandwidth vibrometers. It was noted that modal models with an elastostatic Hertzian contact law are sufficiently accurate and computationally inexpensive when the colliding bodies have purely elastic material behavior. Bhattacharjee and Chatterjee [9,38,39] considered the impact of a ball on a beam. The equations of motion were obtained from modal expansion and consistent prediction of the coefficient of restitution required a sufficiently high number of assumed modes and some level of modal damping. The coefficient of restitution was determined for impacts at many different locations on the beam and for a set of different boundary and initial conditions.

This paper investigates the spatial variation of the coefficient of restitution for impacts of a point mass at different locations along the length of a flexible beam. A continuous impact model [15] is used and energy loss from the rigid body motion of the beam is primarily due to excitement of the vibration modes. Starting from an Euler–Bernoulli beam model which accounts for both internal and external damping, the Galerkin finite element method, followed by modal reduction, is employed to obtain a system of equations amenable to direct time integration. The contact force is modeled using a Hertzian contact law [18], which can be easily modified to include dissipative effects. The finite element method provides flexibility in applying boundary conditions and geometric constraints at locations different from the ends, and does not require prior knowledge of the mode shapes of the beam. Different beam geometries and contact configurations can be accommodated by changing a few model parameters.

In relation to the available literature, the contributions of this work are summarized as follows:

The coefficient of restitution is computed using two definitions, kinematic and kinetic, and their values are compared.

The spatial variation of the coefficient of restitution, the magnitude of the impulse, and the duration of the collision are obtained.

The effect of internal and external damping on the coefficient of restitution is considered.

A simple experimental setup, which uses an encoder instead of a complex motion capture system, is used to obtain the coefficient of restitution of a pinned beam impacting a fixed cylinder at different locations along its length. The experimental results are compared with those obtained through simulations.

The sensitive nature of the spatial variation of the coefficient of restitution and magnitude of impulse, captured through simulations and corroborated through experiments, will be useful for addressing the challenging control problem of non-prehensile manipulation of extended objects using impulsive forces.

## 2 Problem Description

*xy*plane. The cylinder is assumed to have a length that is much smaller than

*ℓ*, and its dynamics is described by that of a point mass. However, the material properties and radius of the cylinder influence the contact force between the colliding masses. The collision occurs at

*x*

_{c}∈ (0,

*ℓ*) and the relative velocity between the cylinder and the beam in the

*y*direction is given by

*w*associated with the rigid body motion. It is assumed that the relative velocity in the

*x*direction is zero and therefore frictional effects can be neglected. For continuous impact modeling of the collision, we define

*t*

_{0}as the instant of first contact between the bodies, and

*t*

_{f}as the instant when the bodies separate such that they do not come into contact again. For a collision to occur,

*V*

_{rel}(

*t*

_{0}) > 0. The collision duration (

*t*

_{f}−

*t*

_{0}) is divided into a compression phase [

*t*

_{0},

*t*

_{c}] marked by decreasing relative velocity with

*V*

_{rel}(

*t*

_{c}) = 0, and a restitution phase (

*t*

_{c},

*t*

_{f}] marked by the separation of the bodies:

*V*

_{rel}(

*t*) < 0 for

*t*>

*t*

_{c}. Since the collision duration is small, the effect of gravity forces can be neglected.

The collision results in a loss of kinetic energy associated with the rigid body motion of the beam, quantified by the dimensionless coefficient of restitution $\epsilon $. The energy lost from the rigid body mode is transferred to the vibration modes of the beam, where it is dissipated due to internal (material) damping. The objective of the present study is to determine the value of $\epsilon $ for a range of *x*_{c} ∈ (0, *ℓ*), different values of *V*_{rel}(*t*_{0}) and different boundary conditions of the beam.

## 3 Mathematical Model

### 3.1 Governing Equations.

*P*≡

*P*(

*t*) at point

*x*

_{c}is given by [40]

*δ*(.) is the Dirac delta function. The governing equation (2) is second order in

*t*, and fourth order in

*x*. It can be solved subject to two initial conditions and four boundary conditions. Assuming that the beam is undeformed before the collision, the initial conditions can be expressed as

*V*

_{0}is the velocity of the center-of-mass of the beam in the

*y*direction and Ω

_{0}is the angular velocity of the beam, defined to be positive along the

*z*-axis, at

*t*=

*t*

_{0}. For a free end, the boundary conditions correspond to zero bending moment and zero shear force, which are described by the relations

*x*=

*x*

_{c}. Its dynamics is described by

*S*

_{0}is the initial velocity of the cylinder in the

*y*direction.

### 3.2 Contact Force.

*x*=

*x*

_{c}and

*r** is the effective

*Gaussian radius of curvature*of the surface. Before and after the impact, the value of

*P*is necessarily zero since ϱ is zero. The potential energy of the elastic deformation at the contact point is given by

The contact law in Eq. (8) does not account for dissipative effects or permanent deformation at the contact point. For the range of relative velocities of impact considered in this work, the losses due to these effects are expected to be negligible compared to the loss due to energy transfer to the vibration modes of the beam.

### 3.3 Nondimensionalization.

*ψ*

_{e},

*ψ*

_{i}denote the nondimensional external and internal damping per unit length of the beam, and $P~$ is the nondimensional contact force; they are given by the relations

## 4 Finite Element Analysis of Collision

### 4.1 Finite Element Model.

*R*≡

*R*(

*ξ*,

*τ*) can be expressed as

*W*≡

*W*(

*ξ*) is the weight function. Substituting Eq. (20) in Eq. (21), we obtainThe term in Eq. (22) is integrated by parts twice to obtainThe elemental weak form of Eq. (22) is given by

*H*

_{01},

*H*

_{11},

*H*

_{02},

*H*

_{12}are the Hermite interpolation polynomials, provided in Appendix.

*W*to be the shape functions in Eq. (A1), we obtain the elemental equation of motion

**M**

^{e},

**C**

^{e}, and

**K**

^{e}are the elemental mass, damping, and stiffness matrices, and

**f**

^{e}is the elemental forcing vector; their expressions are given in Appendix. The global equation of motion of the beam, comprised of

*n*nodes, is obtained by assembling the elemental equations of motion in Eq. (25)

**v**∈

*R*

^{2n}is the vector of degrees-of-freedom of the beam;

**M**,

**C**,

**K**∈

*R*

^{2n×2n}are the mass, damping, and stiffness matrices, which are constant for a fixed number of nodes; and

**f**∈

*R*

^{2n}is the forcing vector.

It must be noted that, in writing Eq. (26), any non-zero boundary terms arising from and are included in the forcing vector **f** [42]. It follows from Eq. (13) that the boundary terms associated with a free end, arising from and , are identically zero.

**K**in Eq. (26) with

**K**

_{pen}

*γ*≫ tr(

**K**) is the penalty term and

**v**

_{con}∈

*R*

^{2n}identifies the constrained degrees-of-freedom.

### 4.2 Transformation Into Modal Coordinates.

*p*modes, we transform Eq. (26) into modal coordinates. We express

**v**as

**u**is the vector of the first

*p*modal displacements and

*ϕ*

_{j}∈

*R*

^{2n},

*j*= 1, 2, …,

*p*, are the first

*p*mass-normalized eigenvectors of the model in Eq. (26), obtained by solving the eigenvalue problem

**Φ**

^{T}, we obtain

**Φ**

^{T}

**MΦ**=

**I**, with

**I**defined to be the identity matrix, and

**Ψ**=

*ψ*

_{e}

**I**+

*ψ*

_{i}

**Λ**; Eq. (30) then represents a system of

*p*decoupled equations. In general,

**Ψ**cannot be expressed as a linear combination of

**I**and

**Λ**, and therefore Eq. (30) represents

*p*coupled equations.

### 4.3 Time Integration.

We adapted the Newmark algorithm [41] to perform direct time integration of the system in Eq. (30) to investigate the collision problem—see Algorithm 1 in Appendix. The accuracy and stability of the integration scheme are governed by parameters *δ* and *α*, with unconditional stability and maximum accuracy guaranteed by the constant average acceleration scheme, where *δ* = 1/2 and *α* = 1/4.

The accuracy of the finite element analysis depends on the number of nodes *n*, the number of modes *p* retained in obtaining Eq. (30), and the choice of the time-step Δ*τ* in Algorithm 1. More accurate results are obtained for larger values of *n* and *p* and smaller values of Δ*τ*, at the cost of computational efficiency. The value of *p* must be chosen to include all vibration modes to which a significant amount of energy is transferred. Similarly, Δ*τ* must be sufficiently small to accurately capture the dynamics on the fast time scale.

### 4.4 Coefficient of Restitution.

*τ*

_{c}and

*τ*

_{f}. Using Eq. (1) and the discussion in Sec. 2,

*τ*

_{c}is identified as the instant when

*τ*

_{f}is identified as the instant of separation of the bodies such that they do not come into contact again, i.e.,

Following the literature [15], we consider the following two definitions of the coefficient of restitution $\epsilon $:

- Kinematic coefficient of restitution(32)$\epsilon N=\u2212V~rel(\tau f)V~rel(\tau 0)$
- Kinetic coefficient of restitution(33)$\epsilon P=\u222b\tau c\tau fP~d\tau \u222b\tau 0\tau cP~d\tau $

## 5 Simulation Results

### 5.1 Geometric and Material Properties.

Simulation results are presented for an aluminum beam colliding with a fixed aluminum cylinder with geometric and material properties identical to that in experiments—see Tables 1 and 2. The center-of-mass of the beam is taken to be at the geometric center of the beam, i.e., $x\xaf=\u2113/2$. To simulate an inertially fixed cylinder, like that in experiments, the value of *m*_{c} is chosen to be arbitrarily large compared to the mass of the beam.

### 5.2 Boundary Conditions.

We simulate collisions of a beam that is pinned at *x* = *x*_{p}, *x*_{p} = 1.905 × 10^{−2} m, i.e., *ξ* = 0.0375, to be consistent with our experiments.^{2} Depending on the value of the offset, the number of elements is chosen such that a node coincides with the location of the pinned joint. If the *j*th node coincides with the pinned joint, the entries of **v**_{con} are all zeros except for the (2*j* − 1)th entry, which is unity.

### 5.3 Finite Element Setup.

The beam is discretized into 80 elements, each of length *η* = 0.0125. This results in *n* = 81 nodes along the length of the beam; **M**, **C**, **K** matrices of dimension 162 × 162; and the pinned joint to coincide with the fourth node, i.e., *j* = 4. The stiffness matrix **K** is modified using Eq. (27); all entries of **v**_{con} are zero, except the 7th entry, which is unity; the value of *γ* was chosen to be tr(**K**) × 10^{3}. In transforming to modal coordinates using Eq. (28), the value of *p* is chosen to be 25; this comprises the rigid body mode and 24 vibration (bending) modes. It was found that using larger values of *p* did not affect the results obtained significantly. The time-step for numerical integration is chosen to be 0.1 *μ*s, which corresponds to Δ*τ* = 6.1765 × 10^{−6}. The small value of the time-step allows the dynamics on the fast time scale to be captured in detail.

### 5.4 Results Without Damping.

When internal and external dampings are zero, the damping matrix **C** is zero, and Eq. (30) represents a system of 25 decoupled equations. Since the cylinder is stationary, *S*_{0} = 0. Three initial angular velocities of the beam are considered: Ω_{0} = {− 4, − 6, − 8}; rad/s. Since the beam is pinned, the initial velocity of the center-of-mass $V0=\Omega 0(x\xaf\u2212xp)$. For all three initial angular velocities, collisions are simulated for 400 points of impact between *ξ* = 0.2 and *ξ* = 1. It should be noted that a single collision may involve multiple sub-impacts, i.e., contact between the beam and the cylinder is made and broken several times over the collision duration.

For all three values of Ω_{0}, the variation of the coefficient of restitution $\epsilon $ with the nondimensional location of impact *ξ*_{c} are shown in Fig. 2. For each Ω_{0}, the plot of $\epsilon $ is actually a superimposition of the plots of $\epsilon N$ and $\epsilon P$ obtained using Eqs. (32) and (33); these plots are indistinguishable from one another in conformity with results in the literature [43]. The plots in Fig. 2 indicate that the value of $\epsilon $ does not depend significantly on the initial velocity of Ω_{0} but strongly depends on the location of impact *ξ*_{c}; this is in conformity with the results in Ref. [9]. In particular, $\epsilon $ assumes high values for low values of *ξ*_{c} and then varies somewhat erratically with increasing *ξ*_{c}, reaching a local minima at *ξ*_{c} ≈ 0.64 and again at *ξ*_{c} ≈ 0.85.

A low value of $\epsilon $ implies that a large fraction of the kinetic energy at *τ* = *τ*_{0} is transferred to the vibration modes of the beam at *τ* = *τ*_{f}. For Ω_{0} = −6 rad/s, the fractions of the kinetic energy transferred to the first four vibration modes of the beam individually and all vibration modes of the beam cumulatively are shown in Fig. 3 as a function of the location of impact *ξ*_{c}. The largest fractions of the initial kinetic energy transmitted to the vibration modes correspond to the local minima of $\epsilon $ at *ξ*_{c} ≈ 0.64 and *ξ*_{c} ≈ 0.85. Summing the kinetic energy of the rigid body mode with the energy transferred to the vibration modes at *τ* = *τ*_{f} gives the initial kinetic energy of the beam to within $\xb11.2%$ due to errors associated with numerical time integration.

For Ω_{0} = −6 rad/s, the magnitude of the nondimensional impulse $\u222b\tau 0\tau fP~d\tau $ is shown in Fig. 4 as a function of *ξ*_{c}. It can be observed that closer proximity of the impact location to the pinned support is associated with larger magnitudes of the impulse. The nondimensional duration of collision (*τ*_{f} − *τ*_{0}) is also shown in Fig. 4 for Ω_{0} = −6 rad/s. Similar to the plot of nondimensional impulse, the duration of collision shows a decreasing trend with increase in the value of *ξ*_{c} until *ξ*_{c} ≈ 0.85; the duration of collision then increases sharply and remains higher for values of *ξ*_{c} exceeding 0.85. For a given value of *ξ*_{c}, increase in the magnitude of Ω_{0} results in an increase in the magnitude of the impulse but no significant change in the collision duration.

### 5.5 Results With Damping.

We consider the effect of internal damping on the coefficient of restitution. For Ω_{0} = −6 rad/s, Fig. 5 shows the variation of $\epsilon $ with *ξ*_{c} for *c*_{i} = 10^{6} Ns/m^{2}, which corresponds to *ψ*_{i} = 9.0432 × 10^{−4}. As in Fig. 2, the plot is comprised of the plots of both $\epsilon N$ and $\epsilon P$, which remain indistinguishable with the inclusion of internal damping. For reference, the plot of $\epsilon $ with *ξ*_{c} in the absence of damping is shown in gray. It can be seen that with inclusion of internal damping, $\epsilon $ varies more smoothly with *ξ*_{c} and assumes lower values for almost all impact locations. The locations of the local minima of $\epsilon $ remain unchanged at *ξ*_{c} ≈ 0.64 and at *ξ*_{c} ≈ 0.85. The magnitude of the nondimensional impulse and the nondimensional duration of collision also vary more smoothly with inclusion of internal damping but there is no significant change in their values; these plots are not provided here.

The inclusion of external damping, of the same magnitude observed in our experimental hardware, produced no perceptible change in the value of $\epsilon $. However, for large values of external damping, the values of $\epsilon P$ and $\epsilon N$ were no longer indistinguishable; the values of $\epsilon P$ matched the values of $\epsilon $ in the absence of damping for all *ξ*_{c} but the values of $\epsilon N$ were lower, with the difference especially prominent for smaller values of *ξ*_{c}.

### 5.6 Mechanics of a Collision With Sub-Impacts.

For a single collision, we present the temporal variation of nondimensional quantities that describe the behavior of the system during the collision. We consider the point of impact to be *ξ*_{c} = 0.68, the initial angular velocity of the beam to be Ω_{0} = −6 rad/s, and no damping. We choose *τ*_{0} = 0. The instant when $V~rel(\tau )$ equals zero, i.e., the end of the compression phase, is *τ*_{c} = 0.0540. This collision ends at *τ*_{f} = 0.0981.^{3}

*v*(

*ξ*

_{c},

*τ*) is shown in Fig. 6(a). Since there are five distinct time intervals during which

*v*(

*ξ*

_{c},

*τ*) < 0 (denoted by ①, ②, ③, ④, and ⑤), the collision consists of five sub-impacts. The nondimensional relative velocity $V~rel(\tau )$ is shown in Fig. 6(b). Since $V~rel(\tau 0)=0.0624$ and $V~rel(\tau f)=\u22120.0397$, the kinematic coefficient of restitution $\epsilon N=0.6362$. The nondimensional contact force $P~(\tau )$ is shown in Fig. 6(c); the five time intervals where $P~(\tau )>0$ are identical to those denoted by ①, ②, ③, ④, ⑤ in Fig. 6(a) and indicate times when the beam and cylinder are in contact. The nondimensional impulses over the compression and restitution phases were found to be

*τ*∈ [

*τ*

_{0},

*τ*

_{f}], these quantities sum to unity within an error of $\xb10.4%$, implying conservation of energy.

We now consider a collision at the same location, with the same initial velocity, but with internal damping *c*_{i} = 10^{6} Ns/m^{2}. The compression phase now ends at *τ*_{c} = 0.0497 and the collision ends at *τ*_{f} = 0.0927. For this collision, only two sub-impacts occur. The kinematic and kinetic coefficients of restitution were found to be $\epsilon N=\epsilon P=0.4950$; these are lower than the numbers obtained in the absence of damping. Due to the presence of damping, $39.09%$ of the initial energy is lost during the collision.

## 6 Experimental Validation

### 6.1 Experimental Setup.

Experiments were performed with the aluminum beam and aluminum cylinder whose geometric and material properties are provided in Tables 1 and 2. The beam is rigidly connected to a vertical shaft at *x*_{p} = 1.905 × 10^{−2} m and is able to rotate freely in the horizontal plane. The vertical shaft is mounted perpendicular to the horizontal TMC MICRO-g optical breadboard (see Fig. 7) by means of a ball-bearing, not visible in the diagram. The axis of the shaft coincides with one of the hole locations on the optical breadboard; the center point of this hole is taken to be the origin of the inertial *XY* frame fixed to the table. The *xy* frame is fixed to the beam, as described in Fig. 1.

*X*and

*Y*directions, thus allowing the cylinder to be placed at (

*X*,

*Y*) = (0.0254

*I*, 0.0254

*J*), where

*I*,

*J*= 1, 2, …. From the geometry of the setup, it can be shown that the location of impact

*x*

_{c}is given by

*I*and

*J*combinations, experiments were conducted with a set of 125 unique values of

*x*

_{c}for which

*ξ*

_{c}was in the range between 0.2 and 1.0. For each value of

*ξ*

_{c}, multiple collisions of the beam were recorded with varying angular velocities.

### 6.2 Data Acquisition and Processing.

The time history of the angular position of the beam for each collision was collected using a Dynapar E23 miniature encoder. The encoder shaft was coupled to the free end of the shaft to which the beam is attached. The encoder has a resolution of 2000 pulses/revolution and the angular position signal was acquired at 2000 Hz using a dSpace DS1104 board in the matlab/simulink environment.

*t*

_{0}) and Ω(

*t*

_{f}), associated with the rigid body motion of the beam were obtained by finding the slopes of the straight lines fitted to the residual angular position data before and after the impact. The linear fit was performed on unfiltered data, as it was observed that using a lowpass Butterworth filter to remove the oscillatory components of the post-collision data did not affect the slope of the fitted lines. The kinematic coefficient of restitution for a collision is given by

### 6.3 Results.

The experimental results for the 125 impact locations with *ξ*_{c} ∈ [0.2, 1.0] are shown in Fig. 8 using circles and error bars indicating their standard deviation. For each impact location, five collisions were recorded on average. The average initial angular velocity was ≈ −3 rad/s. Simulation results for Ω_{0} = −3 rad/s,^{4} with and without damping, are also presented in Fig. 8. For the case with damping, *c*_{i} = 10^{6} Ns/m^{2} and *c*_{e} = 1.37 × 10^{−2} Ns/m^{2}. The value of *c*_{i} matches that chosen in Sec. 5.5; the value of *c*_{e} corresponds to the experimental hardware.^{5} Figure 8 shows a fairly good match between simulation and experimental results. In particular:

The decreasing trend in the value of $\epsilon $ for

*ξ*_{c}∈ [0.2, 0.64], predicted in simulations, is captured in experiments.For

*ξ*_{c}∈ [0.64, 0.8], the experimental results agreed especially well with the simulation results with damping.In the vicinity of the minima at

*ξ*_{c}≈ 0.85, there was some mismatch between experiment and simulation.The increasing trend in the value of $\epsilon $ for

*ξ*_{c}∈ [0.9, 1.0], predicted by simulations, is captured in experiment.

The experimental results showed excellent repeatability, with minor variability for *ξ*_{c} > 0.8, where the simulation results also show more fluctuation.

## 7 Conclusions

This paper considers the problem of frictionless collisions between a flexible circular beam and a cylinder modeled as a point mass. The focus is on the computation of the coefficient of restitution for the collisions, which quantifies the energy lost from the rigid body motion of the beam and cylinder. For this collision configuration, this loss is primarily due to energy transfer to vibration modes of the beam. The discretized equations of motion are obtained using the finite element method and geometric constraints are enforced using the penalty method. Following modal reduction, the equations are solved in forward time using an adaptation of the Newmark algorithm for chosen impact locations. The simulations are capable of capturing the response of the system on a fast time scale and show phenomena such as multiple sub-impacts during a single collision. The coefficient of restitution is computed using the kinematic and kinetic definitions; they yield almost identical results in the absence of significant external damping. A large number of simulations were carried out for collisions between a pinned beam and a fixed cylinder; the coefficient of restitution was found to be sensitive to the location of the point of impact. In the absence of damping, the energy lost from the rigid body motion is completely transferred to the vibration modes of the beam; the nature of the spatial variation of the coefficient of restitution can therefore be attributed to the mode shapes of the beam. The initial relative velocity of impact has negligible effect on the coefficient of restitution and the duration of the collision. It must be noted that different beam geometries and collision configurations can be accommodated by changing only a few model parameters, and geometric constraints can be enforced at arbitrary locations on the beam.

Experiments were conducted using a simple setup and data acquisition scheme that focuses on the capture of the rigid body motion of the beam to validate the simulation results. A fairly good match between simulation and experimental results provides confidence in the observed spatial variation of the coefficient of restitution. The current model does not account for dissipation and permanent deformation at the contact point; the investigation of these effects lie in the scope of our future work. Future work will also focus on obtaining reliable predictions of the coefficient of restitution in the presence of parameter uncertainty. The results obtained in this paper will be useful for realizing the control problem of nonprehensile robotic manipulation of extended objects using impulsive forces. In particular, the results will be used to realize stick-juggling using a robot manipulator, the control design for which has appeared in the literature [3,4].

## Footnotes

In our experiments, described in Sec. 6, the pinned joint is slightly offset from *ξ* = 0 to facilitate the mounting of the beam on a shaft connected to a bearing.

The simulation was carried out for a much longer duration than *τ*_{f} to ensure that the beam and the cylinder did not come into contact again.

Note that the results in Sec. 5.4 indicate that the choice of Ω_{0} does not affect the simulation results significantly.

Note that the results in Sec. 5.5 indicate that small values of *c*_{e} had no perceptible effect on the value of $\epsilon $.

## Acknowledgment

The authors acknowledge the support provided by the National Science Foundation, Grant CMMI-2043464.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

## Nomenclature

*r*=radius of beam cross section (m)

*s*=displacement of cylinder along the

*y*-axis (m)*t*=time (s)

*v*=nondimensional displacement of a point on the beam along the

*y*-axis*w*=displacement of a point on the beam along the

*y*-axis (m)- $x\xaf$ =
*x*coordinate of the center-of-mass of the beam (m)*c*_{e}=external damping per unit length of beam (Ns/m

^{2})*c*_{i}=internal (material) damping per unit length of beam (Ns/m

^{2})*m*_{c}=mass of cylinder (kg)

*r*_{c}=radius of cylinder (m)

*x*_{c}=*x*coordinate of the point of application of*P*(m)*E*=Young’s modulus of beam (Pa)

*I*=area moment of inertia of beam cross section (m

^{4})*P*=magnitude of contact force (N)

*E*_{c}=Young’s modulus of cylinder (Pa)

*xy*=coordinate system fixed to beam; the

*x*-axis coincides with the undeformed axis of the beam and the origin lies at one end of the beam—see Fig. 1*ℓ*=length of beam (m)

- $\epsilon $ =
coefficient of restitution

*μ*=mass per unit length of beam (kg/m)

*ν*=Poisson’s ratio of beam

*ν*_{c}=Poisson’s ratio of cylinder

*ξ*=nondimensional coordinate along the

*x*-axis*σ*=nondimensional displacement of cylinder along the

*y*-axis*τ*=nondimensional time

## Appendix: Finite Element Modeling

*a*)

*b*)

*c*)

*d*)

*H*(·) is the Heaviside unit step function;

*H*(0) ≜ 0.

#### Direct time integration of the collision dynamics

**Initial calculations**

Initialize $\tau =0$ and $u(0)$, $u\u02d9(0)$, and $u\xa8(0)$

Select nondimensional time-step $\Delta \tau $

Choose $\delta =1/2$, $\alpha =1/4$

**for**$\tau =\tau +\Delta \tau do$

Estimate transverse displacement $v(\xi c,\tau )$ using Eq. (24)

Compute nondimensional contact force $P~(\tau +\Delta \tau )$ using $\u03f1~(\tau +\Delta \tau )$ in Eq. (17)

Form nondimensional forcing vector $b(\tau +\Delta \tau )$ using $P~(\tau +\Delta \tau )$ in Eq. (31)

**end for**