Motion of double pendulum colliding with an obstacle of rough surface

The externally excited and damped vibration of the double pendulum in the vertical plane are considered. The pendulum can collide many times during the motion with a motionless obstacle having the rough surface. The double pendulum colliding with this object has been modeled as a piecewise smooth system constrained by the unilateral constraint. In the relatively long time between the collisions, the differential equations govern the motion of the pendulum. When the contact with the obstacle appears, the pendulum exhibits a discontinuous behavior. An important element of the solving algorithm is aimed on the continuous tracking of the position of the pendulum in order to detect the collision with the unilateral constraints and to determine the state vector of the pendulum at the impact time instant. A single collision is described by the Euler’s laws of motion in the integral form. The equations are supplemented by the Poisson’s hypothesis and Coulomb’s law of friction. The friction law is formulated for the instantaneous values of the reaction forces. The values of their impulses depend on the existence of a slip between the contacting bodies. Therefore, the Coulomb law cannot be generalized for the linear impulses of the forces in a simple way. We have applied the Routh method in order to solve the problem. The method has a simple geometrical interpretation in the impulse space. The angular velocities of both pendulum parts as well as the reaction forces at the joints of the system, which change in stepwise manner, have been presented in the paper.


Introduction
It is well known that in many mechanical engineering systems, a contact between individual elements is not continuous. This may be an intended effect or a result of gradual destruction of the mechanical parts. Regardless of the cause, the collisions between various elements entail significant increase in the dynamic loads what leads to the greater destruction of the colliding elements.
In many papers dealing with collisions in multi-body systems consisting of rigid bodies, the influence of friction forces on the kinematic state of the bodies after impact is omitted. An approach to three-dimensional analysis of the collision of rigid bodies with friction was originally formulated by Routh [1]. In recent years, many researchers rediscover the Routh model [2], analyzing it in detail [4] and also applying it in modeling of relatively simple problems concerning collisions of free bodies [5] as well as collisions in the multi-body systems with impacts [6]. On the base of Routh's model, some more advanced models have been recently developed [3,7,8].
In this paper, the collisions between a double pendulum and a motionless obstacle having a rough surface are investigated. The friction forces accompanying the collision are taken into account, and the Routh model is employed in order to modeling the single collision. This paper is continuation of paper [6]. Mathematical model of the motion of the pendulum which many times collides with the obstacle has been extended here, particularly with regard to the idea of modeling as a piecewise smooth system. Moreover, the exact formulas for coefficients of line of sticking and line of maximum compression which are unique for the considered multi-body system are presented in the paper. Simulations illustrating the application of the proposed model have been carried out to demonstrate the influence of friction forces on the collision what, in effect, has a significant impact on the pendulum behavior in longer time.

Model of the system
The behavior of a double pendulum colliding with an obstacle is considered (see Fig. 1a). The system is constrained to move in a fixed vertical plane. The pendulum, suspended at the point O, consists of two rigid bodies linked by a joint A. One of the bodies is a rod of mass m 1 and length l 1 . The second part of the pendulum is composed of a rod of length l 2 and a disk of radius r . The rod and the disk are connected to each other in a rigid manner. Their common mass is equal m 2 . C 1 and C 2 are the mass centers of the both parts. I 1 and I 2 denote the moments of inertia of the parts with respect to the axes which pass through the points C 1 and C 2 , respectively, and are perpendicular to the plane of motion. There is viscous damping at the both joints with a priori known coefficients c 1 and c 2 . Moreover, two torques M 1 (t) = M 01 cos( 1 t) and M 2 (t) = M 02 cos( 2 t) are taken into consideration.
The angles ψ 1 and ψ 2 are chosen as the general coordinates. The obstacle is motionless, and its surface is plane and rough. We assume that the distance D between the obstacle and the point O fulfills the following conditions So, the motion of the pendulum is constrained by unilateral constraint of the form and only the disk placed at the end of the pendulum may come into contact with the rough obstacle.

