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

Paulo Flores, Remco Leine, Christoph Glocker

► To cite this version:

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. Multibody System Dynamics, 2010, 23, pp.165 - 190. 10.1007/s11044-009-9178-y. hal-01435625

HAL Id: hal-01435625

https://hal.science/hal-01435625

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
# Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach

Paulo Flores · Remco Leine · Christoph Glocker

**Abstract** The main purpose of this paper is to present and discuss a methodology for a dynamic modeling and analysis of rigid multibody systems with translational clearance joints. The methodology is based on the non-smooth dynamics approach, in which the interaction of the elements that constitute a translational clearance joint is modeled with multiple frictional unilateral constraints. In the following, the most fundamental issues of the non-smooth dynamics theory are revised. 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 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 a linear complementarity problem, which is embedded in the Moreau time-stepping method. Finally, the classical slider-crank mechanism is considered as a demonstrative application example and numerical results are presented. The results obtained show that the existence of clearance joints in the modeling of multibody systems influences their dynamics response.

**Keywords** Non-smooth dynamics · Clearance joints · Contact-impact mechanics · Linear complementarity problem · Moreau time-stepping method · Rigid multibody systems
# 1 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 [1–8].

Indeed, a number of theoretical and experimental works devoted to the research on multi-body 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. [2]. 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 [9–11]. Moreover, each contact point may be in stick or in slip phase, which greatly enlarges the number of contact configurations. The conditions for switching from one case to another depend on the system’s dynamic response.

Farahanchi and Shaw [12] 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 [13] 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 [14] 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 set-valued 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.
## 2 Basic set-valued elements

### 2.1 The linear complementarity problem

A linear complementarity problem (LCP) is a set of linear equations that can be written as [15, 16]

$$ y = Ax + b \quad (1) $$

subjected to the inequality complementarity conditions

