Yu Wang

Dynamic Modeling and Stability Analysis of Mechanical Systems with Time-Varying Topologies

A model is developed for analyzing mechanical systems with a pair of bodies with topological changes in their kinematic constraints. It is built upon the concept of a Poincare map rather than following the traditional methods of differential equations. The model provides a set of well-defined and naturally-discrete equations of motion and is capable of giving physical insights of dynamic characteristics of deadbeat convergence of multiple collisions and periodic or chaotic responses. The devel opment of a dynamic model and a local stability analysis are presented.

; . G .

1 Introduction Topological changes in a mechanical system often occur when kinematic constraints of the system vary with time. The changes are bound to produce jump discontinuities in the mo tions of bodies in the system, leaving to impulsive reaction and control forces. Examples can be found in the fields of robotics and machine dynamics. Clearances in machine joints and link ages may lead to highly nonlinear or chaotic processes. For a robotic manipulation system for object acquisition, a robot usually interacts through physical contact with an object being manipulated. The object is free to move in a certain range, before being fully controlled by the robot. In both cases the dynamic system being considered undergoes changes in its de grees of freedom, often accompanied by changes in its topol ogy. This paper is concerned with dynamic characteristics of systems with time-varying topologies. Particularly, it focuses on the dynamics of repeated collisions in a system of two planar objects, one of which is under active control.

Since kinematic constraints of a mechanical system generally couple dynamics of its individual bodies, a topological change in the constraints causes changes in the dynamic equations of motion. The dynamic equations are generally of the following form Mij =f(q ,q ,t) + ¢J(q)"A ¢ (q,t ) �0 where q = [qi> Q2, . . • , Qn] T consists of n generalized coor dinates and "A is the Lagrange multiplier vector. The inequality equations describe the kinematic constraints and the differ ential equations are of second order. Exact closed form so lutions for these equations cannot be found, except for simple systems. Hence, numerical methods must be employed.

Investigations of these systems have generally used so called "fitting" methods, in which the equations of motion are ex plicitly expressed between times at which discontinuous events, usually collisions, occur. The discontinuities are treated as boundary conditions to fit into the continuous equations of motion to account for the changes [START_REF] Wehage | Dynamic Analysis of Mechanical Systems With Intermittent Motion[END_REF][START_REF] Haug | Dynamics of Mechanical Systems With Coulomb Friction, Stiction, Impact and Constraint Addition Deletion-Part 1, Theory[END_REF]. These methods require two tasks, automatically to predict the constraint changes and numerically to integrate the differential equations [START_REF] Gilmore | Simulation of Planar Dynamic Mechanical Systems With Changing Topologies: Part !-Characterization and Prediction of the Kinematic Constraint Changes; Part 2-lmplementation Strat egy and Simulation Results for Example Dynamic Systems[END_REF]. Unfortunately, both tasks have many difficulties and the numerical computation is time-consuming. Furthermore, these methods can only be used to show the system dynamics for a specific initial condition. They are not able to anticipate the general dynamic characteristics of the system such as sta bility and regions of stability in a state space.

