Paulo Flores

Remco Leine

Christoph Glocker

Application of the nonsmooth dynamics approach to model and analysis of the contact-impact events in cam-follower systems

Keywords: Nonsmooth dynamics, Augmented Lagrangian approach, Moreau time-stepping method, Contact-impact analysis, Multibody systems

de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The field of multibody system dynamics has its root in classical and analytical methods of dynamics to meet the growing demands in modeling and simulation of complex and advanced mechanical systems in industry and engineering. Multibody systems are ubiquitous in many fields of application, such as aerospace, automotive systems, circuit breakers industry, bipedal locomotion, robotics, biological engineering, and computer graphics [START_REF] Wittenberg | Dynamics of Systems of Rigid Bodies[END_REF][START_REF] Nikravesh | Computer-Aided Analysis of Mechanical Systems[END_REF][START_REF] Haug | Computer-Aided Kinematics and Dynamics of Mechanical Systems-Volume I: Basic Methods[END_REF][START_REF] Shabana | Dynamics of Multibody Systems[END_REF][START_REF] Huston | Multibody Dynamics[END_REF]. Their numerical simulation has become a crucial step not only for the virtual prototyping process in industry, but also in academic fields like global behavior of complex systems, control and stability, in which it is impossible to push forward the studies without reliable simulation software packages [START_REF] Rulka | SIMPACK-a computer program for simulation of large motion multibody systems[END_REF][START_REF] Ryan | ADAMS-Multibody System Analysis Software[END_REF][START_REF] Jiménez | COMPAMM-A Simple and Efficient Code for Kinematic and Dynamic Numerical Simulation of 3-D Multibody System with Realistic Graphics[END_REF]. Numerical simulation must in turn rely on suitable mathematical models. In particular, several points of unilateral contact usually exist in such systems, and impact phenomena and friction are extremely important features in most of multibody systems. As a consequence, multiple impacts may occur quite frequently and become a key point for the numerical simulation of multibody systems. Hurmuzlu and Marghitu [START_REF] Hurmuzlu | Rigid body collision of planar kinematic chain with multiple contact points[END_REF] studied the contact problem in multibody systems, where a planar rigid-body kinematic chain undergoes an external impact and an arbitrary number of internal impacts. Based on the Keller's work [START_REF] Keller | Impact with friction[END_REF], they developed a differential-integral approach and used different models for coefficient of friction. Han and Gilmore [START_REF] Han | Multi body impact motion with friction analysis, simulation, and validation[END_REF] proposed a similar approach, using an algebraic formulation of the equations of motion, the Poisson's model of restitution and the Coulomb's law to define the tangential motion. Different conditions that characterize the motion (slipping, sticking, and reverse sliding) were detected by analyzing velocities and accelerations at the contact points. Han and Gilmore confirmed their simulation results with experiments for two-body and threebody impacts. Based on a canonical form of the equations of motion, Pereira and Nikravesh [START_REF] Pereira | Impact dynamics of multibody systems with frictional contact using joint coordinates and canonical equations of motion[END_REF] presented a methodology that solves this problem in the context of multibody dynamics impact. Haug et al. [START_REF] Haug | Dynamics of mechanical systems with coulomb friction, stiction, impact and constraint addition deletion-I theory[END_REF] solved directly the differential equations of motion by using the Lagrange multiplier technique. Newton's model was used for impact while the Coulomb's law was used for friction [START_REF] Flores | Modeling and analysis of rigid multibody systems with translational clearance joints based on the nonsmooth dynamics approach[END_REF].

Mechanical systems are often modeled as multibody systems with some degree of non-smoothness. Typical examples in mechanics are the noise and vibration produced in railway brakes, impact print hammers, percussion drilling machines, or chattering of machine tools. These effects are due to the nonsmooth characteristics such as clearances, impacts, intermittent contacts, dry friction, or a combination of these effects [START_REF] Tian | A new model for dry and lubricated cylindrical joints with clearance in spatial flexible multibody systems[END_REF][START_REF] Erkaya | Investigation on effect of joint clearance on dynamics of four-bar mechanism[END_REF][START_REF] Ishida | Nonlinear resonances and self-excited oscillations of a rotor caused by radial clearance and collision[END_REF][START_REF] Flores | Spatial rigid-multi-body systems with lubricated spherical clearance joints: modeling and simulation[END_REF][START_REF] Machado | Development of a planar multi-body model of the human knee joint[END_REF][START_REF] Brutti | Modeling 3D revolute joint with clearance and contact stiffness[END_REF]. In nonsmooth systems, the time evolution of the displacement and the velocities is not requested to be smooth. Due to the possible impacts, the velocities are even allowed to undergo jumps at certain time instances in order to fulfill the kinematical restrictions.

