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

HAL
open science

Unilateral problems of dynamics

F. Pfeiffer

► To cite this version:

F. Pfeiffer. Unilateral problems of dynamics. Archive of Applied Mechanics, 1999, 69 (8), pp.503 - 527. 10.1007/s004190050240. hal-01395964

HAL Id: hal-01395964

https://hal.science/hal-01395964

Submitted on 13 Nov 2016

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Unilateral problems of dynamics

F. Pfeiffer

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

**Key words** Unilateral contact, multibody dynamics, complementarity, friction, impact

1

**Introduction**

A large variety of contact configurations is encountered in mechanical and biomechanical systems. Contacts may be closed, and the colliding bodies may detach again. Within a closed contact we may have sliding or stiction, both features connected with local friction. If two bodies come into contact, they usually penetrate into each other which leads to local deformations. If contacts are accompanied by tangential forces and by tangential relative velocities within the contact plane, we obtain not only normal but also tangential deformations. De- pending on the dynamical (or statical) environment, contacts may change their state: from closure to detachment, from sliding to stiction, and vice versa. We call a contact active if it is closed or if there is stiction, otherwise we call a contact passive.

Active contacts always exhibit contact forces which, in the general case of normal and tangential deformations, follow from the local material properties of the colliding bodies, from the external (with respect to the contact) dynamics (statics) and from external forces. Con- sidering contacts from this perspective leads to complicated problems of continuum me- chanics which, as a rule, must be solved by numerical algorithms like FEM or BEM. This approach is too costly for the treatment of dynamical problems and, in many cases, also inadequate. The rigid-body approach gives quicker results, and it applies better to large dy- namical systems.

Under the rigid-body approach we understand a contact behavior characterized, at least in the local contact zone, by lack of deformations and, thus, by rigid-body properties. The contact process is then governed by certain contact laws, like those of Newton and Poisson, and appropriate extensions of them. A fundamental law with respect to the rigid-body model is the complementarity rule, sometimes called corner law or Signorini's law. It states that in contact dynamics either relative kinematical quantities are zero, and their accompanying constraint forces (or constraint force combinations) are not zero, or vice versa. For a closed contact, the relative normal distance and normal velocity of the colliding bodies are zero, and the constraint force in normal direction is not zero, or vice versa. For stiction, the tangential relative velocity is zero, and the constraint force is located within the friction cone, which means that the difference of the static friction force and constraint force is not zero, or vice versa. The resulting inequalities are indispensable for an evaluation of the transitions between the various contact states.
From these properties follows a well-defined indicator yielding the transition phases. For normal passive contacts, the normal relative distance (or velocity) of the colliding bodies indicates the contact state. If it becomes zero, the contact will be active; the indicator “relative distance” becomes a constraint accompanied by a constraint force. The end of an active contact is then indicated by the constraint force. If it changes sign, indicating a change from pressure to tension, we get detachment, and the contact again transits to a passive state.

In the case of friction, we have in the passive state a nonvanishing tangential relative velocity as an indicator. If it becomes zero, there might be a transition from sliding to stiction, depending on the force balance within or on the friction cone. The indicator “relative tangential velocity” becomes then a constraint, leading to tangential constraint forces. The contact remains active as long as the maximum static friction force is larger than the constraint force, which means that there is a force balance within the friction cone. If this “friction saturation” becomes zero, the contact again might go into a passive state with nonzero tangential relative velocity.

If we deal with multibody systems including unilateral constraints there arises the problem of multiple contacts and their interdependences. A straightforward solution of these processes would result in a combinatorial problem of a huge dimension. Therefore, a formulation applying complementarity rules and the resulting inequalities is a must. The mathematical methods developed in this area in the last twenty years assure nearly always an unambiguous solution, even in the cases of mutually dependent contacts.

The field of multibody dynamics with unilateral contacts should be seen before the background of various research activities: in the area of variational inequalities, in connection with convex energy functions on the one side, [20], [21], [27], [36], and in the area of hemivariational inequalities, in connection with non-convex energy functions, [28]. Both areas are young, the first one being some 35 years and the second only 15 years old. The two scientists who are the most connected with this development were Moreau in Montpellier, France, and Panagiotopoulos in Thessaloniki, Greece. Moreau started already in the seventies to formulate the discontinuous properties of nonsmooth mechanics, [21–23]; not much later, Panagiotopoulos considered inequality problems, [27], which led him consistently to the development of hemivariational inequalities, [28]. Most nonsmooth mechanical problems possess nonconvex features. Nearly all applications, at that time, were of statical or quasi-statical nature.

In Sweden, Lötstedt considered a series of practical examples, [15], [16], and established a school, which in more recent time has been continued very successfully by Klarbring and his group, [11], [12]. Klarbring also treats statical and quasistatistical problems with applications to FEM.

The sphere of influence of these scientists extends, since the significance of the area has been noticed. Especially [10] gives an excellent contribution to dry friction, also in connection with frictional impacts. The book [20] presents a proof with respect to the existence of Jean-Moreau impact theory. A very good survey on this topic may be found in [5]. In the meantime, the literature on the subject is increasing considerably, therefore only a few examples are given below.

Papers [13], [14] deal with two-dimensional contacts including dry friction. The resulting linear complementarity problem is solved by a modified simplex-algorithm. The authors of [18], [25] study self-excitations of frictional vibrations. Frictional constraints are described by variational inequalities and evaluated on the basis of measured friction characteristics. The paper [1] discusses, among other topics, also the existence of unambiguous solutions for generalized accelerations concerning frictional contact problems, a question, by the way, which also has been regarded in [15], [16].

Nonsmooth mechanics allows a general theoretical description, but relevant problems must be solved numerically. In spite of many valuable contributions to the numerics of discontinuous systems, [24], the existing algorithms are still extremely time-consuming. As it was indicated above, for plane contacts we get a linear and for spatial contacts a nonlinear complementarity problem. For linear complementarity problems there exist solution procedures, [24], on the basis of extended linear programming theories, and the existence of a solution is assured. This is not the case for nonlinear complementarity. At the time being, various approximations have been suggested, [11], [12], [17], and applied to practical problems, [38] through [44]. The algorithms work, in most cases, convincingly, but the computing times explode because of the iterative character of nearly all solution methods. Therefore, numerical solution of all kinds of complementarity problems is still a topic of current research.

The author's Institute pursues the research in applied and engineering mechanics, and has been involved in unilateral problems since more than fifteen years, exclusively in connection
with multibody dynamics. The fundamentals of the area may be found in books [3], [4], which present the multibody theory on the basis of the projection method, yielding a most efficient theory for multibody applications. In the course of the years, and very much due to practical requirements, this type of multibody theory has been combined with inequality formulations to describe the unilateral behavior of multiple contacts, [29–31]. Especially the complementarity properties, in combination with the kinematics and kinetics of the appropriate indicators, have become the main key to a consistent description of multibody systems with multiple unilateral contacts, [31]. An extension to impact, with the background of Moreau’s work, became then straightforward, [8], [9]. The corresponding theory has recently and successfully been confirmed by extensive measurements, [2].

Based on these fundamentals, a variety of industrial applications has been considered. The paper [6] and the dissertation [7] deal with roller chains used often in the automotive industry. Roller chains typically possess some hundred degrees of freedom and hundreds of chain-wheel contacts, where also stick-slip processes take place. Continuous Variable Transmissions (CVT) chains do not have such large amount of degrees of freedom, but they are particularly difficult due to spatial contacts within the conic wheel discs. Their efficiencies have been evaluated, with an excellent agreement between theory and measurements, [33], [39], [40].

Another area of application is the field of manufacturing assembly processes. Here, during a mating process, contacts are closed and detached, and stick-slip phenomena up to jamming occur. Efforts have been made to describe, in a general way, the relative kinematics of colliding bodies in terms of differential geometry, [19]. The results are convincing, and they are used now in other applications. Another example of manufacturing concerns vibration conveyors, often used to transport, to orient and to select small parts: screws, bolts, electronic components, and the like. Without the theory of multibody systems with unilateral contacts, a treatment of such a system would have been impossible. In the dissertation [44] a theoretical solution of this problem is given and experimentally proven.

In the following, we shall give a survey of the theory including impacts with friction, and we shall present some of its typical applications.

## 2

### Models

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

#### 2.1

##### Multibody dynamics

The theory of multibody dynamics of interconnected rigid and/or elastic bodies is nowadays well established, and commercialized to a large extent. Most of the existing computer codes are based on the so-called projection method, which considers free motion within hyper-spaces tangential to the constraint surfaces. The Jacobians performing these projections are derivations with respect to the constraint equations, usually formulated on a velocity level, [3], [4], [31]. Since equations of motion for multibody systems represent the basis for all further unilateral considerations, we shall give their short survey here.

We consider a multibody system with $f$ degrees of freedom, which we later shall subdivide into “rigid” and “elastic” degrees of freedom. Applying the principle of virtual power (Jourdain) we get

$$ \sum_{m_i} \int_{m_i} \left( \frac{\partial_{R} \mathbf{v}}{\partial \dot{\mathbf{q}}} \right)_i^T (R \mathbf{a} - R \mathbf{f})_i dm_i = 0, \quad (1) $$

where $R\mathbf{v}_i \in \mathbb{R}^3$, $R\mathbf{a}_i \in \mathbb{R}^3$ are the absolute velocity, acceleration of body $i$ in a body-fixed reference system, $R\mathbf{f}_i \in \mathbb{R}^3$ are all forces applied to body $i$ and $\dot{\mathbf{q}} \in \mathbb{R}^f$ are the generalized
velocities; $m_i$ is the mass of body i. The velocity and acceleration vectors must be evaluated in a body-fixed frame which yields, see Fig. 1,

$$ 
\begin{aligned}
R\mathbf{v} &= R\mathbf{v}_A + R\tilde{\omega}_{IR}(R\mathbf{x}_0 + R\bar{\mathbf{r}}_{el}) + R\dot{\bar{\mathbf{r}}}_{el}, \\
R\mathbf{a} &= R\mathbf{a}_A + R\tilde{\omega}_{IRR}\tilde{\omega}_{IRR}\mathbf{x}_0 + R\dot{\tilde{\omega}}_{IRR}\mathbf{x}_0 + R\tilde{\omega}_{IRR}\tilde{\omega}_{IRR}\bar{\mathbf{r}}_{el} \\
&\quad + R\dot{\tilde{\omega}}_{IRR}\bar{\mathbf{r}}_{el} + 2R\tilde{\omega}_{IRR}\dot{\bar{\mathbf{r}}}_{el} + R\ddot{\bar{\mathbf{r}}}_{el}.
\end{aligned} 
\quad (2) $$

The following abbreviations are used (index i always omitted): $\mathbf{v} = R\mathbf{v}_A$, $\mathbf{a} = R\mathbf{a}_A$ – absolute velocity, acceleration of the reference base $R$, expressed in the $R$-system, $\omega = R\omega_{IR}$ – angular velocity of the reference system expressed in the $R$-system, $\mathbf{x}_0 = R\mathbf{x}_0$ – vector from $R$ to the mass element $dm$, in the undeformed configuration, $\bar{\mathbf{r}} = R\bar{\mathbf{r}}_{el}$ – displacement vector, $\mathbf{J}_T = ((\partial_R\mathbf{v}_A)/\partial\dot{\mathbf{q}})$ – Jacobian of translation, $\mathbf{J}_T \in \mathbb{R}^{3f}$, $\mathbf{J}_R = ((\partial_R\omega_{IR})/\partial\dot{\mathbf{q}})$ – Jacobian of rotation, $\mathbf{J}_R \in \mathbb{R}^{3f}$, $\tilde{\mathbf{r}} = \omega \times \mathbf{r}$ – definition of tilde tensor.

Combining Eqs. (1), (2) and the above abbreviations, we come out with the following set of equations of motion:

$$ 
\begin{aligned}
& \sum \int_{m_i} \left[ \mathbf{J}_T^T + \mathbf{J}_R^T (\tilde{\mathbf{x}}_0 + \tilde{\mathbf{r}}) + \left( \frac{\partial \dot{\mathbf{r}}}{\partial \dot{\mathbf{q}}} \right)^T \right] \\
& \times (\mathbf{a} + \tilde{\omega}\tilde{\omega}\mathbf{x}_0 + \dot{\tilde{\omega}}\mathbf{x}_0 + \tilde{\omega}\tilde{\omega}\bar{\mathbf{r}} + \dot{\tilde{\omega}}\bar{\mathbf{r}} + 2\dot{\tilde{\omega}}\dot{\bar{\mathbf{r}}} + \ddot{\bar{\mathbf{r}}} - \mathbf{f}) dm = 0,
\end{aligned} 
\quad (3) $$

where elastic influences are regarded as shown in Fig. 1. The radius-vector $\mathbf{r}$ from an inertial frame (I) to a mass element $dm_i$ of the deformed body $B_i$ can be written as

$$ \mathbf{r} = \mathbf{r}_{IR} + \mathbf{x}_0 + \bar{\mathbf{r}}_{el}, \quad (4) $$

