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 [24] 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:

  • Energy transfer to the flexible modes [79],

  • Dissipation at the contact point while the colliding bodies are in contact [10,11], and

  • Plastic deformation of the bodies due to high stresses at the contact point [1214].

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,1921], 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 [2230]. 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,3237] 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

Consider the collision between an Euler–Bernoulli beam of circular cross section and a short cylinder whose axis is perpendicular to that of the beam—see Fig. 1. The motion of the beam and the cylinder is confined to the 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 xc ∈ (0, ) and the relative velocity between the cylinder and the beam in the y direction is given by
(1)
where w¯ is the component of 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 t0 as the instant of first contact between the bodies, and tf as the instant when the bodies separate such that they do not come into contact again. For a collision to occur, Vrel(t0) > 0. The collision duration (tft0) is divided into a compression phase [t0, tc] marked by decreasing relative velocity with Vrel(tc) = 0, and a restitution phase (tc, tf] marked by the separation of the bodies: Vrel(t) < 0 for t > tc. Since the collision duration is small, the effect of gravity forces can be neglected.
Fig. 1
A circular beam and a cylinder at the instant of collision
Fig. 1
A circular beam and a cylinder at the instant of collision
Close modal

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 ε. 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 ε for a range of xc ∈ (0, ), different values of Vrel(t0) and different boundary conditions of the beam.

3 Mathematical Model

3.1 Governing Equations.

The dynamics of an Euler–Bernoulli beam, subjected to an external transverse load PP(t) at point xc is given by [40]
(2)
where δ(.) 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
(3)
where V0 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 = t0. For a free end, the boundary conditions correspond to zero bending moment and zero shear force, which are described by the relations
(4)
For a pinned end, the boundary conditions correspond to zero deflection and zero bending moment, and are given by
(5)
The cylinder is modeled as a point mass; it is subject to a force equal in magnitude and opposite in direction to the load acting on the beam at x = xc. Its dynamics is described by
(6)
which can be solved subject to the two initial conditions
(7)
where S0 is the initial velocity of the cylinder in the y direction.

3.2 Contact Force.

Assuming purely elastic deformation of the colliding bodies at the point of impact, a Hertzian contact law [18] is used to model the contact force
(8)
(9)
(9)
where ϱ is the instantaneous local deformation of the bodies at x = xc 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
(10)
Remark 1

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.

Introducing the change of variables
in Eq. (2), we get the nondimensional governing equation
(11)
where ψ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
The initial conditions in Eq. (3) are expressed in nondimensional form as
(12)
where ξ¯x¯/. The nondimensional boundary conditions for free and pinned ends are obtained from Eqs. (4) and (5) as
(13)
(14)
The dynamics of the cylinder (6) in nondimensional form becomes
(15)
where m~c is the nondimensional mass of the cylinder. The above equation is subject to the initial conditions
(16)
which follow from Eq. (7).
The contact law in Eqs. (8) and (9) in nondimensional form is
(17)
(18)
From Eq. (10), the nondimensional potential energy is obtained as
(19)

4 Finite Element Analysis of Collision

4.1 Finite Element Model.

The finite element discretization of Eq. (11) is carried out using the Galerkin method [41]. The domain residual RR(ξ, τ) can be expressed as
(20)
The weighted residual form is
(21)
where WW(ξ) is the weight function. Substituting Eq. (20) in Eq. (21), we obtain
(22)
The term in Eq. (22) is integrated by parts twice to obtain
The elemental weak form of Eq. (22) is given by
(23)
For a beam element with two nodes, each node having two degrees-of-freedom, the nondimensional displacement within the element can be approximated by
(24)
where the shape functions H01, H11, H02, H12 are the Hermite interpolation polynomials, provided in Appendix.
Substituting Eq. (24) into Eq. (23) and choosing the weight function W to be the shape functions in Eq. (A1), we obtain the elemental equation of motion
(25)
where Me, Ce, and Ke are the elemental mass, damping, and stiffness matrices, and fe 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)
(26)
where vR2n is the vector of degrees-of-freedom of the beam; M, C, KR2n×2n are the mass, damping, and stiffness matrices, which are constant for a fixed number of nodes; and fR2n 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.

We impose geometric boundary conditions using the penalty method. We replace the stiffness matrix K in Eq. (26) with Kpen
(27)
where γ ≫ tr(K) is the penalty term and vconR2n identifies the constrained degrees-of-freedom.