Contact-impact behavior strongly depends on the material properties of the colliding surfaces, the nature of the contact-impact problem, and the level of contact-impact forces/impulses produced. Therefore, the investigation on the contact-impact field is one of the most challenging and demanding issues the engineering. In addition, contact-impact events can fre-quently occur in multibody systems and in many engineering applications the function of mechanical systems is based on them. Common examples can be described by the contact between tire and road in vehicles, wheel and rail in railway systems, contact in robotics and grasping machines, cam and follower mechanisms, just to mention a few [START_REF] Kim | A two-dimensional tire model on uneven roads for vehicle dynamic simulation[END_REF][START_REF] Pombo | Application of a wheel-rail contact model to railway dynamics in small radius curved tracks[END_REF][START_REF] Miller | From robotic hands to human hands: a visualization and simula-tion engine for grasping research[END_REF][START_REF] Ciocarlie | Grasp analysis using deformable fingers[END_REF][START_REF] Miller | Implementation of multirigid-body dynamics within a robotic grasping simulator[END_REF]. As a result of an impact, the values of the system state variables change very fast, eventually looking like discontinuities in the system velocities. The knowledge of the peak forces developed in the impact process is very important for the dynamic analysis of multibody systems and has consequences in the design process. Thus, the selection of the most adequate contactimpact method used to describe the process correctly plays a key role in the accurate design and analysis of these kinds of systems [START_REF] Flores | On the continuous contact force models for soft materials in multibody dynamics[END_REF][START_REF] Flores | Modeling and simulation of wear in revolute clearance joints in multibody systems[END_REF][START_REF] Mukras | Analysis of planar multibody systems with revolute joint wear[END_REF]. In a broad sense, the different methods to solve the impact problem in multibody mechanical systems are continuous and discontinuous approaches [START_REF] Lankarani | A contact force model with hysteresis damping for impact analysis of multibody systems[END_REF][START_REF] Flores | On the contact detection for contact-impact analysis in multibody systems[END_REF][START_REF] Flores | A parametric study on the dynamic response of planar multibody systems with multiple clearance joints[END_REF]. Within the continuous approach, the methods commonly used are the continuous force model, which is in fact a penalty method, and the unilateral constraint methodology, based on complementarity approaches [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF][START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth Systems[END_REF][START_REF] Potra | A linearly implicit trapezoidal method for integrating stiff multibody dynamics with contact, joints, and friction[END_REF][START_REF] Dimitrakopolous | Analysis of a frictional oblique impact observed in skew bridges[END_REF]. The compliant continuous contact force models, commonly referred as penalty methods, gained significant importance in the context of multibody systems with contacts thanks to their computational simplicity and efficiency. In these models, the contact force is expressed as a continuous function of penetration between contacting bodies. However, one of the main drawbacks associated with these force models is the difficulty to choose the contact parameters such as the equivalent stiffness or the degree of nonlinearity of the penetration, especially for complex contact scenarios.

The complementarity formulations associated with the Moreau time-stepping algorithm for contact modeling in multibody systems have attracted the attention of many researchers [START_REF] Glocker | Dynamical systems with unilateral contacts[END_REF][START_REF] Glocker | Multiple impacts with friction in rigid multibody systems[END_REF][START_REF] Moreau | Numerical aspects of the sweeping process[END_REF][START_REF] Pfeiffer | Unilateral problems of dynamics[END_REF][START_REF] Bauchau | Modeling of unilateral contact conditions with application to aerospace systems involving backlash, freeplay and friction[END_REF][START_REF] Stewart | Finite-dimensional contact mechanics[END_REF][START_REF] Pfeiffer | The idea of complementarity in multibody dynamics[END_REF][START_REF] Förg | Simulation of unilateral constrained systems with many bodies[END_REF][START_REF] Glocker | Formulation and preparation for numerical evaluation of linear complementarity systems in dynamics[END_REF][START_REF] Slavic | Non-linearity and non-smoothness in multi-body dynamics: application to woodpecker toy[END_REF][START_REF] Leine | Stability properties of equilibrium sets of non-linear mechanical systems with dry friction and impact[END_REF][START_REF] Studer | Step size adjustment and extrapolation for time-stepping schemes in non-smooth dynamics[END_REF]. Assuming that the contacting bodies are truly rigid, as opposed to locally deformable or penetrable as in the penalty approaches, the complementarity formulations resolve the contact dynamics problem by using the unilateral constraints to compute contact forces or impulses to prevent penetration from occurring. Thus, at the core of the complementarity approach is an explicit formulation of the unilateral constraints between the contacting rigid bodies [START_REF] Brogliato | Numerical simulations of finite dimensional multibody nonsmooth mechanical systems[END_REF]. One of the main features of unilateral constraints is the impenetrability, which means that points candidates for contact must not cross the boundaries of antagonist bodies. This can be expressed by writing that the distance between contacting bodies or the gap is nonnegative. Also, it is assumed that bodies are not attracting each other, that is, the reaction force is nonnegative, and this reaction force vanishes when the contact is not active [START_REF] Leine | Dynamics and Bifurcations of Non-Smooth Mechanical Systems[END_REF].

When dealing with the formulation of frictional unilateral constraints, it is possible to distinguish active and passive set-valued force laws. An active setvalued force law is always associated with a closed unilateral contact or a frictional contact, while a passive set-valued force law is related to open unilateral contacts. An active set-valued force law can be described at the velocity level by an inclusion. It should be highlighted that a system with active unilateral constraints has variable degrees of freedom, being in general, not known which degree of freedom is removed. This problem is usually solved by looking at all possible solutions and finding the one that is physically consistent. It is obvious that the search for a physical consistent is time-consuming. In addition, from the numerical simulation point of view, it is quite unsuitable to change the number of the minimum generalized coordinates during each time step. Thus, the augmented Lagrangian approach is quite an elegant way to solve this type of problem, being constant the number of generalized coordinates at all instants of time. The number of generalized coordinates is always equal to the number of degrees of freedom of the system without unilateral constraints [START_REF] Anitescu | Formulating dynamic multirigid-body contact problems with friction as solvable linear complementarity problems[END_REF][START_REF] Brogliato | Some perspectives on the analysis and control of complementarity systems[END_REF][START_REF] Flores | Contact-impact analysis in multibody systems based on the nonsmooth dynamics approach[END_REF].

The dynamic modeling and analysis of planar multibody systems that experience contact-impact events is presented and discussed throughout this work. The methodology is based on the nonsmooth dynamics approach, in which the interaction of the colliding bodies is modeled with multiple frictional unilateral constraints. The model of a rigid multibody system is stated as an equality of measures, which are formulated at the velocity-impulse level. The equations of motion are complemented with constitutive laws for the forces and impulses the normal and tangential directions. The unilateral constraints are described by a set-valued force law of the type of Signorini's condition, while the frictional contacts are characterized by a set-valued force law of the type of Coulomb's law for dry friction. The formulation of the generalized contact-impact kinematics in the normal and tangential directions can be performed by obtaining a geometric relation for the gaps of the candidate contact points. The gaps are expressed as functions of the generalized coordinates. The candidate contact points are modeled as hard contacts, being the normal and tangential contact laws formulated as setvalued force laws for frictional unilateral constraints. Furthermore, when a system includes frictional unilateral constraints, the occurring contact forces should be taken into account in the equations of motion. In the present work, due to its simplicity and robustness, this problem is solved using the Moreau time-stepping method, combined with the contact-impact formulated as an augmented Lagrangian approach. Finally, results for a planar rigid multibody system are presented and used to discuss the main assumptions and procedures adopted in this work.

Set-valued force laws for frictional unilateral contacts

In the present work, the normal contact between rigid bodies is characterized by a set-valued force law called Signorini's condition [START_REF] Signorini | Atti della Societa Italian per il Progresso della Scienza[END_REF]. In this section, the Signorini's law, which is an elementary set-valued force law, is formulated in three different ways, namely, as complementarity conditions, as subdifferential equations, and as normal cone of the admissible contact forces [START_REF] Flores | Contact-impact analysis in multibody systems based on the nonsmooth dynamics approach[END_REF]. Figure 1 shows two convex rigid bodies apart from each other by a relative normal gap or distance denoted by g N . This relative normal gap is uniquely defined for convex surfaces, being parallel the tangent planes at the contact points 1 and 2. The relative normal gap is nonnegative due to bodies' impenetrability condition, being the two bodies in contact with each other when g N = 0. On the other hand, the normal contact force λ N is also nonnegative because the bodies cannot attract each other. The normal contact force vanishes when there is no contact, i.e., g N > 0, and can only be positive when contact happens, that is, g N = 0. Thus, under the assumption of impenetrability between the bodies, expressed by g N ≥ 0, only two situations can occur, namely

g N = 0 ∧ λ N ≥ 0 (closed contact) (1) g N > 0 ∧ λ N = 0 (open contact) (2)
Equations ( 1) and (2) represent a complementarity behavior, being always zero the product of the relative normal gap and normal contact force, that is, Thus, the relation between the normal gap and normal contact force can be described by

g N λ N = 0 ( 3 )
g N ≥ 0, λ N ≥ 0, g N λ N = 0 ( 4 )
which represents the inequality complementarity condition between g N and λ N , the so-called Signorini's condition. The inequality complementarity behavior of the normal contact law is depicted in Fig. 2a, which shows a set-valued graph or a corner of admissible combinations between g N and λ N [START_REF] Leine | Stability properties of equilibrium sets of non-linear mechanical systems with dry friction and impact[END_REF]. The magnitude of contact force is denoted by λ N and the direction of the contact force is normal to the bodies' surface, i.e., along the lines 1-2. When two rigid bodies are contacting, the Signorini's condition given by Eq. ( 4) needs to be complemented with an impact law, such as the well known Newton's kinematical law that relates the pre and postimpact velocities to the bodies' normal coefficient of restitution, ε N [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth Systems[END_REF].

The normal contact law can also be expressed as the subdifferential of the indicator function Ψ C N to convex set C N [START_REF] Leine | Stability properties of equilibrium sets of non-linear mechanical systems with dry friction and impact[END_REF]

g N ∈ ∂Ψ C N (-λ N ) (5)
where C N is the admissible set of negative contact forces -λ N as

C N = -λ N ∈ R λ N ≥ 0 = R - (6)
Since the subdifferential of the indicator function is the inverse of the subdifferential of the support function, then Eq. ( 5) can be rewritten using the conjugate as

-λ N ∈ ∂Ψ * C N (g N ) (7)
Alternatively, the normal contact law can be formulated in a compact form by means of the normal cone of C N [START_REF] Leine | Stability properties of equilibrium sets of non-linear mechanical systems with dry friction and impact[END_REF] g

N ∈ N C N (-λ N ) ( 8
)
where N C N (λ N ) denotes the normal cone of the convex set C N at λ N , and C N is the set of admissible contact forces [START_REF] Rulka | SIMPACK-a computer program for simulation of large motion multibody systems[END_REF]. It should be noted that this formulation is identical to the complementarity given by Eq. ( 4).

For closed contacts with g N = 0, the Signorini's law, which is a set-valued law for normal contact at the displacement level can be expressed at the velocity level as

γ N ∈ N C N (-λ N ), g N = 0 ( 9 )
where γ N is the relative normal velocity. It is also common to express the set-valued force law of Eq. ( 9) in terms of impulsive forces instead of contact forces [START_REF] Glocker | Formulation and preparation for numerical evaluation of linear complementarity systems in dynamics[END_REF].

The classical Coulomb's friction law is another typical example that can be considered as a set-valued force law [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth Systems[END_REF]. The friction law of Coulomb states that the sliding friction is proportional to the normal force of a contact. The amount of the static friction force is less than or equal to the maximum static friction force which is also proportional to the normal contact force. Furthermore, the sliding force has the opposite direction of the relative velocity of the frictional contact. Consider again the two contacting rigid bodies depicted in Fig. 1, in which Coulomb friction is present at the contact points 1 and 2. The relative velocity of point 1 with respect to point 2 along their tangent plane is denoted by γ T . If contact between the two bodies takes place, i.e., g N = 0, then the friction phenomenon imposes a tangential force λ T as it is illustrated in Fig. 1b. If the bodies are sliding over each other, then the friction force λ T has the magnitude μλ N and acts in the direction oppose to the relative tangential velocity, that is,

-λ T = μλ N Sgn(γ T ) γ T = 0 ( 10
)
where μ is the friction coefficient and λ N is the normal contact force. If the relative tangential velocity vanishes, i.e. γ T = 0, then the bodies purely roll over each other without slip. Pure rolling, or slip for locally flat objects, is denoted by stick. Thus, if the bodies stick, then the friction force must lie in the interval -μλ N ≤ λ T ≤ μλ N . For unidirectional friction, that is, for planar contact problems, three different scenarios can occur, namely

γ T = 0 ⇒ |λ T | ≤ μλ N (sticking) (11) γ T < 0 ⇒ λ T = +μλ N (negative sliding) (12) γ T > 0 ⇒ λ T = -μλ N (positive sliding) (13)
These three scenarios can be summarized by a setvalued force law as [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth Systems[END_REF]

λ T ∈ ⎧ ⎨ ⎩ -μλ N γ T > 0 [-1, 1]μλ N γ T = 0 μλ N γ T < 0 ( 14
)
Figure 2b shows the Coulomb's friction law as a set-valued force law. The admissible values of the negative tangential force λ T form a convex set C T that is bounded by the values of the normal force [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF]

C T (λ N ) = -λ T -μλ N ≤ λ T ≤ +μλ N (15)
Thus, the Coulomb's law can be expressed with the aid of the indicator function of C T as

γ T ∈ ∂Ψ C T (λ N ) (-λ T ) (16) or γ T ∈ N C T (λ N ) (-λ T ) (17)
Alternatively, Eqs. ( 16) and ( 17) can be written as [START_REF] Leine | Dynamics and Bifurcations of Non-Smooth Mechanical Systems[END_REF]

-λ T ∈ ∂Ψ * C T (λ N ) (γ T ) (18)
In short, the tangential forces are limited by a maximal friction force dependent on μλ N in any direction inside the tangential plane, and where μ represents the friction coefficient. Lower forces correlate with sticking contacts and, therefore, γ T = 0. Forces reaching the boundary of the friction cone may indicate sliding and, therefore, γ T = 0. The full description of the spatial Coulomb's friction law as a set-valued force law can be found in Leine and Glocker [START_REF] Leine | A set-valued force law for spatial Coulomb-Contensou friction[END_REF].

Multibody systems with frictional unilateral constraints

It is known that impacts and frictional phenomena are characterized by unilateral constraints, which usually lead to unsteady dynamical behaviors. Thus, appropriate methodologies and procedures to deal with this class of mechanical systems are required, being the main purpose of the present work. From classical mechanics, the Newton-Euler equations of motion of a multibody system with f degrees of freedom and with only frictionless bilateral constraints can be written as [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF] M uh = 0

q = u ∀t (19)
where M = M(q, t) ∈ R f ×f is the positive definite and symmetric mass matrix, h = h(q, u, t) ∈ R f denotes, in the present work, the vector that contains the differentiable forces (both conservative and non-conservative), such as spring forces, gravitational forces, q = q(t) ∈ R f is the f -dimensional vector of generalized coordinates, u = u(t) ∈ R f addresses the system generalized velocities and u = u(t) ∈ R f is the vector that contains the system accelerations. Joint reaction forces of the bilateral constraints do not appear in the equations of motion [START_REF] Machado | Development of a planar multi-body model of the human knee joint[END_REF] because the coordinates q are minimal Lagrangian coordinates with respect to bilateral constraints, i.e., the vector q represents a set of coordinates that defines uniquely the positions of all bodies in the system when all unilateral contacts are open. The dependence of the system matrices on q, u and t has been omitted in Eq. ( 19) for brevity. The terms M and h can be derived in a straightforward manner, by taking q as a set of classical generalized system coordinates and evaluating Lagrange's equations of second type or the associated virtual work expressions [START_REF] Greenwood | Principles of Dynamics[END_REF].

Adding the contact forces to Eq. ( 19), the dynamic equations of motion of a multibody system with normal and tangential contact forces during an impact can be written at the acceleration level as [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF]

M u -h -w N λ N -w T λ T = 0 (21) q = u ∀t (22)
where w N = w N (q, t) ∈ R f and w T = w T (q, t) ∈ R f represent the generalized normal and tangential force directions, respectively. The normal and tangential contact forces have magnitudes λ Ni and λ Ti for each contact point i. The dual variables to the normal contact forces λ N are the variations of normal gap distances g N , while the dual variables to the generalized friction or tangential forces λ T are the variations of the generalized sliding velocities γ T . The remaining terms of Eq. ( 21) have the same meaning as described above.

It is important to note that Eq. ( 21) requires the existence of the velocities u and accelerations u both being meaningless for the event of an impact. Therefore, it is more adequate to talk about the left and right limit of the velocity at the impact, that is, the pre and postimpact velocity, but never about the velocity at the impact itself, a meaningless term already from the physical point of view. For the case of impacts in multibody systems, Eq. ( 21) have to be substituted by a more suitable formulation, which consists of use the equalities of measures firstly introduced by Moreau [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] and that constitute the general framework for nonsmooth rigid multibody dynamics.

Moreover, motion without impulses implies that λ N (t) is (locally) bounded and time-continuous. The friction force λ T (t) is discontinuous when a slip-stick transition takes place or when the relative sliding velocity of a frictional contact reverses its sign. The acceleration u is not defined when λ T (t) is discontinuous. The set of time instances for which λ T (t) is discontinuous is of measure zero and Eq. ( 21), therefore, holds for almost all t , that is, q = u does not hold at single time instants at which impacts take place. Thus, due to the presence of impulsive forces, a nonsmooth system cannot be described solely by the equations of motion [START_REF] Kim | A two-dimensional tire model on uneven roads for vehicle dynamic simulation[END_REF]. Equalities of measures provide an elegant way to obtain a valid comprehensive description of a nonsmooth system including the impact case. When the equations of motion for the impact case are integrated over a singleton in time yields

M u + -u --w N N -w T T = 0 (23) q = u a.e. ( 24
)
in which u -and u + represent the pre and post impact velocities, N and T denote the normal and tangential impulsive forces, being the remaining terms defined as previously. Note that contact forces are replaced by the impulsive forces, which are well defined in the case of an impact. Furthermore, finite forces, such as gravity or reaction forces from springs and dampers, do not contribute to the integral and, therefore, they are not considered in Eq. ( 23).

The equations of motion without impacts given by Eq. ( 21) and the equations of motion for the impact case ( 23) cannot be used together in the present form, because the former is written at the acceleration level and the second one is developed at the velocity level. Therefore, Eqs. ( 21) and ( 23) should be considered simultaneously. Multiplying Eqs. ( 21) and ( 23) dt and dη, respectively, yields

M u dt -h dt -w N λ N dt -w T λ T dt = 0 (25)
M u + -u -dη -w N N dη -w T T dη = 0 (26)
Adding now Eqs. ( 25) and ( 26) results in

M u dt + u + -u -dη -h dt -w N (λ N dt + N dη) -w T (λ T dt + T dη) = 0 (27)
or more briefly

M du -h dt -w N dP N -w T dP T = 0 ( 28
)
where the Lebesgue measure is represented by dt and dη represents the sum of the Dirac impulse measures at the impact times [START_REF] Leine | Dynamics and Bifurcations of Non-Smooth Mechanical Systems[END_REF].

The measure for the velocities du = u dt + (u +u -) dη is split in Lebesgue measurable part u dt , which is continuous, and the atomic parts which occur at the discontinuity points with the left and right limits u -and u + and the Dirac point measure dη. For impact free motion it holds that du = u dt. Similarly, the measure for the so-called percussions corresponds to a Lagrangian multiplier which gathers both finite contact forces λ and impulsive contact forces , that is, dP = λ dt + dη [START_REF] Förg | Simulation of unilateral constrained systems with many bodies[END_REF][START_REF] Glocker | Formulation and preparation for numerical evaluation of linear complementarity systems in dynamics[END_REF]. In the case of nonimpulsive motion, all measures dη vanish and a formal division by dt yields the classical Newton-Euler equations of motion given by [START_REF] Kim | A two-dimensional tire model on uneven roads for vehicle dynamic simulation[END_REF].

Since the impenetrability condition between colliding bodies is required, let us consider that a multibody system has a total n of frictional unilateral constraints, which can be represented by n inequalities as

g Ni (q, t) ≥ 0, i = 1, . . . , n (29)
where the quantities g Ni are the normal gap functions of the frictional contacts. They are formulated such that, g Ni > 0 indicates an open or positive contact with an Euclidian distance of the contact points given by the value of g Ni , g Ni = 0 corresponds to a closed or active contact, and g Ni < 0 indicates the forbidden overlapping or interpenetration between rigid bodies.

A good treatment of the definition of these inequalities, under the framework of multibody systems formulation, is discussed by Pfeiffer and Glocker [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF] and Glocker [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth Systems[END_REF].

The set of active contacts in the present work is stated as

H (t) = i g Ni (q, t) = 0 ( 30
)
which singles out the contact at which contact-impact forces may occur. In fact, ideally when g Ni (q, t) = 0 the contact is active. However, due to the computation round-off errors accumulation, a tolerance must be introduced in order to accommodate for inaccuracies in the numerical results. Therefore, when the first penetration is within the penetration tolerance it is assumed that such is the moment of the impact and the position and relative velocity of the contact points and the direction of the plane of collision are recorded. This approach has implications with the integration scheme and time steps, being the reader interested in the details on this particular topic referred to the work by Studer et al. [START_REF] Studer | Step size adjustment and extrapolation for time-stepping schemes in non-smooth dynamics[END_REF].

In theory, the constraints are active when the gap vanishes and the constraint forces are such that the gap does not become negative (no penetration). In the numerical scheme the constraint is set active if the gap function becomes zero or negative (i.e., is no longer strictly positive). We therefore allow in the numerical scheme for small interpenetrations, being a numerical approximation of g N = 0. However, the contact law is evaluated on velocity level and it might therefore happen that the constraint drifts resulting in an unwanted penetration of noticeable size. Such a constraint drift can be remedied with a projection to a negative threshold (such that the constraint remains active) if the penetration becomes too large. Furthermore, the time step needs to be small enough to describe hi-frequency bounces. However, if the exact solution has multiple bounces within a time step of the numerical solution, then the numerical solution is still a "good" approximation of the exact solution in the sense that the numerical solution converges to the exact solution for decreasing step size.

In order to define the constitutive force laws which relate the contact-impact impulse measures to the system's kinematics q and u, let first introduce the normal and tangential relative velocities at the contacts as [START_REF] Glocker | On frictionless impact models in rigid-body systems[END_REF]

γ Ni = w T Ni u + wNi (31)
γ Ti = w T Ti u + wTi (32)
where w Ni and w Ti represent the generalized normal and tangential force directions, respectively, and wNi and wTi are the partial derivatives of the of the gap functions with respect to time. The equations of motion [START_REF] Mukras | Analysis of planar multibody systems with revolute joint wear[END_REF] can now be complemented with constitutive laws for normal and tangential contact-impact forces. In the present study, a unilateral version of the Newton's impact law is considered for the normal direction with local coefficient of restitution ε Ni ∈ [0, 1]. The Coulomb's friction law is used for the tangential direction with coefficient of friction μ i , which is complemented by a tangential coefficient of restitution ε Ti ∈ [0, 1].

It is important to note that for the Newton's impact law, the impact, which causes the sudden change in the relative velocity, is accompanied by a normal contact impulse dP N > 0. Suppose that, for any reason, the contact does not participate in the impact, that is, that value of the normal contact impulse is zero, although the contact is closed. This situation happens normally for multiple contact scenarios. Therefore, for this case, it is allowed that the post impact relative velocity to be higher than the value prescribed by Newton's impact law, with the intent to express that the contact is superfluous and could be removed without changing the contact-impact process. Thus, in order to account for these possibilities, two parameters are defined as [START_REF] Glocker | On frictionless impact models in rigid-body systems[END_REF] ξ

Ni := γ + Ni + ε Ni γ - Ni ( 33
)
ξ Ti := γ + Ti + ε Ti γ - Ti (34)
where (γ ± Ni , γ ± Ti ) := (γ Ni , γ Ti )(u ± ). Thus, normal and tangential impact laws can be stated as

-dP Ni ∈ ∂Ψ * C Ni (ξ Ni ) ( 35
)
-dP Ti ∈ ∂Ψ * C Ti (P Ni ) (ξ Ti ) (36)
or, in terms of the normal cones

ξ Ni ∈ N C Ni (-dP Ni ) ( 37
)
ξ Ti ∈ N C Ti (P Ni ) (-dP Ti ) (38)
Finally, the complete description of the dynamics of nonsmooth system, which accounts for both impact and impact-free phases, is given by Eqs. ( 28)- [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF].

Augmented Lagrangian approach to solve contact-impact events

In this section, the augmented Lagrangian approach to solve the contact-impact problem of multibody systems with frictional unilateral constraints is presented. In a simple way, this approach consists of transforming the inclusions in normal and tangential directions into equivalent equations using proximal point of convex analysis, then the problem is solved iteratively as a proximal point formulation [START_REF] Förg | Simulation of unilateral constrained systems with many bodies[END_REF][START_REF] Leine | A set-valued force law for spatial Coulomb-Contensou friction[END_REF][START_REF] Pfeiffer | Numerical aspects of non-smooth multibody dynamics[END_REF]. After discretization, the velocity-impulse equations of motion for nonsmooth mechanical systems can be written as [START_REF] Flores | Modeling and analysis of rigid multibody systems with translational clearance joints based on the nonsmooth dynamics approach[END_REF] M M uh M t -W NM P N -W TM P T = 0 [START_REF] Pfeiffer | Unilateral problems of dynamics[END_REF] in which the subscript M denotes the mid point of the integration time step. Equation ( 39) is used together with the set-valued impulsive force laws,

-P N ∈ ∂Ψ * C N (ξ N ) ( 40
)
-P T ∈ ∂Ψ * C T (P N ) (ξ T ) ( 41
)
where Ψ * C N and Ψ * C T (P N ) denote the support functions of the indicator functions Ψ C N and Ψ C T (P N ) , respectively, being N C N and N C T (P N ) the normal cones associated with the convex sets C N and C T . These two convex sets represent the admissible normal and tangential impulsive forces given by [START_REF] Förg | Simulation of unilateral constrained systems with many bodies[END_REF][START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF]]

C N = -P N ∈ R n P N ≥ 0 (42) C Ti (P Ni ) = -P Ti P Ti ≤ μ i P Ni ; i = 1, . . . , n (43)
The ξ N and ξ T parameters are expressed as [START_REF] Glocker | Formulation and preparation for numerical evaluation of linear complementarity systems in dynamics[END_REF] ξ N = γ NE + ε N γ NA [START_REF] Glocker | Formulation and preparation for numerical evaluation of linear complementarity systems in dynamics[END_REF] ξ T = γ TE + ε T γ TA [START_REF] Slavic | Non-linearity and non-smoothness in multi-body dynamics: application to woodpecker toy[END_REF] in which the actual normal and tangential velocities are given by [START_REF] Flores | Contact-impact analysis in multibody systems based on the nonsmooth dynamics approach[END_REF]

γ NE = W T NM u E + wNM (46)
γ TE = W T TM u E + wTM ( 47
)
γ NA = W T NM u A + wNM ( 48
)
γ TA = W T TM u A + wTM (49)
In short, in each time step the equations of motion ( 39) and the set-valued force laws [START_REF] Bauchau | Modeling of unilateral contact conditions with application to aerospace systems involving backlash, freeplay and friction[END_REF] and ( 41) have to be solved for u, P N and P T . The augmented Lagrangian approach is an elegant way to solve this problem, which transforms the set of algebraic inclusions to a constrained optimization problem [START_REF] Rockafellar | Augmented lagrangians and applications of the proximal point algorithm in convex programming[END_REF][START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF][START_REF] Studer | Solving normal cone inclusion problems in contact mechanics by iterative methods[END_REF]. By definition of proximal point to a convex set C yields that

prox C (x) = argmin ∀x * ∈C x -x * (50)
which represents the closet point in C to its argument. Thus, based on this concept, the normal and tangential impulses are stated as two equalities:

P N = prox C N (P N -rξ N ) (51)
P T = prox C T (P N ) (P T -rξ T ) (52)
with the convex sets C N and C T given by Eqs. ( 42) and [START_REF] Förg | Simulation of unilateral constrained systems with many bodies[END_REF]. It should be mentioned that this approach is dependent on a non-negative arbitrary parameter r, which represents the slope of the regularization function. The value of r should be taken large enough to make the problem well conditioned in the constrained region, but not too high in order to prevent ill-conditioning. The interested reader on the issue of the r-factor strategies for the augmented Lagrangian approach is referred to the work by Foerg et al. [START_REF] Foerg | r-Factor strategies for the augmented Lagrangian approach in multibody contact mechanics[END_REF]. Finally, the contact-impact problem of nonsmooth systems based on the augmented Lagrangian approach can be summarized by the following mathematical relations:

M M u -h M t -W NM Π N (P N , ξ N ) -W TM Π T (P N , P T , ξ T ) = 0 (53) Π N (P N , ξ N ) = prox C N (P N -rξ N ) (54) Π T (P N , P T , ξ T ) = prox C T (P N ) (P T -rξ T ) (55)
Thus, this set of algebraic equations can be easily included in the Moreau time-stepping method, being the saddle point of the augmented Lagrangian found by employing, for instance, the modified Newton method, also called Newton-Raphson method [START_REF] Förg | Simulation of unilateral constrained systems with many bodies[END_REF][START_REF] Leine | A set-valued force law for spatial Coulomb-Contensou friction[END_REF][START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF]. Besides some possible divergences, the Newton root-finding algorithm is considered in the present work due to its simplicity and computational efficiency. In this process, the computed impulsive forces are used as initial guesses in the next step. The Newton approach can fail if, for instance, the contact bodies stick together during the dynamic analysis, because the Newton iteration becomes singular. In order to overcome this numerical difficulty, the LUdecomposition and the forward/backward-substitution must be performed. The reader interested in the details on this particular issue is referred to the work by Förg et al. [START_REF] Förg | Simulation of unilateral constrained systems with many bodies[END_REF].

The computational strategy of the modified Newton method is illustrated in the flowchart of Fig. 3 and is summarized in the basic following steps:

(i) Initialize the counter variable k = 1 and specify initial guesses for the impulses P k N and P k T ; (ii) Solve u k+1 from the following equation:

M M u k+1 = h M t + W NM P k N + W TM P k T (iii) Evaluate ξ k+1 NE ξ k+1 TE from equations ξ k+1 N = W T NM u k A + u k+1 + wNM ξ k+1 T = W T TM u k A + u k+1 + wTM (iv) Project P k+1
N and P k+1 T as

P k+1 N = Π N P k N , ξ k+1 N = prox C N P k N -rξ k+1 N P k+1 T = Π T P k+1 N , P k T , ξ k+1 T = prox C T (P N ) P k T -rξ k+1 T (v)
Compute the tolerance error as

error = P k+1 N -P k N + P k+1 T -P k T (vi)
If the error is lesser than a specified tolerance (tol), then the process has converged and should, therefore, end; otherwise, increment counter k = k + 1 and go to step (ii) to proceed with the process of a new iteration step.

Since the Moreau time-stepping method with an augmented Lagrangian approach involves a good deal of mathematical manipulation, it is convenient to summarize the main steps in an appropriate algorithm. This algorithm, presented in the flowchart of Fig. 4, is developed under the framework of multibody systems formulation and can be condensed in the following steps:

(i) Specify the initial conditions of the problem at hand, t A , t F , t , q A , and u A ; (ii) Define the geometrical, inertial, and material functions, g Ni , M, h, ε Ni , ε Ti , μ i , w Ni , w Ti , wNi , and wTi ; (iii) Compute the mid point state variables:

t M = t A + 1 2 t q M = q A + 1 2 tu A M M = M(q M , t M ) h M = h(q M , u A , t M ) g Ni = g Ni (q M , t M ) H M = i g Ni (q M , t M ) ≤ 0 n i = length(H M )
(iv) For every i ∈ H M evaluate:

W NM = mat w Ni (q M , t M ) ,
W TM = mat w Ti (q M , t M ) , wNM = col wNi (q M , t M ) , wTM = col wTi (q M , t M ) ,

γ NA = col(γ NAi ), γ TA = col(γ T Ai ), ε N = diag(ε Ni ), ε T = diag(ε Ti ), μ = diag(μ i ), k = 1, P k N = P N (H M ), P k T = P T (H M ) (v)
Compute while the Newton method does not converge, i.e., while error > tol: (vi) Compute the positions at the end of the integration time step:

u E = u A + M -1 M h M t + M -1 M W NM P k N + M -1 M W TM P k T γ NE = W T NM u E + wNM γ TE = W T TM u E + wTM ξ k+1 N = γ NE + ε N γ NA ξ k+1 T = γ TE + ε T γ TA P k+1 N = Π N P k N , ξ k+1 N = prox C N P k N -rξ k+1 N P k+1 T = Π T P k+1 N , P k T , ξ k+1 T = prox C T (P N ) P k T -rξ k+1 T error = P k+1 N -P k N + P k+1 T -P k T k = k + 1
q E = q M + 1 2 tu E
(vii) Increment time step:

t A = t A + t
(viii) Update system states' variables q A = q E and u A = u E . Go to step (iii) and proceed with the process for the new time step. These steps must be performed until the final time of analysis is reached.

Analysis of the cam-follower mechanism

This section contains a demonstrative example of application in which the contact-impact events are mod-eled and analyzed under the nonsmooth dynamics approach. This example deals with a cam-follower mechanism with of an industrial application of a cutting file machine [START_REF] Seabra | Theoretical and experimental analysis of an industrial cutting file machine using multibody systems methodology[END_REF][START_REF] Seabra | Kinematics and dynamics study of the cam-follower mechanism of the cutting file machine[END_REF][START_REF] Seabra | Re-design of a cam-follower mechanism of an industrial cutting file chine[END_REF]. Figure 5 shows the overall view and the schematic representation of this machine-tool. The file teeth are produced by impact of the cutting beater (system composed by follower, cylinder, and chisel), with a reciprocate movement. To generate this movement, the cutting bench has a wheel with six rebounds (cam) whose rotation forces the pin to move up. This will lift up the cylinder, to which the chisel is attached, which immediately falls down, when reaching the updead-point, impelled by the spring and its own weight. The impact energy of the chisel depends on the relationship between the spring force and the maximum distance between the chisel and the file (adjusted by a presser foot). The chisel describes a reciprocating motion that always reaches the same up dead point (maximum distance between the chisel and the file), while the pin, rigidly attached with the cylinder which moves the chisel, which always passes by the tops of cam. On the other hand, the down dead point of the chisel is variable, and depends on the impact energy absorbed by the file body. The chisel impact energy depends on the relationship of the regulations of the spring preload and on the maximum distance between the chisel and file. During the file manufacture, in order to obtain a tooth with the appropriate geometry (depth of the penetration), it is required, not only the impact energy should be adequately adjusted, but also the maximum distance between the chisel and the file resulting from the regulation of the presser foot should have a value that allows the chisel to pass above the last produced tooth.

When the machine operates correctly, the kinetic energy produced during the descending chisel movement is totally absorbed by the base body of the file. For that purpose, the presser foot must be adjusted in order to prevent impacts between cam and pin during the descent of the chisel movement. It means that the pin should never collide with cam. When this situation does not happen, as consequence of incorrect positioning of the presser foot, it can be observed that the cutting operation produces a hard and increasing noise. The noise is due to the impact of the pin on the cam, and strongly depends on the spring force. This clash is undesirable for two main reasons; first, because it accelerates the cam and follower wear, and secondly, because it decreases the kinetic energy available for the cutting operation, since part of the energy is be absorbed by that impact. Hence, the file quality is significantly penalized [START_REF] Seabra | Kinematics and dynamics study of the cam-follower mechanism of the cutting file machine[END_REF].

The multibody system of the cutting file machine is made of three rigid bodies (cam-the driver, follower-the driven element, and the ground or frame), one revolute joint, and one translational joint. Figure 6a depicts the kinematic configuration of the cam follower mechanism. It is known that for n b rigid body system with n c independent constraint equations, the mobility or degrees of freedom (DOF) is given by [

2] DOF = 6 × n b -n c ( 56
)
This mathematical expression, usually called as Grüebler equation, can be used to determine the mobility of multibody system. Thus, from Eq. ( 56), the DOF of the cam-follower mechanism is equal to 1, implying one, and only one, motion generator. Since the fol-lower can not rotate about its own axis, and the follower curvature radius is very large when compared to its own dimensions, the follower can be considered flat faced. The flat faced follower has the advantage of a zero degree pressure angle throughout its motion, which is an important feature, since most of camfollower mechanisms are designed with pressure angles as small as possible [START_REF] Chen | Mechanics and Design of Cam Mechanisms[END_REF][START_REF] Norton | Cam Design and Manufacturing Handbook[END_REF]. Figure 6b schematically illustrates the experimental data relative to the follower displacement diagram corresponding to a sixth part of the cam angle rotation, since the cam has six rebounds and the camfollower motion repeats itself six times in each complete cam rotation [START_REF] Seabra | Re-design of a cam-follower mechanism of an industrial cutting file chine[END_REF]. In Fig. 6b, point A represents the maximum follower displacement, point B defines the instant of impact between the follower and file body, point C corresponds to the minimum follower displacement, that is, the maximum penetration/deformation of the body file, and finally, and point D represents the recontact between the cam and follower after the rebound effect. Observing Fig. 6b, it is evident that the follower motion can be divided into two main phases, namely, the fall and the rise movements. In turn, these two phases can be analyzed into two different parts. Starting from maximum follower elevation, point A, the follower motion can be described and summarized by the following steps:

(i) Fall #1-from point A to point B: the follower motion is influenced by three main factors, the gravity effect, the spring action and friction phenomenon that exists between the follower and guide. At point A, the follower is pushed down by preloaded spring and gravity action; (ii) Fall #2-from point B to point C: point B represents the initial instant of impact between the follower (chisel) and the file body. The maximum penetration depth, which corresponds to the edge height, is represented by the distance between points B and C. Point C corresponds to the end of follower fall motion; (iii) Rise #1-from point C to point D: this phase represents the rebound effect caused by the accumulated energy during the contact-impact process between the follower and file body. In this process, there is no contact between the follower and cam due to rebound effect and cam speed; (iv) Rise #2-from point D to point A: the follower is in permanent contact with the cam surface, hence, the follower is rising and the spring is preloading in this process.

The simulation parameters of the cam follower system are listed in Table 1. The system is considered to be frictionless.

In order to keep the analysis simple, the follower motion is considered to be of sinusoidal type, being the displacement expression given by

y = s θ β - 1 2π sin 2πθ β ( 57
)
where s is maximum stroke of the follower, θ represents the angle of the cam rotation corresponding to displacement of the follower y and β is the angle of cam rotation to reach the stroke s. Since the cam follower has one degree of freedom, the variables necessary to define the problem are as follows:

q = (y) (58) u = ( ẏ) (59
Figure 7 shows the behavior of the follower for the data presented above and for a full cam rotation, that is, the follower displacement and the follower velocity. From these two plots, the different contact scenarios between the cam and follower are well visible, namely the continuous or permanent contact and the impact followed be rebounds due to the impacts that take place. The global results are in line with those offered in the literature [START_REF] Seabra | Theoretical and experimental analysis of an industrial cutting file machine using multibody systems methodology[END_REF]. Figure 8 shows an animation sequence of the simulation of the cam follower movement during the first instants after the follower reaches the up dead point. A general and comprehensive methodology for modeling and analyzing contact-impact events in multibody systems has been presented in this work. In the process, the main issues of the nonsmooth dynamics approach have been revised in face of their suitability to represent the interactions developed between the colliding bodies. The unilateral constraints were described by a set-valued force law of the type of Signorini's condition, while the frictional contacts were characterized by a set-valued force law of the type of Coulomb's law for dry friction. The resulting contactimpact problem was formulated and solved as an augmented Lagrangian approach, which was incorporated in the Moreau time-stepping method.

The proposed methodology has been exemplified through the application to a cam-follower mechanism of an industrial cutting-file machine. The proposed methodology is able to capture different phenomena involved in the dynamics of multibody system with multiple contacts, such as impacts followed by rebounds, sliding friction, and stick phenomenon. This type of response strongly depends on the system's dynamic behavior. An important result from this research work is that multibody systems with multiple impacts can have a predictable nonlinear performance.

Fig. 1 Fig. 2 (

 12 Fig. 1 (a) Relative normal gap; (b) Normal and tangential contact forces

Fig. 3

 3 Fig. 3 Computational strategy for the iterative Newton method

Fig. 4

 4 Fig. 4 Flowchart of the Moreau time-stepping algorithm with augmented Lagrangian approach

Fig. 6 1 Follower

 61 Fig. 6 (a) Kinematic configuration of the cam follower mechanism; (b) Follower displacement

h

  = (-mg -K s y) (61) g N = y -R bs (62) W N = (1)(63)wN = (-ẏ)

Fig. 7 Fig. 8

 78 Fig. 7 (a) Follower displacement; (b) follower velocity

Acknowledgements

This work is supported by the Portuguese Foundation for the Science and Technology under the research project BIOJOINTS (PTDC/EME-PME/099764/ 2008). The first author expresses his gratitude to the Portuguese Foundation for the Science and Technology for the postdoctoral scholarship (SFRH/BPD/40067/2007). This research was conducted during a post-doctoral stay of the first author with Professor Christoph Glocker at the Center of Mechanics, ETH Zurich.