where $\mathbf{r}_{IR}$ is the vector from (I) to (R), $\mathbf{x}_0$ – the vector to the mass element in the undeformed configuration and $\bar{\mathbf{r}}_{el}$ – the displacement vector. As in most cases of technical relevancy, we assume that the elastic displacements of the components are very small compared to their geometric dimensions.

This allows the introduction of a Ritz approach for the elastic deformation, [4],

$$ \bar{\mathbf{r}}_i(\mathbf{x}_0, t) = \mathbf{W}_i(\mathbf{x}_0)\mathbf{q}_{el,i}(t), \quad \mathbf{W}_i = (\mathbf{w}_1, \dots, \mathbf{w}_j, \dots)_i = \begin{pmatrix} \bar{\mathbf{r}}_x^T \\ \bar{\mathbf{r}}_y^T \\ \bar{\mathbf{r}}_z^T \end{pmatrix}_i, \quad (5) $$

with $(\mathbf{q}_{el} \in \mathbb{R}^{n_{el}})_i$ and $(\mathbf{W} \in \mathbb{R}^{3,n_{el}})_i$. This well-known superposition of ansatz or shape functions $\mathbf{w}_{ji}(\mathbf{x}_0)$ requires their completeness, [4]. The shape functions may be evaluated by measurements, by FEM-calculations, by analytic approach, or they may be approximated by cubic-spline systems. In any case, they should be chosen in accordance with the elastic behavior of the system.

Combining Eqs. (3–5) we have to consider the dependencies $\mathbf{a}(\ddot{\mathbf{q}})$, $\dot{\omega}(\ddot{\mathbf{q}})$, $\ddot{\mathbf{q}}_{el,i}(\ddot{\mathbf{q}})$, which says that all absolute and elastic accelerations depend on the generalized accelerations $\ddot{\mathbf{q}}$. This property can be expressed by

<figure><img src="image_3.png" /><figcaption>Fig. 1. Coordinates for a deformed body $B_i$</figcaption></figure>
$$
\begin{align}
\mathbf{a} &= \mathbf{J}_T \ddot{\mathbf{q}} + \underline{\mathbf{a}}(\mathbf{q}, \dot{\mathbf{q}}, t), \nonumber \\
\dot{\mathbf{q}} &= \mathbf{J}_R \ddot{\mathbf{q}} + \underline{\dot{\mathbf{q}}}(\mathbf{q}, \dot{\mathbf{q}}, t). \tag{6}
\end{align}
$$

$$
\ddot{\mathbf{q}}_{el,i} = \frac{\partial \mathbf{q}_{el,i}}{\partial \mathbf{q}} \ddot{\mathbf{q}} = \frac{\partial \dot{\mathbf{q}}_{el,i}}{\partial \dot{\mathbf{q}}} \ddot{\mathbf{q}} = \mathbf{J}_E \ddot{\mathbf{q}} . \quad (7)
$$

Equations (3) to (7) can now be put into the form

$$
\sum_i \left\{ \int_{m_i} J_T^T [\mathbf{J}_T - (\tilde{\mathbf{x}}\tilde{\mathbf{r}})\mathbf{J}_R + \mathbf{W}\mathbf{J}_E] dm + \int_{m_i} J_R^T (\tilde{\mathbf{x}}_0 + \tilde{\mathbf{r}}) [\mathbf{J}_T - (\tilde{\mathbf{x}}_0 + \tilde{\mathbf{r}})\mathbf{J}_R + \mathbf{W}\mathbf{J}_E] dm \right. \\
\left. + \int_{m_i} J_E^T \mathbf{W}^T [\mathbf{J}_T - (\tilde{\mathbf{x}}_0 + \tilde{\mathbf{r}})\mathbf{J}_R + \mathbf{W}\mathbf{J}_E] dm \right\} \ddot{\mathbf{q}} + \int_{m_i} [J_T^T + J_R^T (\tilde{\mathbf{x}}_0 + \tilde{\mathbf{r}}) + J_E^T \mathbf{W}^T] \\
\times [\underline{\mathbf{a}} + \tilde{\mathbf{\omega}}\tilde{\mathbf{\omega}}(\tilde{\mathbf{x}}_0 + \tilde{\mathbf{r}}) + \underline{\dot{\mathbf{\omega}}}(\tilde{\mathbf{x}}_0 + \tilde{\mathbf{r}}) + 2\tilde{\mathbf{\omega}}\dot{\tilde{\mathbf{r}}} - \mathbf{f}] dm = 0. \quad (8)
$$

After lengthy and tedious calculations, it is always possible to bring Eqs. (8) into a standard
form which we need for all further considerations. It writes

$$
\mathbf{M}(\mathbf{q}, t)\ddot{\mathbf{q}} - \mathbf{h}(\mathbf{q}, \dot{\mathbf{q}}, t) = 0, \quad \mathbf{q} \in \mathbb{R}^f, \mathbf{h} \in \mathbb{R}^f, \mathbf{M} \in \mathbb{R}^{f \times f}. \qquad (9)
$$

Equations (9) include all bilateral constraints. They represent the maximum number of
minimal (generalized) coordinates, and from this all unilateral contacts are in a state that they
do not block any further degree of freedom. If some of the unilateral constraints become active,
then the remaining number of degrees of freedom is smaller than *f*.

2.2

**Contact kinematics**

Geometry and kinematics are the fundaments in establishing models of dynamical systems. In the case of unilateral contacts, this is especially important because magnitudes of relative kinematics serve as indicators for passive contacts and as constraints for active contacts, see Sec. 1. As most of the applications require more or less arbitrary body contours, it makes sense to derive the kinematical contact equations in a general form applying well-known rules of the differential geometry of surfaces, [8], [19], [31].

In practice we find two types of contacts, two- and three-dimensional ones. For the two-
dimensional case, the contacting bodies lie in a plane thus defining a contact line with given
direction. Only the sense of direction has to be determined. Problems of that kind are con-
nected with linear complementarity.

For the three-dimensional case, the contacting bodies have a spatial form, and the contact takes place in a plane allowing two tangential directions. The resulting direction for the contact process is not known beforehand, and usually has to be determined iteratively. Problems of that kind are connected with nonlinear complementarity.

2.2.1

**Plane contact kinematics**

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

$$
\begin{align}
& {}_B\mathbf{t} = {}_B\mathbf{r}'_{P\Sigma}; \quad \kappa_B \mathbf{n} = {}_B\mathbf{r}''_{P\Sigma}; \quad (\cdot)' = \frac{d}{ds}; \notag \\
& {}_B\mathbf{n} = {}_B\mathbf{b} \times {}_B\mathbf{t}; \quad {}_B\mathbf{b} = {}_B\mathbf{t} \times {}_B\mathbf{n}; \quad {}_B\mathbf{t} = {}_B\mathbf{n} \times {}_B\mathbf{b}. \tag{10} \\
& {}_B\mathbf{n}' = {}_B\mathbf{b} \times {}_B\mathbf{t}' = {}_B\mathbf{b} \times {}_B\mathbf{n}\kappa = -\kappa_B \mathbf{t}, \quad {}_B\mathbf{t}' = \kappa_B \mathbf{n}. \notag
\end{align}
$$

For planar contours the binormal $\mathbf{b}$ is constant. Therefore

$$
{}_B\dot{\mathbf{n}} = {}_B\mathbf{n}'\dot{s} = -\kappa \dot{s} {}_B\mathbf{t}, \quad {}_B\dot{\mathbf{t}} = {}_B\mathbf{t}'\dot{s} = +\kappa \dot{s} {}_B\mathbf{n}. \qquad (11)
$$
<figure><img src="image_4.png" /><figcaption>Fig. 2. Planar contour geometry</figcaption></figure>

On the other hand, the absolute changes of **n** and **t** are given by the Coriolis equation

$$B(\dot{\mathbf{n}}) = B\dot{\mathbf{n}} + B\tilde{\Omega}_B\mathbf{n}, \quad B(\dot{\mathbf{t}}) = B\dot{\mathbf{t}} + B\tilde{\Omega}_B\mathbf{t}, \qquad (12)$$

where we must keep in mind that $B\omega_{IB} = B\Omega$ for body-fixed frames B. Putting (11) into (12) we get a coordinate-free representation of the overall changes $\dot{\mathbf{n}}, \dot{\mathbf{t}}$

$$\dot{\mathbf{n}} = \tilde{\Omega}\mathbf{n} - \kappa\dot{s}\mathbf{t}, \quad \dot{\mathbf{t}} = \tilde{\Omega}\mathbf{t} + \kappa\dot{s}\mathbf{n}, \qquad (13)$$

which we can evaluate in any basis. The main advantage of (13) consists of the eliminated, frame-dependent differentiations $B\dot{\mathbf{n}}$ and $B\dot{\mathbf{t}}$.

In the same manner we proceed with the contour vector $\mathbf{r}_{P\Sigma}$. According to (11), (12) we write

$$B\dot{\mathbf{r}}_{P\Sigma} = B\mathbf{r}'_{P\Sigma}\dot{s} = \dot{s}_B\mathbf{t}, \quad B(\dot{\mathbf{r}}_{P\Sigma}) = B\dot{\mathbf{r}}_{P\Sigma} + B\tilde{\Omega}_B\mathbf{r}_{P\Sigma}, \qquad (14)$$

and eliminate $B\dot{\mathbf{r}}_{P\Sigma}$. Then we get the absolute changes of $\mathbf{r}_{P\Sigma}$,

$$\dot{\mathbf{r}}_{P\Sigma} = \tilde{\Omega}\mathbf{r}_{P\Sigma} + \dot{s}\mathbf{t}. \qquad (15)$$

Due to $\mathbf{v}_{\Sigma} = \mathbf{v}_P + \dot{\mathbf{r}}_{P\Sigma}$, the absolute velocity of the moving contour point is given by

$$\mathbf{v}_{\Sigma} = \mathbf{v}_P + \tilde{\Omega}\mathbf{r}_{P\Sigma} + \dot{s}\mathbf{t}, \qquad (16)$$

where

$$\mathbf{v}_C := \mathbf{v}_P + \tilde{\Omega}\mathbf{r}_{P\Sigma}. \qquad (17)$$

The velocity $\mathbf{v}_C$ results from rigid-body kinematics, and corresponds to the velocity of a body-fixed point at the contour. From (16) and (17) we see that

$$\mathbf{v}_{\Sigma} = \mathbf{v}_C + \dot{s}\mathbf{t}. \qquad (18)$$

Next, we want to derive the absolute acceleration of C by differentiating (17) with respect to time

$$\dot{\mathbf{v}}_C = \dot{\mathbf{v}}_P + \dot{\tilde{\Omega}}\mathbf{r}_{P\Sigma} + \tilde{\Omega}\dot{\mathbf{r}}_{P\Sigma}. \qquad (19)$$

With $\dot{\mathbf{v}}_C = \mathbf{a}_C$, $\dot{\mathbf{v}}_P = \mathbf{a}_P$ and $\dot{\mathbf{r}}_{P\Sigma}$, from Eq. (14) we get

$$\mathbf{a}_C = \mathbf{a}_P + \dot{\tilde{\Omega}}\mathbf{r}_{P\Sigma} + \tilde{\Omega}\tilde{\Omega}\mathbf{r}_{P\Sigma} + \tilde{\Omega}\dot{s}\mathbf{t}, \qquad (20)$$

which is not the acceleration of a body-fixed point on the contour. Only the part

$$\mathbf{a}_Q := \mathbf{a}_P + \dot{\tilde{\Omega}}\mathbf{r}_{P\Sigma} + \tilde{\Omega}\tilde{\Omega}\mathbf{r}_{P\Sigma}, \qquad (21)$$
corresponds to such an acceleration, so we can write

$$ \mathbf{a}_C = \mathbf{a}_Q + \tilde{\Omega}\mathbf{t}\dot{s}. \qquad (22) $$

Later we have to determine the relative velocities of contact points in the normal and tangential directions and their time derivatives. For this purpose we introduce the scalars

$$ v_n = \mathbf{n}^T \mathbf{v}_C; \quad v_t = \mathbf{t}^T \mathbf{v}_C \qquad (23) $$

and state their derivatives as

$$ \dot{v}_n = \dot{\mathbf{n}}^T \mathbf{v}_C + \mathbf{n}^T \dot{\mathbf{v}}_C, \quad \dot{v}_t = \dot{\mathbf{t}}^T \mathbf{v}_C + \mathbf{t}^T \dot{\mathbf{v}}_C . \qquad (24) $$

With $\dot{\mathbf{n}}, \dot{\mathbf{t}}$ from (13), $\dot{\mathbf{v}}_C = \mathbf{a}_C$ from (22), and noting $\mathbf{n}^T \tilde{\Omega} \mathbf{t} = \mathbf{b}^T \Omega$, $\mathbf{t}^T \tilde{\Omega} \mathbf{t} = 0$, we derive

