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
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.
3.2 Contact Force.
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.
4 Finite Element Analysis of Collision
4.1 Finite Element Model.
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.
4.2 Transformation Into Modal Coordinates.
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.
Following the literature [15], we consider the following two definitions of the coefficient of restitution :
- Kinematic coefficient of restitution(32)
- Kinetic coefficient of restitution(33)
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., . 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.
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 . 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 and 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.
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 due to errors associated with numerical time integration.
For Ω0 = −6 rad/s, the magnitude of the nondimensional impulse 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 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 and , 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.
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 and were no longer indistinguishable; the values of matched the values of in the absence of damping for all ξc but the values of 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 equals zero, i.e., the end of the compression phase, is τc = 0.0540. This collision ends at τf = 0.0981.3
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 ; these are lower than the numbers obtained in the absence of damping. Due to the presence of damping, 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.
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.
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.
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 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 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
Direct time integration of the collision dynamics
Initial calculations
Initialize and , , and
Select nondimensional time-step
Choose ,
Estimate transverse displacement using Eq. (24)
Compute nondimensional contact force using in Eq. (17)
Form nondimensional forcing vector using in Eq. (31)
end for