4.2 Transformation Into Modal Coordinates.

To investigate the transfer of energy into the first p modes, we transform Eq. (26) into modal coordinates. We express v as
(28)
where u is the vector of the first p modal displacements and ϕjR2n, j = 1, 2, …, p, are the first p mass-normalized eigenvectors of the model in Eq. (26), obtained by solving the eigenvalue problem
(29)
Using Eq. (28) in Eq. (26), and pre-multiplying both sides by ΦT, we obtain
(30)
where we use the property that ΦT = I, with I defined to be the identity matrix, and
(31)
For a free-free beam, it can be shown that the damping matrix Ψ = ψeI + ψ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.

Remark 2

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.

To compute the coefficient of restitution ε, we identify the nondimensional time instants τc and τf. Using Eq. (1) and the discussion in Sec. 2, τc is identified as the instant when
and τ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 ε:

  1. Kinematic coefficient of restitution
    (32)
  2. Kinetic coefficient of restitution
    (33)
The coefficient of restitution ε will be computed using both definitions in the next section.

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¯=/2. To simulate an inertially fixed cylinder, like that in experiments, the value of mc is chosen to be arbitrarily large compared to the mass of the beam.

Table 1

Properties of aluminum beam

ParameterValue (Unit)
0.5080 (m)
r6.35 × 10−3 (m)
μ0.3433 (kg/m)
I = πr4/41.277 × 10−9 (m4)
E68.3 × 109 (Pa)
ν0.34
ParameterValue (Unit)
0.5080 (m)
r6.35 × 10−3 (m)
μ0.3433 (kg/m)
I = πr4/41.277 × 10−9 (m4)
E68.3 × 109 (Pa)
ν0.34
Table 2

Properties of aluminum cylinder

ParameterValue (Unit)
rc6.35 × 10−3 (m)
mc1000 (kg)
Ec68.3 × 109 (Pa)
νc0.34
ParameterValue (Unit)
rc6.35 × 10−3 (m)
mc1000 (kg)
Ec68.3 × 109 (Pa)
νc0.34

5.2 Boundary Conditions.