$$ \begin{aligned} \dot{v}_n &= \mathbf{n}^T (\mathbf{a}_Q - \tilde{\Omega} \mathbf{v}_C) - \kappa \dot{s} \mathbf{t}^T \mathbf{v}_C + \dot{s} \mathbf{b}^T \Omega, \\ \dot{v}_t &= \mathbf{t}^T (\mathbf{a}_Q - \tilde{\Omega} \mathbf{v}_C) + \kappa \dot{s} \mathbf{n}^T \mathbf{v}_C . \end{aligned} \qquad (25) $$

With this basis, we are able to derive the relative kinematics of two bodies such as relative distances, relative velocities and accelerations. Without going into details and refering especially to [31], we summarize the following relationships, see Fig. 3.

Potential contact points can be characterized by

$$ \mathbf{n}_1^T(s_1) \cdot \mathbf{t}_2(s_2) = 0 \Leftrightarrow \mathbf{n}_2^T(s_2) \cdot \mathbf{t}_1(s_1) = 0, \qquad (26) $$

$$ \mathbf{r}_D^T(s_1, s_2) \cdot \mathbf{t}_1(s_1) = 0, \quad \mathbf{r}_D^T(s_1, s_2) \cdot \mathbf{t}_2(s_2) = 0. \qquad (27) $$

From each set we need only one equation. The relative distance $g_N$ is

$$ g_N(\mathbf{q}, t) = \mathbf{r}_D^T \mathbf{n}_2 = -\mathbf{r}_D^T \mathbf{n}_1 . \qquad (28) $$

Since the normal vectors always point inward, $g_N$ is positive for separation and negative for overlapping. Therefore, a changing sign of $g_N$ from positive to negative indicates a transition from initially separated bodies to contact.

With these equations, and considering Fig. 3, we derive the relative velocities in normal and tangential direction

$$ \dot{g}_N = \mathbf{n}_1^T \mathbf{v}_{C1} + \mathbf{n}_2^T \mathbf{v}_{C2}; \quad \dot{g}_T = \mathbf{t}_1^T \mathbf{v}_{C1} + \mathbf{t}_2^T \mathbf{v}_{C2}, \qquad (29) $$

where $\mathbf{v}_{C1}, \mathbf{v}_{C2}$ are the absolute velocities of the potential contact points $C_1, C_2$. These velocities may be expressed by the generalized velocities $\dot{\mathbf{q}}$ and some Jacobians $\mathbf{J}_{C1}, \mathbf{J}_{C2}$ to give

$$ \mathbf{v}_{C1} = \mathbf{J}_{C1} \dot{\mathbf{q}} + \tilde{\mathbf{j}}_{C1}; \quad \mathbf{v}_{C2} = \mathbf{J}_{C2} \dot{\mathbf{q}} + \tilde{\mathbf{j}}_{C2}. \qquad (30) $$

<figure><img src="image_5.png" /><figcaption>Fig. 3. General orientation of two bodies</figcaption></figure>
Putting (29) into (30) yields

$$ \dot{g}_N = \mathbf{w}_N^T \dot{\mathbf{q}} + \tilde{w}_N; \quad \dot{g}_T = \mathbf{w}_T^T \dot{\mathbf{q}} + \tilde{w}_T, \qquad (31) $$

with

$$ \mathbf{w}_N = \mathbf{J}_{C1}^T \mathbf{n}_1 + \mathbf{J}_{C2}^T \mathbf{n}_2, \quad \mathbf{w}_T = \mathbf{J}_{C1}^T \mathbf{t}_1 + \mathbf{J}_{C2}^T \mathbf{t}_2, \\ \tilde{w}_N = \tilde{\mathbf{j}}_{C1}^T \mathbf{n}_1 + \tilde{\mathbf{j}}_{C2}^T \mathbf{n}_2, \quad \tilde{w}_T = \tilde{\mathbf{j}}_{C1}^T \mathbf{t}_1 + \tilde{\mathbf{j}}_{C2}^T \mathbf{t}_2, \qquad (32) $$

which we use in the following as a representation of the relative velocities. It may be noticed here that a negative value of $\dot{g}_N$ corresponds to an approaching process of the bodies and coincides at vanishing distance $g_N = 0$ with the relative velocity in the normal direction before an impact. In the case of a continuous contact ($g_N = \dot{g}_N = 0$), the term $\dot{g}_T$ shows the relative sliding velocity of the bodies, which we can use to determine the time points of transitions from sliding ($\dot{g}_T \neq 0$) to sticking or rolling ($\dot{g}_T = 0$).

The relevant accelerations follow from a further time differentiation. We get

$$ \ddot{g}_N = \mathbf{w}_N^T \ddot{\mathbf{q}} + \bar{w}_N; \quad \ddot{g}_T = \mathbf{w}_T^T \ddot{\mathbf{q}} + \bar{w}_T, \qquad (33) $$

where $\mathbf{w}_N, \mathbf{w}_T$ are given by (32), and $\bar{w}_N, \bar{w}_T$ are

$$ \bar{w}_N = \mathbf{n}_1^T (\bar{\mathbf{j}}_{Q1} - \tilde{\mathbf{\Omega}}_1 \mathbf{v}_{C1}) - \kappa_1 \dot{s}_1 \mathbf{t}_1^T \mathbf{v}_{C1} + \dot{s}_1 \mathbf{b}_{12}^T \mathbf{\Omega}_1 + \mathbf{n}_2^T (\bar{\mathbf{j}}_{Q2} - \tilde{\mathbf{\Omega}}_2 \mathbf{v}_{C2}) - \kappa_2 \dot{s}_2 \mathbf{t}_2^T \mathbf{v}_{C2} + \dot{s}_2 \mathbf{b}_{12}^T \mathbf{\Omega}_2, \\ \bar{w}_T = \mathbf{t}_1^T (\bar{\mathbf{j}}_{Q1} - \tilde{\mathbf{\Omega}}_1 \mathbf{v}_{C1}) + \kappa_1 \dot{s}_1 \mathbf{n}_1^T \mathbf{v}_{C1} + \mathbf{t}_2^T (\bar{\mathbf{j}}_{Q2} - \tilde{\mathbf{\Omega}}_2 \mathbf{v}_{C2}) + \kappa_2 \dot{s}_2 \mathbf{n}_2^T \mathbf{v}_{C2}, \qquad (34) $$

with

$$ \dot{s}_1 = \frac{\kappa_2 \mathbf{t}_1^T (\mathbf{v}_{C2} - \mathbf{v}_{C1}) - \kappa_2 g_N \mathbf{b}_{12}^T \mathbf{\Omega}_1 + \mathbf{b}_{12}^T (\mathbf{\Omega}_2 - \mathbf{\Omega}_1)}{\kappa_1 + \kappa_2 + g_N \kappa_1 \kappa_2}, \\ \dot{s}_2 = \frac{\kappa_1 \mathbf{t}_1^T (\mathbf{v}_{C2} - \mathbf{v}_{C1}) - \kappa_1 g_N \mathbf{b}_{12}^T \mathbf{\Omega}_2 - \mathbf{b}_{12}^T (\mathbf{\Omega}_2 - \mathbf{\Omega}_1)}{\kappa_1 + \kappa_2 + g_N \kappa_1 \kappa_2}. \qquad (35) $$

The angular velocities $\mathbf{\Omega}_1, \mathbf{\Omega}_2$ relate to the two contacting bodies, Figs. 2, 3.

### 2.2.2 Spatial contact kinematics

Spatial contact kinematics has been presented in the dissertation [19], and since then also has been applied to many practical problems. The situation in this case is, of course, more complex. We still assume that the two approaching bodies are convex, Fig. 4, at least in

<figure><img src="image_6.png" /><figcaption>Fig. 4. Contact geometry of two surfaces</figcaption></figure>
that area where contact points might occur. The two bodies are moving with $\mathbf{v}_i$, $\mathbf{\Omega}_i$, $i = 1, 2$. For the description of a surface $\Sigma$ we need two parameters $s$ and $t$: $\mathbf{r}_\Sigma = \mathbf{r}_\Sigma(s, t)$. The tangents $\mathbf{s}$ and $\mathbf{t}$, which span the tangent plane at a point of the surface, are defined as

$$ \mathbf{s} = \frac{\partial \mathbf{r}_\Sigma}{\partial s}, \quad \mathbf{t} = \frac{\partial \mathbf{r}_\Sigma}{\partial t}. \qquad (36) $$

From these basic vectors, the fundamental magnitudes of the first order are calculated

$$ E = \mathbf{s}^T \mathbf{s}, \quad F = \mathbf{s}^T \mathbf{t}, \quad G = \mathbf{t}^T \mathbf{t}. \qquad (37) $$

The normalized normal vector $\mathbf{n}$ is perpendicular to the tangent plane and pointing outwards

$$ \mathbf{n} = \frac{\mathbf{s} \times \mathbf{t}}{\sqrt{EG - F^2}}. \qquad (38) $$

We further need the fundamental magnitudes of the second order

$$ L = \mathbf{n}^T \frac{\partial^2 \mathbf{r}_\Sigma}{\partial s^2}, \quad M = \mathbf{n}^T \frac{\partial^2 \mathbf{r}_\Sigma}{\partial s \partial t}, \quad N = \mathbf{n}^T \frac{\partial^2 \mathbf{r}_\Sigma}{\partial t^2}. \qquad (39) $$

For a contact point we demand that the normal vector of body 1 ($\mathbf{n}_1$) and the distance vector $\mathbf{r}_D$ are perpendicular to the tangent vectors of body 2 ($\mathbf{s}_2$ and $\mathbf{t}_2$). Thus we obtain four nonlinear equations

$$ \mathbf{n}_1^T \mathbf{s}_2 = 0, \quad \mathbf{r}_D^T \mathbf{s}_2 = 0, \quad \mathbf{n}_1^T \mathbf{t}_2 = 0, \quad \mathbf{r}_D^T \mathbf{t}_2 = 0. \qquad (40) $$

This nonlinear problem has to be solved at every time step of the numerical integration. After the solution is found, the distance $g_N$ between the possible contact points can be calculated as

$$ g_N = \mathbf{n}_1^T \mathbf{r}_D = -\mathbf{n}_2^T \mathbf{r}_D. \qquad (41) $$

Here, $g_N$ is used as an indicator for the contact state. Its value is positive for ‘no contact’, and negative for penetration. The constraints are again formulated on the velocity level, where in the spatial case we have three of them, one in the normal direction $\dot{\mathbf{g}}_N$ and two in the tangential directions $\dot{\mathbf{g}}_S, \dot{\mathbf{g}}_T$

$$ \begin{aligned} \dot{\mathbf{g}}_N(\mathbf{q}, \dot{\mathbf{q}}, t) &= \mathbf{n}_1^T (\mathbf{v}_{\Sigma 2} - \mathbf{v}_{\Sigma 1}), \\ \dot{\mathbf{g}}_S(\mathbf{q}, \dot{\mathbf{q}}, t) &= \mathbf{s}_1^T (\mathbf{v}_{\Sigma 2} - \mathbf{v}_{\Sigma 1}), \\ \dot{\mathbf{g}}_T(\mathbf{q}, \dot{\mathbf{q}}, t) &= \mathbf{t}_1^T (\mathbf{v}_{\Sigma 2} - \mathbf{v}_{\Sigma 1}), \end{aligned} \qquad (42) $$

with $\mathbf{v}_{\Sigma 1}$ and $\mathbf{v}_{\Sigma 2}$ being defined in an analogous way as in (30). Differentiating these equations with respect to time leads to the constraints on the acceleration level

$$ \begin{aligned} \ddot{\mathbf{g}}_N &= \mathbf{n}_1^T (\dot{\mathbf{v}}_{\Sigma 2} - \dot{\mathbf{v}}_{\Sigma 1}) + \dot{\mathbf{n}}_1^T (\mathbf{v}_{\Sigma 2} - \mathbf{v}_{\Sigma 1}), \\ \ddot{\mathbf{g}}_S &= \mathbf{s}_1^T (\dot{\mathbf{v}}_{\Sigma 2} - \dot{\mathbf{v}}_{\Sigma 1}) + \dot{\mathbf{s}}_1^T (\mathbf{v}_{\Sigma 2} - \mathbf{v}_{\Sigma 1}), \\ \ddot{\mathbf{g}}_T &= \mathbf{t}_1^T (\dot{\mathbf{v}}_{\Sigma 2} - \dot{\mathbf{v}}_{\Sigma 1}) + \dot{\mathbf{t}}_1^T (\mathbf{v}_{\Sigma 2} - \mathbf{v}_{\Sigma 1}). \end{aligned} \qquad (43) $$

The time derivatives of the contact point velocities $\mathbf{v}_{\Sigma 1}$ and $\mathbf{v}_{\Sigma 2}$ can be written in the form

$$ \dot{\mathbf{v}}_{\Sigma 1} = \mathbf{J}_{\Sigma 1}(\mathbf{q}, t)\ddot{\mathbf{q}} + \bar{\mathbf{j}}_{\Sigma 1}(\dot{\mathbf{q}}, \mathbf{q}, t), \quad \dot{\mathbf{v}}_{\Sigma 2} = \mathbf{J}_{\Sigma 2}(\mathbf{q}, t)\ddot{\mathbf{q}} + \bar{\mathbf{j}}_{\Sigma 2}(\dot{\mathbf{q}}, \mathbf{q}, t). \qquad (44) $$

