Paulo Flores

Remco Leine

Christoph Glocker

Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach

Keywords: Non-smooth dynamics, Clearance joints, Contact-impact mechanics, Linear complementarity problem, Moreau time-stepping method, Rigid multibody systems 1

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

Manufacturing tolerances, wear and material deformation lead to imperfect joints and, therefore, clearances. These clearances modify the dynamic response of the system, justify the deviations between the numerical predictions and the experimental measurements and eventually lead to important deviations between the projected behavior of the mechanisms and their real outcome. The presence of clearance in joints is a complex and important issue in the realistic modeling of multibody systems. This aspect gains paramount importance due to the demand for the proper design of the real joints in many industrial applications. Over the last few years, extensive work has been done to study the dynamic effect of the revolute joints with clearance in multibody systems. However, translational joints with clearance have received less attention [START_REF] Flores | Revolute joints with clearance in multibody systems[END_REF][START_REF] Flores | Kinematics and Dynamics of Multibody Systems with Imperfect Joints: Models and Case Studies[END_REF][START_REF] Qiang | Dynamics of spatial flexible multibody systems with clearance and lubricated spherical joints[END_REF][START_REF] Liu | The FEM analysis and approximate model for cylindrical joints with clearances[END_REF][START_REF] Crowthera | Impulsive response of an automatic transmission system with multiple clearances: Formulation, simulation and experiment[END_REF][START_REF] Flores | Dynamic analysis for planar multibody mechanical systems with lubricated joints[END_REF][START_REF] Erkaya | A neural-genetic (NN-GA) approach for optimising mechanisms having joints with clearance[END_REF][START_REF] Flores | Modeling and simulation of wear in revolute clearance joints in multibody systems[END_REF].

Indeed, a number of theoretical and experimental works devoted to the research on multibody mechanical systems with realistic joints has been published recently. However, most of these works focus on revolute joints with and without lubrication effects. An extensive literature review on the issue of modeling and simulation of multibody systems with revolute and spherical clearance joints can be found in the work by Flores et al. [START_REF] Flores | Kinematics and Dynamics of Multibody Systems with Imperfect Joints: Models and Case Studies[END_REF]. In contrast to the revolute and spherical clearance joints, not much work has been done to model translational joints with clearance because in this case several different configurations between the joints elements can take place. In fact, the contact configurations of slider and guide include: (i) no contact between the two elements; (ii) one corner of the slider is in contact with the guide surface; (iii) two adjacent slider corners are in contact with the guide surface, which corresponds to have a face of the slider in contact with the guide surface; (iv) two opposite slider corners are in contact with the guide surface [START_REF] Flores | Translational joints with clearance in rigid multibody systems[END_REF][START_REF] Awrejcewicz | Dynamics investigation of three coupled rods with horizontal barrier[END_REF][START_REF] Stoenescu | Dynamic analysis of a planar rigid-link mechanism with rotating slider joint and clearance[END_REF]. Moreover, each contact point may be in stick or in slip phase, which greatly enlarges the number of contact configuration. The conditions for switching from one case to another depend on the system's dynamic response.

Farahanchi and Shaw [START_REF] Farahanchi | Chaotic and periodic dynamics of a slider crank mechanism with slider clearance[END_REF] studied the dynamic response of a planar slider-crank mechanism with slider clearance. They demonstrated how complex the system's response is, which can be chaotic or periodic. More recently, Thümmel and Funk [START_REF] Thümmel | Multibody modelling of linkage mechanisms including friction, clearance and impact[END_REF] used the complementarity approach to model impact and friction in a slider-crank mechanism with both revolute and translational clearance joints. With the purpose to analyze the slider-crank mechanism, Wilson and Fawcett [START_REF] Wilson | Dynamics of slider-crank mechanism with clearance in the sliding bearing[END_REF] derived the equations of motion for all different possible configurations of the slider motion inside the guide, resulting in a total of 40 equations. They also showed how the slider motion in a translational clearance joint depends on the geometry, speed and mass distribution.

In the present work, the non-smooth dynamics approach is used to model the type of multibody systems, due to its simplicity and ability to deal with all possible different configurations in a unified manner. The methodology is based on the non-smooth dynamics approach, in which the interaction of the colliding bodies is modeled with multiple frictional unilateral constraints. The dynamics of rigid multibody systems are 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. In this work, 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 setvalued force law of the type of Coulomb's law for dry friction. The resulting contact-impact problem is formulated and solved as a linear complementarity problem, which is embedded in the Moreau time-stepping method.

Basic set-valued elements

The linear complementarity problem

A linear complementarity problem (LCP) is a set of linear equations that can be written as [START_REF] Cottle | Complementary pivot theory of mathematical programming[END_REF][START_REF] Cottle | The Linear Complementarity Problem[END_REF]

y = Ax + b (1)
subjected to the inequality complementarity conditions

y ≥ 0, x ≥ 0, y T x = 0 ( 2 )
for which the vectors x and y have to be evaluated for given A and b. In other words, the LCP is the problem of finding solutions x ∈ R n and y ∈ R n of ( 1) and ( 2), where b is an n-dimensional constant column and A is a given square matrix of dimension n. The inequality complementarity conditions expressed by ( 2) are often written in the form

0 ≤ y⊥ x ≥ 0 (3)
where y⊥x denotes y T x = 0. An LCP can have a unique solution, multiple solutions or no solution at all [START_REF] Rohn | A note on solvability of a class of linear complementarity problems[END_REF][START_REF] Leine | Periodic motion and bifurcations induced by the Painlevé paradox[END_REF]. All existing solutions can be found using enumerative methods, which treat the problem by a combinatorial evolution of the complementarity condition x i y i = 0. From the complementarity condition it follows that when x i > 0, then y i = 0, and vice versa. An LCP of dimension n provides 2 n different combinations of n variables, which are allowed to be greater than zero at the same time. For large dimensions, enumerative methods become numerically expensive since 2 n grows rapidly. A more efficient algorithm is the complementarity pivot algorithm, usually referred to as Lemke's algorithm [START_REF] Pang | Complementarity formulations and existence of solutions of dynamic multi-rigidbody contact problems with Coulomb friction[END_REF][START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and Coulomb friction[END_REF][START_REF] Leine | Dynamics and Bifurcations of Non-Smooth Mechanical Systems[END_REF].

