Unilateral problems of dynamics

Summary Contact processes may be described by local discretizations, by rigid representation or by mixed methods incorporating both ideas. A rigid body approach is proposed for the dynamics of mechanical systems, achieving good results also for multiple-contact problems. Contacts in multibody systems are mainly considered, with the corresponding contact constraints varying with time, thus generating structure-variant systems. The equations of motion for dynamical systems with such unilateral behavior are discussed, solution methods and applications are presented.


Introduction
A large variety of contact con®gurations is encountered in mechanical and biomechanical systems. Contacts may be closed, and the colliding bodies may detach again. Within a closed contact we may have sliding or stiction, both features connected with local friction. If two bodies come into contact, they usually penetrate into each other which leads to local deformations. If contacts are accompanied by tangential forces and by tangential relative velocities within the contact plane, we obtain not only normal but also tangential deformations. De-pending on the dynamical (or statical) environment, contacts may change their state: from closure to detachment, from sliding to stiction, and vice versa. We call a contact active if it is closed or if there is stiction, otherwise we call a contact passive.
Active contacts always exhibit contact forces which, in the general case of normal and tangential deformations, follow from the local material properties of the colliding bodies, from the external (with respect to the contact) dynamics (statics) and from external forces. Considering contacts from this perspective leads to complicated problems of continuum mechanics which, as a rule, must be solved by numerical algorithms like FEM or BEM. This approach is too costly for the treatment of dynamical problems and, in many cases, also inadequate. The rigid-body approach gives quicker results, and it applies better to large dynamical systems.
Under the rigid-body approach we understand a contact behavior characterized, at least in the local contact zone, by lack of deformations and, thus, by rigid-body properties. The contact process is then governed by certain contact laws, like those of Newton and Poisson, and appropriate extensions of them. A fundamental law with respect to the rigid-body model is the complementarity rule, sometimes called corner law or Signorini's law. It states that in contact dynamics either relative kinematical quantities are zero, and their accompanying constraint forces (or constraint force combinations) are not zero, or vice versa. For a closed contact, the relative normal distance and normal velocity of the colliding bodies are zero, and the constraint force in normal direction is not zero, or vice versa. For stiction, the tangential relative velocity is zero, and the constraint force is located within the friction cone, which means that the difference of the static friction force and constraint force is not zero, or vice versa. The resulting inequalities are indispensable for an evaluation of the transitions between the various contact states.
From these properties follows a well-de®ned indicator yielding the transition phases. For normal passive contacts, the normal relative distance (or velocity) of the colliding bodies indicates the contact state. If it becomes zero, the contact will be active; the indicator ``relative distance'' becomes a constraint accompanied by a constraint force. The end of an active contact is then indicated by the constraint force. If it changes sign, indicating a change from pressure to tension, we get detachment, and the contact again transits to a passive state.
In the case of friction, we have in the passive state a nonvanishing tangential relative velocity as an indicator. If it becomes zero, there might be a transition from sliding to stiction, depending on the force balance within or on the friction cone. The indicator``relative tangential velocity'' becomes then a constraint, leading to tangential constraint forces. The contact remains active as long as the maximum static friction force is larger than the constraint force, which means that there is a force balance within the friction cone. If this``friction saturation'' becomes zero, the contact again might go into a passive state with nonzero tangential relative velocity.
If we deal with multibody systems including unilateral constraints there arises the problem of multiple contacts and their interdependences. A straightforward solution of these processes would result in a combinatiorial problem of a huge dimension. Therefore, a formulation applying complementarity rules and the resulting inequalities is a must. The mathematical methods developed in this area in the last twenty years assure nearly always an unambiguous solution, even in the cases of mutually dependent contacts.
The ®eld of multibody dynamics with unilateral contacts should be seen before the background of various research activities: in the area of variational inequalities, in connection with convex energy functions on the one side, [20], [21], [27], [36], and in the area of hemivariational inequalities, in connection with non-convex energy functions, [28]. Both areas are young, the ®rst one being some 35 years and the second only 15 years old. The two scientists who are the most connected with this development were Moreau in Montpellier, France, and Panagiotopoulos in Thessaloniki, Greece. Moreau started already in the seventies to formulate the discontinuous properties of nonsmooth mechanics, [21±23]; not much later, Panagiotopoulos considered inequality problems, [27], which led him consistently to the development of hemivariational inequalities, [28]. Most nonsmooth mechanical problems possess nonconvex features. Nearly all applications, at that time, were of statical or quasi-statical nature.
In Sweden, Lo Ètstedt considered a series of practical examples, [15], [16], and established a school, which in more recent time has been continued very successfully by Klarbring and his group, [11], [12]. Klarbring also treats statical and quasistatical problems with applications to FEM.
The sphere of in¯uence of these scientists extends, since the signi®cance of the area has been noticed. Especially [10] gives an excellent contribution to dry friction, also in connection with frictional impacts. The book [20] presents a proof with respect to the existence of Jean±Moreau impact theory. A very good survey on this topic may be found in [5]. In the meantime, the literature on the subject is increasing considerably, therefore only a few examples are given below.
Papers [13], [14] deal with two-dimensional contacts including dry friction. The resulting linear complementarity problem is solved by a modi®ed simplex-algorithm. The authors of [18], [25] study self-excitations of frictional vibrations. Frictional constraints are described by variational inequalities and evaluated on the basis of measured friction characteristics. The paper [1] discusses, among other topics, also the existence of unambiguous solutions for generalized accelerations concerning frictional contact problems, a question, by the way, which also has been regarded in [15], [16].
Nonsmooth mechanics allows a general theoretical description, but relevant problems must be solved numerically. In spite of many valuable contributions to the numerics of discontinuous systems, [24], the existing algorithms are still extremely time-consuming. As it was indicated above, for plane contacts we get a linear and for spatial contacts a nonlinear complementarity problem. For linear complementarity problems there exist solution procedures, [24], on the basis of extended linear programming theories, and the existence of a solution is assured. This is not the case for nonlinear complementarity. At the time being, various approximations have been suggested, [11], [12], [17], and applied to practical problems, [38] through [44]. The algorithms work, in most cases, convincingly, but the computing times explode because of the iterative character of nearly all solution methods. Therefore, numerical solution of all kinds of complementarity problems is still a topic of current research.
The author's Institute pursues the research in applied and engineering mechanics, and has been involved in unilateral problems since more than ®fteen years, exclusively in connection with multibody dynamics. The fundamentals of the area may be found in books [3], [4], which present the multibody theory on the basis of the projection method, yielding a most ef®cient theory for multibody applications. In the course of the years, and very much due to practical requirements, this type of multibody theory has been combined with inequality formulations to describe the unilateral behavior of multiple contacts, [29±31]. Especially the complementarity properties, in combination with the kinematics and kinetics of the appropriate indicators, have become the main key to a consistent description of multibody systems with multiple unilateral contacts, [31]. An extension to impact, with the background of Moreau's work, became then straightforward, [8], [9]. The corresponding theory has recently and successfully been con®rmed by extensive measurements, [2].
Based on these fundamentals, a variety of industrial applications has been considered. The paper [6] and the dissertation [7] deal with roller chains used often in the automotive industry. Roller chains typically possess some hundred degrees of freedom and hundreds of chain-wheel contacts, where also stick-slip processes take place. Continuous Variable Transmissions (CVT) chains do not have such large amount of degrees of freedom, but they are particularly dif®cult due to spatial contacts within the conic wheel discs. Their ef®ciencies have been evaluated, with an excellent agreement between theory and measurements, [33], [39], [40].
Another area of application is the ®eld of manufacturing assembly processes. Here, during a mating process, contacts are closed and detached, and stick-slip phenomena up to jamming occur. Efforts have been made to describe, in a general way, the relative kinematics of colliding bodies in terms of differential geometry, [19]. The results are convincing, and they are used now in other applications. Another example of manufacturing concerns vibration conveyors, often used to transport, to orient and to select small parts: screws, bolts, electronic components, and the like. Without the theory of multibody systems with unilateral contacts, a treatment of such a system would have been impossible. In the dissertation [44] a theoretical solution of this problem is given and experimentally proven.
In the following, we shall give a survey of the theory including impacts with friction, and we shall present some of its typical applications.