For systems with small relative distance between the bodies of constraint change, specific methods have been developed to analyze system dynamics [START_REF] Haines | Survey: 2-dimensional Motion and Impact at Revolute Joints[END_REF]. An impact pair model was used to determine specific design parameters of linkages with clearance and intermittent motion mec h anisms (Dubowsky and Freudenstein, 197 1;[START_REF] Lee | On the Dynamics of Intermittent Motion Mechanisms-Part I, Dynamic Model and Response[END_REF]. It is also used to show dynamic response of systems of harmonic excitation [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[END_REF][START_REF] Heiman | Periodic Motions and Bifurcations in Dynamics of an Inclined Impact Pair[END_REF]. However, as opposed to a general anal ysis, these methods are intended for their particular applica tions.

A general approach presented in [START_REF] Zhuravlev | Equations of Motion of Mechanical Systems With Ideal One-Sided Links[END_REF] has also used the impact damper model. It uses a technique of non smooth transformation of generalized coordinates to eliminate the discontinuities in the velocities of bodies with constraint changes. The approach is suited for ideal conservative systems and for some special dissipative systems. It has been expressed in the canonical forms and applied to the analysis of motion of a solid of smooth boundary bouncing on a frictionless plane [START_REF] Ivanov | The Dynamics of Systems with Unilateral Constraints[END_REF]. It is demonstrated that the method is effective to show general properties of system dynamics, such as stability of harmonic and subharmonic motions. Un fortunately, it is difficult to generalize this approach to general mechanical systems with damping and friction, since the energy conservation principle essentially underlines the nonsmooth transformation technique used in the method. In an asymptotic analysis presented in Nagaev ( 197 1), the dynamic description is focused on the events of the discontinuities and equations . . ::

Fig. 1 The physical system of an impact pair of motion are formulated as return maps of the modern theory of dynamical systems. However, the method is applicable to approximate asymptotic behavior only. Both approaches uti lize a key property originally observed by [START_REF] Appell | Traite de Mecanique Ratione/le[END_REF] con cerning changes of system momentum. The approach presented in this paper follows the same observation.

Of particular interest and concern of the present work is the dynamic study of systems in which the topological changes are due to physical contact and collision among a pair of bodies. One reason attributing to this is the need for a general method which could overcome the problem of discontinuity in rep resenting the system dynamics and is efficient for numerical computation. Another motivation for the effort is the impli cation of the method for applications in the field of robotics. Many robotic manipulation tasks require physical interaction of a robot end-effector and the manipulated object. In the presence of collision, the system intermittently changes its char acteristics in response to the actions of the robot. These changes result in difficulties to predict the motions of the object and the collisions raise a large bandwidth of the dynamics to be observed by a sensory-driven feedback control system [START_REF] Buhler | Planning and Control of Robotic Juggling Tasks[END_REF]. A more general method is needed for planning strategies to control the system dynam ics. This is essentially important for challenging robotic tasks requiring catching, throwing, and juggling. Such a method is also of general interest in the practice of automatic design of mechanical part feeders and assembly systems [START_REF] Boothroyd | Statistical Dis tributions of Natural Resting Aspects of Parts for Automatic Handling[END_REF][START_REF] Mason | Mechanics and Planning of Manipulator Pushing Operations[END_REF]Peshkin and Sand erson, 1988).

The purpose of the investigation is to develop a method capable of elucidating the motion characteristics that are due to time-varying topologies. It is hoped that the result of the investigation would provide not only a general dynamic model and efficient computational schemes, but also a deeper un derstanding of the dynamic characteristics. In the author's opinion, the deeper physical insight has not been well reported in the literature and is much needed for system analysis and simulation.

2 Tbe Model of Physical System

The mechanical system to be studied is shown in Fig. 1, which consists of two bodies: a body with strictly convex shape moving in the vertical plane under applied force F and torque T, and a massive plane plate with a specified vertical displace ment. The plate is in theory equivalent to a solid half-space. The model is called an impact pair. It is a simplified version of many typical mechanical connections and robotic manip ulation systems. Because of its simplicity, it has been taken as one of the most frequently used basic models for the inves tigation of mechanical systems with clearances [START_REF] Haines | Survey: 2-dimensional Motion and Impact at Revolute Joints[END_REF]. In spite of its simplicity, the model exhibits the typical behavior which is found in systems containing time-varying kinematic constraints. By studying this model, we shall present a new approach which would shed light on the overall system be havior and gain more physical insight into features of system dynamics. We study the system without contact friction first. Then we discuss the generalization of our approach to frictional systems.

The moving body of mass, m, is subject to applied forces (e.g., gravity) and comes into contact with the massive object.

In the classical theory of rigid body mechanics, its linear and angular velocities undergo jump discontinuities upon colliding the massive body. Energy loss resulting from any impact is modeled by the coefficient of restitution, e, so that the com ponent of relative velocity along the contact normal 1 before and after impact are related by v + = -ev-This is known as Newton's law.

The physical contacts may repeat many times, leading to a finite or an infinite number of collisions. With a contact the system admits one topology; without a contact it exhibits a different topology. The topology varies as a sequence of col lisions take place. We derive analytical equations to describe the system dynamics in the following section.

Discrete Time Formulation

In order to characterize the system dynamics it is important to recognize that the topological changes may cease at a mo ment in a period of time of interest. The object stops bouncing relative to the massive object and the two objects remain con tinuous contact afterwards. When such a situation occurs, the object retains a continuous motion and the velocity discontin uities no longer exist. Hence, the system dynamics can be described by differential equations with algebraic constraints and studied using the methods developed for this class of systems [START_REF] Gear | Ode Methods for the Solution of Differential/ Algebraic Systems[END_REF]. We refer to the phase of motion before the qualitative change as the intermittent phase, and to the phase after as continuous phase.

Therefore, it is desirable that a method for systems with topological changes is not only able to provide a well-suited characterization of motions with multiple collisions, but is also able to identify the conditions for a convergence of the to pological changes. This is of great practical interest in control and planning of robotic manipulation where initially a se quence of collisions are often inevitable but are desired to be eliminated in the course of manipulation. The goal of this section is to develop such a method.

3.1 An Example: A Dropped Tennis Ball. Before car rying out our approach, let's examine a simple example. Con sider a system where a tennis ball is dropped onto a floor. The example appears trivial and may be found in any dynamics textbook; however, it provides a nice illustration of the features of interest. The tennis ball loses mechanical energy at every collision with the floor. At collision n, the impact is governed by Newton's law v,i = -ev;; where Vn is the vertical velocity of the tennis ball at the collision, the superscripts -and + denote the instants just before and just after the collision, and e is the coefficient of restitution. At the next collision (n + 1), t h e velocity of the ball is given as v;; + 1 = -v,i, since mechanical energy is conserved between collisions. Hence, we obtain a difference equation for the consecutive collisions:

( 1) and the elapsed time between the two successive collisions is Tn = -2 ev;; !g where g is the gravity acceleration. These equa tions show that the velocity exponentially decreases as the number of bounces increases, and the process of bouncing ceases after a finite time when the number of bounces n ap proaches infinity. The total elapsed time T from first impact to the end is the sum of the intervals, Fig. 2 The system of the object of smooth geometry and the plane surface where V is the magnitude of the initial approaching velocity v! . After the interval T, the ball will be at rest on the floor.

From this example we see that the system dynamics is de scribed by the discrete representation of the linear difference Eq. ( 1). The intermittent and continuous phases are also quite distinguishable. The behavior of the system in the continuous phase, for this example, is very simple and is governed by a static force balance relation. Similar to this example, the ve locities of the convex body in the general system considered suffer discontinuities at collisions, foiling conventional dif ferential equation approaches. Instead of looking for the con tinuous time history of motion of the system, we shall seek the system state at discrete times, more specifically at a se quence of collision events. In what follows, we shall carry out a theoretical development of the difference equation approach and show that this approach enables us to overcome the draw backs of discontinuities.

3.2 Poincare Maps. The discrete time formulation is ob tained through the mathematical theory of dynamical systems originally proposed by Poincare [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[END_REF]. For an Nth-order continuous-time system, the theory replaces the original system with an (Nl)th-order discrete time system called Poincare map. The replacement reduces the system order and naturally bridges the gap between continuous time and discrete-time systems.

In practice, the Poincare map is often derived by setting a hypersurface of (N -1) dimension in a state space. If the hypersurface is chosen properly, then the state trajectories will repeatedly intersect it. A Poincare map is defined by evolution of the (N -1) variables x at the intersection. Instead of using the continuous trajectory, one can study the system dynamics using the evolution equation

x(n + 1)=f(x(n) ), nan integer
One can also go back to the differential equations to obtain the full continuous time history of the original system if needed.

In order to construct a hypersurface to obtain a correspond ing Poincare map f, one must have a general knowledge of the system behavior beforehand so that the state trajectories indeed cross the hypersurface. Unfortunately, there is no sys tematic method for this construction, although some computer programs provide helpful tools [START_REF] Parker | Insite-A Software Toolkit for the Analysis of Nonlinear Dynamical Systems[END_REF]. In general, for a given system of differential equations, the func tion f is difficult to determine or it can only be determined approximately [START_REF] Guckenheimer | Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields[END_REF]). For our system, however, we can obtain an exact function f, relying on the insight of system behavior. Note that a collision takes place when the two bodies come into a contact and that the bodies will move apart from each other after the collision, provided the convexity that we assumed. If the collision re peatedly occurs, the evolution of collision events could be used to characterize system dynamics. If the collisions cease to exist, the bodies either keep in constant contact or are separated, 0 Fig. 3 An example Poincare surface eliminating the changes in system topology. In this case, the Poincare map for the evolution no longer exists. Therefore, by studying the existence of the Poincare map, we are able to obtain the conditions for qualitative changes in system dy namics. Hence, the Poincare map is a natural choice.

3.3 Quasi-Coordinates and Velocity Discontinuities. In order to seek the state of system at a sequence of collisions, we need to express the geometric conditions of collision ex plicitly. We will limit ourselves to the case of single point contact between the object and the plane surface to avoid the difficulties of multiple point contact. The assumption is sat isfied if the object has a strictly convex shape. The geometric condition for collision can be obtained through a quasi-co ordinate u1 which is the minimum distance between a pair of points on the geometric boundaries of the object and the plane surface, shown in Fig. 2. Note that these points are not fixed relative to the objects themselves if the bouncing object is in motion. This coordinate explicitly determines whether or not the object is in contact with the surface. When they are in contact, u1 = 0; otherwise u1 > 0. In the (x, y) coordinate system embedded in the plane surface, the quasi-coordinate is defined by

uJ(t) =y(t) -h( 8(t)) (2)
where, for the case of the plane surface, y is the coordinate of the center of mass of the object, and h(8) is the projection of the distance from the center of mass to the point of the object boundary on they axis. Note that h(8) is a function of the orientation 8 of the object relative to the surface only, and is independent of the displacements of the object in both x and y directions in this case.

The coordinate u1 can be directly used for defining the Poin care map as U1 = 0 defines the hypersurface y-h(8)=0

(3) For a system shown in Fig. 3, the hypersurface is plotted in the space of coordinates (8, y J. This space is also called the configuration space as it was introduced by Lozano-Perez ( 1983) in the domain of motion planning for robotics. In fact, the hypersurface defined here is the configuration space bound ary of the two bodies. Conceptually, the relative motion of the bodies is described in the configuration space by the motion of a point, which represents their relative kinematic configu rations, relatively bouncing over the configuration space boundary. This provides a different point of view to look at the time-varying topologies.

We must also take into account the velocity discontinuities. For the rigid-body model (Section 2), the velocity discontin uities take place instantaneously and the object displacements remain unchanged during a collision. In other words, in the state space these discontinuities always occur when a trajectory intersects the hypersurface, and the state jumps from one point to another point both on the hypersurface. This leaves two sets of state of intersection on the hypersurface: one for pre collision (t;;, x (t, �)) and the other for post-collision ( t, i, x (t, i )), with n being an integer and x as independent state variable vector, shown in Fig. 4. Therefore, two different Poincare maps can be defined: (1) C governs the discrete trajectory of pre-collision state variables

x -(n+ 1)=f-(x -(n))
where xis an augmented pre-collision state vector, possibly including variable t,, and (2) r + of post-collision state vari ables In this paper we arbitrarily choose the post-collision map, and the set of state vector and the Poincare map are simply denoted as x and f unless they are particularly specified otherwise.

The concept of quasi-coordinate is well known for its special properties in dynamic analysis [START_REF] Appell | Traite de Mecanique Ratione/le[END_REF][START_REF] Neimark | Dynamics of Nonholonomic Systems[END_REF]. The time derivative of a quasi-coordinate is usually a linear function of the derivatives of generalized co ordinates of Lagrange. Our choice of the quasi-coordinate u1 has one particular advantage. It allows us to directly apply Newton's law of impact to determine the velocity discontinuity in u1 as ( 4)

where e is the coefficient of the restitution (Section 2) and u 1 is given as ul(t)=.YUJ-h'(e)iJU) <h'=dhlde)

(5)

If we assume smooth geometric boundaries of the objects, then u 1 (t), which is also called quasi-velocity, has the same order of continuity as y (t) and () (t).

In addition to the quasi-coordinate ul> we take the angular orientation of the object () as the second independent coor dinate. For the frictionless system, there is no tangential impact force, and the collision has no changes in relative motion in the tangential direction. �herefore, we need to express an ad ditional velocity jump in (). In the frictionless case, its expres sion is given as [START_REF] Goldsmith | Impact: The Theory and Physical Behavior of Colliding Solids[END_REF], .

+ . h'(())

.

+ .

e u )-e (r)=h' 2 (())+p 2 (ul (t )-u� (t ) )