A drawback of Lemke's algorithm is that it is not guaranteed to find a solution for arbitrary A (convergence is guaranteed when A is a P -matrix). Other efficient algorithms to solve LCP can be found in the work by Cottle et al. [START_REF] Cottle | The Linear Complementarity Problem[END_REF].

The unilateral primitive

One of the most important multifunctions (or set-valued maps) related to complementarity is the unilateral primitive, denoted by Upr. The unilateral primitive is a maximal monotone set-valued map on R + 0 defined as [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth System[END_REF][START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF] Upr(x) :=

⎧ ⎪ ⎨ ⎪ ⎩ {0} x > 0 (-∞, 0] x = 0 ∅ x < 0 (4)
The graph of the unilateral primitive map is depicted in Fig. 1(a). Thus, each complementarity condition of an LCP can be expressed as one Upr inclusion,

-y ∈ Upr(x) ⇔ y ≥ 0, x ≥ 0, xy = 0 ( 5 )
Unilateral primitives are used in mechanics at the displacement level and at the velocity level to model unilateral geometric and kinematic constraints, such as free plays with stops, sprag clutches among others. The associated set-valued force laws are conveniently stated as inclusions of [START_REF] Crowthera | Impulsive response of an automatic transmission system with multiple clearances: Formulation, simulation and experiment[END_REF].

The Sgn-multifunction

A second maximal monotone set-valued map, frequently used in complementarity problems, is the filled-in relay function Sgn-multifunction, which is defined by [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth System[END_REF][START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF],

Sgn(x) := ⎧ ⎪ ⎨ ⎪ ⎩ {+1} x > 0 [-1, +1] x = 0 {-1} x < 0 (6)
It is important to highlight that, while the classical sgn-function is defined with sgn(0) = 0, the Sgn-multifunction is set-valued at x = 0. The graph of the Sgn-multifunction is shown in Fig. 1(b). An inclusion in the Sgn-multifunction can always be represented by two inclusions involving the unilateral primitive. The decomposition can be written as

-y ∈ Sgn(x) ⇔ ∃x R , x L s.t. ⎧ ⎪ ⎨ ⎪ ⎩ -y ∈ +Upr(x R ) + 1 -y ∈ -Upr(x L ) -1 x = x R -x L (7)
Using [START_REF] Crowthera | Impulsive response of an automatic transmission system with multiple clearances: Formulation, simulation and experiment[END_REF], the (7) can be rewritten in terms of complementarities,

-y ∈ Sgn(x) ⇔ ∃x R , x L s.t. ⎧ ⎪ ⎨ ⎪ ⎩ 1 + y ≥ 0, x R ≥ 0, (1 + y)x R = 0 1 -y ≥ 0, x L ≥ 0, (1 -y)x L = 0 x = x R -x L (8)
This representation has to be used when a problem involving Sgn-multifunctions is formulated as an LCP in its standard form [START_REF] Glocker | Formulation and preparation for numerical evaluation of linear complementarity systems in dynamics[END_REF]. In mechanics, relay functions at the velocity level are used to represent any kind of dry friction. In turn, when expressed at the displacement level, they describe the behavior of pre-stressed springs. More details on this decomposition can be found in the work by Glocker [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth System[END_REF].

Set-valued force laws for frictional unilateral contacts

Set-valued normal contact law

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]. Figure 2 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 perpendicular to the tangent planes at the contact points 1 and 2. The relative normal gap is non-negative due to the impenetrability condition of the bodies. The two bodies in contact with each other when g N = 0. In fact, one of the main features of unilateral contact is the impenetrability condition, which means that the candidate bodies for contact must not cross the boundaries of antagonist bodies. On the other hand, the normal contact force λ N is also non-negative, because the bodies cannot attract each other, that is, the constraint is unilateral. 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) (9) g N > 0 ∧ λ N = 0 (open contact) (10)
Equations ( 9) and ( 10) represent an inequality complementarity behavior, for which the product of the relative normal gap and normal contact force is always zero, that is,

g N λ N = 0 (11)
Thus, the relation between the normal gap and normal contact force can be described by

g N ≥ 0, λ N ≥ 0, g N λ N = 0 ( 12
)
which represents the inequality complementarity condition between g N and λ N , the socalled Signorini's condition.

The inequality complementarity behavior of the normal contact law is depicted in Fig. 3(a) and 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]. When two rigid bodies are contacting, the Signorini's condition given by [START_REF] Farahanchi | Chaotic and periodic dynamics of a slider crank mechanism with slider clearance[END_REF] needs to be complemented with an impact law, such as the well-known Newton's kinematical law that relates the pre-and post-impact velocities to the bodies' normal coefficient of restitution, ε N .

It should be highlighted that use of the Newton's impact law in combination with Coulomb friction can, under circumstances, lead to an (unphysical) energy increase. This typically occurs when there is a wide spread in normal and tangential restitution coefficients. Therefore, alternative methods for the definition of the coefficient of restitution, such as the Poisson's or Stronge's definition, can be considered. Sufficient conditions for energy decrease with Newton's impact law can be found in Leine and van de Wouw [START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF]. These issues have recently been discussed by other authors [START_REF] Djerassi | Collision with friction; Part A: Newton's hypothesis[END_REF][START_REF] Djerassi | Collision with friction; Part B: Poisson's and Stornge's hypotheses[END_REF][START_REF] Bowling | Energetically consistent simulation of simultaneous impacts and contacts in multibody systems with friction[END_REF].