We simulate collisions of a beam that is pinned at x = xp, xp = 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 jth node coincides with the pinned joint, the entries of vcon are all zeros except for the (2j − 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 vcon are zero, except the 7th entry, which is unity; the value of γ was chosen to be tr(K) × 103. 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, S0 = 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=Ω0(x¯xp). 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 ε with the nondimensional location of impact ξc are shown in Fig. 2. For each Ω0, the plot of ε is actually a superimposition of the plots of εN and ε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 ε 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, ε 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.

Fig. 2
Variation of the coefficient of restitution ε with nondimensional location of impact ξc for three initial angular velocities: Ω0 = { − 4, − 6, − 8} rad/s; the plots indicate that ε does not strongly depend on the value of Ω0
Fig. 2
Variation of the coefficient of restitution ε with nondimensional location of impact ξc for three initial angular velocities: Ω0 = { − 4, − 6, − 8} rad/s; the plots indicate that ε does not strongly depend on the value of Ω0
Close modal

A low value of ε 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 ε 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 ±1.2% due to errors associated with numerical time integration.

Fig. 3
Fraction of kinetic energy at τ = τ0 transferred to the vibration modes of the beam at τ = τf for Ω0 = −6 rad/s
Fig. 3
Fraction of kinetic energy at τ = τ0 transferred to the vibration modes of the beam at τ = τf for Ω0 = −6 rad/s
Close modal

For Ω0 = −6 rad/s, the magnitude of the nondimensional impulse τ0τfP~dτ 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.

Fig. 4
Variation of nondimensional impulse and duration of collision with nondimensional location of impact ξc for Ω0 = −6 rad/s
Fig. 4
Variation of nondimensional impulse and duration of collision with nondimensional location of impact ξc for Ω0 = −6 rad/s
Close modal

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 ε with ξc for ci = 106 Ns/m2, which corresponds to ψi = 9.0432 × 10−4. As in Fig. 2, the plot is comprised of the plots of both εN and εP, which remain indistinguishable with the inclusion of internal damping. For reference, the plot of ε with ξc in the absence of damping is shown in gray. It can be seen that with inclusion of internal damping, ε varies more smoothly with ξc and assumes lower values for almost all impact locations. The locations of the local minima of ε 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.

Fig. 5
Effect of internal damping on the spatial variation of the coefficient of restitution ε for Ω0 = −6 rad/s; the plot for ci = 0 is taken from Fig. 2
Fig. 5
Effect of internal damping on the spatial variation of the coefficient of restitution ε for Ω0 = −6 rad/s; the plot for ci = 0 is taken from Fig. 2
Close modal

The inclusion of external damping, of the same magnitude observed in our experimental hardware, produced no perceptible change in the value of ε. However, for large values of external damping, the values of εP and εN were no longer indistinguishable; the values of εP matched the values of ε in the absence of damping for all ξc but the values of ε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(τ) equals zero, i.e., the end of the compression phase, is τc = 0.0540. This collision ends at τf = 0.0981.3

The nondimensional displacement of the contact point 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(τ) is shown in Fig. 6(b). Since V~rel(τ0)=0.0624 and V~rel(τf)=0.0397, the kinematic coefficient of restitution εN=0.6362. The nondimensional contact force P~(τ) is shown in Fig. 6(c); the five time intervals where P~(τ)>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
resulting in a kinetic coefficient of restitution εP=0.6365. Clearly, εN and εP are almost identical. The temporal evolution of the kinetic energy associated with the rigid body motion of the beam, the energy contained in the vibration modes, and the potential energy of elastic deformation at the contact point U~, all expressed as fractions of the initial kinetic energy of the beam, are shown in Fig. 6(d). For all τ ∈ [τ0, τf], these quantities sum to unity within an error of ±0.4%, implying conservation of energy.
Fig. 6
Simulation of a collision of the beam with the cylinder at ξc = 0.68 with Ω0 = −6 rad/s, in the absence of damping
Fig. 6
Simulation of a collision of the beam with the cylinder at ξc = 0.68 with Ω0 = −6 rad/s, in the absence of damping
Close modal

We now consider a collision at the same location, with the same initial velocity, but with internal damping ci = 106 Ns/m2. 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 εN=ε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 xp = 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.

Fig. 7
Experimental setup for collision between an aluminum beam and an aluminum cylinder. The encoder is connected to a dSpace data acquisition system, not shown in the figure.
Fig. 7
Experimental setup for collision between an aluminum beam and an aluminum cylinder. The encoder is connected to a dSpace data acquisition system, not shown in the figure.
Close modal
The aluminum cylinder was screwed at different hole locations on the optical breadboard, aligning the axis of the cylinder each time with the vertical axis. The holes on the optical breadboard are 1.0 in. (0.0254 m) apart in the 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 xc is given by
(34)
By choosing different I and J combinations, experiments were conducted with a set of 125 unique values of xc 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.

From the angular position data, we obtained the angular velocities of the beam before and after the collision as follows. First, the instant of impact was identified based on the knowledge of the angular position of the beam when it comes in contact with the cylinder. Following this, 0.3 s of angular position data on either side of the instant of impact was extracted and separated. It was observed that the pre-collision data were smooth whereas the post-collision data were oscillatory. However, the high-frequency oscillations decayed very rapidly. To eliminate the noise very close to the instant of impact, 0.003 s of data immediately before and after the impact was removed. The angular velocities, Ω(t0) and Ω(tf), 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
(35)

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, ci = 106 Ns/m2 and ce = 1.37 × 10−2 Ns/m2. The value of ci matches that chosen in Sec. 5.5; the value of ce 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 ε 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 ε for ξc ∈ [0.9, 1.0], predicted by simulations, is captured in experiment.

Fig. 8
Experimental and simulation results for ξc ∈ [0.2, 1.0] are compared using four plots. The experimental results are shown using circles and errors bars corresponding to the standard deviation. Simulation results, both in the presence and absence of damping, are presented.
Fig. 8
Experimental and simulation results for ξc ∈ [0.2, 1.0] are compared using four plots. The experimental results are shown using circles and errors bars corresponding to the standard deviation. Simulation results, both in the presence and absence of damping, are presented.
Close modal

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

2

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.

3

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.

4

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

5

Note that the results in Sec. 5.5 indicate that small values of ce had no perceptible effect on the value of ε.

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¯ =

x coordinate of the center-of-mass of the beam (m)

ce =

external damping per unit length of beam (Ns/m2)

ci =

internal (material) damping per unit length of beam (Ns/m2)

mc =

mass of cylinder (kg)

rc =

radius of cylinder (m)

xc =

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 (m4)

P =

magnitude of contact force (N)

Ec =

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)