(6)

Recall that for our system the Poincare map is a function f to express the state variables right after the n + 1 collision, in terms of their values right after the previous collision n.

This evolution can be decomposed as two parts: a mapping f 1 from the post-collision state at t, i to the next pre-collision state at t, � + 1 and then a mapping f2 from the pre-collision state to the post-collision state at t;i+ I> namely, (7)

using the notation for a composite function. The function fz in fact represents the evaluation of the velocity discontinuities at collision n + 1, which h&ve been described above. Function f1 remains to be derived. where jy is the impulsive force at y direction, Fy is the com ponent of acting force in y direction, Tis the acting torque, and y1 and Yz are the absolute displacements of the object and the plane surface, separately. Note that the impulsive force f Y exists only when an impact occurs. When there is no impact, jy = 0. Since we assume the second object of plane boundary to have a relative large mass, y2 and y2 are considered to be continuous at a collision. We also have the following kinematic relationships:

ii 1 <tl =j i(t) -h'OUl -h" iJ 2 (t)

. Y <t> = u1 u> + h '<e> iJ u>

(10) (11)
We now examine the function fl> which is obtained by in tegrating system equations from t;i to t, �+ 1• Let's consider the motion of body at time t between two successive collisions at f11 and tn+ 1• The displacement and velocity of the coordinates u1 and e are given as, taking those at t;i as initial conditions,