The vectors $\dot{\mathbf{n}}_1, \dot{\mathbf{s}}_1$ and $\dot{\mathbf{t}}_1$ are determined by the formulas of Weingarten and Gauss, which express the derivatives of the normal vector and of the tangent vectors in terms of the basic vectors

$$ \dot{\mathbf{n}}_1 = \mathbf{\Omega}_1 \times \mathbf{n}_1 + \frac{\partial \mathbf{n}_1}{\partial s_1} \dot{\mathbf{s}}_1 + \frac{\partial \mathbf{n}_1}{\partial t_1} \dot{\mathbf{t}}_1, \qquad (45) $$
where

$$ \frac{\partial \mathbf{n}_1}{\partial s_1} = \underbrace{\frac{M_1 F_1 - L_1 G_1}{E_1 G_1 - F_1^2}}_{\alpha_1} \mathbf{s}_1 + \underbrace{\frac{L_1 F_1 - M_1 E_1}{E_1 G_1 - F_1^2}}_{\beta_1} \mathbf{t}_1, \quad \frac{\partial \mathbf{n}_1}{\partial t_1} = \underbrace{\frac{N_1 F_1 - M_1 G_1}{E_1 G_1 - F_1^2}}_{\alpha'_1} \mathbf{s}_1 + \underbrace{\frac{M_1 F_1 - N_1 E_1}{E_1 G_1 - F_1^2}}_{\beta'_1} \mathbf{t}_1, $$

$$ \dot{\mathbf{s}}_1 = \mathbf{\Omega}_1 \times \mathbf{s}_1 + \frac{\partial \mathbf{s}_1}{\partial s_1} \dot{s}_1 + \frac{\partial \mathbf{s}_1}{\partial t_1} \dot{t}_1, \quad (46) $$

where

$$ \frac{\partial \mathbf{s}_1}{\partial s_1} = \Gamma_{11,1}^1 \mathbf{s}_1 + \Gamma_{11,1}^2 \mathbf{t}_1 + L_1 \mathbf{n}_1, \quad \frac{\partial \mathbf{s}_1}{\partial t_1} = \Gamma_{12,1}^1 \mathbf{s}_1 + \Gamma_{12,1}^2 \mathbf{t}_1 + M_1 \mathbf{n}_1, $$

$$ \dot{\mathbf{t}}_1 = \mathbf{\Omega}_1 \times \mathbf{t}_1 + \frac{\partial \mathbf{t}_1}{\partial s_1} \dot{s}_1 + \frac{\partial \mathbf{t}_1}{\partial t_1} \dot{t}_1, \quad (47) $$

where

$$ \frac{\partial \mathbf{t}_1}{\partial s_1} = \Gamma_{12,1}^1 \mathbf{s}_1 + \Gamma_{12,1}^2 \mathbf{t}_1 + M_1 \mathbf{n}_1, \quad \frac{\partial \mathbf{t}_1}{\partial t_1} = \Gamma_{22,1}^1 \mathbf{s}_1 + \Gamma_{22,1}^2 \mathbf{t}_1 + N_1 \mathbf{n}_1. $$

The definition of the Christoffel symbols $\Gamma^\Sigma_{\alpha\beta}$, $\alpha, \beta, \sigma = 1, 2$, can be found in standard textbooks. Inserting Eqs. (45)–(47) in (43) yields the constraint equations

$$ \begin{aligned} \ddot{\mathbf{g}}_N &= \mathbf{n}_1^T (\mathbf{J}_{\Sigma2} - \mathbf{J}_{\Sigma1}) \ddot{\mathbf{q}} + \mathbf{n}_1^T (\bar{\mathbf{j}}_{\Sigma2} - \bar{\mathbf{j}}_{\Sigma1}) + (\mathbf{v}_{\Sigma2} - \mathbf{v}_{\Sigma1})^T (\mathbf{\Omega}_1 \times \mathbf{n}_1) + (\mathbf{v}_{\Sigma2} - \mathbf{v}_{\Sigma1})^T \\ & \times [(\alpha_1 \mathbf{s}_1 + \beta_1 \mathbf{t}_1) \dot{s}_1 + (\alpha'_1 \mathbf{s}_1 + \beta'_1 \mathbf{t}_1) \dot{t}_1], \\ \ddot{\mathbf{g}}_S &= \mathbf{s}_1^T (\mathbf{J}_{\Sigma2} - \mathbf{J}_{\Sigma1}) \ddot{\mathbf{q}} + \mathbf{s}_1^T (\bar{\mathbf{j}}_{\Sigma2} - \bar{\mathbf{j}}_{\Sigma1}) + (\mathbf{v}_{\Sigma2} - \mathbf{v}_{\Sigma1})^T (\mathbf{\Omega}_1 \times \mathbf{s}_1) + (\mathbf{v}_{\Sigma2} - \mathbf{v}_{\Sigma1})^T \\ & \times [(\Gamma_{11,1}^1 \mathbf{s}_1 + \Gamma_{11,1}^2 \mathbf{t}_1 + L_1 \mathbf{n}_1) \dot{s}_1 + (\Gamma_{12,1}^1 \mathbf{s}_1 + \Gamma_{12,1}^2 \mathbf{t}_1 + M_1 \mathbf{n}_1) \dot{t}_1], \end{aligned} \quad (48) $$

$$ \begin{aligned} \ddot{\mathbf{g}}_T &= \mathbf{t}_1^T (\mathbf{J}_{\Sigma2} - \mathbf{J}_{\Sigma1}) \ddot{\mathbf{q}} + \mathbf{t}_1^T (\bar{\mathbf{j}}_{\Sigma2} - \bar{\mathbf{j}}_{\Sigma1}) + (\mathbf{v}_{\Sigma2} - \mathbf{v}_{\Sigma1})^T (\mathbf{\Omega}_1 \times \mathbf{t}_1) + (\mathbf{v}_{\Sigma2} - \mathbf{v}_{\Sigma1})^T \\ & \times [(\Gamma_{12,1}^1 \mathbf{s}_1 + \Gamma_{12,1}^2 \mathbf{t}_1 + M_1 \mathbf{n}_1) \dot{s}_1 + (\Gamma_{22,1}^1 \mathbf{s}_1 + \Gamma_{22,1}^2 \mathbf{t}_1 + N_1 \mathbf{n}_1) \dot{t}_1]. \end{aligned} $$

As can be seen, these equations depend only on the Jacobians with respect to the contact points $\mathbf{J}_{\Sigma1}, \mathbf{J}_{\Sigma2}$, the basic vectors of the surfaces and the time derivatives of the contour parameters $\dot{s}_1, \dot{t}_1, \dot{s}_2, \dot{t}_2$. The Jacobians are known from the rigid body algorithm, the basic vectors from the surface description. The time derivatives of the contour parameters can be calculated by deriving Eq. (40) with respect to time

$$ (\mathbf{n}_1^T \mathbf{s}_2)^\bullet = 0, \quad (\mathbf{r}_D^T \mathbf{s}_2)^\bullet = 0, \quad (\mathbf{n}_1^T \mathbf{t}_2)^\bullet = 0, \quad (\mathbf{r}_D^T \mathbf{t}_2)^\bullet = 0, \quad (49) $$

which means, that the conditions for the contact point should not change while the two bodies are moving. Evaluating Eq. (49), we obtain a system of equations which are linear in the derivatives of the contour parameters

$$ \begin{pmatrix} \mathbf{s}_2^T (\alpha_1 \mathbf{s}_1 + \beta_1 \mathbf{t}_1) & \mathbf{s}_2^T (\alpha'_1 \mathbf{s}_1 + \beta'_1 \mathbf{t}_1) & L_2 & M_2 \\ \mathbf{t}_2^T (\alpha_1 \mathbf{s}_1 + \beta_1 \mathbf{t}_1) & \mathbf{t}_2^T (\alpha'_1 \mathbf{s}_1 + \beta'_1 \mathbf{t}_1) & M_2 & N_2 \\ -\mathbf{s}_1^T \mathbf{s}_2 & -\mathbf{s}_1^T \mathbf{s}_2 & \mathbf{s}_2^T \mathbf{s}_2 & \mathbf{s}_2^T \mathbf{t}_2 \\ -\mathbf{s}_1^T \mathbf{t}_2 & -\mathbf{s}_1^T \mathbf{t}_2 & \mathbf{s}_2^T \mathbf{t}_2 & \mathbf{t}_2^T \mathbf{t}_2 \end{pmatrix} \cdot \begin{pmatrix} \dot{s}_1 \\ \dot{t}_1 \\ \dot{s}_2 \\ \dot{t}_2 \end{pmatrix} = \begin{pmatrix} (\mathbf{s}_2 \times \mathbf{n}_1)^T (\mathbf{\Omega}_2 - \mathbf{\Omega}_1) \\ (\mathbf{t}_2 \times \mathbf{n}_1)^T (\mathbf{\Omega}_2 - \mathbf{\Omega}_1) \\ \mathbf{s}_2^T (\mathbf{v}_{\Sigma1} - \mathbf{v}_{\Sigma2}) \\ \mathbf{t}_2^T (\mathbf{v}_{\Sigma1} - \mathbf{v}_{\Sigma2}) \end{pmatrix} \quad (50) $$

This linear problem has to be solved at every time step of numerical integration. Let us summarize the constraint equations in the well known form, by rewriting Eq. (48) as

$$ \ddot{\mathbf{g}}_N = \mathbf{w}_N^T \ddot{\mathbf{q}} + \bar{\mathbf{w}}_N, \quad \ddot{\mathbf{g}}_S = \mathbf{w}_S^T \ddot{\mathbf{q}} + \bar{\mathbf{w}}_S, \quad \ddot{\mathbf{g}}_T = \mathbf{w}_T^T \ddot{\mathbf{q}} + \bar{\mathbf{w}}_T. \quad (51) $$
The terms in Eq. (48), which are linearly dependent on $\ddot{q}$, are collected in the constraint vectors $w_N$, $w_S$ and $w_T$, all the rest is included in the scalars $\bar{w}_N$, $\bar{w}_S$, $\bar{w}_T$.

## 2.3 Multiple unilateral contacts

Multiple contacts in multibody systems include a combinatorial problem of large dimensions. If the state in one contact changes, for example from contact to detachment or from slip to stick, all other contacts are also influenced, which makes a search for a new set of contact configurations necessary. In order to not pursue the combinatorial process, we need extended contact laws which describe unambiguously the transitions for the possible contact states, and which generate only consistent contact configurations. In a first step, we define all contact sets, which can be found in a multibody system, [31]

$$ 
\begin{aligned}
I_A &= \{1, 2, \dots, n_A\} && \text{with } n_A \text{ elements,} \\
I_C &= \{i \in I_A : g_{Ni} = 0\} && \text{with } n_C \text{ elements,} \\
I_N &= \{i \in I_C : \dot{g}_{Ni} = 0\} && \text{with } n_N \text{ elements,} \\
I_T &= \{i \in I_N : |\dot{g}_{Ti}| = 0\} && \text{with } n_T \text{ elements.}
\end{aligned} 
\qquad (52) $$

These sets describe the kinematic state of each contact point. The set $I_A$ consists of $n_A$ indices of all contact points. As an example we consider Eq. (9). It belongs to the set combination $I_A \setminus \{I_C, I_N, I_T\}$. The elements of the set $I_C$ are $n_C$ indices of the unilateral constraints with vanishing normal distance $g_{Ni} = 0$, but an arbitrary relative velocity in the normal direction. In the index set $I_N$ there are $n_N$ indices of the potentially active normal constraints, which fulfill the necessary conditions for a continuous contact (vanishing normal distance $g_{Ni} = 0$, and no relative velocity $\dot{g}_{Ni}$ in the normal direction). The index set $I_N$ includes, for example, all contact states with slipping. The $n_T$ elements of the set $I_T$ are the indices of the potentially active tangential constraints. The corresponding normal constraints are closed, and the relative velocities $\dot{g}_{Ti}$ in the tangential direction are zero. The numbers of elements of the index sets $I_C, I_N$ and $I_T$ are not constant because there are variable states of constraints due to separation and stick-slip phenomena.

As a next step, we must organize all transitions from contact to detachment and from stick to slip and the corresponding reversed transitions. In the normal direction of a contact we find the following situation, [31]:

* passive contact $i$

$$ g_{Ni}(\mathbf{g}, t) \geq 0, \quad \lambda_{Ni} = 0, \quad \text{indicator } g_{Ni}, $$

* transition to contact

$$ g_{Ni}(\mathbf{g}, t) = 0, \quad \lambda_{Ni} \geq 0, \qquad (53) $$

* active contact $i$

$$ g_{Ni}(\mathbf{g}, t) = 0, \quad \lambda_{Ni} > 0, \quad \text{indicator } \lambda_{Ni}, \quad \text{constraint } g_{Ni} = 0, $$

* transition to detachment

