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

Paulo Flores, Remco Leine, Christoph Glocker

► To cite this version:

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. Nonlinear Dynamics, 2012, 69, pp.2117 - 2133. 10.1007/s11071-012-0413-3. hal-01435627

HAL Id: hal-01435627

https://hal.science/hal-01435627

Submitted on 14 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion 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.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Application of the nonsmooth dynamics approach to model and analysis of the contact-impact events in cam-follower systems

Paulo Flores · Remco Leine · Christoph Glocker

**Abstract** The dynamic modeling and analysis of planar rigid 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. 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 in 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 set-valued force law of the type of Coulomb's law for dry friction. The resulting contact-impact problem is formulated and solved as an augmented Lagrangian approach, which is embedded in

the Moreau time-stepping method. The effectiveness of the methodologies presented in this work is demonstrated throughout the dynamic simulation of a cam-follower system of an industrial cutting file machine.

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

## 1 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 [1–5]. 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 [6–8]. 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 ex-
tremely important features in most of multibody sys-
tems. As a consequence, multiple impacts may oc-
cur quite frequently and become a key point for the
numerical simulation of multibody systems. Hurmu-
zlu and Marghitu [9] studied the contact problem in
multibody systems, where a planar rigid-body kine-
matic chain undergoes an external impact and an arbi-
trary number of internal impacts. Based on the Keller's
work [10], they developed a differential-integral ap-
proach and used different models for coefficient of
friction. Han and Gilmore [11] proposed a similar ap-
proach, using an algebraic formulation of the equa-
tions of motion, the Poisson's model of restitution
and the Coulomb's law to define the tangential mo-
tion. Different conditions that characterize the motion
(slipping, sticking, and reverse sliding) were detected
by analyzing velocities and accelerations at the con-
tact points. Han and Gilmore confirmed their simula-
tion results with experiments for two-body and three-
body impacts. Based on a canonical form of the equa-
tions of motion, Pereira and Nikravesh [12] presented
a methodology that solves this problem in the context
of multibody dynamics impact. Haug et al. [13] 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 [14].

Mechanical systems are often modeled as multi-
body systems with some degree of non-smoothness.
Typical examples in mechanics are the noise and vi-
bration produced in railway brakes, impact print ham-
mers, percussion drilling machines, or chattering of
machine tools. These effects are due to the nonsmooth
characteristics such as clearances, impacts, intermit-
tent contacts, dry friction, or a combination of these
effects [15–20]. In nonsmooth systems, the time evo-
lution of the displacement and the velocities is not re-
quested to be smooth. Due to the possible impacts, the
velocities are even allowed to undergo jumps at cer-
tain time instances in order to fulfill the kinematical
restrictions.

Contact-impact behavior strongly depends on the
material properties of the colliding surfaces, the na-
ture 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 en-
gineering. In addition, contact-impact events can fre-

quently occur in multibody systems and in many en-
gineering applications the function of mechanical sys-
tems is based on them. Common examples can be de-
scribed by the contact between tire and road in ve-
hicles, wheel and rail in railway systems, contact in
robotics and grasping machines, cam and follower
mechanisms, just to mention a few [21–25]. As a re-
sult of an impact, the values of the system state vari-
ables change very fast, eventually looking like dis-
continuities in the system velocities. The knowledge
of the peak forces developed in the impact process
is very important for the dynamic analysis of multi-
body systems and has consequences in the design pro-
cess. Thus, the selection of the most adequate contact-
impact method used to describe the process correctly
plays a key role in the accurate design and analy-
sis of these kinds of systems [26–28]. In a broad
sense, the different methods to solve the impact prob-
lem in multibody mechanical systems are continu-
ous and discontinuous approaches [29–31]. Within
the continuous approach, the methods commonly used
are the continuous force model, which is in fact a
penalty method, and the unilateral constraint method-
ology, based on complementarity approaches [32–35].
The compliant continuous contact force models, com-
monly 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 ex-
pressed as a continuous function of penetration be-
tween 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 mod-
eling in multibody systems have attracted the attention
of many researchers [36–47]. Assuming that the con-
tacting bodies are truly rigid, as opposed to locally de-
formable 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 pen-
etration from occurring. Thus, at the core of the com-
plementarity approach is an explicit formulation of
the unilateral constraints between the contacting rigid
bodies [48]. One of the main features of unilateral con-
straints 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 [49].

When dealing with the formulation of frictional unilateral constraints, it is possible to distinguish active and passive set-valued force laws. An active set-valued 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 [50–52].

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 set-valued 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.

## 2 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 [53]. 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 [52]. 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 $\lambda_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 \ge 0$, only two situations can occur, namely

$$g_N = 0 \land \lambda_N \ge 0 \quad \text{(closed contact)} \quad (1)$$