U 1 (t) = U 1 (t, i) + [' U1 (t)dt (
)

Ut(t)=y(t)-h(()(t)) = y ( t, i) + y (t;i) (t-t;i) 8(t )=8(t: )+0(t: )(t-t;i )+ L (L Q2(t)dt) dt ll ll
where =eu: > + e u,; > u-t:

> + L ( L Q2(t )dt) dt n n T(t ) Q2(t ) =--2 mp (17) (18) (19 )
We evaluate these equations at t = t, � +l and, for clarity, we drop superscripts + and -for the continuous quantities, t, u" e, Y2, and Y2• We then obtain u1 (t;;

+ d = u1 (t,i )-(h' (BUn+ 1 ) ) (22) (23)
Up to this point, we have accomplished the derivations for the function f 1, expressing the state at t,� + 1i n terms of the state at t;;-. Using the expressions for the velocity discontinuities in Eqs. ( 4) and ( 6), the entire Poincare function is ready to be completed. Noting that at the collisions the kinematic con straint requires that U1 (t1 1 ) = U1 Un+ 1 ) = 0 and using a convention that variable n representing the time instant right after collision n at t;;and that r(n ) representing the time interval between the collisions n and n + 1, we obtain the system equations in terms of the variables u" u" e' e' t' and r, for the case e > 0, u1(n+1)=0

(24)