Models
By a model, we mean the description of the mechanical, especially dynamical, behavior of the system under consideration, expressed by laws of kinematics and kinetics. It leads to a mathematical description in the form of (nonlinear) differential equations with inequality constraints, which, as a rule, have to be solved numerically. As we treat multibody systems with unilateral contacts, we shall start with a short presentation of elastic multibody dynamics, proceed to contact kinematics in terms of differential geometry, and explain special features of unilateral contact behavior. Then, equations of motion will be formulated, including the nonsmoothness of the system, and, ®nally, a survey on impact with friction will be given described by a newly developed theory.

Multibody dynamics
The theory of multibody dynamics of interconnected rigid and/or elastic bodies is nowadays well established, and commercialized to a large extend. Most of the existing computer codes are based on the so-called projection method, which considers free motion within hyper-spaces tangential to the constraint surfaces. The Jacobians performing these projections are derivations with respect to the constraint equations, usually formulated on a velocity level, [3], [4], [31]. Since equations of motion for multibody systems represent the basis for all further unilateral considerations, we shall give their short survey here.
We consider a multibody system with f degrees of freedom, which we later shall subdivide into ``rigid'' and ``elastic'' degrees of freedom. Applying the principle of virtual power (Jourdain) we get are the absolute velocity, acceleration of body i in a body-®xed reference system, R f i P R 3 are all forces applied to body i and _q P R f are the generalized velocities; m i is the mass of body i. The velocity and acceleration vectors must be evaluated in a body-®xed frame which yields, see Fig. 1, The following abbreviations are used (index i always omitted): v R v A , a R a A ± absolute velocity, acceleration of the reference base R, expressed in the R-system, x R x IR ± angular velocity of the reference system expressed in the R-system, x 0 R x 0 ± vector from R to the mass element dm, in the undeformed con®guration, r R r el ± displacement vector, Combining Eqs. (1), (2) and the above abbreviations, we come out with the following set of equations of motion: where elastic in¯uences are regarded as shown in Fig. 1. The radius-vector r from an inertial frame I to a mass element dm i of the deformed body B i can be written as where r IR is the vector from I to R, x 0 ± the vector to the mass element in the undeformed con®guration and r el ± the displacement vector. As in most cases of technical relevancy, we assume that the elastic displacements of the components are very small compared to their geometric dimensions. This allows the introduction of a Ritz approach for the elastic deformation, [4], with q el P R n el i and W P R 3;n el i . This well-known superposition of ansatz or shape functions w ji x 0 requires their completeness, [4]. The shape functions may be evaluated by measurements, by FEM-calculations, by analytic approach, or they may be approximated by cubicspline systems. In any case, they should be chosen in accordance with the elastic behavior of the system. Combining Eqs. (3±5) we have to consider the dependencies a q, _ x q, q el;i q, which says that all absolute and elastic accelerations depend on the generalized accelerations q. This property can be expressed by Equations (3) to (7) can now be put into the form After lengthy and tedious calculations, it is always possible to bring Eqs. (8) into a standard form which we need for all further considerations. It writes Equations (9) include all bilateral constraints. They represent the maximum number of minimal (generalized) coordinates, and from this all unilateral contacts are in a state that they do not block any further degree of freedom. If some of the unilateral constraints become active, then the remaining number of degrees of freedom is smaller than f .