$$g_N > 0 \land \lambda_N = 0 \quad \text{(open contact)} \quad (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,

$$g_N \lambda_N = 0 \quad (3)$$
<figure><img src="image_2.png" /><figcaption>Fig. 1 (a) Relative normal gap; (b) Normal and tangential contact forces</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2 (a) Signorini's normal contact law; (b) Coulomb's friction law</figcaption></figure>

Thus, the relation between the normal gap and normal contact force can be described by

$$g_N \ge 0, \quad \lambda_N \ge 0, \quad g_N \lambda_N = 0 \quad (4)$$

which represents the inequality complementarity condition between $g_N$ and $\lambda_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 $\lambda_N$ [46].

The magnitude of contact force is denoted by $\lambda_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, $\epsilon_N$ [33].

The normal contact law can also be expressed as the subdifferential of the indicator function $\Psi_{C_N}$ to convex set $C_N$ [46]

$$g_N \in \partial \Psi_{C_N}(-\lambda_N) \quad (5)$$

where $C_N$ is the admissible set of negative contact forces $-\lambda_N$ as

$$C_N = \{-\lambda_N \in \mathbb{R} \mid \lambda_N \ge 0\} = \mathbb{R}^{-} \quad (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

$$-\lambda_N \in \partial \Psi_{C_N}^*(g_N) \quad (7)$$

Alternatively, the normal contact law can be formulated in a compact form by means of the normal cone of $C_N$ [46]

$$g_N \in N_{C_N}(-\lambda_N) \quad (8)$$

where $N_{C_N}(\lambda_N)$ denotes the normal cone of the convex set $C_N$ at $\lambda_N$, and $C_N$ is the set of admissible contact forces (6). 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

$$\gamma_N \in N_{C_N}(-\lambda_N), \quad g_N = 0 \quad (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 [44].

The classical Coulomb's friction law is another typical example that can be considered as a set-valued force law [33]. 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 opposite to the relative tangential velocity, that is,

$$-\lambda_T = \mu\lambda_N \operatorname{Sgn}(\gamma_T) \quad \gamma_T \neq 0 \qquad (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 \quad \Rightarrow \quad |λ_T| ≤ μλ_N \quad (sticking) \qquad (11)$$

$$γ_T < 0 \quad \Rightarrow \quad λ_T = +μλ_N \quad (negative sliding) \qquad (12)$$

$$γ_T > 0 \quad \Rightarrow \quad λ_T = -μλ_N \quad (positive sliding) \qquad (13)$$

These three scenarios can be summarized by a set-valued force law as [33]

$$λ_T \in \begin{cases} -μλ_N & γ_T > 0 \\ [-1, 1]μλ_N & γ_T = 0 \\ μλ_N & γ_T < 0 \end{cases} \qquad (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 [32]

$$C_T(λ_N) = \{-λ_T | -μλ_N ≤ λ_T ≤ +μλ_N\} \qquad (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) \qquad (16)$$

or

$$γ_T ∈ N_{C_T(λ_N)}(-λ_T) \qquad (17)$$

Alternatively, Eqs. (16) and (17) can be written as [49]

$$-λ_T ∈ ∂Ψ_{C_T(λ_N)}^*(γ_T) \qquad (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 [54].

### 3 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 [32]

$$M\dot{u} - h = 0 \qquad (19)$$

$$\dot{q} = u \quad ∀t \qquad (20)$$

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 $\dot{u} = \dot{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 (19) 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 uni- lateral contacts are open. The dependence of the sys- tem 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 classi- cal generalized system coordinates and evaluating La- grange's equations of second type or the associated virtual work expressions [55].

Adding the contact forces to Eq. (19), the dynamic equations of motion of a multibody system with nor- mal and tangential contact forces during an impact can be written at the acceleration level as [32]

$$
\mathbf{M}\dot{\mathbf{u}} - \mathbf{h} - \mathbf{w}_N\boldsymbol{\lambda}_N - \mathbf{w}_T\boldsymbol{\lambda}_T = \mathbf{0} \quad (21)
$$

$$
\dot{\mathbf{q}} = \mathbf{u} \quad \forall t \tag{22}
$$

where $\mathbf{w}_N = \mathbf{w}_N(\mathbf{q}, t) \in \mathbb{R}^f$ and $\mathbf{w}_T = \mathbf{w}_T(\mathbf{q}, t) \in \mathbb{R}^f$ represent the generalized normal and tangential force directions, respectively. The normal and tangen- tial contact forces have magnitudes $\lambda_{Ni}$ and $\lambda_{Ti}$ for each contact point $i$. The dual variables to the normal contact forces $\lambda_N$ are the variations of normal gap dis- tances $g_N$, while the dual variables to the generalized friction or tangential forces $\lambda_T$ are the variations of the generalized sliding velocities $\gamma_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. There-
fore, 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 ve-
locity at the impact itself, a meaningless term already
from the physical point of view. For the case of im-
pacts in multibody systems, Eq. (21) have to be sub-
stituted by a more suitable formulation, which consists
of use the equalities of measures firstly introduced by
Moreau [56] 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 ve-
locity of a frictional contact reverses its sign. The ac-
celeration dot u is not defined when λ_T(t) is discontin-
ous. The set of time instances for which λ_T(t) is dis-
continuous is of measure zero and Eq. (21), therefore,

holds for almost all *t*, that is, $\dot{\mathbf{q}} = \mathbf{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 (21). 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 inte-
grated over a singleton in time yields

$$
\mathbf{M}(\mathbf{u}^{+} - \mathbf{u}^{-}) - \mathbf{w}_{N}\mathbf{\Lambda}_{N} - \mathbf{w}_{T}\mathbf{\Lambda}_{T} = \mathbf{0} \quad (23)
$$

$$
\dot{\mathbf{q}} = \mathbf{u} \quad \text{a.e.} \tag{24}
$$

in which **u**⁻ and **u**⁺ represent the pre and post im-
pact velocities, **Λ**<sub>N</sub> and **Λ**<sub>T</sub> denote the normal and
tangential impulsive forces, being the remaining terms
defined as previously. Note that contact forces are re-
placed 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, there-
fore, 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

$$
\mathbf{M} \dot{\mathbf{u}} \mathrm{d}t - \mathbf{h} \mathrm{d}t - \mathbf{w}_N \boldsymbol{\lambda}_N \mathrm{d}t - \mathbf{w}_T \boldsymbol{\lambda}_T \mathrm{d}t = \mathbf{0} \quad (25)
$$

$$
\mathbf{M}(\mathbf{u}^{+} - \mathbf{u}^{-}) \mathrm{d}\eta - \mathbf{w}_N \mathbf{\Lambda}_N \mathrm{d}\eta - \mathbf{w}_T \mathbf{\Lambda}_T \mathrm{d}\eta = \mathbf{0} \quad (26)
$$

Adding now Eqs. (25) and (26) results in

$$
\begin{equation}
\begin{split}
& \mathbf{M}[\dot{\mathbf{u}} \mathrm{d}t + (\mathbf{u}^+ - \mathbf{u}^-) \mathrm{d}\eta] - \mathbf{h} \mathrm{d}t \\
& \quad - \mathbf{w}_N (\boldsymbol{\lambda}_N \mathrm{d}t + \boldsymbol{\Lambda}_N \mathrm{d}\eta) \\
& \quad - \mathbf{w}_T (\boldsymbol{\lambda}_T \mathrm{d}t + \boldsymbol{\Lambda}_T \mathrm{d}\eta) = \mathbf{0}
\end{split}
\tag{27}
\end{equation}
$$

or more briefly

$$
\mathbf{M} \, \mathrm{d}\mathbf{u} - \mathbf{h} \, \mathrm{d}t - \mathbf{w}_N \, \mathrm{d}\mathbf{P}_N - \mathbf{w}_T \, \mathrm{d}\mathbf{P}_T = \mathbf{0} \quad (28)
$$

where the Lebesgue measure is represented by *dt* and
*dη* represents the sum of the Dirac impulse measures
at the impact times [49].

The measure for the velocities du =dot ut + (u+ - u-)dη is split in Lebesgue measurable part u dot t,
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 = ṁ 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η [43, 44]. 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 (21).

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}(\mathbf{q}, t) \ge 0, \quad i = 1, \dots, n \qquad (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 [32] and Glocker [33].

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

$$H(t) = \{i | g_{Ni}(\mathbf{q}, t) = 0\} \qquad (30)$$

which singles out the contact at which contact-impact forces may occur. In fact, ideally when $g_{Ni}(\mathbf{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. [47].

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 [57]

$$\gamma_{Ni} = \mathbf{w}_{Ni}^{T}\mathbf{u} + \tilde{w}_{Ni} \qquad (31)$$

$$\gamma_{Ti} = \mathbf{w}_{Ti}^{T}\mathbf{u} + \tilde{w}_{Ti} \qquad (32)$$

where $\mathbf{w}_{Ni}$ and $\mathbf{w}_{Ti}$ represent the generalized normal and tangential force directions, respectively, and $\tilde{w}_{Ni}$ and $\tilde{w}_{Ti}$ are the partial derivatives of the of the gap functions with respect to time.

The equations of motion (28) 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 $\epsilon_{Ni} \in [0, 1]$. The Coulomb’s friction law is used for the tangential direction with coefficient of friction $\mu_i$, which is complemented by a tangential coefficient of restitution $\epsilon_{Ti} \in [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 [57]

$$
\xi_{Ni} := \gamma_{Ni}^{+} + \varepsilon_{Ni} \gamma_{Ni}^{-} \quad (33)
$$

$$
\xi_{Ti} := \gamma_{Ti}^{+} + \varepsilon_{Ti} \gamma_{Ti}^{-} \quad (34)
$$

where $(\gamma_{Ni}^{\pm}, \gamma_{Ti}^{\pm}) := (\gamma_{Ni}, \gamma_{Ti})(\mathbf{u}^{\pm})$.

Thus, normal and tangential impact laws can be
stated as

$$
-d\mathbf{P}_{Ni} \in \partial \Psi_{C_{Ni}}^* (\xi_{Ni}) \quad (35)
$$

$$
-d\mathbf{P}_{Ti} \in \partial \Psi_{C_{Ti}(\mathbf{P}_{Ni})}^* (\xi_{Ti}) \quad (36)
$$

or, in terms of the normal cones

$$
\xi_{Ni} \in N_{C_{Ni}}(-d\mathbf{P}_{Ni}) \quad (37)
$$

$$
\xi_{Ti} \in N_{C_{Ti}(\mathbf{P}_{Ni})}(-d\mathbf{P}_{Ti}) \quad (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)–(38).

4 Augmented Lagrangian approach to solve
contact-impact events

In this section, the augmented Lagrangian approach
to solve the contact-impact problem of multibody sys-
tems with frictional unilateral constraints is presented.
In a simple way, this approach consists of transform-
ing the inclusions in normal and tangential directions
into equivalent equations using proximal point of con-
 vex analysis, then the problem is solved iteratively
as a proximal point formulation [43, 54, 58]. After
discretization, the velocity-impulse equations of mo-
tion for nonsmooth mechanical systems can be written
as [14]

$$
\mathbf{M}_M \Delta \mathbf{u} - \mathbf{h}_M \Delta t - \mathbf{W}_{NM} \mathbf{P}_N - \mathbf{W}_{TM} \mathbf{P}_T = \mathbf{0} \quad (39)
$$

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,

$$
-d\mathbf{P}_N \in \partial \Psi_{C_N}^* (\xi_N) \quad (40)
$$

$$
-d\mathbf{P}_T \in \partial \Psi_{CT}^* (\mathbf{P}_N)(\xi_T) \quad (41)
$$

where $\Psi_{C_N}^*$ and $\Psi_{CT(\mathbf{P}_N)}^*$ denote the support functions
of the indicator functions $\Psi_{C_N}$ and $\Psi_{CT(\mathbf{P}_N)}$, respec-
tively, being $N_{C_N}$ and $N_{CT(\mathbf{P}_N)}$ the normal cones as-
sociated with the convex sets $C_N$ and $C_T$. These two
convex sets represent the admissible normal and tan-
gential impulsive forces given by [43, 59]

$$
C_N = \{-\mathbf{P}_N \in \mathbb{R}^n | \mathbf{P}_N \ge \mathbf{0}\} \quad (42)
$$

$$
C_{Ti}(\mathbf{P}_{Ni}) = \{-\mathbf{P}_{Ti} \| \mathbf{P}_{Ti} \| \le \mu_i P_{Ni}; i = 1, \dots, n\} \quad (43)
$$

The $\xi_N$ and $\xi_T$ parameters are expressed as [44]

$$
\xi_N = \gamma_{NE} + \varepsilon_N \gamma_{NA} \quad (44)
$$

$$
\xi_T = \gamma_{TE} + \varepsilon_T \gamma_{TA} \quad (45)
$$

in which the actual normal and tangential velocities
are given by [52]

$$
\gamma_{NE} = \mathbf{W}_{NM}^{\mathrm{T}} \mathbf{u}_{E} + \tilde{\mathbf{w}}_{NM} \quad (46)
$$

$$
\gamma_{TE} = \mathbf{W}_{TM}^{\mathrm{T}} \mathbf{u}_{E} + \tilde{\mathbf{w}}_{TM} \quad (47)
$$

$$
\gamma_{NA} = \mathbf{W}_{NM}^{\mathrm{T}} \mathbf{u}_{A} + \tilde{\mathbf{w}}_{NM} \quad (48)
$$

$$
\gamma_{TA} = \mathbf{W}_{TM}^{\mathrm{T}} \mathbf{u}_{A} + \tilde{\mathbf{w}}_{TM} \quad (49)
$$

In short, in each time step the equations of motion (39) and the set-valued force laws (40) and (41) have to be solved for Δ**u**, **P**<sub>N</sub> and **P**<sub>T</sub>. The augmented Lagrangian approach is an elegant way to solve this problem, which transforms the set of algebraic inclusions to a constrained optimization problem [60–62]. By definition of proximal point to a convex set **C** yields that

$$
\operatorname{prox}_C(x) = \operatorname*{argmin}_{\forall x^* \in C} \|x - x^*\| \quad (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:

$$
\mathbf{P}_N = \operatorname{prox}_{C_N}(\mathbf{P}_N - r\xi_N) \quad (51)
$$

$$
\mathbf{P}_T = \operatorname{prox}_{C_T(\mathbf{P}_N)}(\mathbf{P}_T - r\xi_T) \quad (52)
$$

with the convex sets $C_N$ and $C_T$ given by Eqs. (42) and (43). 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. [63].

Finally, the contact-impact problem of nonsmooth
systems based on the augmented Lagrangian approach
can be summarized by the following mathematical re-
lations:

$$
\begin{equation}
\begin{split}
\mathbf{M}_M \Delta \mathbf{u} - \mathbf{h}_M \Delta t - \mathbf{W}_{NM} \Pi_N(\mathbf{P}_N, \boldsymbol{\xi}_N) \\
\qquad - \mathbf{W}_{TM} \Pi_T(\mathbf{P}_N, \mathbf{P}_T, \boldsymbol{\xi}_T) = \mathbf{0}
\end{split}
\tag{53}
\end{equation}
$$

$$
\Pi_N(\mathbf{P}_N, \xi_N) = \operatorname{prox}_{C_N}(\mathbf{P}_N - r\xi_N) \quad (54)
$$

$$
\Pi_T(\mathbf{P}_N, \mathbf{P}_T, \xi_T) = \operatorname{prox}_{C_T(\mathbf{P}_N)}(\mathbf{P}_T - r\xi_T) \quad (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 [43, 54, 61]. 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 LU-decomposition 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. [43].

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 $\mathbf{P}_N^k$ and $\mathbf{P}_T^k$;

(ii) Solve $\Delta\mathbf{u}^{k+1}$ from the following equation:

$$ \mathbf{M}_M \Delta \mathbf{u}^{k+1} = \mathbf{h}_M \Delta t + \mathbf{W}_{NM} \mathbf{P}_N^k + \mathbf{W}_{TM} \mathbf{P}_T^k $$

(iii) Evaluate $\xi_{NE}^{k+1}$ and $\xi_{TE}^{k+1}$ from equations

$$ \xi_N^{k+1} = \mathbf{W}_{NM}^T (\mathbf{u}_A^k + \Delta \mathbf{u}^{k+1}) + \tilde{\mathbf{w}}_{NM} $$

$$ \xi_T^{k+1} = \mathbf{W}_{TM}^T (\mathbf{u}_A^k + \Delta \mathbf{u}^{k+1}) + \tilde{\mathbf{w}}_{TM} $$

(iv) Project $\mathbf{P}_N^{k+1}$ and $\mathbf{P}_T^{k+1}$ as

$$ \mathbf{P}_N^{k+1} = \Pi_N (\mathbf{P}_N^k, \xi_N^{k+1}) = \text{prox}_{C_N} (\mathbf{P}_N^k - r\xi_N^{k+1}) $$

$$ \mathbf{P}_T^{k+1} = \Pi_T (\mathbf{P}_N^{k+1}, \mathbf{P}_T^k, \xi_T^{k+1}) $$

$$ = \text{prox}_{C_T(\mathbf{P}_N)} (\mathbf{P}_T^k - r\xi_T^{k+1}) $$

(v) Compute the tolerance error as

$$ \text{error} = \| \mathbf{P}_N^{k+1} - \mathbf{P}_N^k \| + \| \mathbf{P}_T^{k+1} - \mathbf{P}_T^k \| $$

(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, \Delta t, \mathbf{q}_A$, and $\mathbf{u}_A$;

(ii) Define the geometrical, inertial, and material functions, $g_{Ni}, \mathbf{M}, \mathbf{h}, \varepsilon_{Ni}, \varepsilon_{Ti}, \mu_i, \mathbf{w}_{Ni}, \mathbf{w}_{Ti}, \tilde{w}_{Ni}$, and $\tilde{w}_{Ti}$;

(iii) Compute the mid point state variables:

$$ t_M = t_A + \frac{1}{2} \Delta t $$

$$ \mathbf{q}_M = \mathbf{q}_A + \frac{1}{2} \Delta t \mathbf{u}_A $$

$$ \mathbf{M}_M = \mathbf{M}(\mathbf{q}_M, t_M) $$

$$ \mathbf{h}_M = \mathbf{h}(\mathbf{q}_M, \mathbf{u}_A, t_M) $$

$$ g_{Ni} = g_{Ni}(\mathbf{q}_M, t_M) $$

$$ H_M = \{i | g_{Ni}(\mathbf{q}_M, t_M) \le 0\} $$

$$ n_i = \text{length}(H_M) $$

(iv) For every $i \in H_M$ evaluate:

$$ \mathbf{W}_{NM} = \text{mat}(\mathbf{w}_{Ni}(\mathbf{q}_M, t_M)), $$

$$ \mathbf{W}_{TM} = \text{mat}(\mathbf{w}_{Ti}(\mathbf{q}_M, t_M)), $$

$$ \tilde{\mathbf{w}}_{NM} = \text{col}(\tilde{w}_{Ni}(\mathbf{q}_M, t_M)), $$

$$ \tilde{\mathbf{w}}_{TM} = \text{col}(\tilde{w}_{Ti}(\mathbf{q}_M, t_M)), $$

$$ \gamma_{NA} = \text{col}(\gamma_{NAi}), $$

$$ \gamma_{TA} = \text{col}(\gamma_{TAi}), $$

$$ \varepsilon_N = \text{diag}(\varepsilon_{Ni}), $$

$$ \varepsilon_T = \text{diag}(\varepsilon_{Ti}), $$

$$ \mu = \text{diag}(\mu_i), \quad k = 1, $$

$$ \mathbf{P}_N^k = \mathbf{P}_N(H_M), $$

$$ \mathbf{P}_T^k = \mathbf{P}_T(H_M) $$

(v) Compute while the Newton method does not converge, i.e., while error > tol:

$$ \mathbf{u}_E = \mathbf{u}_A + \mathbf{M}_M^{-1} \mathbf{h}_M \Delta t + \mathbf{M}_M^{-1} \mathbf{W}_{NM} \mathbf{P}_N^k + \mathbf{M}_M^{-1} \mathbf{W}_{TM} \mathbf{P}_T^k $$

$$ \gamma_{NE} = \mathbf{W}_{NM}^T \mathbf{u}_E + \tilde{\mathbf{w}}_{NM} $$

$$ \gamma_{TE} = \mathbf{W}_{TM}^T \mathbf{u}_E + \tilde{\mathbf{w}}_{TM} $$

$$ \xi_N^{k+1} = \gamma_{NE} + \varepsilon_N \gamma_{NA} $$

$$ \xi_T^{k+1} = \gamma_{TE} + \varepsilon_T \gamma_{TA} $$

$$ \mathbf{P}_N^{k+1} = \Pi_N (\mathbf{P}_N^k, \xi_N^{k+1}) = \text{prox}_{C_N} (\mathbf{P}_N^k - r\xi_N^{k+1}) $$

$$ \mathbf{P}_T^{k+1} = \Pi_T (\mathbf{P}_N^{k+1}, \mathbf{P}_T^k, \xi_T^{k+1}) = \text{prox}_{C_T(\mathbf{P}_N)} (\mathbf{P}_T^k - r\xi_T^{k+1}) $$

$$ \text{error} = \| \mathbf{P}_N^{k+1} - \mathbf{P}_N^k \| + \| \mathbf{P}_T^{k+1} - \mathbf{P}_T^k \| $$

$$ k = k + 1 $$
<figure><img src="image_4.png" /><figcaption>Fig. 3 Computational strategy for the iterative Newton method</figcaption></figure>

(vi) Compute the positions at the end of the integration time step:

$$ \mathbf{q}_E = \mathbf{q}_M + \frac{1}{2} \Delta t \mathbf{u}_E $$

(vii) Increment time step:

$$ t_A = t_A + \Delta t $$

(viii) Update system states' variables $\mathbf{q}_A = \mathbf{q}_E$ and $\mathbf{u}_A = \mathbf{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.

## 5 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 [64–66]. 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 up-dead-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
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Flowchart of the Moreau time-stepping algorithm with augmented Lagrangian approach</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5** (a) Overall view of the cutting file machine; (b) Schematic representation of the corresponding mechanical system</figcaption></figure>
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 pre-load 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 [65].

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 \times n_b - n_c \qquad (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 cam-follower mechanisms are designed with pressure angles as small as possible [67, 68].

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 cam-follower motion repeats itself six times in each complete cam rotation [66]. 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;
<figure><img src="image_7.png" /><figcaption>Fig. 6 (a) Kinematic configuration of the cam follower mechanism; (b) Follower displacement</figcaption></figure>

<table><caption>**Table 1** Parameters used in the dynamic simulation of the cam follower system</caption><tr><td>Follower mass—<i>m</i></td><td>1.0 kg</td></tr><tr><td>Cam speed—<i>ω</i></td><td>20.94 rad/s</td></tr><tr><td>Cam base radius—<i>R</i><sub>b</sub></td><td>0.003 m</td></tr><tr><td>Maximum follower stroke—<i>h</i></td><td>0.017 m</td></tr><tr><td>Pre-load spring</td><td>0.0 N</td></tr><tr><td>Spring stiffness—<i>K</i><sub>s</sub></td><td>240 N/mm</td></tr><tr><td>Initial position—<i>y</i><sub>0</sub></td><td>0.017 m</td></tr><tr><td>Initial velocity—<i>v</i><sub>0</sub></td><td>0.0 m/s</td></tr><tr><td>Gravity acceleration—<i>g</i></td><td>9.81 m/s<sup>2</sup></td></tr><tr><td>Coefficient of restitution—<i>ε</i><sub>N</sub></td><td>0.4</td></tr><tr><td>Integration time step—Δ<i>t</i></td><td>0.0001 s</td></tr></table>

(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 \left( \frac{\theta}{\beta} - \frac{1}{2\pi} \sin \frac{2\pi\theta}{\beta} \right) \quad (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:

$$\mathbf{q} = (y) \quad (58)$$

$$\mathbf{u} = (\dot{y}) \quad (59)$$

$$\mathbf{M} = (m) \quad (60)$$

$$\mathbf{h} = (-mg - K_s y) \quad (61)$$

$$g_N = y - R_b - s \quad (62)$$

$$\mathbf{W}_N = (1) \quad (63)$$

$$\tilde{\mathbf{w}}_N = (-\dot{y}) \quad (64)$$

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 by rebounds due to the impacts that take place. The global results are in line with those offered in the literature [64]. 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.
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 7 (a) Follower displacement; (b) follower velocity</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8 Animation sequence of the virtual simulation of the cam follower movement during the first instants after the follower reaches the up dead point</figcaption></figure>
# 6 Conclusions

A general and comprehensive methodology for modeling and analyzing contact-impact events in multi-body 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 contact-impact 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 re-bounds, 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.

**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.

# References

1. Wittenberg, J.: Dynamics of Systems of Rigid Bodies. Teubner, Stuttgart (1977)

2. Nikravesh, P.E.: Computer-Aided Analysis of Mechanical Systems. Prentice Hall, Englewood Cliffs (1988)

3. Haug, E.J.: Computer-Aided Kinematics and Dynamics of Mechanical Systems—Volume I: Basic Methods. Allyn & Bacon, Boston (1989)

4. Shabana, A.A.: Dynamics of Multibody Systems. Wiley, New York (1989)

5. Huston, R.L.: Multibody Dynamics. Butterworth-Heinemann, Boston (1990)

6. Rulka, W.: SIMPACK—a computer program for simulation of large motion multibody systems. In: Schiehlen, W. (ed.) Multibody Systems Handbook, pp. 265–284. Springer, Berlin (1990)

7. Ryan, R.R.: ADAMS—Multibody System Analysis Software. Multibody Systems Handbook. Springer, Berlin (1990)

8. Jiménez, J.M., Avello, A., García-Alonso, A., and Jalón, J.G.: COMPAMM—A Simple and Efficient Code for Kinematic and Dynamic Numerical Simulation of 3-D Multi-body System with Realistic Graphics. Multibody Systems Handbook. Springer, Berlin (1990)

9. Hurmuzlu, Y., Marghitu, D.B.: Rigid body collision of planar kinematic chain with multiple contact points. Int. J. Robot. Res. 13, 82–89 (1994)

10. Keller, J.B.: Impact with friction. Journal of Applied Mechanics 531(4) (1986)

11. Han, I., Gilmore, B.J.: Multi body impact motion with friction analysis, simulation, and validation. J. Mech. Des. 115, 412–422 (1993)

12. Pereira, M., Nikravesh, P.: Impact dynamics of multibody systems with frictional contact using joint coordinates and canonical equations of motion. Nonlinear Dyn. 9, 53–71 (1996)

13. Haug, E.J., Wu, S.C., Yang, S.M.: Dynamics of mechanical systems with coulomb friction, stiction, impact and constraint addition deletion—I theory. Mech. Mach. Theory 21, 401–406 (1986)

14. Flores, P., Leine, R., Glocker, C.: Modeling and analysis of rigid multibody systems with translational clearance joints based on the nonsmooth dynamics approach. Multibody Syst. Dyn. 23(2), 165–190 (2010)

15. Tian, Q., Liu, C., Machado, M., Flores, P.: A new model for dry and lubricated cylindrical joints with clearance in spatial flexible multibody systems. Nonlinear Dyn. 64(1–2), 25–47 (2011)

16. Erkaya, S., Uzmay, I.: Investigation on effect of joint clearance on dynamics of four-bar mechanism. Nonlinear Dyn. 58(1–2), 179–198 (2009)

17. Ishida, Y., Inagaki, M., Ejima, R., Hayashi, A.: Nonlinear resonances and self-excited oscillations of a rotor caused by radial clearance and collision. Nonlinear Dyn. 57(4), 593–605 (2009)

18. Flores, P., Lankarani, H.M.: Spatial rigid-multi-body systems with lubricated spherical clearance joints: modeling and simulation. Nonlinear Dyn. 60(1–2), 99–114 (2010)

19. Machado, M., Flores, P., Claro, J.C.P., Ambrósio, J., Silva, M., Completo, A., Lankarani, H.M.: Development of a planar multi-body model of the human knee joint. Nonlinear Dyn. 60(3), 459–478 (2010)

20. Brutti, C., Coglitore, C., Valentini, P.P.: Modeling 3D revolute joint with clearance and contact stiffness. Nonlinear Dyn. 66(4), 531–548 (2011)

21. Kim, S., Nikravesh, P.E., Gim, G.: A two-dimensional tire model on uneven roads for vehicle dynamic simulation. Veh. Syst. Dyn. 56(10), 913–930 (2008)

22. Pombo, J.C., Ambrósio, J.A.C.: Application of a wheel–rail contact model to railway dynamics in small radius curved tracks. Multibody Syst. Dyn. 19(1–2), 91–114 (2008)

23. Miller, A., Allen, P., Santos, V., Valero-Cuevas, F.: From robotic hands to human hands: a visualization and simulation
tion engine for grasping research. Ind. Robot **32**(1), 55–63
(2005)

24. Ciocarlie, M., Miller, A., Allen, P.: Grasp analysis using deformable fingers. In: Intelligent Robots and Systems, pp. 4122–4128 (2005)

25. Miller, A.T., Christensen, H.I.: Implementation of multi-rigid-body dynamics within a robotic grasping simulator. In: Proceedings of the 2003 IEEE International Conference on Robotics and Automation, Taipei, Taiwan, 14–19 September, pp. 2262–2268 (2003)

26. Flores, P., Machado, M., Silva, M.T., Martins, J.M.: On the continuous contact force models for soft materials in multibody dynamics. Multibody Syst. Dyn. **25**(3), 357–375 (2011)

27. Flores, P.: Modeling and simulation of wear in revolute clearance joints in multibody systems. Mech. Mach. Theory **44**(6), 1211–1222 (2009)

28. Mukras, S., Kim, N.H., Mauntler, N.A., Schmitz, T.L., Sawyer, W.G.: Analysis of planar multibody systems with revolute joint wear. Wear **268**(5–6), 643–652 (2010)

29. Lankarani, H.M., Nikravesh, P.E.: A contact force model with hysteresis damping for impact analysis of multibody systems. J. Mech. Des. **112**, 369–376 (1990)

30. Flores, P., Ambrósio, J.: On the contact detection for contact-impact analysis in multibody systems. Multibody Syst. Dyn. **24**(1), 103–122 (2010)

31. Flores, P.: A parametric study on the dynamic response of planar multibody systems with multiple clearance joints. Nonlinear Dyn. **61**(4), 633–653 (2010)

32. Pfeiffer, F., Glocker, C.: Multibody Dynamics with Unilateral Contacts. Wiley, New York (1996)

33. Glocker, C.: Set-Valued Force Laws: Dynamics of Non-Smooth Systems. Lecture Notes in Applied Mechanics, vol. 1. Springer, Berlin (2001)

34. Potra, F.A., Anitescu, M., Gavrea, B., Trinkle, J.: A linearly implicit trapezoidal method for integrating stiff multibody dynamics with contact, joints, and friction. Int. J. Numer. Methods Eng. **66**(7), 1079–1124 (2006)

35. Dimitrakopolous, E.G.: Analysis of a frictional oblique impact observed in skew bridges. Nonlinear Dyn. **60**(4), 575–595 (2010)

36. Glocker, C., Pfeiffer, F.: Dynamical systems with unilateral contacts. Nonlinear Dyn. **3**(4), 245–259 (1992)

37. Glocker, C., Pfeiffer, F.: Multiple impacts with friction in rigid multibody systems. Nonlinear Dyn. **7**(4), 471–497 (1995)

38. Moreau, J.J.: Numerical aspects of the sweeping process. Comput. Methods Appl. Mech. Eng. **177**(3–4), 329–349 (1999)

39. Pfeiffer, F.: Unilateral problems of dynamics. Arch. Appl. Mech. **69**(8), 503–527 (1999)

40. Bauchau, O.A., Rodriguez, J., Bottasso, C.J.: Modeling of unilateral contact conditions with application to aerospace systems involving backlash, freeplay and friction. Mech. Res. Commun. **28**(5), 571–599 (2001)

41. Stewart, D.E.: Finite-dimensional contact mechanics. Philos. Trans. R. Soc., Math. Phys. Eng. Sci. **359**(1789), 2467–2482 (2001)

42. Pfeiffer, F.: The idea of complementarity in multibody dynamics. Arch. Appl. Mech. **72**(11–12), 807–816 (2003)

43. Förg, M., Pfeiffer, F., Ulbrich, H.: Simulation of unilateral constrained systems with many bodies. Multibody Syst. Dyn. **14**(2), 137–154 (2005)

44. Glocker, C., Studer, C.: Formulation and preparation for numerical evaluation of linear complementarity systems in dynamics. Multibody Syst. Dyn. **13**(4), 447–463 (2005)

45. Slavic, J., Boltezar, M.: Non-linearity and non-smoothness in multi-body dynamics: application to woodpecker toy. Proceedings of the institution of mechanical engineers, part C. J. Mech. Eng. Sci. **220**(3), 285–296 (2006)

46. Leine, R.I., van de Wouw, N.: Stability properties of equilibrium sets of non-linear mechanical systems with dry friction and impact. Nonlinear Dyn. **51**(4), 551–583 (2008)

47. Studer, C., Leine, R.I., Glocker, C.: Step size adjustment and extrapolation for time-stepping schemes in non-smooth dynamics. Int. J. Numer. Methods Eng. **76**(11), 1747–1781 (2008)

48. Brogliato, B., ten Dam, A.A., Paoli, L., Genot, F., Abadie, M.: Numerical simulations of finite dimensional multibody nonsmooth mechanical systems. Appl. Mech. Rev. **55**(2), 107–150 (2002)

49. Leine, R.I., Nijmeijer, H.: Dynamics and Bifurcations of Non-Smooth Mechanical Systems. Lecture Notes in Applied and Computational Mechanics, vol. 18. Springer, Berlin (2004)

50. Anitescu, M., Potra, F.A.: Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems. Nonlinear Dyn. **14**(3), 231–247 (2007)

51. Brogliato, B.: Some perspectives on the analysis and control of complementarity systems. IEEE Trans. Autom. Control **48**(6), 918–935 (2003)

52. Flores, P.: Contact-impact analysis in multibody systems based on the nonsmooth dynamics approach. Post Doctoral Report, ETH-Zurich Switzerland (2009)

53. Signorini, A.: Sopra alcune questioni di elastostatica. Atti della Societa Italian per il Progresso della Scienza (1933)

54. Leine, R.I., Glocker, C.: A set-valued force law for spatial Coulomb-Contensou friction. Eur. J. Mech. A, Solids **22**(2), 193–216 (2003)

55. Greenwood, D.T.: Principles of Dynamics. Prentice Hall, Englewood Cliffs (1965)

56. Moreau, J.J.: Unilateral contact and dry friction in finite freedom dynamics. In: Moreau, J.J., Panagiotopoulos, P.D. (eds.) Non-Smooth Mechanics and Applications. CISM Courses and Lectures, vol. 302, pp. 1–82. Springer, Wien (1988)

57. Glocker, C.: On frictionless impact models in rigid-body systems. Philos. Trans. R. Soc., Math. Phys. Eng. Sci. **359**(1789), 2385–2404 (2001)

58. Pfeiffer, F., Foerg, M., Ulbrich, H.: Numerical aspects of non-smooth multibody dynamics. Comput. Methods Appl. Mech. Eng. **195**(50–51), 6891–6908 (2006)

59. Leine, R.I., van de Wouw, N.: Stability and Convergence of Mechanical Systems with Unilateral Constraints. Lecture Notes in Applied and Computational Mechanics, vol. 36. Springer, Berlin (2008)

60. Rockafellar, R.T.: Augmented lagrangians and applications of the proximal point algorithm in convex programming. Math. Oper. Res. **1**(2), 97–116 (1976)
61. Alart, P., Curnier, A.: A mixed formulation for frictional contact problems prone to newton like solution methods. Comput. Methods Appl. Mech. Eng. 92(3), 353–375 (1991)

62. Studer, C., Glocker, C.: Solving normal cone inclusion problems in contact mechanics by iterative methods. J. Syst. Des. Dyn. 1(3), 458–467 (2007)

63. Foerg, M., Geier, T., Neumann, L., Ulbrich, H.: r-Factor strategies for the augmented Lagrangian approach in multi-body contact mechanics. In: Proceedings of III European Conference on Computational Mechanics, Lisbon, Portugal (2006), 20 p.

64. Seabra, E., Flores, P., Silva, J.F.: Theoretical and experimental analysis of an industrial cutting file machine using multibody systems methodology. In: Proceedings of EC-COMAS Thematic Conference Multibody Dynamics 2007, Milan, 25–28 June, 2007 (2007), 12 p.

65. Seabra, E.A.R., Flores, P., Claro, J.C.P., Silva, J.C.L.: Kinematics and dynamics study of the cam-follower mechanism of the cutting file machine. In: Internationales Wissenschaftliches Kolloquium, Technische Universität Ilme-nau, September 23–26, 2002 (2002), 12 p.

66. Seabra, E.A.R., Flores, P., Silva, J.C.L.F.: Re-design of a cam-follower mechanism of an industrial cutting file machine. In: Third International Conference on Advanced Engineering Design, Prague, Czech Republic, June 1–4, 2003, (2003), 8 p.

67. Chen, F.Y.: Mechanics and Design of Cam Mechanisms. Pergamon, New York (1982)

68. Norton, R.L.: Cam Design and Manufacturing Handbook. Industrial Press, Inc., New York (2002)