u1 (n + 1 ) = -e [u 1 (n )-[h' (8 (n + 1 ))-h' ( e (n ) )J e ( n ) ((n+l )

+ J Ql(t )dt-[y2(t(n+1 ))-y2(t(n ) )] l ( n) (( 1 1 +1 ) J -h' (8(n+ 1)) J Q2(t )dt 1(11) rl ( n+l ) ( B(n+ 1 )=8(n ) + 8 (n )r(n ) + J J Q2(s )dsdt 1(11) 1(11) (25)
(26)

(( I I + I) + J Q2 (t )dt (27)
1( 11) 1( 11) /( 11)

t(n+1)=t(n )+r(n ) (
-fY2 (t(n + 1))-Y2 (t(n ) )] + [u1 (n ) + h' (B(n ) ) iJ (n )]r(n ) -[h(8(n+1))-h(8(n ) )] = O (29)
The difference Eqs. ( 24)-( 28) can be rewritten into explicit difference equations. Moreover, the system is nonautonomous since the variable t ( n ) is explicitly in the equations. We may ignore variable u1 (n ) since Eq. ( 24) always defines the Poincare section and u1 (n ) = 0. Instead of using u1 (n), we introduce t(n ) as a state variable. Then the equations can be rearranged to become an autonomous system of the difference-algebraic equations of t _ he augmented state variable x ( n ) [t(n )u1(n )8(n )8 (n )]T in a vector form:

x (n + 1 ) = f(x (n ) ,r(n ) ,g(x (n ) ,r(n ) ) ) (30) ¢(x(n),r(n),g(x(n ),r(n ) ) ) = 0 (31) where r(n ) is the elapsed time from t(n ) to t(n + 1) and g represents the applied forces on the object and the absolute motions of the plane surface. The vector Eq. (30) represents the difference Eqs. ( 25)-( 28), while the scalar Eq. (31) defines the algebraic Eq. ( 29).

Functions f and¢ are continuous and differentiable. Hence, through the application of the Poincare's method, we have obtained a well-defined discrete formulation. This develop ment will enable us to apply conventional techniques to analyze the system dynamic characteristics without any problems of discontinuity.

Stability Analysis

The difference-algebraic equations derived above are highly nonlinear, and it is difficult to find exact solutions except for simple systems. However, the behavior of solutions for arbi trary initial conditions as n � oo can be characterized by the local stability of equilibrium states (also called fixed points), and the global regions of their stability. The local property shall provide some insights into the dynamic features of the system, whereas the global analysis shows the overall behavior of the nonlinear systems. In this section, we discuss the fixed points and their stability. The global analysis and the numerical simulation are presented in Wang (199 3).

The characteristics of dynamic behavior is related to the properties of function g which represents the applied forces and the motion of the massive plane object. We shall refer to it as the driving junction. Stability of system response to the driving function is the major concern of this section. This is of practical significance, as discussed in Wang ( 1989b) to the development of control strategies for manipulation tasks in some dynamic environments. We study the general case of dissipative system (0 < e < 1). For a perfectly plastic system (e = 0), the problem becomes trivial and after the first collision the object usually continuously contacts the surface; for a perfectly elastic system (e = 1), the collisions become con servative and the system exhibits special dynamic behavior (Wang, 1989a).

4.1 Equilibrium State and Fixed Point. Our difference algebraic equations are in the semi-explicit form with the ex plicit first order vector difference Eq. ( 30) and the scalar Eq. (31). The system is analogous to the system of differential algebraic equations of continuous time model discussed in [START_REF] Gear | Ode Methods for the Solution of Differential/ Algebraic Systems[END_REF]. We may regard the variable 7(n) as an intermediate variable since it does not possess an evolution, and the algebraic Eq. (31) as a constraint on the intermediate variable. Thus, the system is determinate because the number of intermediate variables is equal to the number of constraints [START_REF] Gear | Ode Methods for the Solution of Differential/ Algebraic Systems[END_REF].

A fixed point x* satisfies x* = f(x* ,7* ,g(x* ,7*))

(32) <f>(x*/,g(x*,7*))=0

(33) and the associated time t * = t( oo) as n � oo could be either finite or infinite, depending on whether E117(n) converges.

Let's at this moment exclude possible periodic equilibrium states, for which both harmonic and subharmonic motions must be considered. Then we find the fixed points for this system and correspondingly

* * x1 = t = constant x; = u� =0 x; = o* = constant * . • x4 = 8 =constant 7* =0
These fixed points are not isolated. In a sub-state space of variables x2, x3, and x4, the fixed points are connected and consist of the plane x2 = 0 (Fig. 5). Every point on the plane is a state of equilibrium, representing the state of zero relative normal velocity at the contact point and of arbitrary relative angular position and velocity.

(t ) -Jl(t )-h (X3)Q2(t )-h (x3)X4 b7 = [1 0 0 0] c7 = [0 1 0 0] if 1 * " * I * * II * * 2 d= 2 [QI(t )-Y2(t )-h (x3)Q2(t )-h (x3)x4 J�O
The stability of the solution z (n) = 0 for the linear system is completely determined by the matrix A. The stability criteria are discussed in many books on difference equations, for ex ample, [START_REF] Bernussou | Point Mapping Stability[END_REF]. Here, we summarize some of rel evant results.

(1) The solution for the linear system is asymptotically stable if and only if all the eigenvalues of A are of absolute values less than one.

(2) The solution is not stable if one of the eigenvalues has absolute value greater than one.

•(3) If there are eigenvalues whose absolute values are equal to one and all other eigenvalues have absolute values less than one, then the linear system is stable when the eigenvalues of absolute values equal to one have in dependent eigenvectors.

The third case in which at least one eigenvalue has absolute value equal to one is referred to as the critical case. For the critical case, the local stability for the linearized system cannot be directly carried over to its original nonlinear system. It usually requires higher order terms in the Taylor series to determine the local stability of the nonlinear system (Bemus sou, 1977). However, if the fixed point is not isolated, then we have a degenerate case, for which the stability characteristics for the linearized system may be used to evaluate the stability of the fixed point for the original nonlinear system.

Eigenvalues of A are found by evaluating lt..I -AI =0

The four eigenvalues are real }q = e and A2 = A3 = A4 = 1 and three of them are repeated eigenvalues with unit magni tude, but the multiple eigenvalue has linearly independent ei genvectors. Together with the first independent eigenvector associated with the eigenvalue 1\1 = e, these eigenvectors are

v1=[0 1-e -Q! -{3]7 V2=[1 0 0 0]7 v3= [O o of V4=[0 0 0 1] 7
Noting that the linearized system is in the critical case and that the fixed point is not isolated, we may directly find the solutions for the linearized system by carefully examining the structure of the linear equations:

where

2 c II z1(n) * e ( 1-e) W Z2(n) = Ce 11 2 II t(n) =--. Ce w * * •• * hI * Q ( *) h II ( *) * 2 0 W = Ql (t ) -Y 2 (t ) - (X3 ) 2 t - X3 X4 -;;!:. ( 36
) (37) (38) (39) ( 40
) (41)
and C is a constant to be determined from the initial conditions. In these solutions, the magnitudes of the variables z(n) and r<n) exponentially decay as n -oo. The sum E;:"r<n) of the intervals j(n) converges to a finite value. This means that in a finite time interval the discrete trajectory approaches the fixed point and the repeated collisions come to an end, although an infinite number of collisions take place. Then a continuous contact begins. These properties resemble those of the tennis ball system shown in Eq. ( 1), even though in addition to the bouncing motion z2 this system has angular motions z3 and z4• Since the first eigenvalue satisfies I A.11 = e < 1 and the mul tiple eigenvalue of magnitude equal to one has linearly inde pendent eigenvectors, the connected fixed point is the degenerate critical case. In a strict sense, the system is not asymptotically stable. For a trajectory S1 in the neighborhood of the fixed point, let us perturb the initial state x(O) slightly to its neighbor x(O) + ox(O). Then the new trajectory S2 starting at the neighbor will not converge to the fixed point x* cor responding to the original trajectory. Instead, it converges to another fixed point x* + ox* in its vicinity (Fig. 5). The per turbation ox(O) in the initial state leads to the disc r epancies ox* in the final state variables except for x; which approaches zero as n -oo despite the perturbation.

The degeneracy of the critical fixed point is a result of the physical nature of the system. The system is of second order, but the energy storage and dissipation elements are not com pletely represented if we treat it as a mass-damper-spring sys tem. In fact, there is no spring-like element to resist the motion of the u1 coordinate. This leads to an extra degeneracy that is not reflected in the difference equations since the variable u1 (n) is ignored in the discrete formulation. Similarly, in co ordinate 0 both spring-like and damper-like elements do not exist. The system motion is essentially free in that coordinate, resulting in two more degrees of degeneracy. If we remove the only dissipative element by making the system perfectly elastic, then the fixed points and their local stability will change dra-matically, causing the system to behave in a totally different manner (Wang, 1989a).

In the solutions for the linear system ( 36)-( 40), Z2(n) is the component of the relative velocity at the contact point along the common normal. It is evaluated at the moment of post collision, and must be non-negative, as is the elapsed time r< n). Therefore, it is evident that these solutions are feasible only for a negative w* and a positive constant C. This means that the conditions for successive collisions require that the object accelerate towards the surface in the course of the col lisions, since w* represents the relative acceleration of the object toward the plane surface. This condition is expressed as This condition has an intuitive meaning. If the object is placed on the surface of the massive object and remains in contact, then the object may have a rolling motion of angular velocity iJ {t). In order to maintain the contact, it is necessary that (Wang, 1989a)

Q1 {t)-Y 2 (t) -h 1 [0 (t)]Q2 (t)-h" [0 (t)] 0 2 (t) < 0
This suggests that the condition for local-stable convergence of a sequence of collisions to a continuous contact is identical to the condition for continuous rolling without separation. Of course, this result is not surprising. When W* > 0, the solutions of the linear equations give either negative time intervals T(n) or negative rebounding ve locities z2 ( n). The solutions for the latter attempt to represent the motion of the object coming from the inside of the massive half-space and moving towards the open half space. Once the object touches the boundary, it is bounced into the solid again. This is mirror-symmetric to the realistic situation that the ob ject can move only in the open half-space. The symmetry of the solutions reflects the fact that the equations only represent geometric boundaries, without a notion of solid occupied vol ume. This problem, called viability, can be resolved by simply discarding solutions of nonphysical significance. For general systems, special care may be required to accomplish this (Hei man et al., 1988).

Moreover, there is a special case that we must take care of.

When w* = 0, the second order terms in the Taylor series expansions are needed to determine stability. The following condition is required wu*> >0

However, this condition guarantees that the collision process ceases as n -oo, but the object will not continue to keep contact with the surface afterwards. In fact, contact will break again just when the successive collisions are about to come to an end. This is because the relative acceleration W(t) between the object and the surface becomes positive after the moment t*, and the plane surface slows relative to the object. The contact cannot pull the object, and the object will leave the contact surface right after the moment. In summary, there are three possibilities at an equilibrium state, depending on the value of the relative acceleration w*:

( 1) The equilibrium is unstable ( w* > 0) and it can never be reached; (2) The equilibrium is stable ( w* < 0). In this case, the bouncing motion becomes deadbeat in a finite time interval as an infinite number of collisions converge. The sta bility is called deadbeat stabiliW; (3) The equilibrium is mar ginally stable ( W* = 0 and W* > 0). These three cases are illustrated in Fig. 6. expansions. Let us examine the elapsed time 7(n ). We have found its solution, Eq. ( 40)

7(n ) =Ce"

with a constant C. For a given initial value 7(0), there corre sponds a unique sequence 7(0), 7(1), 7(2), . . . . This sequence can be represented by an iteration process defined by 7(n + 1) = e7(n ) which is described graphically using rectangular coordinates with axes 7(n ) and 7 (n + 1) by two straight lines defined by 7(n + 1) = 7(n ) and 7(n + 1) = eT(n ) as shown in Fig. 7(a). From the initial value 7(0) the sequence eventually ap proaches the origin, and the iteration process converges with the magnitude of the variable 7(n ) reduced by e for each iteration.

For the critical case ( w * = 0), the approximated system must involve the second-order expansions because the first order term w • vanishes. Generally it is cumbersome to derive the second-order or higher-order expansions (Levy and Less man, 1961). In an asymptotic analysis for a closed related problem, [START_REF] Nagaev | The General Problem of Quasi-Plastic Impact[END_REF] proposed the following form of series expansions

"' X; (n ) = � "'x; (n )E"' for i= 1,3,4 m=O "' x2(n ) = � "'xz (n )E"' m� 2 "' 7 (n ) = � "'7(n )E"'
m=O and an addition series expansion for x 2 ( n )

"' x2(n ) = � "'x2 (n )E"'
m=l where E = Af." and A and 'A are constants.

Similarly, after substituting these expressions into the dif ference-algebraic equations and equating coefficients of like powers of e on both sides of the equations, we again obtain a set of first order difference equations for each order m of E"'.

However, unlike the noncritical case, the equations are non linear and their solutions are polynomials of variable E = A f." instead, where the constant 'A is a positive root of the equation

A3 + 2t- 2 -2eA-e= 0
The equation has a single positive root whose value is always greater than the coefficient of restitution e and less than 1 (e < A < 1). The solution for 7(n ) is 7 (n ) = Df." (f.>e )

�(n) � (n) (a) (b)
where D is a constant. The iteration process 7 (n+ l )=A7 (n ) also converges. However, the rate of convergence is slower than that of the noncritical case since A is greater than e, Fig.

7(b)

. We omit the fine details of the analysis and the reader is referred to [START_REF] Levy | Finite Difference Equations[END_REF] and [START_REF] Nagaev | The General Problem of Quasi-Plastic Impact[END_REF] for further discussions.

5 Systems with Contact Friction

In the presence of contact friction between the bodies, the Poincare hypersurface is defined the same as in the frictionless case, but the Poincare map becomes more complicated. During a collision a frictional impulse Px is generated along the com mon tangent of the contact points. If we adopt Coulomb's law of dry friction, the magnitude of the frictional impulse is related to the magnitude of normal impulse Py by a coefficient of friction l.t as follows

The inequality relationship requires further conditions on con tact modes of impact to determine the exact value of the fric tional impulse. The derivations for these conditions are fairly complicated and lengthy, as presented in Wang (1989a); [START_REF] Wang | Two-dimensional Rigid Body Collisions With Friction[END_REF]. Here, we only outline the results. A contact mode of impact describes how the tangential relative motion of the contact points changes during the collision. In the case of Coulomb's model, the contact points may slide relative to each other, or stick, or stick and then reversely slide. Accord ingly, there are three contact modes: sliding, sticking, and re versed sliding. Conditions for these contacts modes are the functions of the pre-collision system state variables. In general, the state space is divided into three regions, called contact regions, each corresponding to a contact mode. Once the con tact mode is identified, a unique solution to the impulse is determined [START_REF] Wang | Two-dimensional Rigid Body Collisions With Friction[END_REF].

One may recall that the Poincare map consists of two parts: f1 and f2, Eq. ( 7). The former describes the motion of the object free from the contact constraint and maps a post-col lision state to a pre-collision state of the next collision. The function f1 is unchanged for the frictional system, except for an additional set of equations for the motion in the x coor dinate. The second function f2 describes the changes in veloc ities at a collision and maps the pre-collision state into the post-collision state. These velocity discontinuities in the fric tional system are substantially different from those for the frictionless system derived in the previous section. Their an alytical expressions for each contact mode are different. There fore, the final difference-algebraic equations are of the form

x(n+l)= lf s(x(n) ,7(n) ,g(x(n) ,7(n))) for sliding contact fc(x(n),7(n),g(x(n),7(n))) for sticking contact f,(x ( n), 7( n) ,g(x ( n), 7( n))) for reversed sliding contact (42) and ct> (x ( n), 7 ( n) ,g(x ( n), 7 ( n)) ) = 0 (43)

where the vector x is the augmented state variable vector as before but including x(n) and x(n). The exact expressions for these functions are described in Wang (l989a) but not presented here due to space limit.

The mapping function f is defined by the functions fs, fc, and f, for each contact mode respectively and f is continuous but piecewise differentiable with respect to the state variables.

In other words, the function has continuous derivatives inside each contact region (i.e., for each contact mode), and the derivatives approach finite values as any point on the boundary of each region is approached from its interior; but the values are different if the boundary point is approached from interiors of different regions which share the boundary.

The piecewise differentiability result in more complexity for the local and global stability analysis. If the fixed point is in the interior of a contact region then the method of local analysis presented in Section 4 can be directly applied, and we would have similar results to those for the frictionless system. How ever, if the fixed point is on a boundary between contact regions, we must determine in which contact region to perform the local analysis. This depends on the interior of the contact region from which the discrete trajectory approaches the fixed point. Consequently, the results for that contact region should be applied. In the global analysis and numerical simulation [START_REF] Wang | Global Analysis and Simulation of Mechanical Systems with Time-Varying Topologies[END_REF], these factors do not cause any problems if the contact regions are borne in mind and carefully coded in nu merical programs.

Conclusions

Mechanical systems with a pair of bodies may establish a sequence of collisions. When the bodies contact and separate, the topology of system changes drastically. A model of impact pair has been developed to describe the system dynamics by a discrete time representation. A function of point mapping or a Poincare map is defined to explicitly describe an evolution process of the multiple collisions. The resulting difference algebraic equations are well-defined and suited for dynamic analysis and numerical simulation. A local stability analysis has shown physical insight to the deadbeat convergence of multiple collisions and the characteristics of the topological changes. The aspects of global analysis and numerical com putation are the subjects of another study [START_REF] Wang | Global Analysis and Simulation of Mechanical Systems with Time-Varying Topologies[END_REF].

This paper is confined in the domain of planar bodies with one being represented by a massive plane surface. The pre sented approach, however, can be applied to a broader class of systems such as a system of objects both with curved bound aries as discussed in (Wang, 1989a). We have emphasized the deadbeat behavior of a series of collisions in the presence of material inelasticity and contact friction. It is possible that the system would exhibit periodical or chaotic responses to an excitation of the massive body. This will be discussed else where.

  of the contact points is a vertex, the contact normal is defined as the normal of the other body's surface. We do not consider the case of two vertices in contact.

Fig. 4

 4 Fig. 4 Velocity jumps on the Poincare hypersurface

3. 4

 4 Difference-Algebraic Equations. The system dy namic equations in the relative coordinate system (y, ()) are given as mji ( t) = m ( Y1 ( t) -Yz ( t)) = Fy ( t) -mji 2 ( t) + jy ( t)(8) mp 2 0(t) = T(t) -h' (())jy(t)

  After substituting the above kinematic and dynamic equations into them, the expressions become utCt>=. Y<tl-h' <eu>>iJu> =y(t, i)+ r (Ql(t)-jiz(t))dt-h'(()(t)) In = u1 (t;i)-(h' (e (t))-h' (e (t; i))) iJ(t, i)

  )dsdt+ y2(t(n ) )r(n )

Fig. 5

 5 Fig. 5 Discrete trajectories and fixed points

4. 2

 2 Steady-State Solutions and Their Local Stabil ity. Having found the fixed point, we now study its local stability characteristics. Here, stability or instability will be taken in the sense of Lyapunov. The local stability of a fixed point x* is represented by that of the linearized system equa tions. Let z = [z1 z2 z3 z4]7 be the perturbation for the state vector, and!; be the perturbation for variable 7(n), respec tively,x(n)=x*+z(n)7(n) =7* + !;(n)When discussing the local stability, the perturbation variables z and !; are assumed to be small, and when substituting these into Eqs. (30) and (31) and expanding their right-hand side functions around the fixed point x *, and 7 *, the nonlinear terms in the Taylor series expansion are usually neglected. Thus, we obtain a system of linear equations: z(n+ 1)=Az(n) +b!;(n) c7z(n) +d!;(n) =0

4. 3 Fig. 6

 36 Fig. 6 Possible cases of repeated collisions

Fig. 7

 7 Fig. 7 Convergence of variable r{n)

Acknowledgments I would like to thank Matt Mason and Randy Brost for their discussions in the development and presentation of this paper. This work was supported in part by the System Development Foundation and by the National Science Foundation under grant DMC-8520475.