$$ g_{Ni}(\mathbf{q}, t) \geq 0, \quad \lambda_{Ni} = 0. $$

The kinematical magnitudes $g_{Ni}, \dot{g}_{Ni}, \ddot{g}_{Ni}$ are given with Eqs. (28), (31), (33) and (41), (42), (51), where $\dot{g}_{Ni}, \ddot{g}_{Ni}$ are needed when we go to a velocity or an acceleration level. The constraint forces $\lambda_{Ni}$ must be compressive forces. If they change sign, we get separation. The properties Eq. (53) establish a complementarity behavior which may be expressed by $n_N$ (set $I_N$) complementarity conditions (put on an acceleration level)

$$ \ddot{g}_N \geq 0; \quad \lambda_N \geq 0; \quad \ddot{g}_N^T \lambda_N = 0. \qquad (54) $$

The variational inequality

$$ -\ddot{g}_N^T (\lambda_N^* - \lambda_N) \leq 0; \quad \lambda_N \in C_N; \quad \forall \lambda_N^* \in C_N, \qquad (55) $$

is equivalent to the complementary conditions (54). The convex set
$$C_N = \{\lambda_N^* : \lambda_N^* \ge 0\} \qquad (56)$$

contains all admissible contact forces $\lambda_N^*$ in the normal direction, [31], [41].

The complementarity problem defined in Eq. (54) may be interpreted as a corner law, which requires for each contact $\ddot{g}_{Ni} \ge 0$, $\lambda_{Ni} \ge 0$, $\ddot{g}_{Ni}\lambda_{Ni} = 0$. Figure 5 illustrates this property.

With respect to the tangential direction of a contact, we shall confine our considerations to the application of Coulomb's friction law, which in no way means a loss of generality. The complementary behavior is a characteristic feature of all contact phenomena, independent of the specific physical law of contact. Furthermore, we assume that within the infinitesimal small time step for a transition from stick to slip, and vice versa, the coefficients of static and sliding friction are the same, which may be expressed by

$$\lim_{\dot{g}_{Ti} \to 0} \mu_i(\dot{g}_{Ti}) = \mu_{0i}. \qquad (57)$$

For $\dot{g}_{Ti} \ne 0$, any friction law may be applied, see Fig. 6. With this property, Coulomb's friction law distinguishes between the two cases

$$\text{stiction: } |\lambda_{Ti}| < \mu_{0i}\lambda_{Ni} \Rightarrow |\mathbf{g}_{Ti}| = 0 \quad (\text{set } I_T), \qquad (58)$$

$$\text{sliding: } |\lambda_{Ti}| = \mu_{0i}\lambda_{Ni} \Rightarrow |\dot{\mathbf{g}}_{Ti}| > 0 \quad (\text{set } I_N \setminus I_T).$$

Equation (58) formulates the mechanical property that we are for a frictional contact within the friction cone, if the relative tangential velocity is zero, and the tangential constraint force $|\lambda_{Ti}|$ is smaller than the maximum static friction force ($\mu_{0i}\lambda_{Ni}$). Then we have stiction. We are on the friction cone, if we slide with $|\dot{\mathbf{g}}_{Ti}| > 0$. At a transition point, the friction force is ($\mu_{0i}\lambda_{Ni}$), see Eq. (57). In addition, we must regard the fact that in the tangential contact plane we may get one or two directions, according to a plane or a spatial contact. From this we summarize in the following way:

* passive contact i (sliding, set $I_N \setminus I_T$)

$$|\dot{\mathbf{g}}_{Ti}| \ge 0, \quad |\mu_{0i}\lambda_{Ni}| - |\lambda_{Ti}| = 0, \quad \text{indicator } |\mathbf{g}_{Ti}|,$$

* transition slip to stick

$$|\dot{\mathbf{g}}_{Ti}| = 0, \quad |\mu_{0i}\lambda_{Ni}| - |\lambda_{Ti}| \ge 0,$$