Mathematical model
Observe that the unilateral constraint (2) does not act in relatively long intervals. From the point of view of the slow time scale which is suitable to observe the pendulum motion, it becomes active only in some time instants. When the contact with the obstacle occurs, the pendulum exhibits discontinuous behavior. The velocities of the both pendulum parts as well as the reaction forces at the joints and in the contact point with the obstacle undergo rapid changes. The pendulum, the motion of which is constrained by inequality (2), will be treated as piecewise smooth system. The main idea of such approach is presented in Fig. 2. On the right side, we can see the time line involving the whole process. On the time axis, the collisions are represented by several instants t * 1 , t * 2 , t * 3 , . . .. The motion starts from the known initial conditions at time t = 0. All the intervals in which kinematic and dynamic quantities change continuously are separated from each other by the instants t * 1 , t * 2 , t * 3 , . . . in which some of them change abruptly. The pendulum configuration given by ψ 1 (t * − i ), ψ 2 (t * − i ) and its kinematic stateψ 1 (t * − i ),ψ 2 (t * − i ) just before the ith collision decides on the course of the collision. In turn, its kinematic stateψ 1 (t * + i ),ψ 2 (t * + i ) just after the ith collision determines the initial values for the next motion phase when the constraints (1) become inactive.

Motion between collisions
In each relatively long term between the collisions, the pendulum motion is determined by differential equations and initial conditions. The equations of the double pendulum motion derived from the Lagrange equations of the second kind and written in matrix notation are as follows where A = a 11 a 12 a 21 a 22 ,¨ = ψ 1 ψ 2 , Q = q 1 (t, ψ 1 , ψ 2 ,ψ 1 ,ψ 2 ) q 2 (t, ψ 1 , ψ 2 ,ψ 1 ,ψ 2 ) , The configuration of the pendulum at time instant t i and its kinematic state after ith collision are assumed as the initial values for the motion in (i + 1)th interval. Taking into account also the conditions at t = 0, we can write The vectors 00 , 00 are known a priori, while both the vectors 0i , 0i (for i > 0) are determined, using the Routh method, for each of the collision event. Solving equation (3) with respect to the vector¨ we obtain where A −1 is the inverse matrix of A.
The determinant of the matrix A is positive for each pendulum configuration. Hence, the matrix A is regular. Each of the initial value problems (5) with appropriate initial conditions (4) is solved numerically using the fourth order Runge-Kutta method. Usefulness of the method as well as the choice of the suitable time step h have been tested numerically.

Kinematic aspect of the collision
An important element of the algorithm used to solving the initial problem (4) and (5) is the continuous tracking of the pendulum position in order to detect the collision with the obstacle. Accordingly to inequality (2) we define the function f The collision occurs exactly when So, every time when function f changes its sign the Runge-Kutta algorithm is discontinued. Let us suppose that for some interval t j , t j+1 , where t j+1 = t j + h, the relationship respectively. Fulfillment of the condition (9) means that in the found interval t j , t j+1 the collision appears. Let it be the ith collision event from the beginning of the movement. The function f is interpolated linearly in the found interval what allows to determine approximately the instant t * The approximate values of both angles ψ 1 (t * − i ), ψ 2 (t * − i ) as well as the angular velocitiesψ 1 (t * − i ),ψ 2 (t * − i ) are calculated using the linear interpolation. For instance, the value approximating the value of the angle at the instant t * − i is assumed as follows The approximate values of are calculated in the same manner.