$$ y \ge 0, \quad x \ge 0, \quad y^T x = 0 \quad (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 \in \mathbb{R}^n$ and $y \in \mathbb{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 \le y \perp x \ge 0 \quad (3) $$

where $y \perp x$ denotes $y^T x = 0$. An LCP can have a unique solution, multiple solutions or no solution at all [17, 18]. 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 [19–21]. 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. [16].

### 2.2 The unilateral primitive

One of the most important multifunctions (or set-valued maps) related to complementarity is the unilateral primitive, denoted by $\text{Upr}$. The unilateral primitive is a maximal monotone set-valued map on $\mathbb{R}_0^+$ defined as [22, 23]

$$ \text{Upr}(x) := \begin{cases} \{0\} & x > 0 \\ (-\infty, 0] & x = 0 \\ \emptyset & x < 0 \end{cases} \quad (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 \in \text{Upr}(x) \Leftrightarrow y \ge 0, x \ge 0, xy = 0 \quad (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 (5).
<figure><img src="image_2.png" /><figcaption>Fig. 1 (a) The map $x \to \text{Upr}(x)$; (b) The map $x \to \text{Sgn}(x)$; (c) The decomposition $\text{Sgn}(x)$ into $\text{Upr}(x)$</figcaption></figure>

## 2.3 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 [22, 23],

$$ \text{Sgn}(x) := \begin{cases} \{+1\} & x > 0 \\ [-1, +1] & x = 0 \\ \{-1\} & x < 0 \end{cases} \qquad (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 \in \text{Sgn}(x) \Leftrightarrow \exists x_R, x_L \text{ s.t. } \begin{cases} -y \in +\text{Upr}(x_R) + 1 \\ -y \in -\text{Upr}(x_L) - 1 \\ x = x_R - x_L \end{cases} \qquad (7) $$

Using (5), the (7) can be rewritten in terms of complementarities,

$$ -y \in \text{Sgn}(x) \Leftrightarrow \exists x_R, x_L \text{ s.t. } \begin{cases} 1 + y \ge 0, & x_R \ge 0, & (1+y)x_R = 0 \\ 1 - y \ge 0, & x_L \ge 0, & (1-y)x_L = 0 \\ x = x_R - x_L \end{cases} \qquad (8) $$

This representation has to be used when a problem involving Sgn-multifunctions is formulated as an LCP in its standard form [24]. 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 [22].

# 3 Set-valued force laws for frictional unilateral contacts

## 3.1 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 [25]. Figure 2 shows two convex rigid bodies apart
<figure><img src="image_3.png" /><figcaption>Fig. 2 (a) Relative normal gap; (b) Normal and tangential contact forces</figcaption></figure>

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 $\lambda_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 \ge 0$, only two situations can occur, namely,

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

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

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 \tag{12}$$

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. 3(a) and shows a set-valued graph or a corner of admissible combinations between $g_N$ and $\lambda_N$ [26]. When two rigid bodies are contacting, the Signorini's condition given by (12) 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, $\epsilon_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
<figure><img src="image_4.png" /><figcaption>**Fig. 3 (a)** Signorini's normal contact law; **(b)** Coulomb's friction law</figcaption></figure>

decrease with Newton's impact law can be found in Leine and van de Wouw [23]. These issues have recently been discussed by other authors [27–29].

## 3.2 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 [22, 30]. 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 [31, 32]. 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,

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

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

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

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

$$ -\lambda_T \in \mu\lambda_N\operatorname{Sgn}(\gamma_T) \qquad (17) $$

Figure 3(b) shows the Coulomb's friction law as a set-valued force law [22].
# 4 Dynamics of non-smooth rigid multibody systems

## 4.1 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 [33]

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

$$ \dot{\mathbf{q}} = \mathbf{u} \quad \forall t \qquad (19) $$

where $M = M(\mathbf{q}, t) \in \mathbb{R}^{f \times f}$ is the positive definite and symmetric mass matrix, $\mathbf{h} = \mathbf{h}(\mathbf{q}, \mathbf{u}, t) \in \mathbb{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 $\mathbf{h}$, $\mathbf{q} = \mathbf{q}(t) \in \mathbb{R}^f$ is the $f$-dimensional vector of generalized coordinates, $\mathbf{u} = \mathbf{u}(t) \in \mathbb{R}^f$ addresses the system generalized velocities and $\dot{\mathbf{u}} = \dot{\mathbf{u}}(t) \in \mathbb{R}^f$ is the vector that contains the system accelerations.

It is clear that (18) 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 [34]. Thus, adding the contact forces to (18), 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 [23, 32]

$$ M\dot{\mathbf{u}} - \mathbf{h} - W_N\lambda_N - W_T\lambda_T = \mathbf{0} \quad \text{a.e.} \qquad (20) $$

$$ \dot{\mathbf{q}} = \mathbf{u} \quad \forall t \qquad (21) $$

where $W_N = W_N(\mathbf{q}, t) \in \mathbb{R}^{f \times n}$ and $W_T = W_T(\mathbf{q}, t) \in \mathbb{R}^{f \times n}$ gather the generalized normal and tangential force directions $w_{Ni}$ and $w_{Ti}$, respectively. The normal and tangential 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 distances $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 (20) have the same meaning as described above. It is important to note that (20) requires the existence of the velocities $\mathbf{u}$ as well as the existence of accelerations $\dot{\mathbf{u}}$. Motion without impulses implies that $\lambda_N(t)$ is (locally) bounded and time-continuous. The velocities $\mathbf{u}(t)$ therefore exist on non-impulsive time-intervals. The friction force $\lambda_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 $\dot{\mathbf{u}}$ is not defined when $\lambda_T(t)$ is discontinuous. The set of time instances for which $\lambda_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(\mathbf{u}^{+} - \mathbf{u}^{-}) - W_N \mathbf{A}_N - W_T \mathbf{A}_T = \mathbf{0} \quad \text{a.e.} \qquad (22) $$

$$ \mathbf{u}^{+}(t) = \lim_{\Delta t \downarrow 0} \frac{\mathbf{q}(t + \Delta t) - \mathbf{q}(t)}{\Delta t}, \quad \mathbf{u}^{-}(t) = \lim_{\Delta t \uparrow 0} \frac{\mathbf{q}(t + \Delta t) - \mathbf{q}(t)}{\Delta t} \qquad (23) $$

which relates the velocity jump to the impulsive forces $\mathbf{A}_N$ and $\mathbf{A}_T$ in normal and tangential direction respectively. We assume that the velocities $\mathbf{u}(t)$ are of locally bounded variation
(without singular part) and denote $\mathbf{u}^{-}(t)$ and $\mathbf{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 (22).

Following Moreau [35] we will cast the non-impulsive dynamics (20) and the impulsive dynamics (22) in a unified description, by using an equality of measures. This constitutes the general framework for non-smooth rigid multibody dynamics [31, 36].

Multiplying the equation of motion (20) with the Lebesgue measure $dt$ and the impact equation (22) with the atomic measure $d\eta$, being the sum of the Dirac point measures at the impact times, yields

$$ \mathbf{M}\dot{\mathbf{u}}dt - \mathbf{h}dt - \mathbf{W}_N\lambda_N dt - \mathbf{W}_T\lambda_T dt = \mathbf{0} \quad (24) $$

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

Addition of (24) and (25) results in

$$ \mathbf{M}[\dot{\mathbf{u}}dt + (\mathbf{u}^{+} - \mathbf{u}^{-})d\eta] - \mathbf{h}dt - \mathbf{W}_N(\lambda_N dt + \mathbf{\Lambda}_N d\eta) - \mathbf{W}_T(\lambda_T dt + \mathbf{\Lambda}_T d\eta) = \mathbf{0} \quad (26) $$

or more briefly,

$$ \mathbf{M}\mathbf{d}\mathbf{u} - \mathbf{h}dt - \mathbf{W}_N d\mathbf{P}_N - \mathbf{W}_T d\mathbf{P}_T = \mathbf{0} \quad (27) $$

The differential measure for the velocities $\mathbf{d}\mathbf{u} = \dot{\mathbf{u}}dt + (\mathbf{u}^{+} - \mathbf{u}^{-})d\eta$ consists of the Lebesgue measurable part $\dot{\mathbf{u}}dt$, which accounts for absolutely continuous motion, and the atomic parts which accounts for impulsive motion. Hence, for impact free motion it holds that $\mathbf{d}\mathbf{u} = \dot{\mathbf{u}}dt$. Similarly, the measure for the so-called percussions corresponds to a Lagrangian multiplier which gathers both finite contact forces $\lambda$ and impulsive contact forces $\mathbf{\Lambda}$, that is, $d\mathbf{P} = \lambda dt + \mathbf{\Lambda} d\eta$ [37].

In the case of non-impulsive motion, all measures $d\eta$ vanish and a formal division by $dt$ yields the classical Newton-Euler equations of motion given by (20). 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 [24].

## 4.2 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 inclu- sions 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 [35, 38] and Glocker [22].

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}(\mathbf{q}, t) \ge 0, \quad i = 1, \dots, n \quad (28) $$

where the quantities $g_{Ni}$ are the normal gap functions of the frictional contacts. They are for- mulated 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 [39] and Glocker [22].

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

$$H(t) = \{i | g_{Ni}(\mathbf{q}, t) = 0\} \qquad (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 $\mathbf{q}$ and $\mathbf{u}$, let us first introduce the normal and tangential relative velocities at the contacts as [40]

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

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

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 rheonomic terms [22].

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 $\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]$. 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 $d\mathbf{P}_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 [40]

$$\xi_{Ni} := \gamma_{Ni}^{+} + \epsilon_{Ni} \gamma_{Ni}^{-} \qquad (32)$$

$$\xi_{Ti} := \gamma_{Ti}^{+} + \epsilon_{Ti} \gamma_{Ti}^{-} \qquad (33)$$

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

Thus, the normal and tangential impact laws can be stated as two inclusions,

$$-d\mathbf{P}_{Ni} \in \mathrm{Upr}(\xi_{Ni}) \qquad (34)$$

$$-d\mathbf{P}_{Ti} \in \mu_i d\mathbf{P}_{Ni} \mathrm{Sgn}(\xi_{Ti}) \qquad (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)–(35). This problem can be solved by using the Moreau time-stepping method, which is presented and discussed in the next section.
# 5 Moreau time-stepping method

## 5.1 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 [37–42]. These methods are widely used due to their simplicity to implement and their robustness. The time-stepping schemes are based on a time-discretization 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 [35]. The equality of measures (27) together with the set-valued force laws (34) 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 (27). The first step of the Moreau approach consists of the time-discretization of the measure differential equation. Thus, integrating (27) over a small finite time interval $\Delta t$, of which initial and end points are denoted by the indices A and E, yields the following terms:

$$ \int_{\Delta t} \mathbf{M} \mathrm{d} \mathbf{u} \approx \mathbf{M}_M \Delta \mathbf{u} = \mathbf{M}_M (\mathbf{u}_E - \mathbf{u}_A), \quad \mathbf{M}_M = \mathbf{M}(\mathbf{q}_M, t_M) \quad (36) $$

$$ \int_{\Delta t} \mathbf{h} \mathrm{d} t = \Delta \mathbf{h} \approx \mathbf{h}_M \Delta t, \quad \mathbf{h}_M = \mathbf{h}(\mathbf{q}_M, \mathbf{u}_A, t_M) \quad (37) $$

$$ \int_{\Delta t} \mathbf{W}_N \mathrm{d} \mathbf{P}_N = \mathbf{W}_{NM} \mathbf{P}_N, \quad \mathbf{W}_{NM} = \mathbf{W}_N(\mathbf{q}_M, t_M) \quad (38) $$

$$ \int_{\Delta t} \mathbf{W}_T \mathrm{d} \mathbf{P}_T = \mathbf{W}_{TM} \mathbf{P}_T, \quad \mathbf{W}_{TM} = \mathbf{W}_T(\mathbf{q}_M, t_M) \quad (39) $$

where $t_M$ is the midpoint time instant of the compact time interval $[t_A, t_E]$ and $\mathbf{q}_M = \mathbf{q}_A + \frac{1}{2}\mathbf{u}_A \Delta t$ is the midpoint system’s position state. It is clear that the midpoint time instant can be evaluated as

$$ t_M = t_A + \frac{1}{2} \Delta t \quad (40) $$

Finally, after the above discretization, the equations of motion expressed at the velocity level can be written as [24]

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

together with the set-valued contact/impact laws,

$$ -\mathbf{P}_N \in \mathrm{Upr}(\xi_N) \Leftrightarrow -\mathbf{P}_N \in N_{C_N}(\xi_N) \quad (42) $$

$$ -\mathbf{P}_T \in \mu \mathbf{P}_N \mathrm{Sgn}(\xi_T) \Leftrightarrow -\mathbf{P}_T \in N_{C_T}(\mathbf{P}_N)(\xi_T) \quad (43) $$

This set of algebraic inclusions can be solved with a linear complementarity problem (LCP) formulation or by an augmented Lagrangian approach (ALA) [21]. The velocity $\mathbf{u}_E$,
at the end of time-step $t_E = t_A + \Delta t$, is subsequently calculated by using (41). Finally, the positions at the end of the time step are calculated by

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

Note that (42) applies only to active set-valued force laws, $i \in 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 (43). Considering unilateral contacts, Moreau midpoint algorithm calculates the contact distances $g_{Ni}$ of all unilateral contacts at the midpoint $\mathbf{q}_M$ in order to evaluate whether these are active ($g_{Ni} \le 0$) or not ($g_{Ni} > 0$). Only active unilateral contacts can be described by inclusion (42). 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:

<figure><img src="image_5.png" /><figcaption>Fig. 4 Flowchart of the computational procedure for the solution of the equations of motion of constrained rigid multibody systems with frictional unilateral constraints</figcaption></figure>
(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 $\Delta t$, together with the given initial positions $\mathbf{q}_A$ and velocities $\mathbf{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 $\mathbf{q}_M$, assemble the midpoint mass matrix $\mathbf{M}_M$ and the gyroscopic and external forces vector $\mathbf{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 $\mathbf{u}_E$, by using (41); 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 $\mathbf{P}_N$ and $\mathbf{P}_T$ required to compute $\mathbf{u}_E$ for the contact-impact case;

(iv) Compute the position's state at the end time $\mathbf{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.

## 5.2 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 [24]. In order to set up the LCP, let us first introduce the following matrix notation:

$$ \mathbf{W}_{NM} := \mathrm{mat}(\mathbf{w}_{Ni}(\mathbf{q}_M, t_M)) \in \mathbb{R}^{f,i}, \quad i \in H \quad (45) $$

$$ \mathbf{W}_{TM} := \mathrm{mat}(\mathbf{w}_{Ti}(\mathbf{q}_M, t_M)) \in \mathbb{R}^{f,i}, \quad i \in H \quad (46) $$

$$ \tilde{\mathbf{w}}_{NM} := \mathrm{col}(\tilde{\mathbf{w}}_{Ni}(\mathbf{q}_M, t_M)) \in \mathbb{R}^i, \quad i \in H \quad (47) $$

$$ \tilde{\mathbf{w}}_{TM} := \mathrm{col}(\tilde{\mathbf{w}}_{Ti}(\mathbf{q}_M, t_M)) \in \mathbb{R}^i, \quad i \in H \quad (48) $$

$$ \mathbf{P}_N := \mathrm{col}(P_{Ni}) \in \mathbb{R}^i, \quad i \in H \quad (49) $$

$$ \mathbf{P}_T := \mathrm{col}(P_{Ti}) \in \mathbb{R}^i, \quad i \in H \quad (50) $$

$$ \boldsymbol{\gamma}_{NE} := \mathrm{col}(\gamma_{NEi}) \in \mathbb{R}^i, \quad i \in H \quad (51) $$

$$ \boldsymbol{\gamma}_{TE} := \mathrm{col}(\gamma_{TEi}) \in \mathbb{R}^i, \quad i \in H \quad (52) $$

$$ \boldsymbol{\gamma}_{NA} := \mathrm{col}(\gamma_{NAi}) \in \mathbb{R}^i, \quad i \in H \quad (53) $$

$$ \boldsymbol{\gamma}_{TA} := \mathrm{col}(\gamma_{TAi}) \in \mathbb{R}^i, \quad i \in H \quad (54) $$

$$ \boldsymbol{\xi}_N := \mathrm{col}(\boldsymbol{\xi}_{Ni}) \in \mathbb{R}^i, \quad i \in H \quad (55) $$

$$ \boldsymbol{\xi}_T := \mathrm{col}(\boldsymbol{\xi}_{Ti}) \in \mathbb{R}^i, \quad i \in H \quad (56) $$

$$ \boldsymbol{\varepsilon}_N := \mathrm{diag}(\varepsilon_{Ni}) \in \mathbb{R}^i, \quad i \in H \quad (57) $$

$$ \boldsymbol{\varepsilon}_T := \mathrm{diag}(\varepsilon_{Ti}) \in \mathbb{R}^i, \quad i \in H \quad (58) $$

$$ \boldsymbol{\mu} := \mathrm{diag}(\mu_i) \in \mathbb{R}^i, \quad i \in H \quad (59) $$
Thus, the contact-impact problem of non-smooth systems can be summarized by the following mathematical relations:

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

$$ \boldsymbol{\gamma}_{NE} = \mathbf{W}_{NM}^{\mathrm{T}} \mathbf{u}_E + \tilde{\mathbf{w}}_{NM} \quad (61) $$

$$ \boldsymbol{\gamma}_{TE} = \mathbf{W}_{TM}^{\mathrm{T}} \mathbf{u}_E + \tilde{\mathbf{w}}_{TM} \quad (62) $$

$$ \boldsymbol{\gamma}_{NA} = \mathbf{W}_{NM}^{\mathrm{T}} \mathbf{u}_A + \tilde{\mathbf{w}}_{NM} \quad (63) $$

$$ \boldsymbol{\gamma}_{TA} = \mathbf{W}_{TM}^{\mathrm{T}} \mathbf{u}_A + \tilde{\mathbf{w}}_{TM} \quad (64) $$

$$ \boldsymbol{\xi}_N = \boldsymbol{\gamma}_{NE} + \boldsymbol{\varepsilon}_N \boldsymbol{\gamma}_{NA} \quad (65) $$

$$ \boldsymbol{\xi}_T = \boldsymbol{\gamma}_{TE} + \boldsymbol{\varepsilon}_T \boldsymbol{\gamma}_{TA} \quad (66) $$

$$ -\mathbf{P}_N \in \mathrm{Upr}(\boldsymbol{\xi}_N) \quad (67) $$

$$ -\mathbf{P}_T \in \mu \mathbf{P}_N \mathrm{Sgn}(\boldsymbol{\xi}_T) \quad (68) $$

The values of $\boldsymbol{\gamma}_{NA}$ and $\boldsymbol{\gamma}_{TA}$ can be evaluated by using (63) and (64), respectively, since the velocities $\mathbf{u}_A$ are known at the left endpoint of the time interval. Introducing now (61) and (62) into (65) and (66) yields

$$ \boldsymbol{\xi}_N = \mathbf{W}_{NM}^{\mathrm{T}} \mathbf{u}_E + (\tilde{\mathbf{w}}_{NM} + \boldsymbol{\varepsilon}_N \boldsymbol{\gamma}_{NA}) \quad (69) $$

$$ \boldsymbol{\xi}_T = \mathbf{W}_{TM}^{\mathrm{T}} \mathbf{u}_E + (\tilde{\mathbf{w}}_{TM} + \boldsymbol{\varepsilon}_T \boldsymbol{\gamma}_{TA}) \quad (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

$$ -\mathbf{P}_N \in \mathrm{Upr}(\boldsymbol{\xi}_N) \quad \Leftrightarrow \quad \mathbf{P}_N \ge \mathbf{0}, \quad \boldsymbol{\xi}_N \ge \mathbf{0}, \quad \mathbf{P}_N^{\mathrm{T}} \boldsymbol{\xi}_N = \mathbf{0} \quad (71) $$

In turn, the relay function (68) has to be decomposed into two Upr's to achieve the desired complementarity conditions. Thus, (68) yields

$$ -\mathbf{P}_T \in \mu \mathbf{P}_N \mathrm{Sgn}(\boldsymbol{\xi}_T) \quad \Leftrightarrow \quad \exists \boldsymbol{\xi}_R, \boldsymbol{\xi}_L \\ \begin{array}{l} \text{s.t.} \quad \left\{ \begin{array}{ll} \boldsymbol{\mu} \mathbf{P}_N + \mathbf{P}_T \ge \mathbf{0}, & \boldsymbol{\xi}_R \ge \mathbf{0}, & (\boldsymbol{\mu} \mathbf{P}_N + \mathbf{P}_T)^{\mathrm{T}} \boldsymbol{\xi}_R = \mathbf{0} \\ \boldsymbol{\mu} \mathbf{P}_N - \mathbf{P}_T \ge \mathbf{0}, & \boldsymbol{\xi}_L \ge \mathbf{0}, & (\boldsymbol{\mu} \mathbf{P}_N - \mathbf{P}_T)^{\mathrm{T}} \boldsymbol{\xi}_L = \mathbf{0} \\ \boldsymbol{\xi}_T = \boldsymbol{\xi}_R - \boldsymbol{\xi}_L \end{array} \right. \end{array} \quad (72) $$

in which the step height is $[-\boldsymbol{\mu}\mathbf{P}_N, +\boldsymbol{\mu}\mathbf{P}_N]$. In addition, to abbreviate the complementarity conditions of (72) the impulsive friction saturations $\mathbf{P}_R$ and $\mathbf{P}_L$ are defined as [22]

$$ \mathbf{P}_R := \boldsymbol{\mu}\mathbf{P}_N + \mathbf{P}_T, \quad \mathbf{P}_R \in \mathbb{R}^i \quad (73) $$

$$ \mathbf{P}_L := \boldsymbol{\mu}\mathbf{P}_N - \mathbf{P}_T, \quad \mathbf{P}_L \in \mathbb{R}^i \quad (74) $$

together with

$$ \boldsymbol{\xi}_T = \boldsymbol{\xi}_R - \boldsymbol{\xi}_L, \quad \boldsymbol{\xi}_R, \boldsymbol{\xi}_L \in \mathbb{R}^i \quad (75) $$

Then, the whole set of complementarity conditions of (72) can be rewritten as

$$ \mathbf{0} \le \begin{pmatrix} \boldsymbol{\xi}_N \\ \boldsymbol{\xi}_R \\ \mathbf{P}_L \end{pmatrix} \perp \begin{pmatrix} \mathbf{P}_N \\ \mathbf{P}_R \\ \boldsymbol{\xi}_L \end{pmatrix} \ge \mathbf{0} \quad (76) $$
The reason for the special arrangement of $\mathbf{P}_L$ and $\xi_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 [22]. Since the variables $\xi_T$, $\mathbf{P}_T$ and $\mathbf{u}_E$ are not included in (28), they have to be eliminated. Thus, combining (60) and (73), yields

$$ \mathbf{M}_M(\mathbf{u}_E - \mathbf{u}_A) - \mathbf{h}_M \Delta t - (\mathbf{W}_{NM} - \mathbf{W}_{TM}\boldsymbol{\mu})\mathbf{P}_N - \mathbf{W}_{TM}\mathbf{P}_R = \mathbf{0} \quad (77) $$

Substituting now (75) into (70) results in

$$ \xi_R = \mathbf{W}_{TM}^\mathrm{T}\mathbf{u}_E + (\tilde{\mathbf{w}}_{TM} + \boldsymbol{\varepsilon}_T \boldsymbol{\gamma}_{TA}) + \xi_L \quad (78) $$

The elimination of variable $\mathbf{P}_T$ can be done through the combination of (73) and (74), which can be written as

$$ \mathbf{P}_L = 2\boldsymbol{\mu}\mathbf{P}_N - \mathbf{P}_R \quad (79) $$

Since the inversion of mass matrix $\mathbf{M}$ is always possible, (77) can be solved for $\mathbf{u}_E$:

$$ \mathbf{u}_E = \mathbf{u}_A + \mathbf{M}_M^{-1}\mathbf{h}_M \Delta t + \mathbf{M}_M^{-1}(\mathbf{W}_{NM} - \mathbf{W}_{TM}\boldsymbol{\mu})\mathbf{P}_N + \mathbf{M}_M^{-1}\mathbf{W}_{TM}\mathbf{P}_R \quad (80) $$

Now, (63) and (64) are used to express $\mathbf{W}_{NM}^\mathrm{T}\mathbf{u}_A$ and $\mathbf{W}_{TM}^\mathrm{T}\mathbf{u}_A$ in terms of $\boldsymbol{\gamma}_{NA}$ and $\boldsymbol{\gamma}_{TA}$,

$$ \mathbf{W}_{NM}^\mathrm{T}\mathbf{u}_A = \boldsymbol{\gamma}_{NA} - \tilde{\mathbf{w}}_{NM} \quad (81) $$

$$ \mathbf{W}_{TM}^\mathrm{T}\mathbf{u}_A = \boldsymbol{\gamma}_{TA} - \tilde{\mathbf{w}}_{TM} \quad (82) $$

Introducing (80)–(82) into (69) and (78), yields

$$ \xi_N = \mathbf{W}_{NM}^\mathrm{T}\mathbf{M}_M^{-1}\mathbf{h}_M \Delta t + \mathbf{W}_{NM}^\mathrm{T}\mathbf{M}_M^{-1}(\mathbf{W}_{NM} - \mathbf{W}_{TM}\boldsymbol{\mu})\mathbf{P}_N + \mathbf{W}_{NM}^\mathrm{T}\mathbf{M}_M^{-1}\mathbf{W}_{TM}\mathbf{P}_R \\ + (\mathbf{I} + \boldsymbol{\varepsilon}_N)\boldsymbol{\gamma}_{NA} \quad (83) $$

$$ \xi_R = \mathbf{W}_{TM}^\mathrm{T}\mathbf{M}_M^{-1}\mathbf{h}_M \Delta t + \mathbf{W}_{TM}^\mathrm{T}\mathbf{M}_M^{-1}(\mathbf{W}_{NM} - \mathbf{W}_{TM}\boldsymbol{\mu})\mathbf{P}_N + \mathbf{W}_{TM}^\mathrm{T}\mathbf{M}_M^{-1}\mathbf{W}_{TM}\mathbf{P}_R \\ + (\mathbf{I} + \boldsymbol{\varepsilon}_T)\boldsymbol{\gamma}_{TA} + \xi_L \quad (84) $$

Thus, (83), (84) and (79) can be written in a matrix form as

$$ \begin{pmatrix} \xi_N \\ \xi_R \\ \mathbf{P}_L \end{pmatrix} = \begin{pmatrix} \mathbf{W}_{NM}^\mathrm{T}\mathbf{M}_M^{-1}(\mathbf{W}_{NM} - \mathbf{W}_{TM}\boldsymbol{\mu}) & \mathbf{W}_{NM}^\mathrm{T}\mathbf{M}_M^{-1}\mathbf{W}_{TM} & \mathbf{0} \\ \mathbf{W}_{TM}^\mathrm{T}\mathbf{M}_M^{-1}(\mathbf{W}_{NM} - \mathbf{W}_{TM}\boldsymbol{\mu}) & \mathbf{W}_{TM}^\mathrm{T}\mathbf{M}_M^{-1}\mathbf{W}_{TM} & -\mathbf{I} \\ 2\boldsymbol{\mu} & -\mathbf{I} & \mathbf{0} \end{pmatrix} \begin{pmatrix} \mathbf{P}_N \\ \mathbf{P}_R \\ \boldsymbol{\xi}_L \end{pmatrix} + \begin{pmatrix} \mathbf{W}_{NM}^\mathrm{T}\mathbf{M}_M^{-1}\mathbf{h}_M \Delta t + (\mathbf{I} + \boldsymbol{\varepsilon}_N)\boldsymbol{\gamma}_{NA} \\ \mathbf{W}_{TM}^\mathrm{T}\mathbf{M}_M^{-1}\mathbf{h}_M \Delta t + (\mathbf{I} + \boldsymbol{\varepsilon}_T)\boldsymbol{\gamma}_{TA} \\ \mathbf{0} \end{pmatrix} \quad (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 $\mathbf{u}_E$ and positions $\mathbf{q}_E$ for the subsequent time steps are obtained from (80) and (44), respectively.
5.3 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 appropri-
ate 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$, $\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}$;

<figure><img src="image_6.png" /><figcaption>Fig. 5 Flowchart of the Moreau
time-stepping algorithm with an
LCP formulation</figcaption></figure>
(iii) compute the midpoint 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))$$

$$\boldsymbol{\gamma}_{NA} = \text{col}(\gamma_{NAi})$$

$$\boldsymbol{\gamma}_{TA} = \text{col}(\gamma_{TAi})$$

$$\boldsymbol{\varepsilon}_N = \text{diag}(\varepsilon_{Ni})$$

$$\boldsymbol{\varepsilon}_T = \text{diag}(\varepsilon_{Ti})$$

$$\boldsymbol{\mu} = \text{diag}(\mu_i)$$

(v) set up the LCP in the standard form $\mathbf{y} = \mathbf{A}\mathbf{x} + \mathbf{b}$:

$$\mathbf{A} = \begin{pmatrix} \mathbf{W}_{NM}^{\mathrm{T}} \mathbf{M}_{M}^{-1} (\mathbf{W}_{NM} - \mathbf{W}_{TM} \boldsymbol{\mu}) & \mathbf{W}_{NM}^{\mathrm{T}} \mathbf{M}_{M}^{-1} \mathbf{W}_{TM} & \mathbf{0} \\ \mathbf{W}_{TM}^{\mathrm{T}} \mathbf{M}_{M}^{-1} (\mathbf{W}_{NM} - \mathbf{W}_{TM} \boldsymbol{\mu}) & \mathbf{W}_{TM}^{\mathrm{T}} \mathbf{M}_{M}^{-1} \mathbf{W}_{TM} & \mathbf{0} \\ 2\boldsymbol{\mu} & -\mathbf{I} & \mathbf{0} \end{pmatrix}$$

$$\mathbf{b} = \begin{pmatrix} \mathbf{W}_{NM}^{\mathrm{T}} \mathbf{M}_{M}^{-1} \mathbf{h}_{M} \Delta t + (\mathbf{I} + \boldsymbol{\varepsilon}_{N}) \boldsymbol{\gamma}_{NA} \\ \mathbf{W}_{TM}^{\mathrm{T}} \mathbf{M}_{M}^{-1} \mathbf{h}_{M} \Delta t + (\mathbf{I} + \boldsymbol{\varepsilon}_{T}) \boldsymbol{\gamma}_{TA} \\ \mathbf{0} \end{pmatrix}$$

(vi) solve the LCP using an appropriate algorithm:

$$(\mathbf{x}, \mathbf{y}) = \text{LCP}(\mathbf{A}, \mathbf{b})$$

(vii) split the LCP solution according to

$$\mathbf{P}_N = \text{col}(\mathbf{x}_i), \quad i = 1, \dots, n_i$$

$$\mathbf{P}_L = \text{col}(\mathbf{x}_i), \quad i = n_i + 1, \dots, 2n_i$$

$$\mathbf{P}_R = \text{col}(\mathbf{y}_i), \quad i = 2n_i + 1, \dots, 3n_i$$
(viii) evaluate the velocity at the end of the integration time step:

$$ \mathbf{u}_E = \mathbf{u}_A + \mathbf{M}_M^{-1} \mathbf{h}_M \Delta t + \mathbf{M}_M^{-1} (\mathbf{W}_{NM} - \mathbf{W}_{TM} \boldsymbol{\mu}) \mathbf{P}_N + \mathbf{M}_M^{-1} \mathbf{W}_{TM} \mathbf{P}_R $$

(ix) 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 $$

(x) increment time step:

$$ t_A = t_A + \Delta t $$

(xi) update the 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.

# 6 Demonstrative application to a slider-crank mechanism

## 6.1 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><img src="image_7.png" /><figcaption>Fig. 6 Slider-crank mechanism with a translational clearance joint</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 8** Different scenarios for the slider and guide interaction: (a) no contact; (b) one corner in contact with the guide; (c) two adjacent corners in contact with guide; (d) two opposite corners in contact with guide</figcaption></figure>

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.

## 6.2 Lagrange's equations

In order for the translational clearance joint to be simulated in the multibody system environment, it is 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 [43]

$$ \frac{d}{dt}\left(\frac{\partial L}{\partial \dot{q}_i}\right) - \frac{\partial L}{\partial q_i} = 0, \quad i = 1, \dots, f \qquad (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

$$ \mathbf{q} = \begin{pmatrix} \theta_1 \\ \theta_2 \\ \theta_3 \end{pmatrix} \qquad (87) $$

$$ \mathbf{u} = \begin{pmatrix} \omega_1 \\ \omega_2 \\ \omega_3 \end{pmatrix}, \quad \text{with } \dot{\mathbf{q}} = \mathbf{u} \text{ a.e.} \qquad (88) $$

Thus, applying the Lagrange's equation to slider-crank mechanism yields [44]

$$ \begin{pmatrix} M_{11} & M_{12} & M_{13} \\ M_{21} & M_{22} & M_{23} \\ M_{31} & M_{32} & M_{33} \end{pmatrix} \begin{pmatrix} \ddot{\theta}_1 \\ \ddot{\theta}_2 \\ \ddot{\theta}_3 \end{pmatrix} = \begin{pmatrix} h_1 \\ h_2 \\ h_3 \end{pmatrix} \qquad (89) $$

in which

$$ M_{11} = J_1 + \left(\frac{1}{4}m_1 + m_2 + m_3\right)l_1^2 \qquad (90) $$

$$ M_{12} = M_{21} = \left(\frac{1}{2}m_2 + m_3\right)l_1 l_2 \cos(\theta_2 - \theta_1) \qquad (91) $$

$$ M_{13} = M_{31} = M_{23} = M_{32} = 0 \qquad (92) $$

$$ M_{22} = J_2 + \left(\frac{1}{4}m_2 + m_3\right)l_2^2 \qquad (93) $$

$$ M_{33} = J_3 \qquad (94) $$

$$ h_1 = \left(\frac{1}{2}m_2 + m_3\right)l_1 l_2 \sin(\theta_2 - \theta_1)\dot{\theta}_2^2 - \left(\frac{1}{2}m_1 + m_2 + m_3\right)gl_1 \cos\theta_1 \qquad (95) $$

$$ h_2 = -\left(\frac{1}{2}m_2 + m_3\right)l_1 l_2 \sin(\theta_2 - \theta_1)\dot{\theta}_1^2 - \left(\frac{1}{2}m_1 + m_2 + m_3\right)gl_2 \cos\theta_2 \qquad (96) $$

$$ h_3 = 0 \qquad (97) $$
6.3 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 [44]

$$g_{N1} = \frac{d}{2} - l_1 \sin \theta_1 - l_2 \sin \theta_2 + a \sin \theta_3 - b \cos \theta_3 \quad (98)$$

$$g_{T1} = l_1 \cos \theta_1 + l_2 \cos \theta_2 - a \cos \theta_3 - b \sin \theta_3 \quad (99)$$

$$g_{N2} = \frac{d}{2} - l_1 \sin \theta_1 - l_2 \sin \theta_2 - a \sin \theta_3 - b \cos \theta_3 \quad (100)$$

$$g_{T2} = l_1 \cos \theta_1 + l_2 \cos \theta_2 + a \cos \theta_3 - b \sin \theta_3 \quad (101)$$

$$g_{N3} = \frac{d}{2} + l_1 \sin \theta_1 + l_2 \sin \theta_2 - a \sin \theta_3 - b \cos \theta_3 \quad (102)$$

$$g_{T3} = l_1 \cos \theta_1 + l_2 \cos \theta_2 - a \cos \theta_3 + b \sin \theta_3 \quad (103)$$

$$g_{N4} = \frac{d}{2} + l_1 \sin \theta_1 + l_2 \sin \theta_2 + a \sin \theta_3 - b \cos \theta_3 \quad (104)$$

$$g_{T4} = l_1 \cos \theta_1 + l_2 \cos \theta_2 + a \cos \theta_3 + b \sin \theta_3 \quad (105)$$

Then, the **w** vectors and of the **\tilde{w}** scalars associated with each contact point can be obtained as

$$\mathbf{w}_{N1} = \frac{\partial g_{N1}}{\partial \mathbf{q}} = \begin{pmatrix} -l_1 \cos \theta_1 \\ -l_2 \cos \theta_2 \\ a \cos \theta_3 + b \sin \theta_3 \end{pmatrix} \quad (106)$$

$$\mathbf{w}_{T1} = \frac{\partial g_{T1}}{\partial \mathbf{q}} = \begin{pmatrix} -l_1 \sin \theta_1 \\ -l_2 \sin \theta_2 \\ a \sin \theta_3 - b \cos \theta_3 \end{pmatrix} \quad (107)$$

$$\mathbf{w}_{N2} = \frac{\partial g_{N2}}{\partial \mathbf{q}} = \begin{pmatrix} -l_1 \cos \theta_1 \\ -l_2 \cos \theta_2 \\ -a \cos \theta_3 + b \sin \theta_3 \end{pmatrix} \quad (108)$$

<figure><img src="image_9.png" /><figcaption>**Fig. 9** Generic position of the slider inside the guide where the distance between guide upper and lower surfaces is exaggerated for illustration purpose</figcaption></figure>
$$ \mathbf{w}_{T2} = \frac{\partial g_{T2}}{\partial \mathbf{q}} = \begin{pmatrix} -l_1 \sin \theta_1 \\ -l_2 \sin \theta_2 \\ -a \sin \theta_3 - b \cos \theta_3 \end{pmatrix} \qquad (109) $$

$$ \mathbf{w}_{N3} = \frac{\partial g_{N3}}{\partial \mathbf{q}} = \begin{pmatrix} l_1 \cos \theta_1 \\ l_2 \cos \theta_2 \\ -a \cos \theta_3 + b \sin \theta_3 \end{pmatrix} \qquad (110) $$

$$ \mathbf{w}_{T3} = \frac{\partial g_{T3}}{\partial \mathbf{q}} = \begin{pmatrix} -l_1 \sin \theta_1 \\ -l_2 \sin \theta_2 \\ a \sin \theta_3 + b \cos \theta_3 \end{pmatrix} \qquad (111) $$

$$ \mathbf{w}_{N4} = \frac{\partial g_{N4}}{\partial \mathbf{q}} = \begin{pmatrix} l_1 \cos \theta_1 \\ l_2 \cos \theta_2 \\ a \cos \theta_3 + b \sin \theta_3 \end{pmatrix} \qquad (112) $$

$$ \mathbf{w}_{T4} = \frac{\partial g_{T4}}{\partial \mathbf{q}} = \begin{pmatrix} -l_1 \sin \theta_1 \\ -l_2 \sin \theta_2 \\ -a \sin \theta_3 + b \cos \theta_3 \end{pmatrix} \qquad (113) $$

$$ \tilde{w}_{N1} = \tilde{w}_{T1} = \tilde{w}_{N2} = \tilde{w}_{T2} = \tilde{w}_{N3} = \tilde{w}_{T3} = \tilde{w}_{N4} = \tilde{w}_{T4} = 0 \qquad (114) $$

## 6.4 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 slider-crank mechanism with a translational clearance joint are listed in Table 1.

<table><caption>**Table 1** Parameters used in the dynamic simulation of the slider-crank mechanism</caption><tr><td>Geometrical characteristics</td><td>l<sub>1</sub> = 0.1530 m<br/>l<sub>2</sub> = 0.3060 m<br/>a = 0.0500 m<br/>b = 0.0250 m<br/>c = 0.0010 m</td></tr><tr><td>Inertial properties</td><td>m<sub>1</sub> = 0.0380 kg<br/>m<sub>2</sub> = 0.0380 kg<br/>m<sub>3</sub> = 0.0760 kg<br/>J<sub>1</sub> = 7.4 &times; 10<sup>-5</sup> kg m<sup>2</sup><br/>J<sub>2</sub> = 5.9 &times; 10<sup>-4</sup> kg m<sup>2</sup><br/>J<sub>3</sub> = 2.7 &times; 10<sup>-6</sup> kg m<sup>2</sup></td></tr><tr><td>Force elements</td><td>g = 9.81 m/s<sup>2</sup></td></tr><tr><td>Contact parameters</td><td>&epsilon;<sub>N1</sub> = &epsilon;<sub>N2</sub> = &epsilon;<sub>N3</sub> = &epsilon;<sub>N4</sub> = 0.4<br/>&epsilon;<sub>T1</sub> = &epsilon;<sub>T2</sub> = &epsilon;<sub>T3</sub> = &epsilon;<sub>T4</sub> = 0.0<br/>&mu;<sub>1</sub> = &mu;<sub>2</sub> = &mu;<sub>3</sub> = &mu;<sub>4</sub> = 0.01</td></tr><tr><td>Initial conditions</td><td>&theta;<sub>10</sub> = 0.0 rad<br/>&theta;<sub>20</sub> = 0.0 rad<br/>&theta;<sub>30</sub> = 0.0 rad<br/>&omega;<sub>10</sub> = 150.0 rad/s<br/>&omega;<sub>20</sub> = -75.0 rad/s<br/>&omega;<sub>30</sub> = 0.0 rad/s</td></tr></table>
<figure><img src="image_10.png" /><figcaption>Fig. 10 Dimensionless motion of the slider corners</figcaption></figure>

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:

$$ \frac{y_i - c}{c}, \quad (i = 1, 2, 3, 4) \qquad (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 guide-slider 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 well-known 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 [23].
<figure><img src="image_11.png" /></figure>

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

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

<figure><img src="image_14.png" /><figcaption>Fig. 11 (a) Crank speed; (b) connecting-rod speed; (c) connecting-rod portrait; (d) slider portrait</figcaption></figure>

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.

## 7 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 contact-impact 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
<figure><img src="image_15.png" /></figure>

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

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

<figure><img src="image_18.png" /><figcaption>Fig. 12 Influence of the restitution coefficient on the dimensionless motion of corner 1: (a) $\varepsilon = 0.1$; (b) $\varepsilon = 0.4$; (c) $\varepsilon = 0.6$; (d) $\varepsilon = 0.9$</figcaption></figure>

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

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

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

<figure><img src="image_22.png" /><figcaption>Fig. 13 Influence of the restitution coefficient on the slider portraits: (a) $\varepsilon = 0.1$; (b) $\varepsilon = 0.4$; (c) $\varepsilon = 0.6$; (d) $\varepsilon = 0.9$</figcaption></figure>
contacts, being the normal and tangential contact laws formulated as set-valued force laws for frictional unilateral constraints.

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.

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

## References

1. Flores, P., Ambrósio, J.: Revolute joints with clearance in multibody systems. Comput. Struct. **82**, 1359–1369 (2004)

2. Flores, P., Ambrósio, J., Claro, J.C.P., Lankarani, H.M.: Kinematics and Dynamics of Multibody Systems with Imperfect Joints: Models and Case Studies. Lecture Notes in Applied and Computational Mechanics, Vol. 34. Springer, Berlin (2008)

3. Qiang, T., Zhang, Y., Chen, L., Flores, P.: Dynamics of spatial flexible multibody systems with clearance and lubricated spherical joints. Comput. Struct. **87**, 913–929 (2009)

4. Liu, C.-S., Zhang, K., Yang, R.: The FEM analysis and approximate model for cylindrical joints with clearances. Mech. Mach. Theory **42**, 183–197 (2007)

5. Crowthera, A.R., Singha, R., Zhangb, N., Chapman, C.: Impulsive response of an automatic transmission system with multiple clearances: Formulation, simulation and experiment. J. Sound Vib. **306**, 444–466 (2007)

6. Flores, P., Ambrósio, J., Claro, J.P.: Dynamic analysis for planar multibody mechanical systems with lubricated joints. Multibody Syst. Dyn. **12**, 47–74 (2004)

7. Erkaya, S., Uzmay, İ.: A neural–genetic (NN–GA) approach for optimising mechanisms having joints with clearance. Multibody Syst. Dyn. **20**(1), 69–83 (2008)

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

9. Flores, P., Ambrósio, J., Claro, J.C.P., Lankarani, H.M.: Translational joints with clearance in rigid multibody systems. J. Comput. Nonlinear Dyn. **3**(1), 0110071–10 (2008)

10. Awrejcewicz, J., Kudra, G., Lamarque, C.-H.: Dynamics investigation of three coupled rods with horizontal barrier. Meccanica **38**, 687–698 (2003)

11. Stoenescu, E.D., Marghitu, D.B.: Dynamic analysis of a planar rigid-link mechanism with rotating slider joint and clearance. J. Sound Vib. **266**, 394–404 (2003)

12. Farahanchi, F., Shaw, S.W.: Chaotic and periodic dynamics of a slider crank mechanism with slider clearance. J. Sound Vib. **177**(3), 307–324 (1994)

13. Thümmel, T., Funk, K.: Multibody modelling of linkage mechanisms including friction, clearance and impact. In: Proceedings of Tenth World Congress on the Theory of Machine and Mechanisms, Oulu University, Finland, vol. 4, pp. 1375–1386 (1999)

14. Wilson, R., Fawcett, J.N.: Dynamics of slider-crank mechanism with clearance in the sliding bearing. Mech. Mach. Theory **9**, 61–80 (1994)

15. Cottle, R.W., Dantzig, G.B.: Complementary pivot theory of mathematical programming. Linear Algebra Appl. **1**, 103–125 (1968)
16. Cottle, R.W., Pang, J.S., Stone, R.E.: The Linear Complementarity Problem. Academic Press, London (1992)

17. Rohn, J.: A note on solvability of a class of linear complementarity problems. Math. Program. **60**(13), 229–231 (1993)

18. Leine, R.I., Brogliato, B., Nijmeijer, H.: Periodic motion and bifurcations induced by the Painlevé paradox. Eur. J. Mech. A, Solids **21**(5), 869–896 (2002)

19. Pang, J., Trinkle, J.C.: Complementarity formulations and existence of solutions of dynamic multi-rigid-body contact problems with Coulomb friction. Math. Program. **73**(2), 199–226 (1996)

20. Stewart, D.E., Trinkle, J.C.: An implicit time-stepping scheme for rigid body dynamics with inelastic collisions and Coulomb friction. Int. J. Numer. Methods Eng. **39**(15), 2673–2691 (1996)

21. 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)

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

23. 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)

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

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

26. 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)

27. Djerassi, S.: Collision with friction; Part A: Newton's hypothesis. Multibody Syst. Dyn. **21**, 37–54 (2009)

28. Djerassi, S.: Collision with friction; Part B: Poisson's and Stornge's hypotheses. Multibody Syst. Dyn. **21**, 55–70 (2009)

29. Bowling, A., Flickinger, D.M., Harmeyer, S.: Energetically consistent simulation of simultaneous impacts and contacts in multibody systems with friction. Multibody Syst. Dyn. **22**, 27–45 (2009)

30. 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)

31. Jean, M.: The non-smooth contact dynamics method. Comput. Methods Appl. Mech. Eng. **177**(3–4), 235–257 (1999)

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

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

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

35. 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)

36. Panagiotopoulos, P.D., Glocker, C.: Inequality constraints with elastic impacts in deformable bodies. The convex case. Arch. Appl. Mech. **70**(5), 349–365 (2000)

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

38. Moreau, J.J.: Bounded Variation in Time. In: Moreau, J.J., Panagiotopoulos, P.D., Strang, G. (eds.) Topics in Non-smooth Mechanics, pp. 1–74. Birkhäuser, Basel (1988)

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

40. Glocker, C.: On frictionless impact models in rigid-body systems. Philosophical Transactions: Math. Phys. Eng. Sci. **359**(1789), 2385–2404 (2001)

41. Anitescu, M., Potra, F.A., Stewart, D.E.: Time-stepping for three-dimensional rigid body dynamics. Comput. Methods Appl. Mech. Eng. **177**(3), 183–197 (1999)

42. 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)

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

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