<figure><img src="image_7.png" /><figcaption>Fig. 5. Corner law for normal contacts (Signorini's law)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 6. Typical friction characteristic</figcaption></figure>
* active contact i (sticking, set $I_T$)

$$|\dot{\mathbf{g}}_{Ti}| = 0, |\mu_{0i}\lambda_{Ni}| - |\lambda_{Ti}| > 0, \text{ indicator } |\mu_{0i}\lambda_{Ni}| - |\lambda_{Ti}|, \text{ constraint } |\dot{\mathbf{g}}_{Ti}| = 0,$$

* transition stick to slip

$$|\dot{\mathbf{g}}_{Ti}| \ge 0, |\mu_{0i}\lambda_{Ni}| - |\lambda_{Ti}| = 0. \qquad (59)$$

From a numerical standpoint of view, we have to check the indicator for a change of sign, which then requires a subsequent interpolation. For a transition from stick to slip one must examine the possible development of a nonzero relative tangential acceleration as a start for sliding.

Equation (58) put on an acceleration level can then be written in a more detailed form

$$\begin{gathered} |\lambda_{Ti}| < \mu_{0i}\lambda_{Ni} \wedge \ddot{g}_{Ti} = 0 \quad (i \in I_T \text{ sticking}), \\ \lambda_{Ti} = +\mu_{0i}\lambda_{Ni} \wedge \ddot{g}_{Ti} \le 0 \quad (i \in I_N \setminus I_T \text{ negative sliding}), \\ \lambda_{Ti} = -\mu_{0i}\lambda_{Ni} \wedge \ddot{g}_{Ti} \ge 0 \quad (i \in I_N \setminus I_T \text{ positive sliding}). \end{gathered} \qquad (60)$$

This contact law may be represented by a double corner law as indicated in Fig. 7. To transform the law (60) for tangential constraints into a complementarity condition, we must decompose the double corner into single ones. A decomposition into four elementary laws is given in [8], a decomposition into two elements in [37], [38]. In any case, we come out with a complementarity problem of the form, [31],

$$\mathbf{y} = \mathbf{A}\mathbf{x} + \mathbf{b}, \quad \mathbf{y} \ge 0, \quad \mathbf{x} \ge 0, \quad \mathbf{y}^T\mathbf{x} = \mathbf{0}, \quad \mathbf{y}, \mathbf{x} \in \mathbb{R}^{n*}, \qquad (61)$$

where $n^* = n_N + 4n_T$ in the case of decomposition into four and $n^* = n_N + 2n_T$ for a decomposition into two elementary corners. The quantity $\mathbf{x}$ includes the contact forces and one part of the decomposed accelerations, the quantity $\mathbf{y}$ the relative accelerations and, in addition, the friction saturation defined as the difference of static friction force and tangential constraint force ($\mu_{0i}\lambda_{Ni} - |\lambda_{Ti}|$). Equation (61) describes a linear complementarity problem thus being adequate for plane contacts. In the case of spatial contacts, the friction saturation contains the geometric sum of two possible friction directions, leading to a nonlinearity which cannot be solved in a straightforward way. Solution procedures are discussed in [11], [12], [38], [42], [44].

Similar as in the normal case, we can represent the contact law Eq. (60) by a variational inequality of the form, [8], [42],

$$\ddot{\mathbf{g}}_{Ti}^T (\lambda_{Ti}^* - \lambda_{Ti}) \ge 0; \quad \lambda_{Ti} \in C_{Ti}, \quad \forall \lambda_{Ti}^* \in C_{Ti}. \qquad (62)$$

<figure><img src="image_9.png" /><figcaption>Fig. 7. Corner law for tangential constraints</figcaption></figure>
The convex set $C_{Ti}$ contains all admissible contact forces $\lambda_{Ti}^*$ in tangential direction

$$ C_{Ti} = \{ \lambda_{Ti}^* || \lambda_{Ti} | \le \mu_0 \lambda_{Ni}; \forall i \in I_T \}. \qquad (63) $$

For a derivation of the equations of motion including unilateral effects, we must combine the multibody equations (9) with the unilateral constraints (33) or (51) and (60). In a first step, we include the constraint forces into Eq. (9) keeping in mind that, in a system with additional unilateral constraints, the number of degrees of freedom is variable. To avoid difficulties with many different sets of minimal coordinates, we take one set of generalized coordinates (Eq. (9) for the combinations $(I_A \setminus I_c)$, $(I_A \setminus I_N)$, $(I_A \setminus I_T)$, see Eq. (52)), and consider the active unilateral constraints as additional constraints which necessarily are accompanied by constraint forces. We include these constraint forces, which are in fact the contact forces, into the equations of motion (9) by a Lagrange multiplier technique.

The constraint vectors $w_{Ni}$ and $w_{Ti}$ in Eqs. (33), (51) are arranged as columns in the constraint matrices

$$ \mathbf{W}_N = [\dots, \mathbf{w}_{Ni}, \dots] \in \mathbb{R}^{f, n_N}; \quad i \in I_N, \\ \mathbf{W}_T = [\dots, \mathbf{w}_{Ti}, \dots] \in \mathbb{R}^{f, 2n_T}; \quad i \in I_T, \qquad (64) $$

for all active constraints. The constraint matrices are transformation matrices from the space of constraints to the configuration space. The transposed matrices are used for the transition from the configuration space to the space of constraints. The contact forces have the amounts $\lambda_{Ni}$ (normal forces) and the components $\lambda_{Ti1}$ and $\lambda_{Ti2}$ (tangential forces). These elements are combined in the vectors of constraint forces

$$ \lambda_N(t) = \begin{pmatrix} \vdots \\ \lambda_{Ni}(t) \\ \vdots \end{pmatrix} \in \mathbb{R}^{n_N}; \quad i \in I_N, \\ \lambda_T(t) = \begin{pmatrix} \vdots \\ \lambda_{Ti}(t) \\ \vdots \end{pmatrix} \in \mathbb{R}^{2n_T}; \quad i \in I_T, \qquad (65) $$

with $\lambda_{Ti}(t) = [\lambda_{Ti1}(t), \lambda_{Ti2}(t)]^T$. In general, the contact forces are time-varying quantities. By the constraint vectors and matrices in Eq. (51), the contact forces can be expressed in the configuration space. These forces are then added to Eq. (9) to give

$$ \mathbf{M}\ddot{\mathbf{q}} - \mathbf{h} - \sum_{i \in I_N} (\mathbf{W}_{Ni}\lambda_N + \mathbf{W}_{Ti}\lambda_{Ti}) = 0. \qquad (66) $$

For the index sets, see Eq. (52). The contact forces $\lambda_{Ti}$ in Eq. (66) can be passive forces of sticking contacts or active forces of sliding contacts. We express the tangential forces of the $n_N - n_T$ sliding contacts by the corresponding normal forces, using Coulomb's friction law, by

$$ \lambda_{Ti} = -\mu_i (|\dot{\mathbf{g}}_{Ti}|) \frac{\dot{\mathbf{g}}_{Ti}}{|\dot{\mathbf{g}}_{Ti}|} \lambda_{Ni}; \quad i \in I_N \setminus I_T, \qquad (67) $$

where the coefficients $\mu_i(|\dot{\mathbf{g}}_{Ti}|)$ of sliding friction may depend on time. The negative sign relates to the opposite direction of relative velocity and friction force. The sliding forces of Eq. (66) in the configuration space are then

$$ \mathbf{W}_{Ti} \lambda_{Ti} = -\mu_i (|\dot{\mathbf{g}}_{Ti}|) \mathbf{W}_{Ti} \frac{\dot{\mathbf{g}}_{Ti}}{|\dot{\mathbf{g}}_{Ti}|} \lambda_{Ni}; \quad i \in I_N \setminus I_T. \qquad (68) $$
A substitution of these forces into Eq. (66) yields the equations of motion

$$ M(\mathbf{q}, t)\ddot{\mathbf{q}}(t) - \mathbf{h}(\mathbf{q}, \dot{\mathbf{q}}, t) - [\mathbf{W}_N + \mathbf{H}_R, \mathbf{W}_T] \begin{pmatrix} \lambda_N(t) \\ \lambda_T(t) \end{pmatrix} = 0, \qquad (69) $$

with the additional contact forces as Lagrange multipliers. The matrices $W_N$ and $W_T$ contain components from Eqs. (33) or (51). The matrix $\mathbf{H}_R \in \mathbb{R}^{f, n_N}$ of the sliding contacts has the same dimension as the constraint matrix $W_N$. For $n_T \le n_N$, $\mathbf{H}_R$ consists of $n_N - n_T$ columns

$$ -\mu_i \mathbf{W}_{Ti} \frac{\dot{\mathbf{g}}_{Ti}}{|\dot{\mathbf{g}}_{Ti}|}; \quad i \in I_N \setminus I_T, $$

while the other $n_T$ columns contain only zero elements.

The relative accelerations of the active normal and tangential constraints in Eqs. (33) or (51) can be combined by means of the constraint matrices (64) in the matrix notation. Together with Eq. (69), we get the system of equations

$$ \begin{aligned} & \mathbf{M}\ddot{\mathbf{q}} - \mathbf{h} - [\mathbf{W}_N + \mathbf{H}_R, \mathbf{W}_T] \begin{pmatrix} \lambda_N \\ \lambda_T \end{pmatrix} = 0 \in \mathbb{R}^f, \\ & \ddot{\mathbf{g}}_N = \mathbf{W}_N^T\ddot{\mathbf{q}} + \bar{\mathbf{w}}_N \in \mathbb{R}^{n_N}, \\ & \ddot{\mathbf{g}}_T = \mathbf{W}_T^T\ddot{\mathbf{q}} + \bar{\mathbf{w}}_T \in \mathbb{R}^{n_T}. \end{aligned} \qquad (70) $$

The unknown quantities are the generalized accelerations $\ddot{\mathbf{q}} \in \mathbb{R}^f$, the contact forces in the normal direction $\lambda_N \in \mathbb{R}^{n_N}$ and in the tangential direction $\lambda_T \in \mathbb{R}^{2n_T}$, as well as the corresponding relative accelerations $\ddot{\mathbf{g}}_N \in \mathbb{R}^{n_N}$ and $\ddot{\mathbf{g}}_T \in \mathbb{R}^{2n_T}$. For the determination of the $f + 2(n_N + 2n_T)$ quantities, we have up to now $f + n_N + 2n_T$ equations. In the following, Eqs. (70) will be completed by including the missing $n_N + 2n_T$ contact laws. In general, the kinematic equations are dependent on each other, if there is more than one contact point per rigid body. The situation results in linearly dependent columns of the constraint matrices $\mathbf{W}_N, \mathbf{W}_T$ in Eqs. (70). Such constraints are called dependent constraints.

Before including the contact laws we shall evaluate Eqs. (70) a bit further. Introducing the magnitudes

$$ \ddot{\mathbf{g}} = \begin{pmatrix} \ddot{\mathbf{g}}_N \\ \ddot{\mathbf{g}}_T \end{pmatrix}; \quad \lambda = \begin{pmatrix} \lambda_N \\ \lambda_T \end{pmatrix}; \quad \bar{\mathbf{w}} = \begin{pmatrix} \bar{\mathbf{w}}_N \\ \bar{\mathbf{w}}_T \end{pmatrix}; \mathbf{W} = (\mathbf{W}_N, \mathbf{W}_T); \quad \mathbf{N}_G = (\mathbf{H}_R, 0), \qquad (71) $$

we can write Eq. (70) as

$$ \mathbf{M}\ddot{\mathbf{g}} - \mathbf{h} - (\mathbf{W} + \mathbf{N}_G)\lambda = 0, \quad \ddot{\mathbf{g}} = \mathbf{W}^T\ddot{\mathbf{q}} + \bar{\mathbf{w}}. \qquad (72) $$

Outside the transition events and, thus, for a not changing contact configuration, the relative accelerations $\ddot{\mathbf{g}}$ are zero. In this case, Eqs. (72) have a solution for $\ddot{\mathbf{g}}$ and $\lambda$, which writes

$$ \lambda = -[\mathbf{W}^T \mathbf{M}^{-1}(\mathbf{W} + \mathbf{N}_G)]^{-1}(\mathbf{W}^T \mathbf{M}^{-1}\mathbf{h} + \bar{\mathbf{w}}), \quad \ddot{\mathbf{q}} = \mathbf{M}^{-1}[\mathbf{h} + (\mathbf{W} + \mathbf{N}_G)\lambda]. \qquad (73) $$

The first equation of (73) may also be expressed in the form $A\lambda + \mathbf{b} = 0$ with A and b following from (73).

We combine now the equations of motion in the form (70) with the contact Eqs. (55) and (62), which results in the following set containing all unilateral processes in the contacts under consideration:

$$ \begin{aligned} & \mathbf{M}\ddot{\mathbf{q}} - [\mathbf{W}_N + \mathbf{H}_R, \mathbf{W}_T] \begin{pmatrix} \lambda_N \\ \lambda_T \end{pmatrix} - \mathbf{h}(\mathbf{q}, \dot{\mathbf{q}}, t) = 0, \\ & \begin{pmatrix} \ddot{\mathbf{g}}_N \\ \ddot{\mathbf{g}}_T \end{pmatrix} = \begin{pmatrix} \mathbf{W}_N^T \\ \mathbf{W}_T^T \end{pmatrix} \ddot{\mathbf{q}} + \begin{pmatrix} \bar{\mathbf{w}}_N \\ \bar{\mathbf{w}}_T \end{pmatrix}, \qquad (74) \\ & -\ddot{\mathbf{g}}_N^T(\lambda_N^* - \lambda_N) \le 0; \quad \lambda_N \in C_N; \quad \forall \lambda_N^* \in C_N; \quad I \in I_N, \\ & -\ddot{\mathbf{g}}_{Ti}^T(\lambda_{Ti}^* - \lambda_{Ti}) \le 0; \quad \lambda_{Ti} \in C_{Ti}(\lambda_{Ni}); \quad \forall \lambda_{Ti}^* \in C_{Ti}(\lambda_{Ni}); \quad i \in I_T. \end{aligned} $$
The system (74) is not solvable in this form. Therefore the variational inequalities are transformed into equalities. From this we get a nonlinear system of equations, which represent linear or nonlinear complementarity problems, depending on the type of contact, i.e. plane or spatial contacts. Linear complementarity problem can be solved by algorithms related to linear programming methods, for example Lemke's algorithm, nonlinear complementarity problems require iterative algorithms, [11], [12], [24].

## 2.4 Impact with friction

Impacts with and without friction play an important role in many machines and mechanisms. They have been a research topic at the author's institute since many years, starting with impacts without friction, [29], and proceeding to a theory which describes impacts with friction in a rather general way, [8]. The new theory has been convincingly verified by a large variety of experiments, [2], and will be shortly presented in the following.

We assume, as usual, that an impact takes place in an infinitesimal short time and without any change of position, orientation and all nonimpulsive forces. Nevertheless, we zoom virtually the impact time, establish the equations for a compression and for an expansion phase, and apply these equations again to an infinitesimal short time interval. The evaluation has to be performed on a velocity level, which we realize by formal integration of the equations of motion, the constraints and the contact laws, Eqs. (74). Denoting the beginning of an impact, the end of compression and the end of expansion by the indices A, C, E, respectively, we get for $\Delta t = t_E - t_A$

$$M(\dot{\mathbf{q}}_C - \dot{\mathbf{q}}_A) - (\mathbf{W}_N \mathbf{W}_T) \begin{pmatrix} \mathbf{\Lambda}_{NC} \\ \mathbf{\Lambda}_{TC} \end{pmatrix} = 0,$$

$$M(\dot{\mathbf{q}}_E - \dot{\mathbf{q}}_C) - (\mathbf{W}_N \mathbf{W}_T) \begin{pmatrix} \mathbf{\Lambda}_{NE} \\ \mathbf{\Lambda}_{TE} \end{pmatrix} = 0, \qquad (75)$$

$$\mathbf{\Lambda}_i = \lim_{t_E \to t_A} \int_{t_A}^{t_E} \lambda_i dt .$$

Here, $\mathbf{\Lambda}_{NC}$, $\mathbf{\Lambda}_{TC}$ are the impulses in the normal and tangential direction which are transferred during compression, and $\mathbf{\Lambda}_{NE}$, $\mathbf{\Lambda}_{TE}$ those of expansion. Defining $\dot{\mathbf{q}}_A = \dot{\mathbf{q}}(t_A)$, $\dot{\mathbf{q}}_C = \dot{\mathbf{q}}(t_C)$, $\dot{\mathbf{q}}_E = \dot{\mathbf{q}}(t_E)$, we express the relative velocities as

$$\begin{pmatrix} \dot{\mathbf{g}}_{NA} \\ \dot{\mathbf{g}}_{TA} \end{pmatrix} = \begin{pmatrix} \mathbf{W}_N^T \\ \mathbf{W}_T^T \end{pmatrix} \dot{\mathbf{q}}_A + \begin{pmatrix} \tilde{\mathbf{w}}_N \\ \tilde{\mathbf{w}}_T \end{pmatrix},$$

$$\begin{pmatrix} \dot{\mathbf{g}}_{NE} \\ \dot{\mathbf{g}}_{TE} \end{pmatrix} = \begin{pmatrix} \mathbf{W}_N^T \\ \mathbf{W}_T^T \end{pmatrix} \dot{\mathbf{q}}_E + \begin{pmatrix} \tilde{\mathbf{w}}_N \\ \tilde{\mathbf{w}}_T \end{pmatrix}, \qquad (76)$$

$$\begin{pmatrix} \dot{\mathbf{g}}_{NC} \\ \dot{\mathbf{g}}_{TC} \end{pmatrix} = \begin{pmatrix} \mathbf{W}_N^T \\ \mathbf{W}_T^T \end{pmatrix} \dot{\mathbf{q}}_C + \begin{pmatrix} \tilde{\mathbf{w}}_N \\ \tilde{\mathbf{w}}_T \end{pmatrix}.$$

Considering in a first step the compression phase, and combining Eqs. (75) and (76), we come out with

$$\begin{pmatrix} \dot{\mathbf{g}}_{NC} \\ \dot{\mathbf{g}}_{TC} \end{pmatrix} = \underbrace{\begin{pmatrix} \mathbf{W}_N^T \\ \mathbf{W}_T^T \end{pmatrix}}_{G} \cdot \mathbf{M}^{-1} \begin{pmatrix} \mathbf{W}_N \\ \mathbf{W}_T \end{pmatrix}^T \cdot \begin{pmatrix} \mathbf{\Lambda}_{NC} \\ \mathbf{\Lambda}_{TC} \end{pmatrix} + \begin{pmatrix} \dot{\mathbf{g}}_{NA} \\ \dot{\mathbf{g}}_{TA} \end{pmatrix}, \qquad (77)$$

where G is called the matrix of projected inertia. It consists of four blocks $G_{NN}...G_{TT}$. Equation (77) allows to calculate the relative velocities $\dot{\mathbf{g}}_{NC}$ and $\dot{\mathbf{g}}_{TC}$ at the end of the compression phase, depending on the velocities at the beginning of the impact $\dot{\mathbf{g}}_{NA}$ and $\dot{\mathbf{g}}_{TA}$ under the influence of the contact impulses $\mathbf{\Lambda}_{NC}$ and $\mathbf{\Lambda}_{TC}$. To calculate these impulses, two impact laws in normal and tangential direction are necessary. As already indicated, magnitudes of relative kinematics and constraint forces (here: impulses) are complementary quantities. In the normal direction these are $\dot{\mathbf{g}}_{NC}$ and $\mathbf{\Lambda}_{NC}$. In the tangential direction, we have the relative tangential velocity vector $\dot{\mathbf{g}}_{TC}$ and the friction saturation $\mathbf{\Lambda}_{TC} - (\mathrm{diag}\mu_i)\mathbf{\Lambda}_{NC}$. Decomposing the tangential behavior we obtain, [2],
$$ \Lambda_{TCV,i} = \Lambda_{TC,i} + \mu_i \Lambda_{TN,i}, $$

$$ \dot{\mathbf{g}}_{TC,i} = \dot{\mathbf{g}}_{TC,i}^{+} - \dot{\mathbf{g}}_{TC,i}^{-}, $$

$$ \Lambda_{TCV,i}^{(+)} = \Lambda_{TCV,i}, \qquad (78) $$

$$ \Lambda_{TCV,i}^{(-)} = -\Lambda_{TCV,i} + 2\mu_i \Lambda_{NC,i}. $$

Together with Eq. (77) this results in a Linear Complementary Problem (LCP) in a standard form $\mathbf{y} = \mathbf{A}\mathbf{x} + \mathbf{b}$ with $\mathbf{x} \ge 0, \mathbf{y} \ge 0$ and $\mathbf{x}^T\mathbf{y} = 0$

$$ \underbrace{\left( \begin{array}{c} \dot{\mathbf{g}}_{NC} \\ \dot{\mathbf{g}}_{TC}^{+} \\ \mathbf{\Lambda}_{TCV}^{(-)} \end{array} \right)}_{\mathbf{y}} = \underbrace{\left( \begin{array}{ccc} \mathbf{G}_{NN} - \mathbf{G}_{NT}\mathbf{\mu} & \mathbf{G}_{NT} & 0 \\ \mathbf{G}_{TN} - \mathbf{G}_{TT}\mathbf{\mu} & \mathbf{G}_{TT} & \mathbf{E} \\ 2\mathbf{\mu} & -\mathbf{E} & 0 \end{array} \right)}_{\mathbf{A}} \underbrace{\left( \begin{array}{c} \mathbf{\Lambda}_{NC} \\ \mathbf{\Lambda}_{TCV}^{(+)} \\ \dot{\mathbf{g}}_{TC}^{-} \end{array} \right)}_{\mathbf{x}} + \underbrace{\left( \begin{array}{c} \dot{\mathbf{g}}_{NA} \\ \dot{\mathbf{g}}_{TA} \\ 0 \end{array} \right)}_{\mathbf{b}}, \qquad (79) $$

where $\mathbf{\mu}$ is a diagonal matrix, containing the friction coefficients of the contacts. The problem can be solved numerically. The velocities $\dot{\mathbf{g}}_{NC}, \dot{\mathbf{g}}_{TC}$ and the impulses $\Lambda_{NC}, \Lambda_{TC}$ are either part of the result or can be obtained by transformation (78) and by $\Lambda_{TC} = \Lambda_{TCV}^{(+)} - \mathbf{\mu}\Lambda_{NC}$.

In the compression phase, the kinetic energy of the colliding bodies is transformed into potential energy. During the expansion phase parts of this energy are transformed back. This recovery process is governed by two coefficients of restitution in normal and tangential directions $\varepsilon_N$ and $\varepsilon_T$, respectively.

In the case of multiple impacts, Poisson's hypothesis does not guarantee impenetrability of bodies, [8]. Therefore, the impact law is enhanced by an additional condition. The normal impulse during the phase of restitution $\Lambda_{NE}$ is, at the minimum, $\varepsilon_E \Lambda_{NC}$ in each contact, but can be arbitrary high to avoid penetration. In this case, the bodies remain in contact after the impact. This impact law with a complementary character is drawn in Fig. 8.

In the tangential direction, the impact law is ruled by three effects: at first, a minimum impulse $\varepsilon_T \Lambda_{TC}$ must be transferred. On the other hand, the impulse must not exceed the friction limit. Between these limits, a tangential relative velocity $\dot{\mathbf{g}}_{TE0}$ is prescribed. This velocity allows a restitution of the stored energy, if during the phase of restitution sticking occurs, [2]. The existence of $\dot{\mathbf{g}}_{TE0}$ was not introduced in [8], and also not in Moreau's impact model [10], but from the measurements, [2], it could be derived as a necessary correction. With

$$ \dot{\mathbf{g}}_{TE0} = \mathbf{G}_{TN}\varepsilon_N\mathbf{\Lambda}_{NC} + \mathbf{G}_{TT}\varepsilon_N\mathbf{\Lambda}_{TC}, \qquad (80) $$

one can calculate $\dot{\mathbf{g}}_{TE0}$ for all contacts, and the tangential impact law looks like the one drawn in Fig. 9. Here, $\varepsilon_N$ and $\varepsilon_T$ are diagonal matrices containing the different coefficients for all contacts.

To formulate the equation for the phase of restitution as an LCP, similar to that for the compression phase, the two matrices

$$ \mathbf{S}^{+} = \mathrm{diag}\{\frac{1}{2}[1 + \mathrm{sign}(\mathbf{\Lambda}_{TC,i})]\}, \quad \mathbf{S}^{-} = \mathrm{diag}\{\frac{1}{2}[1 - \mathrm{sign}(\mathbf{\Lambda}_{TC,i})]\}, \qquad (81) $$

are introduced. After some transformations similar to those of the compression phase, the LCP writes

<figure><img src="image_10.png" /><figcaption>Fig. 8. Impact law for the phase of expansion in normal direction</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 9. Impact law for the phase of restitution in tangential direction</figcaption></figure>

$$ \begin{pmatrix} \dot{\mathbf{g}}_{NE} \\ \dot{\mathbf{g}}_{TEV}^{+} \\ \mathbf{\Lambda}_{TEV}^{(-)} \end{pmatrix} = \begin{pmatrix} \mathbf{G}_{NN} - \mathbf{G}_{NT}\mathbf{S}^{-}\mathbf{\mu} & \mathbf{G}_{NT} & 0 \\ \mathbf{G}_{TN} - \mathbf{G}_{TT}\mathbf{S}^{-}\mathbf{\mu} & \mathbf{G}_{TT} & \mathbf{E} \\ \mathbf{\mu} & -\mathbf{E} & 0 \end{pmatrix} \begin{pmatrix} \mathbf{\Lambda}_{NP} \\ \mathbf{\Lambda}_{TEV}^{(+)} \\ \dot{\mathbf{g}}_{TEV}^{-} \end{pmatrix} + \mathbf{b}, \qquad (82) $$

with

$$ \mathbf{b} = \begin{pmatrix} \mathbf{G}_{NN}\varepsilon_N\mathbf{\Lambda}_{NC} + \mathbf{G}_{NT}\mathbf{S}^{+}\varepsilon_N\varepsilon_T\mathbf{\Lambda}_{TC} - \mathbf{G}_{NT}\mathbf{S}^{-}\mathbf{\mu}\varepsilon_N\mathbf{\Lambda}_{NC} + \dot{\mathbf{g}}_{NC} \\ \mathbf{G}_{TT}(\mathbf{S}^{+} - \mathbf{E})\varepsilon_N\varepsilon_T\mathbf{\Lambda}_{TC} - \mathbf{G}_{TT}\mathbf{S}^{-}\mathbf{\mu}\varepsilon_N\mathbf{\Lambda}_{NC} + \dot{\mathbf{g}}_{TC} \\ \mathbf{\mu}\varepsilon_N\mathbf{\Lambda}_{NC} - \varepsilon_N\varepsilon_T|\mathbf{\Lambda}_{TC}| \end{pmatrix}. \qquad (83) $$

After solution, the velocities $\dot{\mathbf{g}}_{NC}$, $\dot{\mathbf{g}}_{TC}$ and the impulses $\mathbf{\Lambda}_{NC}$, $\mathbf{\Lambda}_{TC}$ are either part of the result or can be obtained by the transformations

$$ \dot{\mathbf{g}}_{TE} = \dot{\mathbf{g}}_{TEV}^{+} - \dot{\mathbf{g}}_{TEV}^{-} + \dot{\mathbf{g}}_{TE0}, \qquad (84) $$

$$ \mathbf{\Lambda}_{NE} = \mathbf{\Lambda}_{NP} + \varepsilon_N \mathbf{\Lambda}_{NC}, \quad \mathbf{\Lambda}_{TE} = \mathbf{\Lambda}_{TEV}^{(+)} + \mathbf{\Lambda}_{TEL} = \mathbf{\Lambda}_{TEV}^{(+)} + \mathbf{S}^{+}\varepsilon_N\varepsilon_T\mathbf{\Lambda}_{TC} - \mathbf{S}^{-}\mathbf{\mu}\mathbf{\Lambda}_{NE}. \qquad (85) $$

If the impulses during the two phases of the impact are known, one can calculate the motion $\dot{\mathbf{q}}_E$ of the multibody system at the end of the impact

$$ \dot{\mathbf{q}}_E = \dot{\mathbf{q}}_A \mathbf{M}^{-1} [\mathbf{W}_N(\mathbf{\Lambda}_{NC} + \mathbf{\Lambda}_{NE}) + \mathbf{W}_T(\mathbf{\Lambda}_{TC} + \mathbf{\Lambda}_{TE})] \qquad (86) $$

The above impact theory has been verified in [2]. More than 600 impact measurements have been performed. Figure 10 depicts a typical result for the material pair PVC/PVC. The comparison with measurements is excellent, as in many other cases.

## 3 Applications

From a large variety of industrial applications, we shall present in the following three typical examples, namely chains, vibration conveyors and chimney dampers. We do not go into details, because all applications have been published elsewhere. Therefore we shall give only a short summary for each case. In all cases, the theoretical methods are based on those of Sec. 2.

<figure><img src="image_12.png" /><figcaption>Fig. 10. Dimensionless tangential relative velocity after versus before the impact</figcaption></figure>
## 3.1 CVT chains, [40]

The dynamics of chains has been investigated in [6], [7], [39], [40]. The dissertation [7] considers roller chains, where the main problem lies in the many degrees of freedom, usually some hundred, on one side, from which follows a large number of contacts with guides and sprocket wheels, on the other side. An additional problem arises due to the dynamical behavior of the tension devices which requires detailed and complicated models.

Another type of chains are the CVT chains applied in continuous variable transmissions, where they transmit power by friction. At the time being, there exist three different chain configurations, the rocker pin chains by PIV in Germany, the metal-pushed V-belt chain by van Doorne in the Netherlands and the Borg-Warner chain in the US. We shall focus here on the PIV chain, which was investigated in all details in [40]. Figure 11 represents a typical CVT configuration. A complete CVT gear consists of two discs with hydraulically shifted sheaves as indicated in Fig. 11. One can change the gear ratio continuously, by opening the discs in one pulley and, at the same time, closing it in the second pulley. Therefore, the continuous variable properties are achieved by axially moving sheaves that are shifted hydraulically. This device also applies to the chain the necessary forces.

The mechanical model for the plane motion of the chain includes altogether 63 links with 189 degrees of freedom. Each pulley possesses two translational and one rotational degree of freedom with respect to rigid body motion. Pitching effects of the axially movable pulley are regarded by additional degrees of freedom. Both pulleys are additionally modeled as elastic bodies, applying a Ritz approach within the context of elastic multibody theory, [3], [4], [40]. Special care has been taken to solve the problem of unilateral contacts in connection with elastically deformed surfaces.

The results have been computed for stationary operation with constant driving speed and an output torque on geared level. Computing time for this case amounts to about eight hours on a SUN workstation. Figure 12 shows the contact forces acting on a pair of rocker pins during one revolution and the tensile force in the related chain link.

As long as the chain link is part of a strand, no contact forces work on its rocker pins. When it comes into contact with one of the pulleys, the pins are pressed between the two sheaves, and hence the normal force increases. Its amplitude depends on the geometry of the sheaves and the transmitted power. The frictional force is a function of this normal force and the relative velocity between the pulley and the pins. It is split into one radial and one circumferential component. The radial contact force coincides with a radial movement of the chain link that equals a power dissipation. In contrast, the circumferential contact force causes the changes of the tensile force in the corresponding chain link, which leads to different tensile force levels in the two strands that agree with the transmitted torque.

<figure><img src="image_13.png" /><figcaption>Fig. 11. Typical rocker pin chain configuration</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 12. Forces acting on the rocker pins and in a chain link</figcaption></figure>

Due to the mechanical model, the simulation provides an integrated value of the tensile force in the plates of a chain link, whereas measurements (by G. Sauer and K.T. Renius, Lehrstuhl für Landmaschinen, Technische Universität München) were performed for the tensile force in a clasp plate, s. left-hand side of Fig. 13. Therefore, it was necessary to determine the distribution of the tensile force on the plates of the chain links. Using the results of the dynamic simulation shown in Fig. 12, and modeling the pair of rocker pins as bending beams and the plates as linear springs, we obtained the graph of the clasp plate, s. in the right-hand side of Fig. 13. The comparison of simulation and measurement confirms the mechanical model.

## 3.2

### Vibration conveyor, [44]

Vibratory feeders are used in automatic assembly to feed small parts. They are capable to store, transport, orient and isolate the parts. An oscillating track with frequencies up to 100 Hz excites the transportation process, which is mainly based on impact and friction phenomena between the parts and the track. Vibratory feeders are applied for a wide variety of parts and for lots of different tasks. In the majority of cases, the parts are available as a sort of bulk material that is stored in a container. The transportation process, starting in this reservoir, is often combined with orienting devices that orient parts, or select only the parts having already the right orientation. Figure 14 shows an example of a vibratory bowl feeder with an orienting device. Each kind of parts, with its special geometry and mechanical properties, requires an individual adaption of the feeder. This individual tuning comprises the development of suitable track and orienting device geometries and the adjustment of the excitation

<figure><img src="image_15.png" /><figcaption>Fig. 13. Tensile force in the clasp plate</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 14. Vibratory bowl feeder and mechanical model</figcaption></figure>

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

parameters frequency and amplitude. Due to the complex mechanics of the feeding process this design is usually done by trial and error, without any theoretical background. A complete dynamical model of the transportation process allows a theoretical investigation and, consequently, an improvement in the properties of the feeder.

Friction and impact phenomena between the parts and the track are the most important mechanical properties of transportation processes. Consequently, the required dynamical model has to deal with unilateral constraints, dry friction and multiple impacts. The mechanical model of the vibratory feeder can be split in two parts: the transportation process and the base device. The dissertation [44] focuses on modeling and simulation of the transportation process. The modeling of the base device can be done with well-known standard techniques for multibody systems. Friction and impact effects have a fundamental importance for the transportation of parts. Changing contact configurations between the parts and the track and also between the parts itself are characteristic for the feeding process. The contacts appear either continuously for a certain time interval, or for a discrete time at impact. A structure-varying multibody system with unilateral constraints with friction is an ideal technique for modeling the feeding process. Its formulation results in a set of differential equations with inequality constraints requiring special mathematical and numerical methods, see Sec. 2.

For the verification of the developed model of the transportation process an experimental vibratory feeder was built, allowing different measurements concerning the impact model and the average transportation rates. Figure 15 shows the principle of the device. The track, fixed on leaf springs, is excited with an electro-magnetic shaker with a frequency of about 50 Hz. The eigenfrequency of the system is at 52 Hz. The resulting vibration amplitude reaches a maximum value of about 2 mm. The track has an inclination angle $α = 3°$, the angle between the track and the direction of the vibration is $β = 15°$. For accurate contactless measurement of the motion of the transported part, six laser distance sensors were applied. An eddy current sensor was used for the vibration measurement of the track.

For a comparison of the theory and the measurements, the averaged transportation rate was used. Figure 16 presents the result, which looks good before the background of the complexity of the problem. An interesting finding is the fact that the averaged transportation velocity does not depend very much on the number of parts and also not on the type of modeling, plane or spatial, [44]. Therefore, the design of vibration conveyors can be carried through considering one part only. For the layout of orienting devices we need, of course, a spatial theory.

<figure><img src="image_18.png" /><figcaption>Fig. 15. Test setup and part measurement</figcaption></figure>

<figure><img src="image_19.png" /></figure>
<figure><img src="image_20.png" /><figcaption>Fig. 16. Simulation and measurement of the average transportation rate with $\varepsilon = 0.84$, $\mu = 0.23$</figcaption></figure>

## 3.3

### Chimney dampers, [34]

Towerlike structures, like steel chimneys, may be excited to vibrations by wind vortex streets. Such a mechanism becomes dangerous if the first eigenfrequency of the structure is small enough to be exited by not-too-large wind velocities which appear quite frequently. Oscillations of that type can be significantly damped by application of a simple idea. A pendulum is attached to the chimney top. Its mass and length are tuned to the chimneys first eigenfrequency by classical analysis, and, additionally, optimal viscous damping is also evaluated by classical formulas, given for example by Den Hartog's method. In a second step, the determined optimal damping behavior is not realized by viscous means but by arranging a package of circular plates, which are moved by the end of the pendulum rod through internal circular holes of the plates. The approximate realization of viscous damping by dry friction is accompanied by impacts of the pendulum rod in the plate holes and by stick-slip transitions between the plates. To achieve a best fit to the optimal viscous case, an optimization of the complete pendulum-plate-system is carried through, applying multibody theory with unilateral contacts.

Figure 17 illustrates the basic principle. At the top of the chimney, a pendulum is arranged which damps the oscillations. To achieve best damping efficiency, the damper is optimized in two steps. In a first step, a classical damper is assumed working with viscous damping. For the resonance frequency area such a damper can be optimized with respect to length, mass and viscous damping of the pendulum. As damping will be realized by dry friction within a plate package, the plates must be selected in such a way that they perform a damping effect coming as nearest as possible to the optimal viscous damping. Therefore, in a second step the plate package will be optimized with respect to damping efficiency.

To verify the theory, two types of experiments have been performed. In a first test, a pendulum-plate configuration has been arranged in a car-like frame with wheels, which could be excited periodically, Fig. 18. The results compare well with the theory.

A second test has been performed with a real steel chimney, which was bent by a steel rope. By releasing suddenly this rope, the chimney starts to vibrate. This process has been measured and compared with the corresponding theory based on the equations of Sec. 2. Figure 19

<figure><img src="image_21.png" /><figcaption>Fig. 17. Principal configuration, viscous and plate model, [34]</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>Fig. 18. Vibrating test-setup and a comparison theory/measurements</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 19. Pendulum vibrations in a real chimney (height 80 m, first bending eigenfrequency 0.76 Hz)</figcaption></figure>

illustrates the comparison between the experiment and simulation, which also confirms the method presented.

If we approximate the optimal viscous system of Fig. 17 (middle) by a package of plates, we can do that in a best way by optimizing the number of plates, their radii and thicknesses. We find that the number of plates affects the damping efficiency only partly. On the left side of Fig. 20, the total mass of the plates was kept constant, also the radius of the top at the ground damper plate. The number of plates was varied from one to eight. When doing so, the distribution of the plate holes was kept linear. In the case of one plate, the diameter of the plate hole was half as big as the largest one in all other cases.

Figure 20 shows that the damping mechanism works best with one plate and worst with two plates. The reason lies in the fact, that in the case of two plates only the upper one is moved effectively. The ground plate is ineffective, and hence its damping is too small. As an important result it should be noticed, that the application of a large number of plates does not make sense. Obviously the damping characteristic cannot be improved by using more than five or six plates.

<figure><img src="image_24.png" /><figcaption>Fig. 20. Variation of plate number (left) and of graduation (right) of the plate hole diameters</figcaption></figure>
A further sensitive parameter is the distribution of the hole diameters over all plates, called graduation. The right graph of Fig. 20 shows with number 1 the behavior of the system with small, and with graph number 5-with large graduation. The other curves correspond to some intermediate graduations. Keeping the total mass of the plates constant and increasing graduation from 1 to 5 leads to a decrease of damping and to rising chimney displacements for the considered configuration.

4

**Summary**

Classical theory of rigid or elastic multibody system dynamics may be characterized by d'Alembert's principle in the form of Lagrange and extended by Jourdain principle, allowing to project the equations of motion of all interconnected bodies into their free directions according to their individual constraints. These constraints are bilateral in the classical theory and do not include contact phenomena.

If we want to consider dynamical problems with unilateral contacts, especially multibody systems with multiple contacts, and if we further want to model such contacts by unilateral constraints, we must enlarge the multibody theory by certain rules describing unilateral features. One fundamental property, sometimes addressed to as Signorini's law, consists in the fact, that for each contact either quantities of relative kinematics are zero, and the corresponding constraint forces are not zero, or vice versa. This establishes a linear or nonlinear complementarity problem for each of the contacts, which has to be added to the classical multibody formalism. In a more mathematical sense, the nonsmooth character of contact dynamics may be formulated in terms of variational or hemivariational inequalities, depending on the convex or non-convex features of the related energies. Representations of that kind allow for elegant mathematical developments, which however have again to be decomposed for practical applications.

**References**

1. Baraff, D.: Issues in computing contact forces for non-penetrating rigid bodies. Algorithmica 10 (1993) 292–352

2. Beitelschmidt, M.: Reibstöße in Mehrkörpersystemen. Dissertation TU-München, Lehrstuhl B für Mechanik, 1998

3. Bremer, H.: Dynamik und Regelung mechanischer Systeme. Teubner Studienbücher, Mechanik, Stuttgart: B.G. Teubner 1988

4. Bremer, H.; Pfeiffer, F.: Elastische Mehrkörpersysteme. Teubner Studienbücher, Mechanik Stuttgart: B.G. Teubner 1992

5. Brogliatto, B.: Nonsmooth Impact Dynamics. London: Springer 1996

6. Fritz, P.; Pfeiffer, F.: Dynamics of High Speed Roller Chain Drives. Proceeding of the 1995 Design Engineering Technical Conference, Vol. 3, Part A, DE-Vol. 84-1, pp. 573–584, Boston, 1995

7. Fritz, P.: Dynamik schnellaufender Kettentriebe. Fortschrittberichte VDI, Reihe 11, Nr. 253, Düsseldorf: VDI-Verlag 1998

8. Glocker, Ch.: Dynamik von Starrkörpersystemen mit Reibung und Stößen. Fortschritt-Berichte VDI, Reihe 18: Mechanik/Bruchmechanik, Nr. 182, Düsseldorf: VDI-Verlag 1995

9. Glocker, Ch.; Pfeiffer, F.: Multiple impacts with friction in multibody systems. Nonlinear Dyn 7 (1995) 471–497

10. Jean, M.; Moreau, J. J.: Unilaterality and Dry Friction in the Dynamics of Rigid Body Collections. In: Curnier, A (ed.) Proceedings Contact Mechanics International Symposium. EPFL, PPUR, pp. 31–48, Lausanne, Switzerland 1992

11. Klarbring, A.: Mathematical Programming and Augmented Lagrangian Methods for Frictional Contact Problems. In: Curnier, A. (ed.) Proceedings Contact Mechanics International Symposium, EPFL, PPUR, pp. 409–422, Lausanne, Switzerland 1992

12. Klarbring, A.: Mathematical Programming in Contact Problems. In: Aliabadi, M.H. (ed.) Computational Methods for Contact Problems, Elsevier, 1994

13. Kwak, B. M.: Complementarity problem formulation of three-dimensional frictional contact. ASME J. Appl. Mech. 58 (1991) 134–140

14. Kwak, B. M.; Park, J. K.: Three-dimensional frictional analysis using the homotopy method. ASME J. Appl. Mech. 61 (1994) 703–709

15. Lötstedt, P.: Coulomb friction in two-dimensional rigid body systems. Z. Angew. Math. Mech. 61 (1981) 605–615

16. Lötstedt, P.: Mechanical systems of rigid bodies subject to unilateral constraints. SIAM J. Appl. Math. 42 (1982) 281–296

17. Mangasarian, O. L.: Equivalence of the complementarity problem to a system of nonlinear equations. SIAM J. Appl. Math. 31(1) (1976) 89–92

18. Martins, J. A. C.; Oden, J. T.; Simoes, F. M. F.: A study of static and kinetic friction. Int. J. Eng. Sci., 28(1) (1990) 29–92
19. Meitinger, Th.: Modellierung und Simulation von Montageprozessen. In: Fortschritt-Berichte VDI, Reihe 2, Nr. 476, Düsseldorf: VDI-Verlag 1998

20. Monteiro Marques, M. D. P.: Differential Inclusions in Nonsmooth Mechanical Problems. Basel: Birkhäuser 1993

21. Moreau, J. J.: Application of Convex Analysis to Some Problems of Dry Friction. In: Trends in Applications of Pure Mathematics to Mechanics, Vol. 2, London, 1979

22. Moreau, J. J.: Une formulation du contact a frottement sec; application au calcul numerique. Technical Report 13, C.R. Acad. Sci. Paris, Serie II, 1986

23. Moreau, J. J.: Unilateral Contact and Dry Friction in Finite Freedom Dynamics. In: Non-Smooth Mechanics and Applications, CISM Courses and Lectures, Vol. 302, Wien: Springer Verlag 1988

24. Murty, K. G.: Linear Complementarity, Linear and Nonlinear Programming. In: White, D.J (ed.) Sigma Series in Applied Mathematics. Berlin: Heldermann Verlag 1988

25. Oden, J. T.; Martins, J. A. C.: Models and computational methods for dynamic friction phenomena. Comput. Methods Appl. Mech. Eng. 52 (1985) 527-634

26. Newton, I.: Principia, Corol. VI, 1687

27. Panagiotopoulos, P. D.: Inequality Problems in Mechanics and Applications. Boston, Basel, Stuttgart: Birkhäuser 1985

28. Panagiotopoulos, P. D.: Hemivariational Inequalities. Berlin, Heidelberg: Springer Verlag 1993

29. Pfeiffer, F.: Mechanische Systeme mit unstetigen Übergängen. Ing. Arch. 54 (1984) 232-240

30. Pfeiffer, F.: Multibody Dynamics with Multiple Unilateral Contacts. Proceedings of the XIXth International Congress of Theoretical and Applied Mechanics, Kyoto, 1996

31. Pfeiffer, F.; Glocker, Ch.: Multibody Dynamics with Unilateral Contacts. In: Wiley Series in Nonlinear Science, New York: J. Wiley 1996

32. Pfeiffer, F.: Assembly processes with robotic systems. Robotics and Autonomous Systems 19 (1996) 151-166

33. Pfeiffer, F.; Fritz, P.; Srnik, J.: Nonlinear vibrations of chains. J. Vib. and Control 3 (1997) 397-410

34. Pfeiffer, F.; Stiegelmeyr, A.: Damping Towerlike Structures by Dry Friction. Proc. of DETC '97, ASME Design Eng. Techn. Conf., 1997

35. Poisson, S. D.: Traité de méchanique. Bachelier, Paris, 1833

36. Rockafellar, R. T.: Convex Analysis. Princeton University Press, 1972

37. Rossmann, Th.; Pfeiffer, F.: Efficient Algorithms for Nonsmooth Dynamics. In: Symposium on Non-linear Dynamics and Stochastic Mechanics 1997, Int. Mechanical Engineering Congress and Exposition, Dallas, 1997

38. Rosmann, Th.: Eine Laufmaschine für Rohre. Fortschrittberichte VDI, Reihe 8, Nr. 732, Düsseldorf: VDI-Verlag 1998

39. Srnik, J.; Pfeiffer, F.: Dynamics of CVT Chain Drives: Mechanical Model and Verification. Proc. of DETC '97, ASME Design Eng. Techn. Conf., 1997

40. Srnik, J.: Dynamik von CVT-Keilkettengetrieben. In: Fortschritt-Berichte VDI, Reihe 12, Nr. 372, Düsseldorf: VDI-Verlag 1999

41. Wö sle, M.; Pfeiffer, F.: Dynamics of multibody systems containing dependent unilateral constraints with friction. J. Vib. and Control 2(2) (1996) 161-192

42. Wö sle, M.: Dynamik von räumlichen strukturvarianten Starrkörpersystemen. Fortschrittberichte VDI, Reihe 19, Nr. 213, Düsseldorf: VDI-Verlag 1997

43. Wolfsteiner, P.; Pfeiffer, F.: The Parts Transportation in a Vibratory Feeder. IUTAM Symposium on Unilateral Multibody Dynamics, München, 1998

44. Wolfsteiner, P.: Dynamik von Vibrationsförderern. Dissertation TU-München, Lehrstuhl B für Mechanik, 1999