Contact kinematics
Geometry and kinematics are the fundaments in establishing models of dynamical systems. In the case of unilateral contacts, this is especially important because magnitudes of relative kinematics serve as indicators for passive contacts and as constraints for active contacts, see Sec. 1. As most of the applications require more or less arbitrary body contours, it makes sense to derive the kinematical contact equations in a general form applying well-known rules of the differential geometry of surfaces, [8], [19], [31]. In practice we ®nd two types of contacts, two-and three-dimensional ones. For the twodimensional case, the contacting bodies lie in a plane thus de®ning a contact line with given direction. Only the sense of direction has to be determined. Problems of that kind are connected with linear complementarity.
For the three-dimensional case, the contacting bodies have a spatial form, and the contact takes place in a plane allowing two tangential directions. The resulting direction for the contact process is not known beforehand, and usually has to be determined iteratively. Problems of that kind are connected with nonlinear complementarity.

Plane contact kinematics
Plane contact kinematics has been presented in the dissertation [8], and since then applied to many practical problems, see also [31]. We start with the geometry of a single body in motion as indicated in Fig. 2. We assume a convex contour, and describe it by a parameter s. Connecting with s a moving trihedral t; n; b and introducing a body-®xed frame B we write

10
For planar contours the binormal B b is constant. Therefore On the other hand, the absolute changes of n and t are given by the Coriolis equation where we must keep in mind that B x IB B X for body-®xed frames B. Putting (11) into (12) we get a coordinate-free representation of the overall changes _ n; _ t _ n Xn À j_ st; _ t Xt j_ sn ; 13 which we can evaluate in any basis. The main advantage of (13) consists of the eliminated, frame-dependent differentiations B _ n and B _ t. In the same manner we proceed with the contour vector r PR . According to (11), (12) we write The velocity v C results from rigid-body kinematics, and corresponds to the velocity of a body-®xed point at the contour. From (16) and (17) we see that Next, we want to derive the absolute acceleration of C by differentiating (17) with respect to time which is not the acceleration of a body-®xed point on the contour. Only the part Planar contour geometry corresponds to such an acceleration, so we can write a C a Q Xt_ s : 22 Later we have to determine the relative velocities of contact points in the normal and tangential directions and their time derivatives. For this purpose we introduce the scalars and state their derivatives as (22), and noting With this basis, we are able to derive the relative kinematics of two bodies such as relative distances, relative velocities and accelerations. Without going into details and refering especially to [31], we summarize the following relationships, see Fig. 3. Potential contact points can be characterized by From each set we need only one equation. The relative distance g N is g N q; t r T D n 2 Àr T D n 1 : 28 Since the normal vectors always point inward, g N is positive for separation and negative for overlapping. Therefore, a changing sign of g N from positive to negative indicates a transition from initially separated bodies to contact. With these equations, and considering Fig. 3, we derive the relative velocities in normal and tangential direction where v C1 ; v C2 are the absolute velocities of the potential contact points C 1 ; C 2 . These velocities may be expressed by the generalized velocities _ q and some Jacobians 32 which we use in the following as a representation of the relative velocities. It may be noticed here that a negative value of _ g N corresponds to an approaching process of the bodies and coincides at vanishing distance g N 0 with the relative velocity in the normal direction before an impact. In the case of a continuous contact g N _ g N 0, the term _ g T shows the relative sliding velocity of the bodies, which we can use to determine the time points of transitions from sliding _ g T T 0 to sticking or rolling _ g T 0. The relevant accelerations follow from a further time differentiation. We get where w N ; w T are given by (32), and w n ; The angular velocities X 1 ; X 2 relate to the two contacting bodies, Figs. 2, 3.