Set-valued tangential contact law

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 System[END_REF][START_REF] Leine | A set-valued force law for spatial Coulomb-Contensou friction[END_REF]. The Coulomb law states that the sliding friction is proportional to the normal force of a contact. The magnitude 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 is in opposite direction to the relative velocity of the frictional contact [START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Pfeiffer | Numerical aspects of non-smooth multibody dynamics[END_REF]. Consider again the two contacting rigid bodies depicted in Fig. 2, 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 is illustrated in Fig. 2(b). If the bodies are sliding over each other, then the friction force λ T has the magnitude μλ N and acts in the direction opposed to the relative tangential velocity, that is,

-λ T = μλ N Sgn(γ T ), γ T = 0 ( 13
)
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 no-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) ( 14
)
γ T < 0 ⇒ λ T = +μλ N (negative sliding) ( 15
)
γ T > 0 ⇒ λ T = -μλ N (positive sliding) (16)
These three scenarios can be summarized by a set-valued force law as

-λ T ∈ μλ N Sgn(γ T ) (17)
Figure 3(b) shows the Coulomb's friction law as a set-valued force law [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth System[END_REF].

Dynamics of non-smooth rigid multibody systems

Equations of motion

From classical mechanics, it is well known that 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] Haug | Computer-Aided Kinematics and Dynamics of Mechanical Systems. Volume I: Basic Methods[END_REF] M uh = 0 (18)

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 represents the vector of all external and gyroscopic forces acting on the system forces originating from springs and dampers are also included in vector h, 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.

It is clear that [START_REF] Leine | Periodic motion and bifurcations induced by the Painlevé paradox[END_REF] represents a classical second-order differential equation that describes the dynamic behavior of a multibody system without any contacts and contact forces. Therefore, when a system includes frictional unilateral constraints, the occurring contact forces should be taken into account in the equations of motion. In general, the magnitudes of the normal and tangential contact forces are added to the equations of motion by using the Lagrange multiplier technique [START_REF] Nikravesh | Computer-Aided Analysis of Mechanical Systems[END_REF]. Thus, adding the contact forces to [START_REF] Leine | Periodic motion and bifurcations induced by the Painlevé paradox[END_REF], the dynamic equations of motion of a rigid multibody system with normal and tangential contact forces can, for non-impulsive motion, be written on the acceleration level as [START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF][START_REF] Pfeiffer | Numerical aspects of non-smooth multibody dynamics[END_REF]

M u -h -W N λ N -W T λ T = 0 a.e. ( 20
) q = u ∀t (21)
where W N = W N (q, t) ∈ R f ×n and W T = W T (q, t) ∈ R f ×n gather the generalized normal and tangential force directions w Ni and w T i , respectively. The normal and tangential contact forces have magnitudes λ Ni and λ T i 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 ( 20) have the same meaning as described above. It is important to note that (20) requires the existence of the velocities u as well as the existence of accelerations u. Motion without impulses implies that λ N (t) is (locally) bounded and timecontinuous. The velocities u(t) therefore exist on non-impulsive time-intervals. 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 (20), therefore, holds for almost all t . Impulsive motion is described by the impact equation,