Collision laws
Modeling the collision, we assume commonly used assumptions. All displacements and rotations during the collision are negligible, therefore The colliding bodies remain rigid and interact on each other only at one point. In time scale adapted for description of pendulum motion the duration of each collision is very short. However, during the collision the forces as well as the velocities and the accelerations change rapidly and these changes require to be described in a different time scale. We introduce the fast scale time τ being more suitable for the collisions modeling. For each collision, the time τ runs from 0 to τ f i . The impulses of the contact forces at the contact point M as well as the forces at joints O and A are bounded when time τ f i tends to zero. The forces whose impulses are taken into account are shown in Fig. 3. The impulses of the gravity forces, the damping forces and both torques may be neglected due to very short duration of collision. The Cartesian coordinate system Mxyz depicted in Fig. 3 is assumed as the collision reference frame. The axis y along which the angular velocities ω 1 , ω 2 of the pendulum parts are directed is perpendicular to the plane xz.
The Euler laws of motion in the integral form govern the behavior of the pendulum during the collision. Applying the above-mentioned assumptions, we can write , v 2z are components of the velocities v 1 and v 2 of the mass centers C 1 and C 2 , ω 1y =ψ 1 , ω 2y =ψ 2 are components of the angular velocities ω 1 and ω 2 . The angles ψ * 1i = ψ 1 (t * − i ) and ψ * 2i = ψ 2 (t * − i ) determine the pendulum configuration at the beginning of the ith collision.  Fig. 3. For instance, the impulse of the normal force FF Equations (17) and (18) The following commonly known relationships hold that reduce the number of the independent kinematic quantities in Eqs. (13)-(18) from six to two. We choose ω 1y (τ ) and ω 2y (τ ) as the mutually independent quantities and then rewrite the system of Eqs. (13)-(18) where ω 1y (τ ) = ω 1y (τ ) − ω 1y (0), ω 2y (τ ) = ω 2y (τ ) − ω 2y (0). The velocity at the contact point M may be decomposed into two mutually perpendicular components in the following way where i and k are unit vectors of the axes of the collision reference frame shown in Fig. 3. The vector s is called the sliding velocity, and c is the closing velocity. According to the Poisson hypothesis, we distinguish two phases of the collision. They are separated from each other by an instant τ m i at which the closing velocity c is equal zero. The magnitude of the normal force F achieves its maximum at τ m i . The impulses of the normal force F in the first and in the second phase are coupled as followsF where k is the coefficient of restitution, and The Coulomb law has been employed to describe the relationship between the friction force T and normal force F. Regarding that character of this relationship depends on the existence of the slip between the contacting bodies, we write where μ is the friction coefficient, and T = T · i is the projection of the vector T on the common tangent at the contact point M. Equations (21)-(26), written for the instant τ f i , together with Eq. (28) and relationships (30) create the set of collision laws. They allow to determine the kinematic state of the double pendulum after ith collision.

Routh's method and collision of double pendulum
All of the collision laws, except the friction laws, have the global nature. In other words they are written in the integral formulation and concern the whole collision. The law given by (30) is formulated for instantaneous values and due to its dual form cannot be generalized for the impulsesF andT in a simple way, i.e., by integrating. The value of the impulseT at the moment τ f i depends on existence of the slip at the contact point M because the magnitude as well as the direction of the force T depend on the sliding velocity s at the contact point M. In what follows we illustrate how to solve this difficulty applying the Routh method. It should be emphasized that Routh's method is an exact one.
In order to apply this method, we should determine the relations between the sliding and the closing velocities on one hand and the impulses of the forces at the contact point Mon the other hand. The assumption the position of the double pendulum does not change during the collision yields Solving the system of Euler's laws (21)-(26) with respect to the following quantities: ω 1y , ω 2y ,H ,V ,R,P and then substituting the solutions for ω 1y (τ ) and ω 2y (τ ) into Eqs. (32) and (33) we obtain the sought dependencies The constants α i , ε i , γ i are defined as follows where the determinant |A| is defined by Eq. (6). The constants α i , γ i ε i depend on the geometrical and inertial parameters of the pendulum but also on its configuration at the beginning of the collision what means that they must be determined separately for each collision.
The Routh method has a clear geometrical interpretation in the impulse space. In the case of planar motion, this space is two-dimensional. According to the nature of a normal force F, its impulse is always nonnegative. Therefore, we will say further about the impulse semi-plane.   Stepwise changes of mechanical energy of the pendulum So, the disk, having changed the direction of motion, still slides on the obstacle. The change of the sign of the vectors s and T causes that on the line of friction appears the discontinuity in slope. Graphs 4d and 4e concern the situation when the line of friction crosses the line of sticking in the second phase of the collision. For each of the graphs shown in Fig. 4 exists its counterpart that presents an arrangement with the line of friction which is symmetric in respect of the horizontal axis.