Spatial contact kinematics
Spatial contact kinematics has been presented in the dissertation [19], and since then also has been applied to many practical problems. The situation in this case is, of course, more complex. We still assume that the two approaching bodies are convex, Fig. 4, at least in Fig. 4. Contact geometry of two surfaces that area where contact points might occur. The two bodies are moving with v i , X i ; i 1; 2. For the description of a surface R we need two parameters s and t : r R r R s; t. The tangents s and t, which span the tangent plane at a point of the surface, are de®ned as s or R os ; t or R ot : 36 From these basic vectors, the fundamental magnitudes of the ®rst order are calculated The normalized normal vector n is perpendicular to the tangent plane and pointing outwards We further need the fundamental magnitudes of the second order For a contact point we demand that the normal vector of body 1 n 1 and the distance vector r D are perpendicular to the tangent vectors of body 2 s 2 and t 2 . Thus we obtain four nonlinear equations This nonlinear problem has to be solved at every time step of the numerical integration. After the solution is found, the distance g N between the possible contact points can be calculated as Here, g N is used as an indicator for the contact state. Its value is positive for`no contact', and negative for penetration. The constraints are again formulated on the velocity level, where in the spatial case we have three of them, one in the normal direction _ g N and two in the with v R1 and v R2 being de®ned in an analoguous way as in (30). Differentiating these equations with respect to time leads to the constraints on the acceleration level The time derivatives of the contact point velocities v R1 and v R2 can be written in the form The vectors _ n 1 ; _ s 1 and _ t 1 are determined by the formulas of Weingarten and Gauss, which express the derivatives of the normal vector and of the tangent vectors in terms of the basic vectors The de®nition of the Christoffel symbols C R ab , a; b; r 1; 2, can be found in standard textbooks. Inserting Eqs. (45)±(47) in (43) yields the constraint equations As can be seen, these equations depend only on the Jacobians with respect to the contact points J R1 , J R2 , the basic vectors of the surfaces and the time derivatives of the contour parameters _ The Jacobians are known from the rigid body algorithm, the basic vectors from the surface description. The time derivatives of the contour parameters can be calculated by deriving Eq. (40) with respect to time which means, that the conditions for the contact point should not change while the two bodies are moving. Evaluating Eq. (49), we obtain a system of equations which are linear in the derivatives of the contour parameters This linear problem has to be solved at every time step of numerical integration. Let us summarize the constraint equations in the well known form, by rewriting Eq. (48) as The terms in Eq. (48), which are linearly dependent on q, are collected in the constraint vectors w N ; w S and w T , all the rest is included in the scalars w N ; w S ; w T .

Multiple unilateral contacts
Multiple contacts in multibody systems include a combinatorial problem of large dimensions.
If the state in one contact changes, for example from contact to detachment or from slip to stick, all other contacts are also in¯uenced, which makes a search for a new set of contact con®gurations necessary. In order to not pursue the combinatorial process, we need extended contact laws which describe unambiguously the transitions for the possible contact states, and which generate only consistent contact con®gurations. In a ®rst step, we de®ne all contact sets, which can be found in a multibody system, [31] I A f1; 2; . . . ; n A g with n A elements; I C fi P I A : g Ni 0g with n C elements; I N fi P I C : _ g Ni 0g with n N elements; I T fi P I N : j _ g Ti j 0g with n T elements :

52
These sets describe the kinematic state of each contact point. The set I A consists of n A indices of all contact points. As an example we consider Eq. (9). It belongs to the set combination I A nfI C ; I N ; I T g. The elements of the set I C are n C indices of the unilateral constraints with vanishing normal distance g Ni 0, but an arbitrary relative velocity in the normal direction. In the index set I N there are n N indices of the potentially active normal constraints, which ful®ll the necessary conditions for a continuous contact (vanishing normal distance g Ni 0, and no relative velocity _ g Ni in the normal direction). The index set I N includes, for example, all contact states with slipping. The n T elements of the set I T are the indices of the potentially active tangential constraints. The corresponding normal constraints are closed, and the relative velocities _ g Ti in the tangential direction are zero. The numbers of elements of the index sets I C ; I N and I T are not constant because there are variable states of constraints due to separation and stick-slip phenomena.
As a next step, we must organize all transitions from contact to detachment and from stick to slip and the corresponding reversed transitions. In the normal direction of a contact we ®nd the following situation, [31]: passive contact i g Ni g; t ! 0; k Ni 0; indicator g Ni ; transition to contact g Ni g; t 0; k Ni ! 0; active contact i g Ni g; t 0; k Ni > 0; indicator k Ni ; constraint g Ni 0; transition to detachment g Ni q; t ! 0; k Ni 0 :

53
The kinematical magnitudes g Ni ; _ g Ni ; g Ni are given with Eqs. (28), (31), (33) and (41), (42), (51), where _ g Ni ; g Ni are needed when we go to a velocity or an acceleration level. The constraint forces k Ni must be compressive forces. If they change sign, we get separation. The properties Eq. (53) establish a complementarity behavior which may be expressed by n N (set I N ) complementarity conditions (put on an acceleration level) The variational inequality is equivalent to the complementary conditions (54). The convex set C N fk Ã N : k Ã N ! 0g 56 contains all admissible contact forces k Ã Ni in the normal direction, [31], [41]. The complementarity problem de®ned in Eq. (54) may be interpreted as a corner law, which requires for each contact g Ni ! 0; k Ni ! 0; g Ni k Ni 0. Figure 5 illustrates this property. With respect to the tangential direction of a contact, we shall con®ne our considerations to the application of Coulomb's friction law, which in no way means a loss of generality. The complementary behavior is a characteristic feature of all contact phenomena, independent of the speci®c physical law of contact. Furthermore, we assume that within the in®nitesimal small time step for a transition from stick to slip, and vice versa, the coef®cients of static and sliding friction are the same, which may be expressed by lim _ g Ti 30 l i _ g Ti l 0i : 57 For _ g Ti T 0, any friction law may be applied, see Fig. 6. With this property, Coulomb's friction law distinguishes between the two cases stiction: jk Ti j < l 0i k Ni A j_ g Ti j 0 set I T ; sliding: jk Ti j l 0i k Ni A j_ g Ti j > 0 set I N nI T : 58 Equation (58) formulates the mechanical property that we are for a frictional contact within the friction cone, if the relative tangential velocity is zero, and the tangential constraint force jk Ti j is smaller than the maximum static friction force l 0i k Ni . Then we have stiction.
We are on the friction cone, if we slide with j _ g Ti j > 0. At a transition point, the friction force is l 0i k Ni , see Eq. (57). In addition, we must regard the fact that in the tangential contact plane we may get one or two directions, according to a plane or a spatial contact. From this we summarize in the following way: passive contact i sliding, set I N nI T j _ g Ti j ! 0; jl 0i k Ni j À jk Ti j 0; indicator j _ g Ti j; transition slip to stick j _ g Ti j 0; jl 0i k Ni j À jk Ti j ! 0;  active contact i sticking, set I T j _ g Ti j 0; jl 0i k Ni j À jk Ti j > 0; indicator jl 0i k Ni j À jk Ti j; constraint j _ g Ti j 0; transition stick to slip j _ g Ti j ! 0; jl 0i k Ni j À jk Ti j 0 : 59 From a numerical standpoint of view, we have to check the indicator for a change of sign, which then requires a subsequent interpolation. For a transition from stick to slip one must examine the possible development of a nonzero relative tangential acceleration as a start for sliding. Equation (58) put on an acceleration level can then be written in a more detailed form jk Ti j < l 0i k Ni g Ti 0 i P I T sticking; k Ti l 0i k Ni g Ti 0 i P I N nI T negative sliding; k Ti Àl 0i k Ni g Ti ! 0 i P I N nI T positive sliding :

60
This contact law may be represented by a double corner law as indicated in Fig. 7. To transform the law (60) for tangential constraints into a complementarity condition, we must decompose the double corner into single ones. A decomposition into four elementary laws is given in [8], a decomposition into two elements in [37], [38]. In any case, we come out with a complementarity problem of the form, [31], y Ax b; y ! 0; x ! 0; y T x 0; y; x P R nÃ ; 61 where n Ã n N 4n T in the case of decomposition into four and n Ã n N 2n T for a decomposition into two elementary corners. The quantity x includes the contact forces and one part of the decomposed accelerations, the quantity y the relative accelerations and, in addition, the friction saturation de®ned as the difference of static friction force and tangential constraint force l 0i k Ni À jk Ti j. Equation (61) describes a linear complementarity problem thus being adequate for plane contacts. In the case of spatial contacts, the friction saturation contains the geometric sum of two possible friction directions, leading to a nonlinearity which cannot be solved in a straightforward way. Solution procedures are discussed in [11], [12], [38], [42], [44]. Similar as in the normal case, we can represent the contact law Eq. (60) by a variational inequality of the form, [8], [ The convex set C Ti contains all admissible contact forces k Ã Ti in tangential direction C Ti k Ã Ti jjk Ti j l 0 k Ni ; V i P I T È É : 63 For a derivation of the equations of motion including unilateral effects, we must combine the multibody equations (9) with the unilateral constraints (33) or (51) and (60). In a ®rst step, we include the constraint forces into Eq. (9) keeping in mind that, in a system with additional unilateral constraints, the number of degrees of freedom is variable. To avoid dif®culties with many different sets of minimal coordinates, we take one set of generalized coordinates (Eq. (9) for the combinations (I A nI c ), (I A nI N ), (I A nI T ), see Eq. (52)), and consider the active unilateral constraints as additional constraints which necessarily are accompanied by constraint forces. We include these constraint forces, which are in fact the contact forces, into the equations of motion (9) by a Lagrange multiplier technique.
The constraint vectors w Ni and w Ti in Eqs. . . . .
65 with k Ti t k Ti1 t; k Ti2 t T . In general, the contact forces are time-varying quantities. By the constraint vectors and matrices in Eq. (51), the contact forces can be expressed in the con-®guration space. These forces are then added to Eq. (9) to give For the index sets, see Eq. (52). The contact forces k Ti in Eq. (66) can be passive forces of sticking contacts or active forces of sliding contacts. We express the tangential forces of the n N À n T sliding contacts by the corresponding normal forces, using Coulomb's friction law, by where the coef®cients l i j _ g Ti j of sliding friction may depend on time. The negative sign relates to the opposite direction of relative velocity and friction force. The sliding forces of Eq. (66) in the con®guration space are then A substitution of these forces into Eq. (66) yields the equations of motion with the additional contact forces as Lagrange multipliers. The matrices W N and W T contain components from Eqs. (33) or (51). The matrix H R P R f ;n N of the sliding contacts has the same dimension as the constraint matrix W N . For n T n N ; H R consists of n N À n T columns while the other n T columns contain only zero elements. The relative accelerations of the active normal and tangential constraints in Eqs. (33) or (51) can be combined by means of the constraint matrices (64) in the matrix notation. Together with Eq. (69), we get the system of equations The unknown quantities are the generalized accelerations q P R f , the contact forces in the normal direction k N P R n N and in the tangential direction k T P R 2n T , as well as the corresponding relative accelerations g N P R n N and g T P R 2n T . For the determination of the f 2n N 2n T quantities, we have up to now f n N 2n T equations. In the following, Eqs. (70) will be completed by including the missing n N 2n T contact laws. In general, the kinematic equations are dependent on each other, if there is more than one contact point per rigid body. The situation results in linearly dependent columns of the constraint matrices W N ; W T in Eqs. (70). Such constraints are called dependent constraints.
Before including the contact laws we shall evaluate Eqs. (70) a bit further. Introducing the magnitudes we can write Eq. (70) as Outside the transition events and, thus, for a not changing contact con®guration, the relative accelerations g are zero. In this case, Eqs. (72) have a solution for g and k, which writes The ®rst equation of (73) may also be expressed in the form Ak b 0 with A and b following from (73).
We combine now the equations of motion in the form (70) with the contact Eqs. (55) and (62), which results in the following set containing all unilateral processes in the contacts under consideration: The system (74) is not solvable in this form. Therefore the variational inequalities are transformed into equalities. From this we get a nonlinear system of equations, which represent linear or nonlinear complementarity problems, depending on the type of contact, i.e. plane or spatial contacts. Linear complementarity problem can be solved by algorithms related to linear programming methods, for example Lemke's algorithm, nonlinear complementarity problems require iterative algorithms, [11], [12], [24].

Impact with friction
Impacts with and without friction play an important role in many machines and mechanisms. They have been a research topic at the author's institute since many years, starting with impacts without friction, [29], and proceeding to a theory which describes impacts with friction in a rather general way, [8]. The new theory has been convincingly veri®ed by a large variety of experiments, [2], and will be shortly presented in the following. We assume, as usual, that an impact takes place in an in®nitesimal short time and without any change of position, orientation and all nonimpulsive forces. Nevertheless, we zoom virtually the impact time, establish the equations for a compression and for an expansion phase, and apply these equations again to an in®nitesimal short time interval. The evaluation has to be performed on a velocity level, which we realize by formal integration of the equations of motion, the constraints and the contact laws, Eqs. (74). Denoting the beginning of an impact, the end of compression and the end of expansion by the indices A; C; E, respectively, we get for Here, K NC ; K TC are the impulses in the normal and tangential direction which are transferred during compression, and K NE ; K TE those of expansion. De®ning _ q A _ qt A ; _ q C _ qt C ; _ q E _ qt E , we express the relative velocities as Considering in a ®rst step the compression phase, and combining Eqs. (75) and (76), we come out with where G is called the matrix of projected inertia. It consists of four blocks G NN . . . G TT . Equation (77) allows to calculate the relative velocities _ g NC and _ g TC at the end of the compression phase, depending on the velocities at the beginning of the impact _ g NA and _ g TA under the in¯uence of the contact impulses K NC and K TC . To calculate these impulses, two impact laws in normal and tangential direction are necessary. As already indicated, magnitudes of relative kinematics and constraint forces (here: impulses) are complementary quantities. In the normal direction these are _ g NC and K NC . In the tangential direction, we have the relative tangential velocity vector _ g TC and the friction saturation K TC À diagl i K NC . Decomposing the tangential behavior we obtain, [2], K TCV;i K TC;i l i K TN;i ; _ g TC;i _ g TC;i À _ g À TC;i ; K TCV;i K TCV;i ; K À TCV;i ÀK TCV;i 2l i K NC;i :

78
Together with Eq. (77) this results in a Linear Complementary Problem (LCP) in a standard form y Ax b with x ! 0; y ! 0 and x T y 0 where l is a diagonal matrix, containing the friction coef®cients of the contacts. The problem can be solved numerically. The velocities _ g NC ; _ g TC and the impulses K NC ; K TC are either part of the result or can be obtained by transformation (78) and by K TC K TCV À lK NC . In the compression phase, the kinetic energy of the colliding bodies is transformed into potential energy. During the expansion phase parts of this energy are transformed back. This recovery process is governed by two coef®cients of restitution in normal and tangential directions e N and e T , respectively.
In the case of multiple impacts, Poisson's hypothesis does not guarantee impenetrability of bodies, [8]. Therefore, the impact law is enhanced by an additional condition. The normal impulse during the phase of restitution K NE is, at the minimum, e E K NC in each contact, but can be arbitrary high to avoid penetration. In this case, the bodies remain in contact after the impact. This impact law with a complementary character is drawn in Fig. 8.
In the tangential direction, the impact law is ruled by three effects: at ®rst, a minimum impulse e T K TC must be transferred. On the other hand, the impulse must not exceed the friction limit. Between these limits, a tangential relative velocity _ g TE0 is prescribed. This velocity allows a restitution of the stored energy, if during the phase of restitution sticking occurs, [2]. The existence of _ g TE0 was not introduced in [8], and also not in Moreau's impact model [10], but from the measurements, [2], it could be derived as a necessary correction. With _ g TE0 G TN e N K NC G TT e N K TC ; 80 one can calculate _ g TE0 for all contacts, and the tangential impact law looks like the one drawn in Fig. 9. Here, e N and e T are diagonal matrices containing the different coef®cients for all contacts.
To formulate the equation for the phase of restitution as an LCP, similar to that for the compression phase, the two matrices are introduced. After some transformations similar to those of the compression phase, the LCP writes After solution, the velocities _ g NC ; _ g TC and the impulses K NC ; K TC are either part of the result or can be obtained by the transformations If the impulses during the two phases of the impact are known, one can calculate the motion _ q E of the multibody system at the end of the impact The above impact theory has been veri®ed in [2]. More than 600 impact measurements have been performed. Figure 10 depicts a typical result for the material pair PVC/PVC. The comparison with measurements is excellent, as in many other cases.

Applications
From a large variety of industrial applications, we shall present in the following three typical examples, namely chains, vibration conveyors and chimney dampers. We do not go into details, because all applications have been published elsewhere. Therefore we shall give only a short summary for each case. In all cases, the theoretical methods are based on those of Sec. 2. . Fig. 10. Dimensionless tangential relative velocity after versus before the impact 3.1 CVT chains, [40] The dynamics of chains has been investigated in [6], [7], [39], [40]. The dissertation [7] considers roller chains, where the main problem lies in the many degrees of freedom, usually some hundred, on one side, from which follows a large number of contacts with guides and sprocket wheels, on the other side. An additional problem arises due to the dynamical behavior of the tension devices which requires detailed and complicated models.

____ ____
Another type of chains are the CVT chains applied in continuous variable transmissions, where they transmit power by friction. At the time being, there exist three different chain con®gurations, the rocker pin chains by PIV in Germany, the metal-pushed V-belt chain by van Doorne in the Netherlands and the Borg-Warner chain in the US. We shall focus here on the PIV chain, which was investigated in all details in [40]. Figure 11 represents a typical CVT con®guration. A complete CVT gear consists of two discs with hydraulically shifted sheaves as indicated in Fig. 11. One can change the gear ratio continuously, by opening the discs in one pulley and, at the same time, closing it in the second pulley. Therefore, the continuous variable properties are achieved by axially moving sheaves that are shifted hydraulically. This device also applies to the chain the necessary forces.
The mechanical model for the plane motion of the chain includes altogether 63 links with 189 degrees of freedom. Each pulley possesses two translational and one rotational degree of freedom with respect to rigid body motion. Pitching effects of the axially movable pulley are regarded by additional degrees of freedom. Both pulleys are additionally modeled as elastic bodies, applying a Ritz approach within the context of elastic multibody theory, [3], [4], [40]. Special care has been taken to solve the problem of unilateral contacts in connection with elastically deformed surfaces.
The results have been computed for stationary operation with constant driving speed and an output torque on geared level. Computing time for this case amounts to about eight hours on a SUN workstation. Figure 12 shows the contact forces acting on a pair of rocker pins during one revolution and the tensile force in the related chain link.
As long as the chain link is part of a strand, no contact forces work on its rocker pins. When it comes into contact with one of the pulleys, the pins are pressed between the two sheaves, and hence the normal force increases. Its amplitude depends on the geometry of the sheaves and the transmitted power. The frictional force is a function of this normal force and the relative velocity between the pulley and the pins. It is split into one radial and one circumferential component. The radial contact force coincides with a radial movement of the chain link that equals a power dissipation. In contrast, the circumferential contact force causes the changes of the tensile force in the corresponding chain link, which leads to different tensile force levels in the two strands that agree with the transmitted torque.  Due to the mechanical model, the simulation provides an integrated value of the tensile force in the plates of a chain link, whereas measurements (by G. Sauer and K.T. Renius, Lehrstuhl fu Èr Landmaschinen, Technische Universita Èt Mu Ènchen) were performed for the tensile force in a clasp plate, s. left-hand side of Fig. 13. Therefore, it was necessary to determine the distribution of the tensile force on the plates of the chain links. Using the results of the dynamic simulation shown in Fig. 12, and modeling the pair of rocker pins as bending beams and the plates as linear springs, we obtained the graph of the clasp plate, s. in the right-hand side of Fig. 13. The comparison of simulation and measurement con®rms the mechanical model.

3.2
Vibration conveyor, [44] Vibratory feeders are used in automatic assembly to feed small parts. They are capable to store, transport, orient and isolate the parts. An oscillating track with frequencies up to 100 Hz excites the transportation process, which is mainly based on impact and friction phenomena between the parts and the track. Vibratory feeders are applied for a wide variety of parts and for lots of different tasks. In the majority of cases, the parts are available as a sort of bulk material that is stored in a container. The transportation process, starting in this reservoir, is often combined with orienting devices that orient parts, or select only the parts having already the right orientation. Figure 14 shows an example of a vibratory bowl feeder with an orienting device. Each kind of parts, with its special geometry and mechanical properties, requires an individual adaption of the feeder. This individual tuning comprises the development of suitable track and orienting device geometries and the adjustment of the excitation parameters frequency and amplitude. Due to the complex mechanics of the feeding process this design is usually done by trial and error, without any theoretical background. A complete dynamical model of the transportation process allows a theoretical investigation and, consequently, an improvement in the properties of the feeder. Friction and impact phenomena between the parts and the track are the most important mechanical properties of transportation processes. Consequently, the required dynamical model has to deal with unilateral constraints, dry friction and multiple impacts. The mechanical model of the vibratory feeder can be split in two parts: the transportation process and the base device. The dissertation [44] focuses on modeling and simulation of the transportation process. The modeling of the base device can be done with well-known standard techniques for multibody systems. Friction and impact effects have a fundamental importance for the transportation of parts. Changing contact con®gurations between the parts and the track and also between the parts itself are characteristic for the feeding process. The contacts appear either continuously for a certain time interval, or for a discrete time at impact. A structurevarying multibody system with unilateral constraints with friction is an ideal technique for modeling the feeding process. Its formulation results in a set of differential equations with inequality constraints requiring special mathematical and numerical methods, see Sec. 2.
For the veri®cation of the developed model of the transportation process an experimental vibratory feeder was built, allowing different measurements concerning the impact model and the average transportation rates. Figure 15 shows the principle of the device. The track, ®xed on leaf springs, is excited with an electro-magnetic shaker with a frequency of about 50 Hz. The eigenfrequency of the system is at 52 Hz. The resulting vibration amplitude reaches a maximum value of about 2 mm. The track has an inclination angle a 3 , the angle between the track and the direction of the vibration is b 15 . For accurate contactless measurement of the motion of the transported part, six laser distance sensors were applied. An eddy current sensor was used for the vibration measurement of the track.
For a comparison of the theory and the measurements, the averaged transportation rate was used. Figure 16 presents the result, which looks good before the background of the complexity of the problem. An interesting ®nding is the fact that the averaged transportation velocity does not depend very much on the number of parts and also not on the type of modeling, plane or spatial, [44]. Therefore, the design of vibration conveyors can be carried through considering one part only. For the layout of orienting devices we need, of course, a spatial theory.  Fig. 15. Test setup and part measurement 3.3 Chimney dampers, [34] Towerlike structures, like steel chimneys, may be excited to vibrations by wind vortex streets. Such a mechanism becomes dangerous if the ®rst eigenfrequency of the structure is small enough to be exited by not-too-large wind velocities which appear quite frequently. Oscillations of that type can be signi®cantly damped by application of a simple idea. A pendulum is attached to the chimney top. Its mass and length are tuned to the chimneys ®rst eigenfrequency by classical analysis, and, additionally, optimal viscous damping is also evaluated by classical formulas, given for example by Den Hartog's method. In a second step, the determined optimal damping behavior is not realized by viscous means but by arranging a package of circular plates, which are moved by the end of the pendulum rod through internal circular holes of the plates. The approximate realization of viscous damping by dry friction is accompanied by impacts of the pendulum rod in the plate holes and by stick-slip transitions between the plates. To achieve a best ®t to the optimal viscous case, an optimization of the complete pendulumplate-system is carried through, applying multibody theory with unilateral contacts. Figure 17 illustrates the basic principle. At the top of the chimney, a pendulum is arranged which damps the oscillations. To achieve best damping ef®ciency, the damper is optimized in two steps. In a ®rst step, a classical damper is assumed working with viscous damping. For the resonance frequency area such a damper can be optimized with respect to length, mass and viscous damping of the pendulum. As damping will be realized by dry friction within a plate package, the plates must be selected in such a way that they perform a damping effect coming as nearest as possible to the optimal viscous damping. Therefore, in a second step the plate package will be optimized with respect to damping ef®ciency.
To verify the theory, two types of experiments have been performed. In a ®rst test, a pendulum-plate con®guration has been arranged in a car-like frame with wheels, which could be excited periodically, Fig. 18. The results compare well with the theory.
A second test has been performed with a real steel chimney, which was bent by a steel rope. By releasing suddenly this rope, the chimney starts to vibrate. This process has been measured and compared with the corresponding theory based on the equations of Sec. 2. Figure 19 Fig. 17. Principal con®guration, viscous and plate model, [34] illustrates the comparison between the experiment and simulation, which also con®rms the method presented.
If we approximate the optimal viscous system of Fig. 17 (middle) by a package of plates, we can do that in a best way by optimizing the number of plates, their radii and thicknesses. We ®nd that the number of plates affects the damping ef®ciency only partly. On the left side of Fig. 20, the total mass of the plates was kept constant, also the radius of the top at the ground damper plate. The number of plates was varied from one to eight. When doing so, the distribution of the plate holes was kept linear. In the case of one plate, the diameter of the plate hole was half as big as the largest one in all other cases. Figure 20 shows that the damping mechanism works best with one plate and worst with two plates. The reason lies in the fact, that in the case of two plates only the upper one is moved effectively. The ground plate is ineffective, and hence its damping is too small. As an important result it should be noticed, that the application of a large number of plates does not make sense. Obviously the damping characteristic cannot be improved by using more than ®ve or six plates.  A further sensitive parameter is the distribution of the hole diameters over all plates, called graduation. The right graph of Fig. 20 shows with number 1 the behavior of the system with small, and with graph number 5-with large graduation. The other curves correspond to some intermediate graduations. Keeping the total mass of the plates constant and increasing graduation from 1 to 5 leads to a decrease of damping and to rising chimney displacements for the considered con®guration.

Summary
Classical theory of rigid or elastic multibody system dynamics may be characterized by d'Alembert's principle in the form of Lagrange and extended by Jourdain principle, allowing to project the equations of motion of all interconnected bodies into their free directions according to their individual constraints. These constraints are bilateral in the classical theory and do not include contact phenomena.
If we want to consider dynamical problems with unilateral contacts, especially multibody systems with multiple contacts, and if we further want to model such contacts by unilateral constraints, we must enlarge the multibody theory by certain rules describing unilateral features. One fundamental property, sometimes addressed to as Signorini's law, consists in the fact, that for each contact either quantities of relative kinematics are zero, and the corresponding constraint forces are not zero, or vice versa. This establishes a linear or nonlinear complementarity problem for each of the contacts, which has to be added to the classical multibody formalism. In a more mathematical sense, the nonsmooth character of contact dynamics may be formulated in terms of variational or hemivariational inequalities, depending on the convex or non-convex features of the related energies. Representations of that kind allow for elegant mathematical developments, which however have again to be decomposed for practical applications.