ε =

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

The Hermite interpolation polynomials are
(A1)
The elemental mass, stiffness, damping matrices, and the forcing vector are
(A2a)
(A2b)
(A2c)
(A2d)
where H(·) is the Heaviside unit step function; H(0) ≜ 0.

Direct time integration of the collision dynamics

Algorithm 1

Initial calculations

Initialize τ=0 and u(0), u˙(0), and u¨(0)

Select nondimensional time-step Δτ

Choose δ=1/2, α=1/4

Calculate integration constants:
Form effective stiffness matrix:
forτ=τ+Δτdo

 Estimate transverse displacement v(ξc,τ) using Eq. (24)

 Estimate nondimensional deformation at contact point ϱ~(τ+Δτ) by modifying Eq. (18) as follows:

 Compute nondimensional contact force P~(τ+Δτ) using ϱ~(τ+Δτ) in Eq. (17)

 Form nondimensional forcing vector b(τ+Δτ) using P~(τ+Δτ) in Eq. (31)

 Form effective forcing vector at τ+Δτ as follows:
 Solve for displacements at τ+Δτ:
 Calculate accelerations and velocities at τ+Δτ:

end for

References

1.
Brogliato
,
B.
,
2016
,
Nonsmooth Mechanics: Models, Dynamics and Control. Communications and Control Engineering
,
Springer International Publishing
,
Cham, Switzerland
.
2.
Brogliato
,
B.
, and
Rio
,
A.
,
2000
, “
On the Control of Complementary-Slackness Juggling Mechanical Systems
,”
IEEE Trans. Automat. Contr.
,
45
(
2
), pp.
235
246
.
3.
Kant
,
N.
, and
Mukherjee
,
R.
,
2021
, “
Non-Prehensile Manipulation of a Devil-Stick: Planar Symmetric Juggling Using Impulsive Forces
,”
Nonlinear Dyn.
,
103
(
3
), pp.
2409
2420
.
4.
Khandelwal
,
A.
,
Kant
,
N.
, and
Mukherjee
,
R.
,
2023
, “
Nonprehensile Manipulation of a Stick Using Impulsive Forces
,”
Nonlinear Dyn.
,
111
(
1
), pp.
113
127
.
5.
Mathis
,
F. B.
, and
Mukherjee
,
R.
,
2016
, “
Apex Height Control of a Two-Mass Robot Hopping on a Rigid Foundation
,”
Mech. Mach. Theory
,
105
, pp.
44
57
.
6.
Allafi
,
A.
, and
Mukherjee
,
R.
,
2019
Apex Height Control of a Two-DOF Ankle-Knee-Hip Robot Hopping on a Rigid Foundation
,”
2019 ASME Dynamic Systems and Control Conference
,
Park City, UT
,
Oct. 8–11
.
7.
Raman.
,
C. V.
,
1920
, “
On Some Applications of Hertz’s Theory of Impact
,”
Phys. Rev.
,
15
(
4
), pp.
277
284
.
8.
Zener
,
C.
, and
Feshbach
,
H.
,
2021
, “
A Method of Calculating Energy Losses During Impact
,”
ASME J. Appl. Mech.
,
6
(
2
), pp.
A67
A70
.
9.
Bhattacharjee
,
A.
,
2019
, “
New Approximations in Vibroimpact Problems
,” Ph.D. thesis,
Department of Mechanical Engineering, Indian Institute of Technology
,
Kanpur, India
.
10.
Hunt
,
K. H.
, and
Crossley
,
F. R. E.
,
1975
, “
Coefficient of Restitution Interpreted as Damping in Vibroimpact
,”
ASME J. Appl. Mech.
,
42
(
2
), pp.
440
445
.
11.
Lankarani
,
H. M.
, and
Nikravesh
,
P. E.
,
1994
, “
Continuous Contact Force Models for Impact Analysis in Multibody Systems
,”
Nonlinear Dyn.
,
5
(
2
), pp.
193
207
.
12.
Thornton
,
C.
,
1997
, “
Coefficient of Restitution for Collinear Collisions of Elastic-Perfectly Plastic Spheres
,”
ASME J. Appl. Mech.
,
64
(
2
), pp.
383
386
.
13.
Jackson
,
R. L.
,
Green
,
I.
, and
Marghitu
,
D. B.
,
2010
, “
Predicting the Coefficient of Restitution of Impacting Elastic-Perfectly Plastic Spheres
,”
Nonlinear Dyn.
,
60
(
3
), pp.
217
229
.
14.
Yigit
,
A. S.
,
Christoforou
,
A. P.
, and
Majeed
,
M. A.
,
2011
, “
A Nonlinear Visco-Elastoplastic Impact Model and the Coefficient of Restitution
,”
Nonlinear Dyn.
,
66
(
4
), pp.
509
521
.
15.
Seifried
,
R.
,
Schiehlen
,
W.
, and
Eberhard
,
P.
,
2010
, “
The Role of the Coefficient of Restitution on Impact Problems in Multi-body Dynamics
,”
Proc. Inst. Mech. Eng., Part K: J. Multi-body Dyn.
,
224
(
3
), pp.
279
306
.
16.
Pfeiffer
,
F.
, and
Glocker
,
C.
,
1996
,
Multibody Dynamics With Unilateral Contacts
, 1st ed.,
John Wiley & Sons, Ltd
,
Germany
.
17.
Brach
,
R. M.
,
1989
, “
Rigid Body Collisions
,”
ASME J. Appl. Mech.
,
56
(
1
), pp.
133
138
.
18.
Popov
,
V. L.
,
2017
,
Contact Mechanics and Friction
,
Springer Berlin Heidelberg
,
Berlin, Heidelberg
.
19.
Gharib
,
M.
, and
Hurmuzlu
,
Y.
,
2012
, “
A New Contact Force Model for Low Coefficient of Restitution Impact
,”
ASME J. Appl. Mech.
,
79
(
6
), p.
064506
.
20.
Zhang
,
J.
,
Li
,
W.
,
Zhao
,
L.
, and
He
,
G.
,
2020
, “
A Continuous Contact Force Model for Impact Analysis in Multibody Dynamics
,”
Mech. Mach. Theory
,
153
, p.
103946
.
21.
Rodrigues Da Silva
,
M.
,
Marques
,
F.
,
Tavares Da Silva
,
M.
, and
Flores
,
P.
,
2022
, “
A Compendium of Contact Force Models Inspired by Hunt and Crossley’s Cornerstone Work
,”
Mech. Mach. Theory
,
167
, p.
104501
.
22.
Khulief
,
Y. A.
, and
Shabana
,
A. A.
,
1986
, “
Impact Responses of Multi-body Systems With Consistent and Lumped Masses
,”
J. Sound Vib.
,
104
(
2
), pp.
187
207
.
23.
Wagg
,
D. J.
, and
Bishop
,
S. R.
,
2002
, “
Application of Non-Smooth Modeling Techniques to the Dynamics of a Flexible Impacting Beam
,”
J. Sound Vib.
,
256
(
5
), pp.
803
820
.
24.
Wagg
,
D. J.
,
2004
, “
A Note on Using the Collocation Method for Modelling the Dynamics of a Flexible Continuous Beam Subject to Impacts
,”
J. Sound Vib.
,
276
(
3
), pp.
1128
1134
.
25.
Wagg
,
D. J.
,
2007
, “
A Note on Coefficient of Restitution Models Including the Effects of Impact Induced Vibration
,”
J. Sound Vib.
,
300
(
3
), pp.
1071
1078
.
26.
Vyasarayani
,
C. P.
,
Sandhu
,
S. S.
, and
McPhee
,
J.
,
2012
, “
Nonsmooth Modeling of Vibro-Impacting Euler-Bernoulli Beam
,”
Adv. Acoust. Vibr.
,
2012
, p.
e268595
.
27.
Vyasarayani
,
C. P.
,
McPhee
,
J.
, and
Birkett
,
S.
,
2009
, “
Modeling Impacts Between a Continuous System and a Rigid Obstacle Using Coefficient of Restitution
,”
ASME J. Appl. Mech.
,
77
(
2
), p.
021008
.
28.
Yigit
,
A. S.
,
Ulsoy
,
A. G.
, and
Scott
,
R. A.
,
1990
, “
Dynamics of a Radially Rotating Beam With Impact, Part 1: Theoretical and Computational Model
,”
ASME J. Vib. Acoust.
,
112
(
1
), pp.
65
70
.
29.
Yigit
,
A. S.
,
Ulsoy
,
A. G.
, and
Scott
,
R. A.
,
1990
, “
Dynamics of a Radially Rotating Beam With Impact, Part 2: Experimental and Simulation Results
,”
ASME J. Vib. Acoust.
,
112
(
1
), pp.
71
77
.
30.
Palas
,
H.
,
Hsu
,
W. C.
, and
Shabana
,
A. A.
,
1992
, “
On the Use of Momentum Balance and the Assumed Modes Method in Transverse Impact Problems
,”
ASME J. Vib. Acoust.
,
114
(
3
), pp.
364
373
.
31.
Stoianovici
,
D.
, and
Hurmuzlu
,
Y.
,
1996
, “
A Critical Study of the Applicability of Rigid-Body Collision Theory
,”
ASME J. Appl. Mech.
,
63
(
2
), pp.
307
316
.
32.
Schiehlen
,
W.
, and
Seifried
,
R.
,
2004
, “
Three Approaches for Elastodynamic Contact in Multibody Systems
,”
Multibody Sys.Dyn.
,
12
(
1
), pp.
1
16
.
33.
Schiehlen
,
W.
,
Hu
,
B.
, and
Seifried
,
R.
,
2005
, “Multiscale Methods for Multibody Systems With Impacts,”
Advances in Computational Multibody Systems
, Computational Methods in Applied Sciences,
Ambrósio
,
J. A.
, ed.,
Springer Netherlands
,
Dordrecht
, pp.
95
124
.
34.
Schiehlen
,
W.
,
Seifried
,
R.
, and
Eberhard
,
P.
,
2006
, “
Elastoplastic Phenomena in Multibody Impact Dynamics
,”
Comput. Methods Appl. Mech. Eng.
,
195
(
50–51
), pp.
6874
6890
.
35.
Seifried
,
R.
,
Schiehlen
,
W.
, and
Eberhard
,
P.
,
2005
, “
Numerical and Experimental Evaluation of the Coefficient of Restitution for Repeated Impacts
,”
Int. J. Impact Eng.
,
32
(
1
), pp.
508
524
.
36.
Seifried
,
R.
,
Minamoto
,
H.
, and
Eberhard
,
P.
,
2010
, “
Viscoplastic Effects Occurring in Impacts of Aluminum and Steel Bodies and Their Influence on the Coefficient of Restitution
,”
ASME J. Appl. Mech.
,
77
(
4
), p.
041008
.
37.
Minamoto
,
H.
,
Seifried
,
R.
,
Eberhard
,
P.
, and
Kawamura
,
S.
,
2011
, “
Analysis of Repeated Impacts on a Steel Rod With Visco-Plastic Material Behavior
,”
Eur. J. Mech. A/Solids
,
30
(
3
), pp.
336
344
.
38.
Bhattacharjee
,
A.
, and
Chatterjee
,
A.
,
2017
, “
Interplay Between Dissipation and Modal Truncation in Ball-Beam Impact
,”
ASME J. Comput. Nonlinear Dyn.
,
12
(
6
), p.
061018
.
39.
Bhattacharjee
,
A.
, and
Chatterjee
,
A.
,
2020
, “
Restitution Modeling in Vibration-Dominated Impacts Using Energy Minimization Under Outward Constraints
,”
Int. J. Mech. Sci.
,
166
, p.
105215
.
40.
Svedholm
,
C.
,
Zangeneh
,
A.
,
Pacoste
,
C.
,
François
,
S.
, and
Karoumi
,
R.
,
2016
, “
Vibration of Damped Uniform Beams With General End Conditions Under Moving Loads
,”
Eng. Struct.
,
126
, pp.
40
52
.
41.
Bathe
,
K.-J.
,
2014
,
Finite Element Procedures
, 2nd ed.,
K. J. Bathe
,
Watertown, MA
.
42.
Mohsen
,
M. F. N.
,
1982
, “
Some Details of the Galerkin Finite Element Method
,”
Appl. Math. Model.
,
6
(
3
), pp.
165
170
.
43.
Stronge
,
W. J.
,
2018
,
Impact Mechanics
, 2nd ed.,
Cambridge University Press
,
Cambridge, MA
.