M(u + -u -) -W N N -W T T = 0 a.e. ( 22
) u + (t) = lim t↓0 q(t + t) -q(t) t , u -(t) = lim t↑0 q(t + t) -q(t) t ( 23
)
which relates the velocity jump to the impulsive forces N and T in normal and tangential direction respectively. We assume that the velocities u(t) are of locally bounded variation (without singular part) and denote u -(t) and u + (t) as the pre-and post-impact velocity respectively. Furthermore, note that finite forces, such as gravity or reaction forces from springs and dampers, are non-impulsive, and do not occur in [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth System[END_REF]. Following Moreau [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] we will cast the non-impulsive dynamics [START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and Coulomb friction[END_REF] and the impulsive dynamics [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth System[END_REF] in a unified description, by using an equality of measures. This constitutes the general framework for non-smooth rigid multibody dynamics [START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Panagiotopoulos | Inequality constraints with elastic impacts in deformable bodies. The convex case[END_REF].

Multiplying the equation of motion [START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and Coulomb friction[END_REF] with the Lebesgue measure dt and the impact equation [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth System[END_REF] with the atomic measure dη, being the sum of the Dirac point measures at the impact times, yields

M udt -hdt -W N λ N dt -W T λ T dt = 0 (24) M u + -u -dη -W N N dη -W T T dη = 0 (25)
Addition of ( 24) and ( 25) results in

M udt + u + -u -dη -hdt -W N (λ N dt + N dη) -W T (λ T dt + T dη) = 0 (26)
or more briefly,

Mdu -hdt -W N dP N -W T dP T = 0 (27)
The differential measure for the velocities du = udt + (u +u -)dη consists of the Lebesgue measurable part udt , which accounts for absolutely continuous motion, and the atomic parts which accounts for impulsive motion. Hence, for impact free motion it holds that du = udt . 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].

In the case of non-impulsive motion, all measures dη vanish and a formal division by dt yields the classical Newton-Euler equations of motion given by [START_REF] Stewart | An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and Coulomb friction[END_REF]. The basic idea of the use of equalities of measures in multibody dynamics with unilateral constraints is to treat impulsive and non-impulsive dynamics in a unified way, i.e. with a single integration process, which opens the possibility to handle both within a single discretization [START_REF] Glocker | Formulation and preparation for numerical evaluation of linear complementarity systems in dynamics[END_REF].

Impact laws

In this paragraph, the resolution of the equations of motion expressed in the form of the equality of measures ( 27) is briefly presented and discussed in a review manner. The inclusions that are necessary to solve the frictional unilateral contact events in an autonomous multibody system, based on the Newton's impact law combined with the Coulomb's friction law, are also stated. In addition, the force laws are related to the systems' kinematics. The reader interested in the detailed description of this formulation is referred to the work of Moreau [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Moreau | Bounded Variation in Time[END_REF] and Glocker [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth System[END_REF].

Since the impenetrability condition between colliding bodies is required, let us consider a MBS with n of frictional unilateral constraints, which can be represented by n inequalities as

g Ni (q, t) ≥ 0, i = 1, . . . , n (28)
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 Euclidean 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 rigorous treatment of the definition of these inequalities, within the framework of multibody systems formulation, is presented and 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 System[END_REF]. The set of active contacts in the present work is stated as

H (t) = i|g Ni (q, t) = 0 ( 29
)
which singles out the contact(s) at which contact-impact forces may occur.

In order to define the constitutive force laws which relate the contact-impact impulse measures to the system's kinematics q and u, let us 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 ( 30
)
γ T i = w T T i u + wT i ( 31
)
where w Ni and w T i represent the generalized normal and tangential force directions, respectively, and wNi and wT i are rheonomic terms [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth System[END_REF]. The equations of motion ( 27) 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 ε T i ∈ [0, 1]. For the case of a completely elastic contact the coefficient of restitution is equal to unity, while for a perfectly inelastic contact the coefficient of restitution assumes the value of zero.

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, we allow the post-impact relative velocity to be higher than the value prescribed by Newton's impact law, with the intention 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 ( 32
)
ξ T i := γ + T i + ε T i γ - T i ( 33
)
where

(γ + Ni , γ - T i ) := (γ Ni , γ T i )(u ±
). Thus, the normal and tangential impact laws can be stated as two inclusions,

-dP Ni ∈ Upr(ξ Ni ) ( 34
)
-dP T i ∈ μ i dP Ni Sgn(ξ T i ) (35)
Finally, the complete description of the dynamics of non-smooth system, which accounts for both impact and impact-free phases, is given by ( 27)- [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]. This problem can be solved by using the Moreau time-stepping method, which is presented and discussed in the next section.

Moreau time-stepping method

Time discretization based on the Moreau midpoint rule

The time-stepping methods provide a discrete numerical scheme suitable for the simulation of non-smooth systems [START_REF] Förg | Simulation of unilateral constrained systems with many bodies[END_REF][START_REF] Moreau | Bounded Variation in Time[END_REF][START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF][START_REF] Glocker | On frictionless impact models in rigid-body systems[END_REF][START_REF] Anitescu | Time-stepping for three-dimensional rigid body dynamics[END_REF][START_REF] Studer | Step size adjustment and extrapolation for time-stepping schemes in non-smooth dynamics[END_REF]. These methods are widely used due to their simplicity to implement and their robustness. The time-stepping schemes are based on a timediscretization of the system dynamics. The whole set of discretized equations and constraints is used to compute the next state of the motion. Among the various time-stepping methods available in the literature, the Moreau midpoint method is one of the most popular and is considered in the present work [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]. The equality of measures [START_REF] Djerassi | Collision with friction; Part A: Newton's hypothesis[END_REF] together with the setvalued force laws [START_REF] Nikravesh | Computer-Aided Analysis of Mechanical Systems[END_REF] and ( 35) form a measure differential inclusion which describes the time evolution of a multibody system with discontinuities in the generalized velocities, that is, a non-smooth dynamical system. A general way to solve this mathematical problem consists of applying the Moreau time-stepping method, which does not make use of the classical equations of motion, which relate the accelerations to forces, but considers the equations of motion at the velocity level [START_REF] Djerassi | Collision with friction; Part A: Newton's hypothesis[END_REF]. The first step of the Moreau approach consists of the timediscretization of the measure differential equation. Thus, integrating [START_REF] Djerassi | Collision with friction; Part A: Newton's hypothesis[END_REF] over a small finite time interval t , of which initial and end points are denoted by the indices A and E, yields the following terms:

t Mdu ≈ M M u = M M (u E -u A ), M M = M(q M , t M ) ( 36
) t hdt = h ≈ h M t, h M = h(q M , u A , t M ) ( 37
) t W N dP N = W NM P N , W NM = W N (q M , t M ) ( 38
) t W T dP T = W T M P T , W T M = W T (q M , t M ) ( 39
)
where t M is the midpoint time instant of the compact time interval [t A , t E ] and q M = q A + 1 2 u A t is the midpoint system's position state. It is clear that the midpoint time instant can be evaluated as

t M = t A + 1 2 t ( 40
)
Finally, after the above discretization, the equations of motion expressed at the velocity level can be written as [START_REF] Glocker | Formulation and preparation for numerical evaluation of linear complementarity systems in dynamics[END_REF] M M (u Eu A )h M t -W NM P N -W T M P T = 0 [START_REF] Anitescu | Time-stepping for three-dimensional rigid body dynamics[END_REF] together with the set-valued contact/impact laws,

-P N ∈ Upr(ξ N ) ⇔ -P N ∈ N C N (ξ N ) ( 42
)
-P T ∈ μP N Sgn(ξ T ) ⇔ -P T ∈ N C T (P N ) (ξ T ) (43)
This set of algebraic inclusions can be solved with a linear complementarity problem (LCP) formulation or by an augmented Lagrangian approach (ALA) [START_REF] Leine | Dynamics and Bifurcations of Non-Smooth Mechanical Systems[END_REF]. The velocity u E , at the end of time-step t E = t A + t , is subsequently calculated by using [START_REF] Anitescu | Time-stepping for three-dimensional rigid body dynamics[END_REF]. Finally, the positions at the end of the time step are calculated by

q E = q M + 1 2 tu E (44)
Note that (42) applies only to active set-valued force laws, i ∈ H (t), i.e. set-valued force laws that can be described at the velocity level. As friction elements are naturally defined at the velocity level, they are always active and can always be described by [START_REF] Greenwood | Principles of Dynamics[END_REF]. Considering unilateral contacts, Moreau midpoint algorithm calculates the contact distances g Ni of all unilateral contacts at the midpoint q M in order to evaluate whether these are active (g Ni ≤ 0) or not (g Ni > 0). Only active unilateral contacts can be described by inclusion [START_REF] Studer | Step size adjustment and extrapolation for time-stepping schemes in non-smooth dynamics[END_REF]. Unilateral contacts that are non-active, thus open, are disregarded because it is assumed that their contact force contribution is equal to zero.

Figure 4 shows the flowchart of the general computational strategy, based on the Moreau time-stepping method, to solve the equations of motion for rigid multibody systems with frictional unilateral constraints, which can be summarized by the following steps: Fig. 4 Flowchart of the computational procedure for the solution of the equations of motion of constrained rigid multibody systems with frictional unilateral constraints (i) Start the analysis by defining the initial conditions of the problem at hand, namely the initial time t A , final time of simulation t F , time step t , together with the given initial positions q A and velocities u A ; (ii) According to the Moreau midpoint rule compute the midpoint time instant t M , the end time of the interval t E , evaluate the position's state at the midpoint instants q M , assemble the midpoint mass matrix M M and the gyroscopic and external forces vector h M , and compute the midpoint states of the potential or candidate contact-impact points H M ; (iii) Check for contact-impact between contacting bodies. If there is not any contact-impact (open contacts), then calculate the velocity at the end time u E , by using [START_REF] Anitescu | Time-stepping for three-dimensional rigid body dynamics[END_REF]; otherwise (at least one closed contact) solve the contact problem (for instance by formulating it as a LCP or by using ALA) in order to obtain the impulsive forces P N and P T required to compute u E for the contact-impact case; (iv) Compute the position's state at the end time q E , by solving (44); (v) Increment the time step. If the current time is smaller than the intended final simulation time, then update the position and velocity variables and go to step (ii) to proceed with the process of a new time step; otherwise stop the simulation.

Formulation of the contact-impact problem as an LCP

In this section, the LCP formulation to solve the contact-impact problem of multibody systems with frictional unilateral constraints is presented, which closely follows the work by Glocker and Studer [START_REF] Glocker | Formulation and preparation for numerical evaluation of linear complementarity systems in dynamics[END_REF]. In order to set up the LCP, let us first introduce the following matrix notation:

W NM := mat(w Ni (q M , t M )) ∈ R f,i , i ∈ H (45) W T M := mat(w T i (q M , t M )) ∈ R f,i , i ∈ H (46) wNM := col( wNi (q M , t M )) ∈ R i , i ∈ H ( 47
)
wT M := col( wT i (q M , t M )) ∈ R i , i ∈ H (48)

P N := col(P Ni ) ∈ R i , i ∈ H ( 49
)
P T := col(P T i ) ∈ R i , i ∈ H (50) γ NE := col(γ NEi ) ∈ R i , i ∈ H (51) γ T E := col(γ T Ei ) ∈ R i , i ∈ H (52) γ NA := col(γ NAi ) ∈ R i , i ∈ H (53) γ T A := col(γ T Ai ) ∈ R i , i ∈ H (54) ξ N := col(ξ Ni ) ∈ R i , i ∈ H (55) ξ T := col(ξ T i ) ∈ R i , i ∈ H (56) ε N := diag(ε Ni ) ∈ R i , i ∈ H (57) ε T := diag(ε T i ) ∈ R i , i ∈ H (58) μ := diag(μ i ) ∈ R i , i ∈ H (59)
Thus, the contact-impact problem of non-smooth systems can be summarized by the following mathematical relations:

M M (u E -u A ) -h M t -W NM P N -W T M P T = 0 (60) γ NE = W T NM u E + wNM (61) γ T E = W T T M u E + wT M (62) γ NA = W T NM u A + wNM (63) γ T A = W T T M u A + wT M ( 64
)
ξ N = γ NE + ε N γ NA ( 65
)
ξ T = γ T E + ε T γ T A ( 66
) -P N ∈ Upr(ξ N ) ( 67
)
-P T ∈ μP N Sgn(ξ T ) (68)
The values of γ NA and γ T A can be evaluated by using ( 63) and ( 64), respectively, since the velocities u A are known at the left endpoint of the time interval. Introducing now ( 61) and ( 62) into ( 65) and (66) yields

ξ N = W T NM u E + ( wNM + ε N γ NA ) (69) ξ T = W T T M u E + ( wT M + ε T γ T A ) (70)
Now, it should be mentioned that the inclusions for the contact-impact force laws need to be formulated as complementarity conditions. The unilateral primitive of (67) results in

-P N ∈ Upr(ξ N ) ⇔ P N ≥ 0, ξ N ≥ 0, P T N ξ N = 0 (71)
In turn, the relay function (68) has to be decomposed into two Upr's to achieve the desired complementarity conditions. Thus, (68) yields

-P T ∈ μP N Sgn(ξ T ) ⇔ ∃ξ R , ξ L s.t. ⎧ ⎪ ⎨ ⎪ ⎩ μP N + P T ≥ 0, ξ R ≥ 0, (μP N + P T ) T ξ R = 0 μP N -P T ≥ 0, ξ L ≥ 0, (μP N -P T ) T ξ L = 0 ξ T = ξ R -ξ L (72) in which the step height is [-μP N , +μP N ].
In addition, to abbreviate the complementarity conditions of (72) the impulsive friction saturations P R and P L are defined as [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth System[END_REF] P R := μP N + P T , P R ∈ R i (73)

P L := μP N -P T , P L ∈ R i (74)
together with

ξ T = ξ R -ξ L , ξ R , ξ L ∈ R i (75)
Then, the whole set of complementarity conditions of (72) can be rewritten as

0 ≤ ⎛ ⎝ ξ N ξ R P L ⎞ ⎠ ⊥ ⎛ ⎝ P N P R ξ L ⎞ ⎠ ≥ 0 (76)
The reason for the special arrangement of P L and ξ L in (76), must be sought in optimization theory. Without this special arrangement, one is not able to be set up the LCP formulation without additional matrix inversion processes [START_REF] Glocker | Set-Valued Force Laws: Dynamics of Non-Smooth System[END_REF]. Since the variables ξ T , P T and u E are not included in [START_REF] Djerassi | Collision with friction; Part B: Poisson's and Stornge's hypotheses[END_REF], they have to be eliminated. Thus, combining (60) and ( 73), yields

M M (u E -u A ) -h M t -(W NM -W T M μ)P N -W T M P R = 0 (77)
Substituting now (75) into (70) results in

ξ R = W T T M u E + ( wT M + ε T γ T A ) + ξ L ( 78
)
The elimination of variable P T can be done through the combination of ( 73) and ( 74), which can be written as

P L = 2μP N -P R (79)
Since the inversion of mass matrix M is always possible, (77) can be solved for u E :

u E = u A + M -1 M h M t + M -1 M (W NM -W T M μ)P N + M -1 M W T M P R (80)
Now, ( 63) and ( 64) are used to express W T NM u A and W T T M u A in terms of γ NA and γ T A ,

W T NM u A = γ NA -wNM (81)
W T T M u A = γ T A -wT M (82)
Introducing ( 80)-( 82) into ( 69) and (78), yields

ξ N = W T NM M -1 M h M t + W T NM M -1 M (W NM -W T M μ)P N + W T NM M -1 M W T M P R + (I + ε N )γ NA (83) ξ R = W T T M M -1 M h M t + W T T M M -1 M (W NM -W T M μ)P N + W T T M M -1 M W T M P R + (I + ε T )γ T A + ξ L (84)
Thus, (83), ( 84) and ( 79) can be written in a matrix form as

⎛ ⎝ ξ N ξ R P L ⎞ ⎠ = ⎛ ⎝ W T NM M -1 M (W NM -W T M μ) W T NM M -1 M W T M 0 W T T M M -1 M (W NM -W T M μ) W T T M M -1 M W T M I 2μ -I 0 ⎞ ⎠ ⎛ ⎝ P N P R ξ L ⎞ ⎠ + ⎛ ⎝ W T NM M -1 M h M t + (I + ε N )γ NA W T T M M -1 M h M t + (I + ε T )γ T A 0 ⎞ ⎠ (85)
Equations ( 85) together with the complementarity conditions (76) form the LCP for the contact-impact analysis of multibody systems with frictional unilateral constraints. The dimension of this LCP is 3n, where n represents the number of active contacts. The LCP (85) is solved in each integration time step. Then, the velocities u E and positions q E for the subsequent time steps are obtained from (80) and [START_REF] Flores | Contact-impact analysis in multibody systems based on the non-smooth dynamics approach[END_REF], respectively.

Computational strategy adopted

Since the Moreau time-stepping method with an LCP formulation 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. 5, is developed under the framework of MBS 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 , ε T i , μ i , w Ni , w T i , wNi and wT i ; (iii) compute the midpoint 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 T M = mat w T i (q M , t M ) wNM = col wNi (q M , t M ) wT M = col wT i (q M , t M ) γ NA = col(γ NAi ) γ T A = col(γ T Ai ) ε N = diag(ε Ni ) ε T = diag(ε T i ) μ = diag(μ i ) (v) set up the LCP in the standard form y = Ax + b: A = ⎛ ⎝ W T NM M -1 M (W NM -W T M μ) W T NM M -1 M W T M 0 W T T M M -1 M (W NM -W T M μ) W T T M M -1 M W T M I 2μ -I 0 ⎞ ⎠ b = ⎛ ⎝ W T NM M -1 M h M t + (I + ε N )γ NA W T T M M -1 M h M t + (I + ε T )γ T A 0 ⎞ ⎠ (
vi) solve the LCP using an appropriate algorithm:

(x, y) = LCP(A, b)
(vii) split the LCP solution according to

P N = col(x i ), i = 1, . . . , n i P L = col(x i ), i = n i + 1, . . . , 2n i P R = col(y i ), i = 2n i + 1, . . . , 3n i
(viii) evaluate the velocity at the end of the integration time step:

u E = u A + M -1 M h M t + M -1 M (W NM -W T M μ)P N + M -1 M W T M P R
(ix) compute the positions at the end of the integration time step:

q E = q M + 1 2 tu E (x) increment time step: t A = t A + t
(xi) update the 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.

6 Demonstrative application to a slider-crank mechanism

System's description

This section deals with the dynamic modeling and analysis of a planar slider-crank mechanism with a translation clearance joint. This multibody mechanical system consists of four rigid bodies, which represent the ground, the crank, the connecting rod and the slider. The body numbers and their center of mass are shown in Fig. 6. The ground, the crank, the connecting rod and the slider are constrained via ideal revolute joints. The center of mass of each body is considered to be located at the mid distance of the bodies' total length. The translational clearance joint is composed by a guide and a slider. This joint has a finite clearance, which is constant along the length of the slider. The crank rotates with a constant angular velocity of 150 rad/s. The initial configuration is taken with the crank and the connecting rod collinear, being the initial positions and velocities necessary to start the dynamic analysis are obtained from kinematic simulation of the slider-crank mechanism with ideal joints only. The system is under the action of gravity force, which is taken to act in the negative Z direction. Figure 7 shows a translational clearance joint. The clearance c is defined as the difference between the distance of the guide and the slider surfaces. The geometric characteristics of the translational clearance joint are the slider length 2a, the slider width 2b, and the distance between the guide surfaces d. In an ideal translational joint the two bodies translate with respect to each other parallel to the line of translation, so that, there is neither rotation between the bodies nor a relative translation motion in the direction perpendicular to the axis of the joint. The existence of a clearance in a translational joint introduces two extra degrees of freedom. Hence, the slider can move 'freely' inside the guide limits, until it reaches the guide surfaces.

The modeling of translational clearance joints is a complex task, due to the several possible contact configurations between the slider and guide. Figure 8 illustrates four different scenarios for the slider configuration relative to guide surface, namely:

(i) no contact between the two elements: the slider is in free flight motion inside the guide; (ii) one corner of the slider is in contact with the guide surface; (iii) two adjacent slider corners are in contact with the guide surface, which implies that a face of slider is in contact with the guide surface; (iv) two opposite slider corners are in contact with the guide surface. The conditions for switching from one case to another depend on the system's dynamic response as well as on the material colliding properties.

Lagrange's equations

In order for the translational clearance joint to be simulated in the multibody system environment, is it first required that the system's equations of motion be derived. In this work the Lagrange's equation of second type is used and it can be written as [START_REF] Greenwood | Principles of Dynamics[END_REF] d dt

∂L ∂ qi - ∂L ∂q i = 0, i = 1, . . . , f ( 86
)
where L is the Lagrangian of the system, that is, the difference between kinetic and potential energies, expressed in terms of the generalized coordinates and their time derivatives. The equations represented by (86) are also called as Euler-Lagrange's equations of motion, because although Lagrange was the first to formulate them specifically as the equations of motion, they were previously derived by Euler as the conditions under which a point passes from one specific place and time to another in such a way that the integral of a given function L with respect to time is stationary. Since the slider-crank mechanism represented in Fig. 6 has three degrees of freedom, three is also the number of generalized coordinates that uniquely represent the system's configuration. Furthermore, the crank, the connecting rod and the slider have masses m i and moments of inertia with respect to the principal central axes perpendicular to the plane of motion J i , where i = 1, 2 and 3. Thus, the vector of generalized coordinates and velocities are defined as

q = ⎛ ⎝ θ 1 θ 2 θ 3 ⎞ ⎠ (87) u = ⎛ ⎝ ω 1 ω 2 ω 3 ⎞ ⎠ , with q = u a.e. ( 88
)
Thus, applying the Lagrange's equation to slider-crank mechanism yields [START_REF] Flores | Contact-impact analysis in multibody systems based on the non-smooth dynamics approach[END_REF] ⎛

⎝ M 11 M 12 M 13 M 21 M 22 M 23 M 31 M 32 M 33 ⎞ ⎠ ⎛ ⎝ θ1 θ2 θ3 ⎞ ⎠ = ⎛ ⎝ h 1 h 2 h 3 ⎞ ⎠ ( 89
)
in which

M 11 = J 1 + 1 4 m 1 + m 2 + m 3 l 2 1 ( 90
)
M 12 = M 21 = 1 2 m 2 + m 3 l 1 l 2 cos(θ 2 -θ 1 ) (91) M 13 = M 31 = M 23 = M 32 = 0 ( 92
)
M 22 = J 2 + 1 4 m 2 + m 3 l 2 2
(93)

M 33 = J 3 (94) h 1 = 1 2 m 2 + m 3 l 1 l 2 sin(θ 2 -θ 1 ) θ 2 2 - 1 2 m 1 + m 2 m 3 gl 1 cos θ 1 ( 95
)
h 2 = - 1 2 m 2 + m 3 l 1 l 2 sin(θ 2 -θ 1 ) θ 2 1 - 1 2 m 2 + m 3 gl 2 cos θ 2 ( 96
)
h 3 = 0 (97)

Gap functions

In order to determine the gap functions let us consider Fig. 9 where a generic position of the slider inside the guide is illustrated with the purpose to represent the closed kinematic chain of each potential contact point. From analysis of Fig. 9 and considering the system kinematics, the mathematical expressions of the gap functions can be written as [START_REF] Flores | Contact-impact analysis in multibody systems based on the non-smooth dynamics approach[END_REF]

g N1 = d 2 -l 1 sin θ 1 -l 2 sin θ 2 + a sin θ 3 -b cos θ 3 (98) g T 1 = l 1 cos θ 1 + l 2 cos θ 2 -a cos θ 3 -b sin θ 3 (99) g N2 = d 2 -l 1 sin θ 1 -l sin θ 2 -a sin θ 3 -b cos θ 3 (100) g T 2 = l cos θ 1 + l 2 cos θ 2 + a cos θ 3 -b sin θ 3 (101) g N3 = d 2 + l 1 sin θ 1 + l 2 sin θ 2 -a sin θ 3 -b cos θ 3 (102) g T 3 = l 1 cos θ 1 + l 2 cos θ 2 -a cos θ 3 + b sin θ 3 (103) g N4 = d 2 + l 1 sin θ 1 + l 2 sin θ 2 + a sin θ 3 -b cos θ 3 (104) g T 4 = l 1 cos θ 1 + l 2 cos θ 2 + a cos θ 3 + b sin θ 3 (105)
Then, the w vectors and of the w scalars associated with each contact point can be obtained as

w N1 = ∂g N1 ∂q = ⎛ ⎝ -l 1 cos θ 1 -l 2 cos θ 2 a cos θ 3 + b sin θ 3 ⎞ ⎠ (106) w T 1 = ∂g T 1 ∂q = ⎛ ⎝ -l 1 sin θ 1 -l 2 sin θ 2 a sin θ 3 -b cos θ 3 ⎞ ⎠ (107) w N2 = ∂g N2 ∂q = ⎛ ⎝ -l 1 cos θ 1 -l 2 cos θ 2 -a cos θ 3 + b sin θ 3 ⎞ ⎠ (108)
Fig. 9 Generic position of the slider inside the guide where the distance between guide upper and lower surfaces is exaggerated for illustration purpose

w T 2 = ∂g T 2 ∂q = ⎛ ⎝ -l 1 sin θ 1 -l 2 sin θ 2 -a sin θ 3 -b cos θ 3 ⎞ ⎠ (109) w N3 = ∂g N3 ∂q = ⎛ ⎝ l 1 cos θ 1 l 2 cos θ 2 -a cos θ 3 + b sin θ 3 ⎞ ⎠ (110) w T 3 = ∂g T 3 ∂q = ⎛ ⎝ -l 1 sin θ 1 -l 2 sin θ 2 a sin θ 3 + b cos θ 3 ⎞ ⎠ (111) w N4 = ∂g N4 ∂q = ⎛ ⎝ l 1 cos θ 1 l 2 cos θ 2 a cos θ 3 + b sin θ 3 ⎞ ⎠ (112) w T 4 = ∂g T 4 ∂q ⎛ ⎝ -l 1 sin θ 1 -l 2 sin θ 2 -a sin θ 3 + b cos 3 ⎞ ⎠ (113)
wN1 = wT 1 = wN2 = wT 2 = wN3 = wT 3 = wN4 = wT 4 = 0 (114)

Results and discussion

The geometrical characteristics, the inertial properties, the force elements, the contact parameters and the initial conditions necessary to perform the dynamic analysis of the slidercrank mechanism with a translational clearance joint are listed in Table 1. Figure 10 shows the corners motion in a dimensionless form for two full crank rotations, in which the free slider motion and contact-impact events can be observed. Figure 11 illustrates the crank speed, the connecting-rod speed and the portraits relative to connecting rod and slider for two complete crank rotations. The normalized slider corner motions are evaluated using the following relation:

Contact parameters

ε N 1 = ε N 2 = ε N 3 = ε N 4 = 0.4 ε T 1 = ε T 2 = ε T 3 = ε T 4 = 0.0 μ 1 = μ 2 = μ 3 = μ 4 =
y i -c c , (i = 1, 2, 3, 4) (115)
where y i represents the y coordinate of the slider corners and c is the clearance size.

The dimensionless slider trajectories are shown in Fig. 10, where the different types of motion between the slider and guide observed are associated with the different guideslider configurations, i.e., no contact, impact followed by rebound and permanent contact between the joint elements. The effects of impact between the slider and guide surfaces are also quite visible in the plots of Figs. 11(b) and 11(c), namely, one can observe the discontinuities in the connecting-rod speed. On the other hand, the smooth changes in the speed indicate that the slider and guide surfaces are in permanent contact for long periods, as is illustrated in the slider portrait of Fig. 11(d). It should be highlighted that some numerical difficulties can arise when the clearance size is very small, which will lead to the wellknown drift problem. In these situations, one possible way to overcome those difficulties consists of a projection technique, in which the excessive penetration between the slider and guide surfaces is eliminated in each time step in order to avoid the further interpretation of the bodies. When this scheme is implemented, special attention should be paid to the conservation of the system's energy, since it can lead to overestimated total system energy associated with the contact-impact phenomena [START_REF] Leine | Stability and Convergence of Mechanical Systems with Unilateral Constraints[END_REF]. Figure 12 shows the influence of the value of the restitution coefficient on the dimensionless motion of the corner 1. For this purpose four different restitution coefficient values were considered, namely, 0.1, 0.4, 0.6 and 0.9. From these plots it can be observed that the methodology is valid for different set of material properties, being the system's response different when restitution coefficient varies. That is, for lower values, the rebounds are fewer and the slider and guide tend to have long periods of permanent or continuous contact, as is illustrated in Fig. 12(a). For higher values of the restitution coefficient, the free flight motion of the slider inside the guide is dominant, as Fig. 12(d) shows. Consequently, the slider portraits phases are also affected by the value of the coefficient of restitution, as it can be observed in the plots of Fig. 13.

Concluding remarks

A comprehensive investigation of contact-impact analysis in multibody systems based on the non-smooth dynamics approach was presented in this work. The methodology was based on the non-smooth dynamics approach, in which the interaction of the colliding bodies is modeled with multiple frictional unilateral constraints. The dynamics of rigid multibody systems were stated as an equality of measures, which were formulated at the velocity-impulse level. The equations of motion were complemented with constitutive laws for the forces and impulses in normal and tangential directions. The formulation of the generalized contactimpact kinematics in the normal and tangential directions was performed by obtaining a geometric relation for the gaps of the candidate contact points. The gaps were expressed as functions of the generalized coordinates. The candidate contact points were modeled as hard In this work, 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 contact-impact problem was formulated and solved as a linear complementarity problem and with the augmented Lagrangian approach, which were embedded in the Moreau time-stepping method. Finally, the effectiveness of the presented methodologies was demonstrated through the study of the slider-crank mechanism with a translational clearance joint. The main results obtained from this research work showed that the effect of the contact-impact phenomena can have a predictable nonlinear behavior. This nonlinearity aspect is more evident when the system includes friction phenomena. With the knowledge of nonlinearities in multibody systems, chaotic behavior may be eliminated with suitable design and/or parameter changes of a mechanical system. This feature plays a crucial role in the dynamics, design and control of general multibody systems of common application.

Fig. 1

 1 Fig. 1 (a) The map x → Upr(x); (b) The map x → Sgn(x); (c) The decomposition Sgn(x) into Upr(x)

Fig. 2

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

Fig. 3

 3 Fig. 3 (a) Signorini's normal contact law; (b) Coulomb's friction law

Fig. 5

 5 Fig. 5 Flowchart of the Moreau time-stepping algorithm with an LCP formulation

Fig. 6 Fig. 7 Fig. 8

 678 Fig. 6 Slider-crank mechanism with a translational clearance joint

Geometrical characteristics l 1 = 1 = 7 . 4 × 10 -5 kg m 2 J 2 = 5 . 9 × 10 -4 kg m 2 J 3 = 2 . 7 × 10 -6 kg m 2 Force elements g = 9 .

 11741022591023271029 0.1530 m l 2 = 0.3060 m a = 0.0500 m b = 0.0250 m c = 0.0010 m Inertial properties m 1 = 0.0380 kg m 2 = 0.0380 kg m 3 = 0.0760 kg J 81 m/s 2

θ 10 =Fig. 10

 1010 Fig. 10 Dimensionless motion of the slider corners

Fig. 11 (

 11 Fig. 11 (a) Crank speed; (b) connecting-rod speed; (c) connecting-rod portrait; (d) slider portrait

Fig. 12 9 Fig. 13

 12913 Fig. 12 Influence of the restitution coefficient on the dimensionless motion of corner 1: (a) ε = 0.1; (b) ε = 0.4; (c) ε = 0.6; (d) ε = 0.9

Table 1

 1 Parameters used in the dynamic simulation of the slider-crank mechanism

Acknowledgements This work is supported by the Portuguese Foundation for the Science and Technology under a postdoctoral scholarship (SFRH/BPD/40067/2007). We would also like to acknowledge the anonymous reviewers for their insightful comments and suggestions on an earlier draft of this paper.