Results
Let us consider the motion of double pendulum whose parameters are as follows: m 1 = 3.2 kg, m 2 = 3.8 kg, l 1 = 0.5 m, l 2 = 0.6 m, r = 0.05 m. The mass center of the second pendulum part is distant from the joint A by l c ≈ 0.3627 m. The moments of inertia of the both pendulum parts with respect to their central principal axes are: I 1 ≈ 6.76 × 10 −2 kg m 2 , I 2 ≈ 0.164 kg m 2 . During the motion, the pendulum collides with a fixed obstacle that is placed at the distance D = 0.9 m from the point O. The value of the coefficient of restitution k was assumed as 0.9. It is assumed that do not act the external torques and also there is no damping at the both joints in order to focus on the influence of the collisions on the pendulum motion. The results of two simulations will be presented farther on. These two simulations differ from each other only in the occurrence of friction force between the obstacle and the disk. One of them concerns the case when the pendulum collides with the obstacle having perfectly smooth surface (μ = 0). In the second simulation, the value of the friction coefficient μ = 0.15. The values describing the initial state of the pendulum were assumed for the both simulations the same. They are listed below: Fig. 5, the loss of mechanical energy versus time of motion is shown for the case with the friction (red line) and without the friction (black line). Under the given assumptions (no torques and no damping), the system may lose the energy only as a result of collisions. Indeed, we can see that between the collisions the total energy is constant. The energy decreases stepwise after each contact with the obstacle. For comparison, the energy of the pendulum at the stable equilibrium position ψ 1 = 0 and ψ 2 = 0 is about of −40.007 J. The both graphs shown in Fig. 5 differ from each other. The numbers of the impacts with the obstacle are different. The total energy of the pendulum after the first collision decreased much more when the friction force was taken into account. The comparison, concerning the first collision, is justified because only then, the kinematic state of the pendulum just before the collision is the same for the both considered cases. Greater global energy loss which is observed in the case when the friction force exists is also a result of the bigger number of collisions. It can be seen that the pendulum collides six times in time of simulation for the case μ = 0, and ten times, when μ = 0.15.
In Figs. 6 and 7, the time histories of the angular velocities of the both part of the pendulum are shown. In these figures, the graphs, drawn in black, concern the case when the friction force does not act during the collisions, whereas the graphs drawn in red-concern collisions with the friction force. When it is assumed that the obstacle is perfectly smooth, the angular velocities of the pendulum parts are subject to the discontinuous changes as a result only of the oblique eccentric collision. Appearance of the friction force at the contact point is an additional circumstance which causes the discontinuities of the angular velocities. The visible differences between the time dependencies of angular velocities in both compared cases involve the shape of the graphs and also the number of the discontinuity points (six and ten, respectively).
The vertical P and horizontal R components of the reaction force at the joint O as functions of time are presented in Fig. 8. As previously, the black lines refer to the case without the friction force, and the red lines concern the case when μ = 0.15. The divergence between the graphs (presented in one figure), which is caused by the influence of the friction force, becomes more and more visible after each subsequent collision.
The time histories of the vertical V and horizontal H components of the internal force at the joint A are shown in Fig. 9. The black lines refer to the case without the friction force, whereas the red lines-the case when the friction coefficient μ = 0.15. The differences observed for the both components deepen with time, i.e., with the growing number of collisions.

Conclusions
Models of collision phenomena taking into account friction forces at contact are more suitable for real problems. The obtained results indicate that the friction forces have a significant impact on the collisions. Possibility of the non-sliding type of contact significantly complicates the solving of the problem, even for the relatively simple models with only one point of contact. In the paper, the Routh method has been applied to solve the problem concerning the collision of multi-body system. Combination of the exact method with numerical calculations in the framework of the piecewise smooth model allows us to investigate multiple collisions of the mechanical system taken under consideration. The proposed model employing the Routh method is relatively simple. This approach does not require neither the formulating nor solving the differential equations which describe the collisions. The process of collision is solved on the base of the algebraic formulation.
The numerical simulations have been done using the original software written in Fortran 95.