- Email: [email protected]
A linear complementarity algorithm for rigid body impact with friction
Eur. J. Mech. A/Solids 18 (1999) 703–717 Elsevier, Paris
A linear complementarity algorithm for rigid body impact with friction Lars Johansson Linköping University, Department of Mechanical Engineering, S-581 83 Linköping, Sweden (Received 10 June 1998; revised and accepted 18 November 1998) Abstract – In this paper an algorithm for the impact of rigid bodies against rigid obstacles is developed. Starting from a normal contact law, Coulomb’s law of friction, the equations of motion and a simple discretization scheme, a linear complementarity problem (LCP) is formulated and used to compute the post-impact velocities. The predictions of the algorithm is compared to an experiment, designed such that the tangential impulses at impact are important, with encouraging results. Elsevier, Paris contact / impact / friction / linear complementarity problem
1. Introduction This paper is concerned with a numerical algorithm for the impact of rigid bodies against rigid obstacles, and with its experimental verification. Some problems of this kind can be treated by specifying the quotient between the relative normal velocity of approach and separation (Newton, 1687), that is, by introducing the classical coefficient of restitution. This is sometimes generalized to the so-called Poisson’s hypothesis (Kilmister and Reeve, 1966), in order to avoid some anomalies arising when the theory is applied to rigid bodies. In many cases however, such as in the experiment described in Section 4 below, it is necessary to take both the normal and the tangential impulses at the impact into account to achieve reasonable agreement with observations. This class of problems has been the object of several recent studies, invariably leading to the introduction of one or more additional constitutive parameters, such as the coefficient of friction. Thus in (Brach, 1989), a quotient between normal and tangential impulses are introduced and several bounds, based on physical assumptions, are derived for this quotient, whereas in (Walton, 1993) a model for spheres with three parameters is used which is, for the case of spheres, equivalent to the model in the present paper. In (Stronge, 1990), the division of the impact process into a compression and an expansion phase is analysed, and the problem is treated using a coefficient of restitution relating energies rather than velocities. In this paper, a contact law and Coulomb’s friction law is formulated in terms of velocities and impulses, in a manner similar to that used for quasistatic contact problems in (Johansson and Klarbring, 1992). These laws are combined with the laws of motion and, after discretization, form a Linear Complementarity Problem (LCP) which is solved to obtain the velocities after impact. From a modelling point of view, the approach taken in the present paper is similar to that taken by (Moreau, 1988), which is much concerned with the mathematical setting of the problem, but where a very elegant algorithm for the case of completely inelastic impacts is also given. See also (Moreau, 1994) where granular materials are simulated with an algorithm involving iteration between the contact laws and the equations of motion, applying the contact laws to a linear combination of the pre and post impact velocities.
704
L. Johansson
Figure 1. Geometry of the contact.
It has, of course, been noted before that contact problems lend themselves to formulation in terms of complementarity conditions. Recent work involving LCP based algorithms for rigid body impact include (Pfeiffer and Glocker, 1996) and (Anitescu et al., 1995).
2. Computational algorithm In this section the computational algorithm is described. The treatment of the contacts is described in detail in Sections 2.1–2.3 and an outline of the overall algorithm is given in Section 2.4. 2.1. Equations of motion The motion of a single body in contact with a rigid wall, shown in figure 1, will now be studied. We are assuming that it has been determined that the body studied is violating a contact condition and the subsequent motion is to be determined. The aim is to compute the position, angular orientation, velocity and angular velocity at time t j +1 , assuming that these quantities are known at an earlier time t j . The body is acted upon by unknown contact impulses Pn and Pt , and by additional forces and a moment which are assumed to be known at times t j and t j +1 , so that a suitable approximation to the corresponding impulses Px , Py and M during the time interval can be calculated. In the case of external forces prescribed as depending on the position or velocity of the body itself or another body, these cannot be known beforehand at time t j +1 . This issue is resolved iteratively as described in Section 2.4.
A linear complementarity algorithm for rigid body impact with friction
705
Denoting the components of the velocity of the contact point at time t j by x˙ j and y˙ j and at time t j +1 by x˙ and y˙ j +1 , the equations of motion can be expressed in terms of these motions as j +1
x˙ j +1 = x˙ j +
1 1 1 r1 Pt + Pn + Px + M, ¯ m1 m Izz I
(1)
1 1 1 r2 Pn + Pt + Py + M. (2) ¯ m2 m I I zz Here 1/m1 = (r12 /Izz + 1/m), 1/m2 = (r22 /Izz + 1/m) and 1/I¯ = r1 r2 /Izz , with m denoting the mass and Izz the mass moment of inertia about the centre of gravity of the body. r1 and r2 are the distances between the contact point and centre of gravity as shown in figure 1, which are assumed constant during the time interval [t j , t j +1 ]. y˙ j +1 = y˙ j +
2.2. Contact conditions The contact law in the normal direction is written in terms of the normal velocity and the normal impulse y˙ > 0,
Pn > 0,
yP ˙ n = 0.
(3)
It cannot be assumed that y˙ is continuous. This is obvious for the case of impact, but in fact other cases of discontinuous velocities can occur, as is shown by the Painleve example (Brogliato, 1996). It is assumed, following (Moreau, 1994), that the right and left limits of y˙ exist, and that the normal contact law can be applied to a linear combination of these limits. Thus we put
en en 1 1 y˙ − + y˙ + > 0, Pn > 0, y˙ − + y˙ + Pn = 0, (4) 1 + en 1 + en 1 + en 1 + en where en is a constitutive parameter, to be chosen according to the situation to be modelled. Aiming for a linear complementarity problem formulation, we follow the approach used for quasistatic elastic contact problems in (Johansson and Klarbring, 1992), and introduce the following quantities relevant to the tangential direction (
ψ1 = −Pt + µPn , ψ2 = Pt + µPn
and
(5)
x˙ − |x| ˙ , χ1 = −
2
˙ χ = x˙ + |x| . 2
(6)
2
These definitions are made such that two complementary pairs are formed. Thus if x˙ is positive then χ1 will be zero and, if Pt obeys Coulomb’s law of friction, ψ1 will be equal to 2µPn , and therefore positive. If, on the other hand, x˙ is negative, then χ1 will be equal to the absolute value of x˙ and ψ1 will be zero. Similarly χ2 and ψ2 form a complementarity pair. We thus have: (
ψ1 > 0, χ1 > 0,
ψ1 χ1 = 0,
ψ2 > 0, χ2 > 0,
ψ2 χ2 = 0.
(7)
706
L. Johansson
As is shown by the Painleve example (Brogliato, 1996), it cannot be assumed that the tangential velocity x˙ is continuous either. It is assumed that the right and left limits of x˙ exist and that the conditions (7) remain valid when (6) is replaced by et x˙ − + x˙ + − |et x˙ − + x˙ + | , χ1 = −
2(1 + et )
e x˙ − + x˙ + + |et x˙ − + x˙ + | χ2 = t ,
(8)
2(1 + et )
where et is a constitutive parameter. It should be noted (Moreau, 1994), that the constitutive parameters en and et introduced above will determine the nature of the contact. Considering the case of impact there is a nonzero normal impulse, Pn > 0, and (4) implies y˙ + = −en y˙ − .
(9)
Thus the parameter en can be interpreted as a restitution coefficient, with en = 0 corresponding to purely plastic impact and en = 1 to purely elastic impact. In the tangential direction there can be either stick or slip. If the coefficient of friction is sufficiently large that there is no slip, we have χ1 = χ2 = 0. (8) then gives x˙ + = −et x˙ −
(10)
and et can be interpreted as a tangential restitution coefficient. It should be noted, however, that the interpretation of et as a tangential restitution coefficient is possible only for large values of µ, since the tangential impulse is limited by Coulomb’s law of friction, and in the sliding case, for a nonspherical body, there is a coupling between normal and tangential contact impulses which complicates things further. With the introduction of these restitution coefficients, there is a total number of three parameters describing the impact process; µ, en and et . This is somewhat problematic, since these parameters must be determined experimentally. In particular there is no obvious simple method for measuring et . This increases the scope for adjusting computations to fit a particular experiment. This might be perceived as an advantage, but since the main reason for developing computational algorithms is to be able to make a priori predictions of experimental outcomes, it is actually more likely to be a disadvantage. To overcome this difficulty it would be convenient to have some rule of thumb for the choice of et . To this end the choices et = 0 and et = en were considered. Of these two possibilities et = en should be preferred since et = 0 could not, by any stretch of the imagination, represent the behaviour observed in Horak’s experiments as described in Section 5 below. Also, if bodies with a nonspherical shape is considered, it is possible to construct examples where energy is gained in an impact if en 6= et , as discussed in Section 3.3. Finally, if the linear combinations introduced in Eqs (4) and (8) are interpreted as specifying an intermediate time between the beginning and end of the impact process at which the contact laws are applied, it seems natural to put en = et . 2.3. Discretization and linear complementarity problem For the discretization, a sequence of times [t 1 , . . . , t j , t j +1 , . . . , t n ] covering the desired time interval is introduced. Next, it is assumed that if there is a violation of a contact condition at time t j then the conditions (4) and (7) can be applied approximately with the right and left limits of the velocities at a particular instant
A linear complementarity algorithm for rigid body impact with friction
707
replaced with the velocities at the beginning and end of the timestep. Thus, Eqs (4) and (8) are written: en 1 φ= y˙ j + y˙ j +1 > 0, Pn > 0, 1 + en 1 + en
en 1 y˙ j + y˙ j +1 Pn = 0, 1 + en 1 + en
(11)
et x˙ j + x˙ j +1 − |et x˙ j + x˙ j +1 | χ = − , 1 2(1 + e ) t
(12)
et x˙ j + x˙ j +1 + |et x˙ j + x˙ j +1 | . χ2 =
2(1 + et )
j +1
Using Eq. (12) to eliminate x from Eq. (1) and the result to eliminate Pt from Eqs (2) and (5) and also using the definition of φ in Eq. (11), three equations are obtained, which together with the complementarity conditions (7) and (11), form the following linear complementarity problem: LCP problem Assuming that x˙ j , y˙ j , Px , Py and M are known, find φ, ψ1 , ψ2 , Pn , χ1 and χ2 such that
ψ1 = ψ2
φ
m1 m1 r1 (1 + et )m1 j 1 1 r2 − M Py + x˙ − Px + (1 + en )m (1 + en )Izz I¯ (1 + en )I¯ (1 + en )I¯m m1 m1 r1 j m1 (1 + et )x˙ + M Px + m Izz m m r 1 1 1 j −m1 (1 + et )x˙ − M Px − m Izz 1 1 m1 (1 + et )m1 (1 + et )m1 − 2 − 1 + en m2 I¯ (1 + en )I¯ (1 + en )I¯ Pn m +µ + 1 χ1 m (1 + e ) −m (1 + e ) 1 t 1 t , I¯ χ2 m1 µ− −m1 (1 + et ) m1 (1 + et ) I¯
j y˙ −
φ > 0, Pn > 0,
ψ1 > 0, χ1 > 0, ψ2 > 0, χ2 > 0,
φPn = 0, ψ1 χ1 = 0, ψ2 χ2 = 0.
This problem is recognized as a linear complementarity problem, or LCP, a class of problems which has been the subject of much study, see, e.g. (Murty, 1988). The solution to the linear complementarity problem w = q + Mz, w > 0,
z > 0,
wt z = 0
can be found using a pivot algorithm. Thus, if all elements in q are greater than or equal to zero, w = q and z = 0 is a solution. If not, the equations are rearranged so that one or more of the unknowns in w are exchanged for unknowns in z. This procedure is repeated in some systematic way until an arrangement is found where all elements in q are greater than or equal to zero. The solution is then obtained by putting w = q and z = 0, where, of course, the unknowns occupying w and z are not the same as before the pivots. Specifically, in the computations for this paper, q was scanned in the order q1 , q2 , q3 and as soon as a qi < 0 was found the corresponding unknowns wi and zi were exchanged, and the procedure was repeated until a solution was found.
708
L. Johansson
Once the LCP problem is solved, Pn , Pt , x˙ j +1 and y˙ j +1 can be extracted and used to compute the velocity and rotation of the centre of gravity at time t j +1 . 2.4. The overall algorithm The aim is to compute the position, angular orientation, velocity and angular velocity at time t j +1 , assuming that these quantities are known at an earlier time t j . The motion of a single body can be computed directly with the equations in Sections 2.1–2.3. In this case, the applied loads can be assumed to be known both at times t j and t j +1 and suitable approximations of the corresponding impulses are easily computed. Thus, it is first checked if the body violates any of the contact conditions. If it is, the velocities at time t j +1 are computed by solving the LCP problem in Section 2.3. This calculation proceeds from the position the body had at the end of the previous timestep, accepting without correction a small violation of the contact conditions. If no contact condition is violated the equations of motion are used directly, without contact impulses. If, instead, the motion of a system of bodies interconnected with springs and dashpots is sought, there will be coupling forces which cannot be assumed to be known at time t j +1 , since, in order to incorporate the forces from the connecting springs and dashpots, they must be assumed to depend on positions and velocities at time t j +1 . This difficulty is resolved iteratively. Thus it is first assumed that the coupling forces at time t j +1 are the same as at time t j . The position and velocities at time t j +1 are then computed, making it possible to compute better values for the coupling loads. This process is repeated until the position of each body is within a prescribed tolerance between two consecutive iterations. 3. The special case of impacting spheres This section is concerned with the detailed analysis of the algorithm for the special case r2 = 0, covering spheres, cylinders at right angles with the plane of motion, and also the case when a rod impacts at right angles against a surface. In Section 3.1 below the LCP is analysed for this special case. Then, in Section 3.2 it is shown that the present model is equivalent to a model by Walton, see (Walton, 1993) in this special case. Section 3.3, finally, is devoted to the question of energy consistency of the algorithm. 3.1. The LCP algorithm for the case r2 = 0 In this section, the case r2 = 0 is studied. This results in 1/I¯ = 0 and m2 = m, and the coupling terms Pn /I¯ and Pt /I¯ between the normal and tangential directions disappear in Eqs (1) and (2). Since one of the aims of this section is to compare the present model with the Walton model, and since this model is concerned with what happens at the instant of an impact, the LCP problem is applied to an infinitesimal time interval [t − , t + ]. This also means that the impulses Px , Py and M from external loads become zero in Eqs (1) and (2), if it is assumed that these impulses are derived from forces which are bounded. This also makes sense when studying the energy consistency of the model, since impulses from external forces is a legitimate way for a body to gain energy. In short, we study the behaviour of the model without making the approximation of replacing [t − , t + ] with [t j , t j +1 ], with external forces assumed bounded and with r2 = 0. With these assumptions the first equation in the LCP problem, which is for the normal direction, becomes independent of the remaining two, and can be written together with the corresponding complementarity conditions as 1 φ = y˙ − + (13) Pn , φ > 0, Pn > 0, φPn = 0. m(1 + en )
A linear complementarity algorithm for rigid body impact with friction
709
For the case y˙ − > 0, Eqs (13), (5) and (7) imply Pn = Pt = 0. This is the case when contact is breaking. For the remainder of this section the case y˙ − < 0, giving nonzero contact impulses, will be studied. Equation (13) then gives the normal contact impulse as Pn = −m(1 + en )y˙ − .
(14)
Using (14), the LCP problem takes the form LCP problem, r2 = 0, y˙ − < 0
m1 (1 + et )x˙ − − µm(1 + en )y˙ − (
ψ1 > 0, χ1 > 0,
ψ1 χ1 = 0,
ψ2 > 0, χ2 > 0,
ψ2 χ2 = 0.
−1 χ1 = + m1 (1 + et ) 1 , −1 1 ψ2 −m1 (1 + et )x˙ − − µm(1 + en )y˙ − χ2
ψ1
(15)
(16)
This LCP admits two different types of solutions: stick and slip. Analysing first the case of stick, this occurs when µ is sufficiently large that both elements in the first column matrix to the right in (16) is positive (with y˙ − < 0). The solution is then obtained by putting χ1 = χ2 = 0 in (15), which gives, together with (5) Pt = −m1 (1 + et )x˙ − .
(17)
Turning next to the slip case, it is first noted that slip can be either in the positive or the negative direction, depending on whether x˙ − > 0 or x˙ − < 0, the case x˙ − = 0 being covered by the stick solution discussed above. These two cases behave in exactly the same way (in the r2 = 0 case presently studied, but not in the general case), the only difference being the slip direction, and we choose to study the x˙ − > 0 case. When x˙ − is increased, the second row element of the first column matrix to the right in Eq. (15) will become negative; the LCP solver will rearrange the equations so that ψ2 is exchanged for χ2 . Equation (15) will then take the form
−2µm(1 + en )y˙ − 0 −1 ψ1 χ1 = + . 1 µm(1 + en ) − − 1 x˙ + y˙ χ2 ψ2 m1 (1 + et ) m1 (1 + et )
(18)
A solution satisfying (18) and (16) is obtained by putting χ1 = ψ2 = 0. This requires that both members of the first column matrix to the right are positive, which is ascertained since we are studying the x˙ − > 0, y˙ − < 0 case and since −m1 (1 + et )x˙ − − µm(1 + en )y˙ − < 0
(19)
is necessary for the stick assumption to be abandoned. From (18) and (5) and since y˙ − < 0 we have Pt = µm(1 + en )y˙ − 6 0.
(20)
−m1 (1 + et )x˙ − 6 Pt 6 0
(21)
Equations (19), (20) and (17) give
710
L. Johansson
covering both the stick and slip in the positive direction cases. The case of slip in the negative direction can be given a completely analogous treatment, giving the same condition but with reversed inequalities. Finally, (18) and (8) gives µm(1 + en ) − y˙ , m1 where, again, the case of slip in the negative direction can be given an analogous treatment. x˙ + = x˙ − +
(22)
3.2. Comparison with the Walton model For the case of spheres, it is stated in (Moreau, 1994) that the present model is equivalent to the model for impacting spheres involving three parameters described in (Walton, 1993), and it can be seen from the equations in Section 3.1 that this is indeed the case. The Walton model is also described in (Foerster et al., 1994) where it is used for comparison with experiments. In fact, Eqs (9), (10), (19) and (22) are direct counterparts in the present notation to the relevant equations in (Foerster et al., 1994). In the Walton model, a coefficient of restitution is introduced for the normal direction, equivalent to Eq. (9). In the tangential direction, both a coefficient of restitution and a coefficient of friction (relating impulses rather than forces) are introduced, equivalent to Eqs (10) and (22). The coefficient of restitution is applied when the tangential relative velocity is small and the coefficient of friction when it is large. The switch between these two equations for the tangential direction is when they predict the same tangential velocity, as determined by Eq. (19), although Walton actually begins with slip assumptions and swiches to stick if necessary. It should be noted that the crucial point in specializing the present model to that of Walton is putting r2 = 0; then the normal and tangential equations of motion (1) and (2), become decoupled. Thus, in the present model, with r2 6= 0 and slip conditions, the normal impulse will depend on the tangential impulse. It should also be noted that the discretization leading to Eqs (11 ) and (12) has no counterpart in Walton’s model, instead requiring the calculation of the exact moment of impact. Finally, it should be noted that Walton’s theory is aimed at impact situations, whereas the present algorithm also aims to take care of conditions of sliding contact persisting over several timesteps. 3.3. Energy considerations In this section energy consistency will be discussed, that is, whether the algorithm has the desirable property that energy is either dissipated or conserved, but not created. For the case en = et this question is answered in the affirmative in (Moreau, 1994). For the case en 6= et and a nonspherical body it is possible to construct examples which violates energy consistency. For example, putting r1 = r > 0, r2 = −r, Izz = mr 2 , µ = 1, x˙ − = v0 > 0, y˙ − = −2v0 and θ˙ − = 0 gives such an example. Here it will be verified that the case en 6= et is energy consistent, with the (rather restrictive) assumption r2 = 0. Writing the kinetic energy T at times t − and t + and expressing the change in velocities in terms of the contact impulses applied under this time period, the change in kinetic energy (for r2 = 0) can be expressed as 1T = x˙ − Pt + y˙ − Pn + It is sufficient for energy consistency that
1 2 1 2 Pt + P . 2m1 2m n
(23)
A linear complementarity algorithm for rigid body impact with friction
711
1 2 − Pt 6 0, x˙ Pt +
2m1
1 2 y˙ − Pn + Pn 6 0.
(24)
2m
Inserting (14) and (21) into (24) it is seen that the energy consistency condition is fulfilled, assuming that 0 6 en 6 1 and 0 6 et 6 1. This leaves only the y˙ − > 0 case, which results in zero contact impulses and therefore zero change in kinetic energy.
4. Experimental verification In order to test the predictions of the algorithm described above, an experiment was performed with a bouncing rubber ball, of the type available in toy stores. This experiment is outlined in this section. Note that, in the description of the experimental setup below, measured values are given with more significant figures than actually warranted by the precision of the measurements, and the last digit should not be trusted. For the coefficient of friction an approximate interval is given, giving some indication of the precision of these measurements. The ball, with mass 0.00859 kg and diameter 0.0262 m, was thrown into the space between two surfaces, bouncing first on the lower surface, then on the upper and then a second time on the lower, before exiting in approximately the same direction as it was entering. The mechanism is that the ball starts to rotate rapidly at the first impact. This causes the relative tangential velocity at the second impact to be large enough so that the ball changes direction and returns again in the same general direction as it was thrown from. The process was photographed in stroboscopic light, to produce photographs of the type shown in figure 2, which is the photograph actually used for comparison with computations later in this section.
Figure 2. Photograph from the experiments.
712
L. Johansson
The experimental setup consisted of two wooden tables which were painted black and placed one on top of the other and weighted down by a steel plate. The distance between the upper and lower surface was 0.530 m and the width of the tables were 0.600 m. The ruler seen on the lower table is 0.3245 m long, which provides a length scale. The plane of motion of the ball was a few centimetres in front of the ruler and the camera was placed at a distance of about 1 metre; to compensate approximately for this, a correction of 3% was applied to the length scale implied by measuring the ruler on the photograph. The flight of the ball was photographed from a distance of about 1 metre with a Nikon SLR camera with a Nikkor 35 mm lens set at aperture f/5.6, using Kodak TMAX 400 film. The ball was thrown by hand after releasing the shutter, which was set to 1 s. The setup was illuminated by a Bruel and Kjaer strobe light. The frequency was set to 30 Hz, but a more accurate time interval of 1t = 0.0309 s between flashes was computed by analysing a photograph of a ball in free fall taken for this purpose. The normal coefficient of restitution was determined by measuring the rebound from a 1 m drop and was found to be en = 0.91. The coefficient of friction was measured by pulling four balls glued to an aluminium weight along the surfaces. These measurements were made independently by three different persons, not involved in the computations or other parts of the experiment, in order to avoid any bias. It was found that the dynamic coefficient of friction was µ = 1.13 ± 0.2 for the lower surface and µ = 0.95 ± 0.1 for the upper. The lower surface showed an appreciable difference between static and dynamic friction, the static value being about 1.6. The upper surface showed no significant difference between static and dynamic coefficients of friction. In the computations, the dynamic values cited above were used. To compare computations with experiment, initial conditions were obtained by measuring the two first positions of the ball in figure 2, which are before the first bounce, knowing the time in between. The initial rotational velocity was set to zero in the computations, which is an approximation since the balls were thrown by hand and some rotation was unavoidable. However, the rotational velocity induced by the impacts was very high, the value obtained in the computer simulation after the first impact being 567 rad/s, and any initial rotational velocity was ignored compared to this. This leaves the tangential coefficient of restitution et as the only remaining numerical value needed for the computations. The time interval used in the computations was taken as 0.000309 s, which is one hundredth of the time interval of the strobe light.
Figure 3. Comparison between experiments and computations using et = e = 0.91.
A linear complementarity algorithm for rigid body impact with friction
713
Figure 4. Comparison between experiments and computations using et = 0.68.
Figures 3 and 4 show comparisons between experiment and computations. In these figures the crosses are the experimental positions of the ball as measured from figure 2, and the filled circles shows the computed positions at the same times. Dotted lines are used to connect the measured and computed positions at corresponding times. The circles are numbered 1 through 12, showing the order of the positions of the ball in order of increasing time. In figure 3, the value for the tangential coefficient of restitution are taken as equal to measured value for normal coefficient of restitution, et = en = 0.91, whereas in figure 4 the value of the tangential coefficient is adjusted to et = 0.68 in order to obtain a better fit.
5. Comparison with published experiments In a recent paper (Foerster et al., 1994) a number of experiments involving glass and acetate spheres colliding against each other and against massive slabs are reported. The results are compared with the Walton model, and since this model is equivalent with the present model for the case of spheres, as seen in Section 3 above, these results also provide corroboration for the present algorithm. Since figures comparing the Walton model with experiments are given in (Foerster et al., 1994) this comparison is not repeated here. Next, calculations with the algorithm outlined above are compared to experimental findings by Horak, taken from a secondary source (Goldsmith, 1960). In these experiments a spinning ball is dropped on to a plane surface. The spin velocity is then retarded due to the friction between the ball and the surface. Less obvious, but experimentally observed, is that there is an interval of the initial rotational velocity for which the spin velocity is not only retarded but reversed; the ball rotates in the opposite direction after the impact.
714
L. Johansson
One of the experiments concerns a rubber ball with mass m = 0.1134 kg and radius 0.0359 m, which was dropped from a height of 0.799 m with various initial angular velocities and allowed to drop onto a marble slab, measuring the angular velocity after impact. The coefficient of friction was given as µ = 0.24 and the coefficient of restitution as en = 0.82. It is not stated in (Goldsmith, 1960) if the ball is massive or hollow, but from comparing the mass and the size of the ball and from the fact that rather similar data are given for a tennis ball, it is concluded that the ball had the form of a thin shell. In figure 5, data from this experiment is compared with various computations with the above algorithm. First, a comparison was made between Horak’s experiments and two computations using Horak’s data for µ and en , one using et = 0 and one using et = en . It is seen in figure 5 that neither of these curves is in very close quantitative agreement with the experimental curve, but the computation using et = en agrees qualitatively in that it has a region of initial angular velocities which gives rise to rotation in the opposite direction after impact. Next, the coefficients of friction and restitution were varied to obtain a closer quantitative fit with the experimental data, with a resonable fit being obtained for µ = 0.08, et = en = 0.92. This fit could possibly be improved further by adjusting et independently. It was, however, not possible to obtain a good fit by using Horak’s values for µ and en and adjusting et only. It is, of course, somewhat unsatisfactory that it was necessary to adjust the experimental coefficients of friction and restitution to obtain a good fit between experiments and computation. A plausible explanation for this discrepancy is that the balls in Horak’s experiments are to soft for the rigid body assumptions in this paper to be accurate. Since the balls were similar to tennis balls it seems possible that point contact is to crude an approximation to use in this case. If contact is spread over an appreciable area, it might be necessary to include an impulse moment at the contact. Indeed (Brach, 1989), mentioning Horak’s experiments in passing, indicates that this is the case.
Figure 5. Comparison between Horak’s experiment and computations.
A linear complementarity algorithm for rigid body impact with friction
715
6. A numerical example The problems treated in Sections 4 and 5 above concerns the motion of a single spherical body, resulting in a simplification of the behaviour due to the fact that there is no coupling between the normal and tangential directions. Also in these problems there are no unknown forces, except at the contacts, and therefore the iterative procedure outlined in Section 2.4 is not needed. In this section an, example is given where the flight and impact of solid bar is compared to that of two point masses connected with a spring, and the iterative procedure is needed to resolve the force in the spring. The solid bar has a mass of 2 kg and a length of 2 m. The moment of inertia is 2 kgm2 about the centre of mass, which means that the mass distribution is not uniform, but the mass is concentrated to the ends. The two point masses each has a mass of 1 kg and are connected by a spring with an unloaded length of 2 m and a stiffness of 1000 N/m. The bar and the particle system is released from a height of 20 with an initial velocity of 10 m/s to the right. A gravitational acceleration of 9.81 m/s2 is acting downwards. After some time the end of the bar and the lower particle hits a plane with a coefficient of friction of µ = 0.1 and inclination of 0.1 rad. The coefficients of restitution were et = en = 1 and the timestep in the computation was 0.0008333333 s. In figure 6 the bodies are plotted every one hundredth timestep, illustrating the difference in motion between the rigid body and the particle–spring system after the impact. This difference is further illustrated by the energy plots, for the bar in figure 7 and for the particle–spring system in figure 8. It is seen how part of the energy is transferred to a vibratory mode. Also, from the curves for total energy, it is seen that there are actually a sequence of two impacts for the particle–spring system, not only a single impact as for the bar.
Figure 6. Comparison of motion of bar and particle–spring system.
716
L. Johansson
Figure 7. Loss of energy at the impacts, and transfer of energy between different forms for bar.
Figure 8. Loss of energy at the impacts, and transfer of energy between different forms for particle-spring system.
A linear complementarity algorithm for rigid body impact with friction
717
7. Discussion In this paper, a numerical algorithm is developed for the treatment for rigid body impact against rigid obstacles aiming, in particular, to a reasonable treatment of the case where the tangential contact impulse is important. The modelling does not include a normal and/or tangential coefficient of restitution at the outset, but arise as a consequence of the way the normal contact law and a tangential law of friction are applied. An experiment has been performed to verify the predictions of the proposed algorithm, with encouraging results as described in Section 4. The comparison with a published experiment concerning hollow rather than massive balls, as described in Section 5, are somewhat less encouraging, although a qualitative agreement is achieved. It seems that point contact is to crude an approximation and it is necessary to include an impulse moment to capture this experiment accurately. In order to fit the approach taken in the present paper, it would be necessary to find a constitutive relation for such a moment impulse formulated on complementarity form. Acknowledgement The author wishes to thank Mr. Bo Skoog for his invaluable help with the experiments. References Anitescu M., Cremer J.F., Potra, F.A., 1995. Formulating 3D contact dynamics problems. Reports on Computational Mathematics, No. 80/1995, Department of Mathematics, The University of Iowa. Brach R.M., 1989. Rigid body collisions. J. Appl. Mech. 56, 133–138. Brogliato B., 1996. Nonsmooth Impact Mechanics: Models, Dynamics and Control. Springer, Berlin. Foerster S.F., Louge M.Y, Chang. H., Allia K., 1994. Measurements of the collision properties of small spheres. Physics of Fluids 6, 1108–1115. Goldsmith E., 1960. Impact. Edward Arnold, London. Johansson L., Klarbring A., 1992. The rigid punch problem with friction using variational inequalities and linear complementarity. Mech. Struct. Mach. 20, 293–319. Kilmister C.W., Reeve, J.E., 1966. Rational Mechanics. Longmans, London. Moreau J.J., 1988. Unilateral Contact and Dry Friction. In: Moreau J.J., Panagiotopoulos P.D. (Eds), Nonsmooth Mechanics and Applications, CISM Courses and Lectures No. 302, Springer, Wien, pp. 1–82. Moreau J.J., 1994. Some numerical methods in multibody dynamics: Application to granular materials. Eur. J. Mech. A/Solids 13, 93–114. Murty K.G., 1988. Linear Complementarity, Linear and Nonlinear Programming. Heldermann, Berlin. Newton I., 1687. Philosophiae Naturalis Principia Mathematica. Swedish translation by C.V.L. Charlier, Gleerups, Lund, 1927–31 and Liber, Malmö, 1986. Pfeiffer F., Glocker Ch., 1996. Multibody Dynamics with Unilateral Contacts. Wiley, New York. Stronge W.J., 1990. Rigid body collisions with friction. Proc. Roy. Soc. Lond. A 431, 169–181. Walton O.R., 1993. Numerical simulation of inelastic, frictional particle–particle interactions. In: Roco M.C. (Ed.), Particulate Two-Phase Flow, Butterworth–Heinemann, Stoneham, pp. 884–911.
A linear complementarity algorithm for rigid body impact with friction Lars Johansson Linköping University, Department of Mechanical Engineering, S-581 83 Linköping, Sweden (Received 10 June 1998; revised and accepted 18 November 1998) Abstract – In this paper an algorithm for the impact of rigid bodies against rigid obstacles is developed. Starting from a normal contact law, Coulomb’s law of friction, the equations of motion and a simple discretization scheme, a linear complementarity problem (LCP) is formulated and used to compute the post-impact velocities. The predictions of the algorithm is compared to an experiment, designed such that the tangential impulses at impact are important, with encouraging results. Elsevier, Paris contact / impact / friction / linear complementarity problem
1. Introduction This paper is concerned with a numerical algorithm for the impact of rigid bodies against rigid obstacles, and with its experimental verification. Some problems of this kind can be treated by specifying the quotient between the relative normal velocity of approach and separation (Newton, 1687), that is, by introducing the classical coefficient of restitution. This is sometimes generalized to the so-called Poisson’s hypothesis (Kilmister and Reeve, 1966), in order to avoid some anomalies arising when the theory is applied to rigid bodies. In many cases however, such as in the experiment described in Section 4 below, it is necessary to take both the normal and the tangential impulses at the impact into account to achieve reasonable agreement with observations. This class of problems has been the object of several recent studies, invariably leading to the introduction of one or more additional constitutive parameters, such as the coefficient of friction. Thus in (Brach, 1989), a quotient between normal and tangential impulses are introduced and several bounds, based on physical assumptions, are derived for this quotient, whereas in (Walton, 1993) a model for spheres with three parameters is used which is, for the case of spheres, equivalent to the model in the present paper. In (Stronge, 1990), the division of the impact process into a compression and an expansion phase is analysed, and the problem is treated using a coefficient of restitution relating energies rather than velocities. In this paper, a contact law and Coulomb’s friction law is formulated in terms of velocities and impulses, in a manner similar to that used for quasistatic contact problems in (Johansson and Klarbring, 1992). These laws are combined with the laws of motion and, after discretization, form a Linear Complementarity Problem (LCP) which is solved to obtain the velocities after impact. From a modelling point of view, the approach taken in the present paper is similar to that taken by (Moreau, 1988), which is much concerned with the mathematical setting of the problem, but where a very elegant algorithm for the case of completely inelastic impacts is also given. See also (Moreau, 1994) where granular materials are simulated with an algorithm involving iteration between the contact laws and the equations of motion, applying the contact laws to a linear combination of the pre and post impact velocities.
704
L. Johansson
Figure 1. Geometry of the contact.
It has, of course, been noted before that contact problems lend themselves to formulation in terms of complementarity conditions. Recent work involving LCP based algorithms for rigid body impact include (Pfeiffer and Glocker, 1996) and (Anitescu et al., 1995).
2. Computational algorithm In this section the computational algorithm is described. The treatment of the contacts is described in detail in Sections 2.1–2.3 and an outline of the overall algorithm is given in Section 2.4. 2.1. Equations of motion The motion of a single body in contact with a rigid wall, shown in figure 1, will now be studied. We are assuming that it has been determined that the body studied is violating a contact condition and the subsequent motion is to be determined. The aim is to compute the position, angular orientation, velocity and angular velocity at time t j +1 , assuming that these quantities are known at an earlier time t j . The body is acted upon by unknown contact impulses Pn and Pt , and by additional forces and a moment which are assumed to be known at times t j and t j +1 , so that a suitable approximation to the corresponding impulses Px , Py and M during the time interval can be calculated. In the case of external forces prescribed as depending on the position or velocity of the body itself or another body, these cannot be known beforehand at time t j +1 . This issue is resolved iteratively as described in Section 2.4.
A linear complementarity algorithm for rigid body impact with friction
705
Denoting the components of the velocity of the contact point at time t j by x˙ j and y˙ j and at time t j +1 by x˙ and y˙ j +1 , the equations of motion can be expressed in terms of these motions as j +1
x˙ j +1 = x˙ j +
1 1 1 r1 Pt + Pn + Px + M, ¯ m1 m Izz I
(1)
1 1 1 r2 Pn + Pt + Py + M. (2) ¯ m2 m I I zz Here 1/m1 = (r12 /Izz + 1/m), 1/m2 = (r22 /Izz + 1/m) and 1/I¯ = r1 r2 /Izz , with m denoting the mass and Izz the mass moment of inertia about the centre of gravity of the body. r1 and r2 are the distances between the contact point and centre of gravity as shown in figure 1, which are assumed constant during the time interval [t j , t j +1 ]. y˙ j +1 = y˙ j +
2.2. Contact conditions The contact law in the normal direction is written in terms of the normal velocity and the normal impulse y˙ > 0,
Pn > 0,
yP ˙ n = 0.
(3)
It cannot be assumed that y˙ is continuous. This is obvious for the case of impact, but in fact other cases of discontinuous velocities can occur, as is shown by the Painleve example (Brogliato, 1996). It is assumed, following (Moreau, 1994), that the right and left limits of y˙ exist, and that the normal contact law can be applied to a linear combination of these limits. Thus we put
en en 1 1 y˙ − + y˙ + > 0, Pn > 0, y˙ − + y˙ + Pn = 0, (4) 1 + en 1 + en 1 + en 1 + en where en is a constitutive parameter, to be chosen according to the situation to be modelled. Aiming for a linear complementarity problem formulation, we follow the approach used for quasistatic elastic contact problems in (Johansson and Klarbring, 1992), and introduce the following quantities relevant to the tangential direction (
ψ1 = −Pt + µPn , ψ2 = Pt + µPn
and
(5)
x˙ − |x| ˙ , χ1 = −
2
˙ χ = x˙ + |x| . 2
(6)
2
These definitions are made such that two complementary pairs are formed. Thus if x˙ is positive then χ1 will be zero and, if Pt obeys Coulomb’s law of friction, ψ1 will be equal to 2µPn , and therefore positive. If, on the other hand, x˙ is negative, then χ1 will be equal to the absolute value of x˙ and ψ1 will be zero. Similarly χ2 and ψ2 form a complementarity pair. We thus have: (
ψ1 > 0, χ1 > 0,
ψ1 χ1 = 0,
ψ2 > 0, χ2 > 0,
ψ2 χ2 = 0.
(7)
706
L. Johansson
As is shown by the Painleve example (Brogliato, 1996), it cannot be assumed that the tangential velocity x˙ is continuous either. It is assumed that the right and left limits of x˙ exist and that the conditions (7) remain valid when (6) is replaced by et x˙ − + x˙ + − |et x˙ − + x˙ + | , χ1 = −
2(1 + et )
e x˙ − + x˙ + + |et x˙ − + x˙ + | χ2 = t ,
(8)
2(1 + et )
where et is a constitutive parameter. It should be noted (Moreau, 1994), that the constitutive parameters en and et introduced above will determine the nature of the contact. Considering the case of impact there is a nonzero normal impulse, Pn > 0, and (4) implies y˙ + = −en y˙ − .
(9)
Thus the parameter en can be interpreted as a restitution coefficient, with en = 0 corresponding to purely plastic impact and en = 1 to purely elastic impact. In the tangential direction there can be either stick or slip. If the coefficient of friction is sufficiently large that there is no slip, we have χ1 = χ2 = 0. (8) then gives x˙ + = −et x˙ −
(10)
and et can be interpreted as a tangential restitution coefficient. It should be noted, however, that the interpretation of et as a tangential restitution coefficient is possible only for large values of µ, since the tangential impulse is limited by Coulomb’s law of friction, and in the sliding case, for a nonspherical body, there is a coupling between normal and tangential contact impulses which complicates things further. With the introduction of these restitution coefficients, there is a total number of three parameters describing the impact process; µ, en and et . This is somewhat problematic, since these parameters must be determined experimentally. In particular there is no obvious simple method for measuring et . This increases the scope for adjusting computations to fit a particular experiment. This might be perceived as an advantage, but since the main reason for developing computational algorithms is to be able to make a priori predictions of experimental outcomes, it is actually more likely to be a disadvantage. To overcome this difficulty it would be convenient to have some rule of thumb for the choice of et . To this end the choices et = 0 and et = en were considered. Of these two possibilities et = en should be preferred since et = 0 could not, by any stretch of the imagination, represent the behaviour observed in Horak’s experiments as described in Section 5 below. Also, if bodies with a nonspherical shape is considered, it is possible to construct examples where energy is gained in an impact if en 6= et , as discussed in Section 3.3. Finally, if the linear combinations introduced in Eqs (4) and (8) are interpreted as specifying an intermediate time between the beginning and end of the impact process at which the contact laws are applied, it seems natural to put en = et . 2.3. Discretization and linear complementarity problem For the discretization, a sequence of times [t 1 , . . . , t j , t j +1 , . . . , t n ] covering the desired time interval is introduced. Next, it is assumed that if there is a violation of a contact condition at time t j then the conditions (4) and (7) can be applied approximately with the right and left limits of the velocities at a particular instant
A linear complementarity algorithm for rigid body impact with friction
707
replaced with the velocities at the beginning and end of the timestep. Thus, Eqs (4) and (8) are written: en 1 φ= y˙ j + y˙ j +1 > 0, Pn > 0, 1 + en 1 + en
en 1 y˙ j + y˙ j +1 Pn = 0, 1 + en 1 + en
(11)
et x˙ j + x˙ j +1 − |et x˙ j + x˙ j +1 | χ = − , 1 2(1 + e ) t
(12)
et x˙ j + x˙ j +1 + |et x˙ j + x˙ j +1 | . χ2 =
2(1 + et )
j +1
Using Eq. (12) to eliminate x from Eq. (1) and the result to eliminate Pt from Eqs (2) and (5) and also using the definition of φ in Eq. (11), three equations are obtained, which together with the complementarity conditions (7) and (11), form the following linear complementarity problem: LCP problem Assuming that x˙ j , y˙ j , Px , Py and M are known, find φ, ψ1 , ψ2 , Pn , χ1 and χ2 such that
ψ1 = ψ2
φ
m1 m1 r1 (1 + et )m1 j 1 1 r2 − M Py + x˙ − Px + (1 + en )m (1 + en )Izz I¯ (1 + en )I¯ (1 + en )I¯m m1 m1 r1 j m1 (1 + et )x˙ + M Px + m Izz m m r 1 1 1 j −m1 (1 + et )x˙ − M Px − m Izz 1 1 m1 (1 + et )m1 (1 + et )m1 − 2 − 1 + en m2 I¯ (1 + en )I¯ (1 + en )I¯ Pn m +µ + 1 χ1 m (1 + e ) −m (1 + e ) 1 t 1 t , I¯ χ2 m1 µ− −m1 (1 + et ) m1 (1 + et ) I¯
j y˙ −
φ > 0, Pn > 0,
ψ1 > 0, χ1 > 0, ψ2 > 0, χ2 > 0,
φPn = 0, ψ1 χ1 = 0, ψ2 χ2 = 0.
This problem is recognized as a linear complementarity problem, or LCP, a class of problems which has been the subject of much study, see, e.g. (Murty, 1988). The solution to the linear complementarity problem w = q + Mz, w > 0,
z > 0,
wt z = 0
can be found using a pivot algorithm. Thus, if all elements in q are greater than or equal to zero, w = q and z = 0 is a solution. If not, the equations are rearranged so that one or more of the unknowns in w are exchanged for unknowns in z. This procedure is repeated in some systematic way until an arrangement is found where all elements in q are greater than or equal to zero. The solution is then obtained by putting w = q and z = 0, where, of course, the unknowns occupying w and z are not the same as before the pivots. Specifically, in the computations for this paper, q was scanned in the order q1 , q2 , q3 and as soon as a qi < 0 was found the corresponding unknowns wi and zi were exchanged, and the procedure was repeated until a solution was found.
708
L. Johansson
Once the LCP problem is solved, Pn , Pt , x˙ j +1 and y˙ j +1 can be extracted and used to compute the velocity and rotation of the centre of gravity at time t j +1 . 2.4. The overall algorithm The aim is to compute the position, angular orientation, velocity and angular velocity at time t j +1 , assuming that these quantities are known at an earlier time t j . The motion of a single body can be computed directly with the equations in Sections 2.1–2.3. In this case, the applied loads can be assumed to be known both at times t j and t j +1 and suitable approximations of the corresponding impulses are easily computed. Thus, it is first checked if the body violates any of the contact conditions. If it is, the velocities at time t j +1 are computed by solving the LCP problem in Section 2.3. This calculation proceeds from the position the body had at the end of the previous timestep, accepting without correction a small violation of the contact conditions. If no contact condition is violated the equations of motion are used directly, without contact impulses. If, instead, the motion of a system of bodies interconnected with springs and dashpots is sought, there will be coupling forces which cannot be assumed to be known at time t j +1 , since, in order to incorporate the forces from the connecting springs and dashpots, they must be assumed to depend on positions and velocities at time t j +1 . This difficulty is resolved iteratively. Thus it is first assumed that the coupling forces at time t j +1 are the same as at time t j . The position and velocities at time t j +1 are then computed, making it possible to compute better values for the coupling loads. This process is repeated until the position of each body is within a prescribed tolerance between two consecutive iterations. 3. The special case of impacting spheres This section is concerned with the detailed analysis of the algorithm for the special case r2 = 0, covering spheres, cylinders at right angles with the plane of motion, and also the case when a rod impacts at right angles against a surface. In Section 3.1 below the LCP is analysed for this special case. Then, in Section 3.2 it is shown that the present model is equivalent to a model by Walton, see (Walton, 1993) in this special case. Section 3.3, finally, is devoted to the question of energy consistency of the algorithm. 3.1. The LCP algorithm for the case r2 = 0 In this section, the case r2 = 0 is studied. This results in 1/I¯ = 0 and m2 = m, and the coupling terms Pn /I¯ and Pt /I¯ between the normal and tangential directions disappear in Eqs (1) and (2). Since one of the aims of this section is to compare the present model with the Walton model, and since this model is concerned with what happens at the instant of an impact, the LCP problem is applied to an infinitesimal time interval [t − , t + ]. This also means that the impulses Px , Py and M from external loads become zero in Eqs (1) and (2), if it is assumed that these impulses are derived from forces which are bounded. This also makes sense when studying the energy consistency of the model, since impulses from external forces is a legitimate way for a body to gain energy. In short, we study the behaviour of the model without making the approximation of replacing [t − , t + ] with [t j , t j +1 ], with external forces assumed bounded and with r2 = 0. With these assumptions the first equation in the LCP problem, which is for the normal direction, becomes independent of the remaining two, and can be written together with the corresponding complementarity conditions as 1 φ = y˙ − + (13) Pn , φ > 0, Pn > 0, φPn = 0. m(1 + en )
A linear complementarity algorithm for rigid body impact with friction
709
For the case y˙ − > 0, Eqs (13), (5) and (7) imply Pn = Pt = 0. This is the case when contact is breaking. For the remainder of this section the case y˙ − < 0, giving nonzero contact impulses, will be studied. Equation (13) then gives the normal contact impulse as Pn = −m(1 + en )y˙ − .
(14)
Using (14), the LCP problem takes the form LCP problem, r2 = 0, y˙ − < 0
m1 (1 + et )x˙ − − µm(1 + en )y˙ − (
ψ1 > 0, χ1 > 0,
ψ1 χ1 = 0,
ψ2 > 0, χ2 > 0,
ψ2 χ2 = 0.
−1 χ1 = + m1 (1 + et ) 1 , −1 1 ψ2 −m1 (1 + et )x˙ − − µm(1 + en )y˙ − χ2
ψ1
(15)
(16)
This LCP admits two different types of solutions: stick and slip. Analysing first the case of stick, this occurs when µ is sufficiently large that both elements in the first column matrix to the right in (16) is positive (with y˙ − < 0). The solution is then obtained by putting χ1 = χ2 = 0 in (15), which gives, together with (5) Pt = −m1 (1 + et )x˙ − .
(17)
Turning next to the slip case, it is first noted that slip can be either in the positive or the negative direction, depending on whether x˙ − > 0 or x˙ − < 0, the case x˙ − = 0 being covered by the stick solution discussed above. These two cases behave in exactly the same way (in the r2 = 0 case presently studied, but not in the general case), the only difference being the slip direction, and we choose to study the x˙ − > 0 case. When x˙ − is increased, the second row element of the first column matrix to the right in Eq. (15) will become negative; the LCP solver will rearrange the equations so that ψ2 is exchanged for χ2 . Equation (15) will then take the form
−2µm(1 + en )y˙ − 0 −1 ψ1 χ1 = + . 1 µm(1 + en ) − − 1 x˙ + y˙ χ2 ψ2 m1 (1 + et ) m1 (1 + et )
(18)
A solution satisfying (18) and (16) is obtained by putting χ1 = ψ2 = 0. This requires that both members of the first column matrix to the right are positive, which is ascertained since we are studying the x˙ − > 0, y˙ − < 0 case and since −m1 (1 + et )x˙ − − µm(1 + en )y˙ − < 0
(19)
is necessary for the stick assumption to be abandoned. From (18) and (5) and since y˙ − < 0 we have Pt = µm(1 + en )y˙ − 6 0.
(20)
−m1 (1 + et )x˙ − 6 Pt 6 0
(21)
Equations (19), (20) and (17) give
710
L. Johansson
covering both the stick and slip in the positive direction cases. The case of slip in the negative direction can be given a completely analogous treatment, giving the same condition but with reversed inequalities. Finally, (18) and (8) gives µm(1 + en ) − y˙ , m1 where, again, the case of slip in the negative direction can be given an analogous treatment. x˙ + = x˙ − +
(22)
3.2. Comparison with the Walton model For the case of spheres, it is stated in (Moreau, 1994) that the present model is equivalent to the model for impacting spheres involving three parameters described in (Walton, 1993), and it can be seen from the equations in Section 3.1 that this is indeed the case. The Walton model is also described in (Foerster et al., 1994) where it is used for comparison with experiments. In fact, Eqs (9), (10), (19) and (22) are direct counterparts in the present notation to the relevant equations in (Foerster et al., 1994). In the Walton model, a coefficient of restitution is introduced for the normal direction, equivalent to Eq. (9). In the tangential direction, both a coefficient of restitution and a coefficient of friction (relating impulses rather than forces) are introduced, equivalent to Eqs (10) and (22). The coefficient of restitution is applied when the tangential relative velocity is small and the coefficient of friction when it is large. The switch between these two equations for the tangential direction is when they predict the same tangential velocity, as determined by Eq. (19), although Walton actually begins with slip assumptions and swiches to stick if necessary. It should be noted that the crucial point in specializing the present model to that of Walton is putting r2 = 0; then the normal and tangential equations of motion (1) and (2), become decoupled. Thus, in the present model, with r2 6= 0 and slip conditions, the normal impulse will depend on the tangential impulse. It should also be noted that the discretization leading to Eqs (11 ) and (12) has no counterpart in Walton’s model, instead requiring the calculation of the exact moment of impact. Finally, it should be noted that Walton’s theory is aimed at impact situations, whereas the present algorithm also aims to take care of conditions of sliding contact persisting over several timesteps. 3.3. Energy considerations In this section energy consistency will be discussed, that is, whether the algorithm has the desirable property that energy is either dissipated or conserved, but not created. For the case en = et this question is answered in the affirmative in (Moreau, 1994). For the case en 6= et and a nonspherical body it is possible to construct examples which violates energy consistency. For example, putting r1 = r > 0, r2 = −r, Izz = mr 2 , µ = 1, x˙ − = v0 > 0, y˙ − = −2v0 and θ˙ − = 0 gives such an example. Here it will be verified that the case en 6= et is energy consistent, with the (rather restrictive) assumption r2 = 0. Writing the kinetic energy T at times t − and t + and expressing the change in velocities in terms of the contact impulses applied under this time period, the change in kinetic energy (for r2 = 0) can be expressed as 1T = x˙ − Pt + y˙ − Pn + It is sufficient for energy consistency that
1 2 1 2 Pt + P . 2m1 2m n
(23)
A linear complementarity algorithm for rigid body impact with friction
711
1 2 − Pt 6 0, x˙ Pt +
2m1
1 2 y˙ − Pn + Pn 6 0.
(24)
2m
Inserting (14) and (21) into (24) it is seen that the energy consistency condition is fulfilled, assuming that 0 6 en 6 1 and 0 6 et 6 1. This leaves only the y˙ − > 0 case, which results in zero contact impulses and therefore zero change in kinetic energy.
4. Experimental verification In order to test the predictions of the algorithm described above, an experiment was performed with a bouncing rubber ball, of the type available in toy stores. This experiment is outlined in this section. Note that, in the description of the experimental setup below, measured values are given with more significant figures than actually warranted by the precision of the measurements, and the last digit should not be trusted. For the coefficient of friction an approximate interval is given, giving some indication of the precision of these measurements. The ball, with mass 0.00859 kg and diameter 0.0262 m, was thrown into the space between two surfaces, bouncing first on the lower surface, then on the upper and then a second time on the lower, before exiting in approximately the same direction as it was entering. The mechanism is that the ball starts to rotate rapidly at the first impact. This causes the relative tangential velocity at the second impact to be large enough so that the ball changes direction and returns again in the same general direction as it was thrown from. The process was photographed in stroboscopic light, to produce photographs of the type shown in figure 2, which is the photograph actually used for comparison with computations later in this section.
Figure 2. Photograph from the experiments.
712
L. Johansson
The experimental setup consisted of two wooden tables which were painted black and placed one on top of the other and weighted down by a steel plate. The distance between the upper and lower surface was 0.530 m and the width of the tables were 0.600 m. The ruler seen on the lower table is 0.3245 m long, which provides a length scale. The plane of motion of the ball was a few centimetres in front of the ruler and the camera was placed at a distance of about 1 metre; to compensate approximately for this, a correction of 3% was applied to the length scale implied by measuring the ruler on the photograph. The flight of the ball was photographed from a distance of about 1 metre with a Nikon SLR camera with a Nikkor 35 mm lens set at aperture f/5.6, using Kodak TMAX 400 film. The ball was thrown by hand after releasing the shutter, which was set to 1 s. The setup was illuminated by a Bruel and Kjaer strobe light. The frequency was set to 30 Hz, but a more accurate time interval of 1t = 0.0309 s between flashes was computed by analysing a photograph of a ball in free fall taken for this purpose. The normal coefficient of restitution was determined by measuring the rebound from a 1 m drop and was found to be en = 0.91. The coefficient of friction was measured by pulling four balls glued to an aluminium weight along the surfaces. These measurements were made independently by three different persons, not involved in the computations or other parts of the experiment, in order to avoid any bias. It was found that the dynamic coefficient of friction was µ = 1.13 ± 0.2 for the lower surface and µ = 0.95 ± 0.1 for the upper. The lower surface showed an appreciable difference between static and dynamic friction, the static value being about 1.6. The upper surface showed no significant difference between static and dynamic coefficients of friction. In the computations, the dynamic values cited above were used. To compare computations with experiment, initial conditions were obtained by measuring the two first positions of the ball in figure 2, which are before the first bounce, knowing the time in between. The initial rotational velocity was set to zero in the computations, which is an approximation since the balls were thrown by hand and some rotation was unavoidable. However, the rotational velocity induced by the impacts was very high, the value obtained in the computer simulation after the first impact being 567 rad/s, and any initial rotational velocity was ignored compared to this. This leaves the tangential coefficient of restitution et as the only remaining numerical value needed for the computations. The time interval used in the computations was taken as 0.000309 s, which is one hundredth of the time interval of the strobe light.
Figure 3. Comparison between experiments and computations using et = e = 0.91.
A linear complementarity algorithm for rigid body impact with friction
713
Figure 4. Comparison between experiments and computations using et = 0.68.
Figures 3 and 4 show comparisons between experiment and computations. In these figures the crosses are the experimental positions of the ball as measured from figure 2, and the filled circles shows the computed positions at the same times. Dotted lines are used to connect the measured and computed positions at corresponding times. The circles are numbered 1 through 12, showing the order of the positions of the ball in order of increasing time. In figure 3, the value for the tangential coefficient of restitution are taken as equal to measured value for normal coefficient of restitution, et = en = 0.91, whereas in figure 4 the value of the tangential coefficient is adjusted to et = 0.68 in order to obtain a better fit.
5. Comparison with published experiments In a recent paper (Foerster et al., 1994) a number of experiments involving glass and acetate spheres colliding against each other and against massive slabs are reported. The results are compared with the Walton model, and since this model is equivalent with the present model for the case of spheres, as seen in Section 3 above, these results also provide corroboration for the present algorithm. Since figures comparing the Walton model with experiments are given in (Foerster et al., 1994) this comparison is not repeated here. Next, calculations with the algorithm outlined above are compared to experimental findings by Horak, taken from a secondary source (Goldsmith, 1960). In these experiments a spinning ball is dropped on to a plane surface. The spin velocity is then retarded due to the friction between the ball and the surface. Less obvious, but experimentally observed, is that there is an interval of the initial rotational velocity for which the spin velocity is not only retarded but reversed; the ball rotates in the opposite direction after the impact.
714
L. Johansson
One of the experiments concerns a rubber ball with mass m = 0.1134 kg and radius 0.0359 m, which was dropped from a height of 0.799 m with various initial angular velocities and allowed to drop onto a marble slab, measuring the angular velocity after impact. The coefficient of friction was given as µ = 0.24 and the coefficient of restitution as en = 0.82. It is not stated in (Goldsmith, 1960) if the ball is massive or hollow, but from comparing the mass and the size of the ball and from the fact that rather similar data are given for a tennis ball, it is concluded that the ball had the form of a thin shell. In figure 5, data from this experiment is compared with various computations with the above algorithm. First, a comparison was made between Horak’s experiments and two computations using Horak’s data for µ and en , one using et = 0 and one using et = en . It is seen in figure 5 that neither of these curves is in very close quantitative agreement with the experimental curve, but the computation using et = en agrees qualitatively in that it has a region of initial angular velocities which gives rise to rotation in the opposite direction after impact. Next, the coefficients of friction and restitution were varied to obtain a closer quantitative fit with the experimental data, with a resonable fit being obtained for µ = 0.08, et = en = 0.92. This fit could possibly be improved further by adjusting et independently. It was, however, not possible to obtain a good fit by using Horak’s values for µ and en and adjusting et only. It is, of course, somewhat unsatisfactory that it was necessary to adjust the experimental coefficients of friction and restitution to obtain a good fit between experiments and computation. A plausible explanation for this discrepancy is that the balls in Horak’s experiments are to soft for the rigid body assumptions in this paper to be accurate. Since the balls were similar to tennis balls it seems possible that point contact is to crude an approximation to use in this case. If contact is spread over an appreciable area, it might be necessary to include an impulse moment at the contact. Indeed (Brach, 1989), mentioning Horak’s experiments in passing, indicates that this is the case.
Figure 5. Comparison between Horak’s experiment and computations.
A linear complementarity algorithm for rigid body impact with friction
715
6. A numerical example The problems treated in Sections 4 and 5 above concerns the motion of a single spherical body, resulting in a simplification of the behaviour due to the fact that there is no coupling between the normal and tangential directions. Also in these problems there are no unknown forces, except at the contacts, and therefore the iterative procedure outlined in Section 2.4 is not needed. In this section an, example is given where the flight and impact of solid bar is compared to that of two point masses connected with a spring, and the iterative procedure is needed to resolve the force in the spring. The solid bar has a mass of 2 kg and a length of 2 m. The moment of inertia is 2 kgm2 about the centre of mass, which means that the mass distribution is not uniform, but the mass is concentrated to the ends. The two point masses each has a mass of 1 kg and are connected by a spring with an unloaded length of 2 m and a stiffness of 1000 N/m. The bar and the particle system is released from a height of 20 with an initial velocity of 10 m/s to the right. A gravitational acceleration of 9.81 m/s2 is acting downwards. After some time the end of the bar and the lower particle hits a plane with a coefficient of friction of µ = 0.1 and inclination of 0.1 rad. The coefficients of restitution were et = en = 1 and the timestep in the computation was 0.0008333333 s. In figure 6 the bodies are plotted every one hundredth timestep, illustrating the difference in motion between the rigid body and the particle–spring system after the impact. This difference is further illustrated by the energy plots, for the bar in figure 7 and for the particle–spring system in figure 8. It is seen how part of the energy is transferred to a vibratory mode. Also, from the curves for total energy, it is seen that there are actually a sequence of two impacts for the particle–spring system, not only a single impact as for the bar.
Figure 6. Comparison of motion of bar and particle–spring system.
716
L. Johansson
Figure 7. Loss of energy at the impacts, and transfer of energy between different forms for bar.
Figure 8. Loss of energy at the impacts, and transfer of energy between different forms for particle-spring system.
A linear complementarity algorithm for rigid body impact with friction
717
7. Discussion In this paper, a numerical algorithm is developed for the treatment for rigid body impact against rigid obstacles aiming, in particular, to a reasonable treatment of the case where the tangential contact impulse is important. The modelling does not include a normal and/or tangential coefficient of restitution at the outset, but arise as a consequence of the way the normal contact law and a tangential law of friction are applied. An experiment has been performed to verify the predictions of the proposed algorithm, with encouraging results as described in Section 4. The comparison with a published experiment concerning hollow rather than massive balls, as described in Section 5, are somewhat less encouraging, although a qualitative agreement is achieved. It seems that point contact is to crude an approximation and it is necessary to include an impulse moment to capture this experiment accurately. In order to fit the approach taken in the present paper, it would be necessary to find a constitutive relation for such a moment impulse formulated on complementarity form. Acknowledgement The author wishes to thank Mr. Bo Skoog for his invaluable help with the experiments. References Anitescu M., Cremer J.F., Potra, F.A., 1995. Formulating 3D contact dynamics problems. Reports on Computational Mathematics, No. 80/1995, Department of Mathematics, The University of Iowa. Brach R.M., 1989. Rigid body collisions. J. Appl. Mech. 56, 133–138. Brogliato B., 1996. Nonsmooth Impact Mechanics: Models, Dynamics and Control. Springer, Berlin. Foerster S.F., Louge M.Y, Chang. H., Allia K., 1994. Measurements of the collision properties of small spheres. Physics of Fluids 6, 1108–1115. Goldsmith E., 1960. Impact. Edward Arnold, London. Johansson L., Klarbring A., 1992. The rigid punch problem with friction using variational inequalities and linear complementarity. Mech. Struct. Mach. 20, 293–319. Kilmister C.W., Reeve, J.E., 1966. Rational Mechanics. Longmans, London. Moreau J.J., 1988. Unilateral Contact and Dry Friction. In: Moreau J.J., Panagiotopoulos P.D. (Eds), Nonsmooth Mechanics and Applications, CISM Courses and Lectures No. 302, Springer, Wien, pp. 1–82. Moreau J.J., 1994. Some numerical methods in multibody dynamics: Application to granular materials. Eur. J. Mech. A/Solids 13, 93–114. Murty K.G., 1988. Linear Complementarity, Linear and Nonlinear Programming. Heldermann, Berlin. Newton I., 1687. Philosophiae Naturalis Principia Mathematica. Swedish translation by C.V.L. Charlier, Gleerups, Lund, 1927–31 and Liber, Malmö, 1986. Pfeiffer F., Glocker Ch., 1996. Multibody Dynamics with Unilateral Contacts. Wiley, New York. Stronge W.J., 1990. Rigid body collisions with friction. Proc. Roy. Soc. Lond. A 431, 169–181. Walton O.R., 1993. Numerical simulation of inelastic, frictional particle–particle interactions. In: Roco M.C. (Ed.), Particulate Two-Phase Flow, Butterworth–Heinemann, Stoneham, pp. 884–911.