# Formulatin of spatial contact situations in rigid multibody systems

Christopher Glocker

► To cite this version:

Christopher Glocker. Formulatin of spatial contact situations in rigid multibody systems. Computer Methods in Applied Mechanics and Engineering, 1999, 177, pp.199 - 214. 10.1016/S0045-7825(98)00381-8. hal-01406289

HAL Id: hal-01406289

https://hal.science/hal-01406289

Submitted on 1 Dec 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_1.png" /></figure>
# Formulation of spatial contact situations in rigid multibody systems

Ch. Glocker

Lehrstuhl B für Mechanik, Technische Universität München, D-85747 Garching, Germany

The paper gives an overview of possible formulations of spatial contact situations as they occur in dynamic systems composed of rigid bodies. We will discuss formulations of the contact problem of sliding and rolling as linear and nonlinear complementarity problems in standard form, where general approximations of convex sets by polytopes are also presented. For completeness, a short treatise on spatial contact kinematics and the structure of the equations of motion is given.

## 1. Introduction

Multibody systems are usually considered to be bilaterally coupled or to possess interconnections in the form of force laws. On the other side, a wide area of mechanical systems and machine applications can be characterized by unilateral behavior including unilateral constraints with impact and friction phenomena. Events of that kind occurring in multiple contacts of multibody systems may depend on each other so that each transition from slip to stick and each impact event changes the complete contact configuration. A description of these processes conveniently starts with a complementarity law of all contact situations which states that either terms of relative kinematics are zero and corresponding force expressions are not zero or vice versa.

In the case of a closed sliding contact between two smooth surfaces the process of formulating the corresponding linear complementarity problem is straightforward and easy because of the symmetric structure of the Signorini-Fichera-conditions for discrete systems. First difficulties usually occur when this problem is accompanied by sticking contacts that might undergo a transition to sliding. This situation requires a formulation of the Coulomb friction law on the level of acceleration, sometimes called the rolling contact problem. The formulation of the complementarity problem that is also linear in the planar case and in the spatial case when friction pyramid approximations are used is more complicated. It may lead to improper structures that forbid a treatment of dependent constraints as they usually occur in multibody systems because of additional inversion procedures. Finally, if the conic friction law is considered, nonlinearities are involved in the formulation of the complementarity conditions, and one cannot expect linear complementarity to be available. In most cases this situation is treated with special algorithms such as in [1,9,8] but no formulation as a standard problem is given.

The main objective of this paper is therefore the formulation of the conic Coulomb friction law as a nonlinear complementarity problem (NCP) in standard form. We will also provide some linear complementarity formulations (LCP) of approximations of the friction cone that are also used by other authors and that may probably turn out to be the most suitable. In our formulation only one contact will be considered which is sufficient to demonstrate the structure of the underlying LCPs and NCPs. The extension to multiple contacts is trivial since only certain matrices have to be enlarged and rearranged. This procedure can be found in [13], for example, but only for the friction pyramid model. Impacts are also not treated. Since the structure of the impact equations is very close to the Newton/Euler equations [3,11], one might immediately transfer the results if the impact laws are chosen as in [8].

The paper is organized as follows: Sections 2–4 give a short overview of the structure of the equations of
motion, when spatial contact situations are treated. In Sections 5–9 the contact laws are introduced, and the complementarity formulations are derived. In Section 2 some basic geometric and kinematic relations of single rigid bodies are put together and are then used in Section 3 to formulate the kinematic equations of spatial contact problems. The distance, the relative velocities, and their changes over time which we call relative accelerations, are defined there. Section 4 starts with Newton’s second law and Euler’s axiom and describes how the contact forces are considered in the equations of motion. In Section 5 the contact laws are introduced, consisting of the impenetrability condition in the normal direction and Coulomb friction in the tangential directions. Both contact laws are first stated in their original form at the level of displacements and velocities, respectively, and then transferred to the acceleration level. Using the results obtained so far, the sliding contact problem is treated in Section 6 and formulated as a linear complementarity problem. Section 7 deals with the rolling contact problem by using friction pyramid approximations of the Coulomb friction law. Based on the approximation techniques presented in [9,5], we pass the LCP formulation process [10] and generalize this approach to arbitrary polytopes. This leads to LCPs with a dimension that is equal to the number of inequalities used for the approximation of the friction cone. In Section 8 and 9 the conic friction law is discussed. Standard NCP formulations are presented for the conic friction law which are minimal with respect to the number of equations in use in Section 8, whereas a continuous and differentiable formulation is given in Section 9.

## 2. Geometry of surfaces

In this section some basic kinematic relations are put together which will later be needed when contact kinematical magnitudes such as relative velocities and distances have to be derived. The approach is straightforward and does not use the tensor calculus notation from differential geometry.

We consider two bodies moving with $v_{P_i}$, $\Omega_i$ ($i = 1, 2$) where $v_{P_i}$ denotes the velocity of the body-fixed point $P$ of body $i$ and $\Omega_i$ its angular velocity. Both bodies are assumed to be convex and smooth, at least in a neighborhood of a potential contact point, and may be represented by a parametrization of their surfaces in the form $\mathbf{r}_{P\Sigma i} = \mathbf{r}_{P\Sigma i}(\mathbf{q}, \xi_i, \eta_i)$, see Fig. 1 [7].

Note that $\mathbf{r}_{P\Sigma}$ only depends on the surface parameters $\xi$ and $\eta$ when stated in the corresponding body-fixed frame whereas an additional dependence on the generalized coordinates $q$ is observed when using arbitrary frames. For each of the bodies one defines an orthonormal basis ($u, v, n$) via the tangents $s := \partial_{\xi} \mathbf{r}_{P\Sigma}$ and $t := \partial_{\eta} \mathbf{r}_{P\Sigma}$ in order to obtain

<figure><img src="image_2.png" /><figcaption>Fig. 1. The spatial contact problem.</figcaption></figure>
$$ \mathbf{n} := \frac{\mathbf{s} \times \mathbf{t}}{|\mathbf{s} \times \mathbf{t}|}, \quad \mathbf{u} := \frac{\mathbf{s}}{|\mathbf{s}|}, \quad \mathbf{v} := \mathbf{n} \times \mathbf{u} \qquad (1) $$

where we agree that $(\xi, \eta)$ are ordered such that $\mathbf{n}$ becomes the outward normal. The (absolute) changes in time of the contour vector $\dot{\mathbf{r}}_{P\Sigma}$ and the trihedral $(\dot{\mathbf{u}}, \dot{\mathbf{v}}, \dot{\mathbf{n}})$ are

$$ \begin{aligned} \dot{\mathbf{r}}_{P\Sigma} &= \mathbf{\Omega} \times \mathbf{r}_{P\Sigma} + \mathbf{R}\zeta \quad \text{with} \quad \mathbf{R} = \partial_{\xi} \mathbf{r}_{P\Sigma} = (\mathbf{s}, \mathbf{t}) \\ \dot{\mathbf{k}} &= \mathbf{\Omega} \times \mathbf{k} + \mathbf{K}\zeta \quad \text{with} \quad \mathbf{K} = \partial_{\xi} \mathbf{k} \end{aligned} \qquad (2) $$

where $\zeta = (\xi, \eta)^T$ and $(\mathbf{k}, \mathbf{K})$ is any of $(\mathbf{u}, \mathbf{U}), (\mathbf{v}, \mathbf{V}), (\mathbf{n}, \mathbf{N})$. The velocity $\mathbf{v}_{\Sigma}$ of a point moving on the surface in terms of the rigid body motion $(\mathbf{v}_{\mathrm{P}}, \mathbf{\Omega})$ is thus given by $\mathbf{v}_{\Sigma} = \mathbf{v}_{\mathrm{P}} + \dot{\mathbf{r}}_{P\Sigma}$ which yields

$$ \mathbf{v}_{\Sigma} = \mathbf{v}_{\mathrm{C}} + \mathbf{R}\zeta \quad \text{with} \quad \mathbf{v}_{\mathrm{C}} = \mathbf{v}_{\mathrm{P}} + \mathbf{\Omega} \times \mathbf{r}_{P\Sigma} \qquad (3) $$

Here, $\mathbf{v}_{\mathrm{C}}$ denotes the rigid body velocity portion of $\mathbf{v}_{\Sigma}$. Analogously, we may find the acceleration $\mathbf{a}_{\mathrm{C}} = \dot{\mathbf{v}}_{\mathrm{C}}$ by differentiation of the second equation in (3),

$$ \mathbf{a}_{\mathrm{C}} = \mathbf{a}_{\mathrm{Q}} + \mathbf{\Omega} \times \mathbf{R}\zeta \quad \text{with} \quad \mathbf{a}_{\mathrm{Q}} = \mathbf{a}_{\mathrm{P}} + \dot{\mathbf{\Omega}} \times \mathbf{r}_{P\Sigma} + \mathbf{\Omega} \times (\mathbf{\Omega} \times \mathbf{r}_{P\Sigma}). \qquad (4) $$

Finally, we express $\mathbf{v}_{\mathrm{C}}$ and $\mathbf{a}_{\mathrm{C}}$ by the generalized magnitudes $\dot{\mathbf{q}}$ and $\ddot{\mathbf{q}}$,

$$ \mathbf{v}_{\mathrm{C}} = \mathbf{J}_{\mathrm{C}}\dot{\mathbf{q}} + \hat{\mathbf{t}}_{\mathrm{C}}, \quad \mathbf{a}_{\mathrm{C}} = \mathbf{J}_{\mathrm{C}}\ddot{\mathbf{q}} + \bar{\mathbf{t}}_{\mathrm{C}} \qquad (5) $$

where $\bar{\mathbf{t}}_{\mathrm{C}}$ is given by $\bar{\mathbf{t}}_{\mathrm{C}} = \dot{\mathbf{J}}_{\mathrm{C}}\dot{\mathbf{q}} + \dot{\hat{\mathbf{t}}}_{\mathrm{C}}$, and $\mathbf{J}_{\mathrm{C}}$ is the Jacobian of the contact point as it will turn out later.

### 3. Contact kinematics

Within the framework as presented in the foregoing section one is able to derive the main equations of contact kinematics in dynamics which are: The definition of the contact points, their distance, their relative velocities with respect to the normal and the tangential directions, and the corresponding relative accelerations. In order to do the first step one introduces a distance vector

$$ \mathbf{r}_{\mathrm{D}} := \mathbf{r}_{O\Sigma 2} - \mathbf{r}_{O\Sigma 1}, \qquad (6) $$

cf. [3], where $\mathbf{r}_{O\Sigma i}$ are vectors pointing from a common inertially fixed point $O$ to the surfaces under consideration.

Next, the two tangent planes spanned by $(\mathbf{u}_i, \mathbf{v}_i)$ are adjusted parallel to each other such that the distance vector becomes parallel to each of the normals. The surface parameters $(\zeta_1^*, \zeta_2^*)$ corresponding with that situation solve the four nonlinear equations

$$ \mathbf{n}_1^\mathrm{T}\mathbf{u}_2 = 0, \quad \mathbf{n}_1^\mathrm{T}\mathbf{v}_2 = 0, \quad \mathbf{r}_\mathrm{D}^\mathrm{T}\mathbf{u}_1 = 0, \quad \mathbf{r}_\mathrm{D}^\mathrm{T}\mathbf{v}_1 = 0 \qquad (7) $$

and are called the contact parameters; the corresponding points at the surfaces given by $\mathbf{r}_{P\Sigma i}(\zeta_i^*)$ are called the contact points. When the bodies are moving under the conditions (7) for the contact points should not change [7,3], thus

$$ \begin{aligned} \dot{\mathbf{n}}_1^\mathrm{T}\mathbf{u}_2 + \mathbf{n}_1^\mathrm{T}\dot{\mathbf{u}}_2 &= 0, \quad \dot{\mathbf{r}}_\mathrm{D}^\mathrm{T}\mathbf{u}_1 + \mathbf{r}_\mathrm{D}^\mathrm{T}\dot{\mathbf{u}}_1 &= 0 \\ \dot{\mathbf{n}}_1^\mathrm{T}\mathbf{v}_2 + \mathbf{n}_1^\mathrm{T}\dot{\mathbf{v}}_2 &= 0, \quad \dot{\mathbf{r}}_\mathrm{D}^\mathrm{T}\mathbf{v}_1 + \mathbf{r}_\mathrm{D}^\mathrm{T}\dot{\mathbf{v}}_1 &= 0 \end{aligned} \qquad (8) $$

must hold which constitutes a system of four linear equations in $\dot{\zeta}_i^*$ when substituting $(\dot{\mathbf{n}}_i, \dot{\mathbf{u}}_i, \dot{\mathbf{v}}_i)$ from (2) and $\mathbf{v}_{\Sigma i}$ from (3) with $\dot{\mathbf{r}}_\mathrm{D} = \mathbf{v}_{\Sigma 2} - \mathbf{v}_{\Sigma 1}$.

After the solution $(\zeta_1^*, \zeta_2^*)$ of (7) has been found the distance $g_N$ of the contact points can be calculated as

$$ g_N = \mathbf{n}_1^\mathrm{T}\mathbf{r}_\mathrm{D} \qquad (9) $$

with values greater than zero for separation and values less than zero for overlapping. The relative velocities of the contact points with respect to the three directions $(\mathbf{u}_1, \mathbf{v}_1, \mathbf{n}_1)$ are then given by
$$ \dot{g}_K = (v_{C2} - v_{C1})^T k_1, \quad (k, k_1) = (u_1, v_1), \quad (u, v) = (u_1, v_1), \quad (n, n_1) \tag{10} $$

where it can be easily verified that differentiation of (9) with respect to time results in (10) for $(k, k_1) = (n, n_1)$. The temporal changes of the relative velocities $\dot{g}_K$ in (10) are given by

$$ \ddot{g}_K = (a_{C2} - a_{C1})^T k_1 + (v_{C2} - v_{C1})^T \dot{k}_1 \tag{11} $$

where $a_{Ci}$ is known from (4), $\dot{k}_1$ may be taken from (2), and $\dot{\xi}_i^*$ involved in $a_{Ci}$ and $\dot{k}_1$ is the solution of (8).

Substituting $v_{Ci}$ and $a_{Ci}$ from (5) into (10), (11) yields a representation of the relative velocities and accelerations in the form

$$ \dot{g}_K = w_K^T \dot{q} + \hat{w}_K, \quad \ddot{g}_K = w_K^T \ddot{q} + \bar{w}_K \tag{12} $$

with $w_K = (J_{C2} - J_{C1})^T k_1$, $\hat{w}_K = (\hat{I}_{C2} - \hat{I}_{C1})^T k_1$, and $\bar{w}_K = \dot{w}_K^T \dot{q} + \dot{w}_K^T \ddot{q}$. Finally we, arrange the two tangential relative velocities in the vector $\dot{g}_T$, i.e. $\dot{g}_T^T := (\dot{g}_U, \dot{g}_V)$ which gives us the sliding direction of the contact points in the common tangent plane. With the abbreviations $W_T := (w_U, w_V)$, $\dot{w}_T^T := (\dot{w}_U, \dot{w}_V)$, $\bar{w}_T^T := (\bar{w}_U, \bar{w}_V)$ we rewrite (12) in the form

$$ \begin{aligned} \dot{g}_N &= w_N^T \dot{q} + \hat{w}_N, & \dot{g}_T &= W_T^T \dot{q} + \bar{w}_T \\ \ddot{g}_N &= w_N^T \ddot{q} + \bar{w}_N, & \ddot{g}_T &= W_T^T \ddot{q} + \bar{w}_T \end{aligned} \tag{13} $$

which will serve as the contact kinematic equations in standard form in the sequel.

## 4. Kinetics

The dynamics of a multibody system with $n$ bodies and one contact is described by the Newton-Euler equations in the form

$$ \sum_{i=1}^{n} [J_{S}^{T}(\dot{\boldsymbol{p}} - \boldsymbol{F}) + J_{R}^{T}(\dot{\boldsymbol{L}} - \boldsymbol{M})]_{i} - J_{C1}^{T}\boldsymbol{F}_{C1} - J_{C2}^{T}\boldsymbol{F}_{C2} = 0, \tag{14} $$

see e.g. [3,11]. The terms $p_i$ and $L_i$ denote the momentum and the moment of momentum of body $i$, and $F_i$ and $M_i$ are the applied forces and moments, respectively. The contact forces acting at the contact points $C_1$ and $C_2$ are denoted by $F_{C1}$ and $F_{C2}$. In order to perform a projection into the configuration space all terms have to be premultiplied by the corresponding Jacobians. In particular, $J_{C1}$ and $J_{C2}$ are the Jacobians of the contact points which have already been introduced in Eq. (5).

Decomposing the contact forces into components corresponding to the trihedral $(u_1, v_1, n_1)$ admits, together with the cutting principle, a unique representation in the form

$$ F_{C2} = -F_{C1} := n_1 \lambda_N + u_1 \lambda_U + v_1 \lambda_V \tag{15} $$

Putting (15) into (14) and using the abbreviations $w_K = (J_{C2} - J_{C1})^T k_1$ and $W_T = (w_U, w_V)$ from Eq. (12) and below we arrive with

$$ M \ddot{q} - h - w_N \lambda_N - W_T \lambda_T = 0 \tag{16} $$

where the vector $\lambda_T$ of the tangential contact forces is defined to be $\lambda_T^T = (\lambda_U, \lambda_V)$. The sum in (14) is represented by the term $(M\ddot{q} - h)$ with $M$ as the symmetric and positive definite mass matrix of the system and $h$ as the vector of the gyroscopical accelerations and the classical applied forces. The contact laws are still missing and have to be chosen according to the physical properties of the contacts.

## 5. Contact laws

Generally, every force occurring in a rigid multibody system may be expressed in terms of some relative kinematic magnitudes such as relative displacements and velocities. Classical applied forces are represented by continuous functions whereas classical bilateral constraints already require a representation via set-valued
mappings. Contact laws as used in the sequel belong to the latter class and constitute mixed-type force characteristics, acting in some areas as classical applied forces, in other areas as constraints.

With $g_N$ being the distance of the contact points and $\lambda_N$ the corresponding normal force one might express impenetrability of rigid bodies by the Signorini–Fichera condition [1,3,5,6,9,10]

$$g_N \ge 0, \quad \lambda_N \ge 0, \quad g_N \lambda_N = 0 \qquad (17)$$

This condition may be transferred to the velocity level and, in a second step, to the acceleration level [3,11,10]

$$ \begin{aligned} g_N > 0 & \Rightarrow \lambda_N = 0 \\ g_N = 0, \quad \dot{g}_N > 0 & \Rightarrow \lambda_N = 0 \\ g_N = 0, \quad \dot{g}_N = 0: \quad \lambda_N \ge 0, \quad \ddot{g}_N \ge 0, \quad \lambda_N \ddot{g}_N = 0 \end{aligned} \qquad (18) $$

describing the detachment of a closed contact and also the behavior of open contacts. Note that the occurrence of a new closed contact is usually described by the conditions $g_N = 0$, $\dot{g}_N < 0$ leading to an impact which is not covered by (18). Impacts are already excluded by the formulation of the equations of motion (16), because a proper formulation of impacts must be done in terms of equalities of measures leading to measure differential equations [8] which are not described by (16). Furthermore, the process of proceeding from (17) to (18) requires certain continuity assumptions on the solution $q(t)$. In this context one has to understand $\dot{q}(t)$ to be the right limit $u^+(t)$ of some velocity function of bounded variations $u(t)$, i.e. the post impact velocity in the case of a collision. For impact-free motion one has $\dot{q}(t) = u^+(t) = u^-(t)$ since $q(t)$ is differentiable. This also applies to $\dot{g}_N(t)$ through (10). Similarly, the accelerations $\ddot{q}(t)$ (and $\ddot{g}_N(t)$ through (11), respectively) have to be understood as the right derivatives of the velocity functions $u(t)$, i.e. the changes in the velocities that describe the further evolution of the system [2].

In the tangential directions we assume Coulomb friction given by the friction law [1,5,6,9,10]

$$ \begin{cases} \dot{g}_T = 0 \quad \Rightarrow \quad |\lambda_T| \le \mu_0 \lambda_N \\ \dot{g}_T \ne 0 \quad \Rightarrow \quad |\lambda_T| = \mu_N \lambda_N, \quad \lambda_T = -\beta \dot{g}_T, \quad \beta \ge 0 \end{cases} \qquad (19) $$

where $\mu$ denotes the coefficient of friction which depends smoothly on $\dot{g}_T$, and $\mu_0 = \mu(\dot{g}_T = 0)$. Similarly we obtain from (19) the friction law on the acceleration level in the form [10,14]

$$ \begin{aligned} \dot{g}_T \ne 0 & \Rightarrow \lambda_T = -e\mu\lambda_N \\ \dot{g}_T = 0: \quad \begin{cases} \ddot{g}_T = 0 \quad \Rightarrow \quad |\lambda_T| \le \mu_0 \lambda_N \\ \ddot{g}_T \ne 0 \quad \Rightarrow \quad |\lambda_T| = \mu_0 \lambda_N, \quad \lambda_T = -\beta \ddot{g}_T, \quad \beta \ge 0 \end{cases} \end{aligned} \qquad (20) $$

with $e$ being the unit vector of the sliding direction, i.e. $e = \dot{g}_T / |\dot{g}_T|$. Note that even in the case of impact free motion, Coulomb friction obviously causes acceleration jumps for transitions from sliding to stiction. Thus, even in this impact free case, the accelerations $\ddot{q}$ are not defined at the transition points but only their left and right limits [2]. The first line in Eq. (20) describes sliding of the contact up to velocities $\dot{g}_T = 0$, whereas the second condition takes into account stiction as well as the transitions from stiction to sliding, and is valid for the right derivatives of the velocity functions.

The Newton/Euler equations (16), the contact laws in the normal (18) and the tangential (20) directions, and the contact kinematic terms from (9) and (13) as far as needed provide a complete set of equations in order to finally solve (16) for the unknowns $\ddot{q}$. In the case of closed contacts this is usually done by two steps [3,10]: First, the accelerations $\ddot{q}$ are eliminated with the help of (16) from the relative accelerations in the second line of (13). Together with the contact laws they constitute a system of inequalities which has to be solved for the relative accelerations $\ddot{g}$ and the contact forces $\lambda$. When this solution has been found, $\lambda$ is re-substituted into the Newton/Euler equations (16) which permits finally the calculation of $\dot{q}$.

In the remaining part of the paper we concentrate on the formulation of the inequality system for sliding and rolling contact situations. We present LCP and NCP formulations in standard form for friction pyramids and for the friction cone, respectively. In order to avoid unessential and cumbersome indexing only one contact (as already done in Eq. (16)) is considered. The blow up to multi-contact configurations is straightforward and
consists of nothing more than a rearrangement of certain matrices, see e.g. [10] for spatial problems and [3,11] for the planar case.

## 6. Sliding contacts

In this section we derive the linear complementarity problem for a closed contact ($g_N = 0, \dot{g}_N = 0$) which is sliding ($\dot{g}_T \neq 0$). This situation is described by the Newton/Euler equations (16), the normal contact law for $g_N = 0, \dot{g}_N = 0$ from the third line in (18), and the tangential contact law for sliding contacts $\dot{g}_T \neq 0$ from the first equation in (20). Additionally as kinematical relation the third equation in (13) is needed, since $\ddot{g}_N$ is involved in the normal contact law. Altogether these equations write

$$ 
\begin{aligned}
& M\ddot{q} - h - w_N\lambda_N - W_T\lambda_T = 0 \\
& \ddot{g}_N = w_N^T\ddot{q} + \bar{w}_N \\
& \lambda_N \ge 0, \quad \ddot{g}_N \ge 0, \quad \lambda_N \ddot{g}_N = 0 \\
& \lambda_T = -e\mu\lambda_N
\end{aligned} 
\quad (21) 
$$

Using the fourth equation $\lambda_T$ is substituted into the first one which is then solved for $\ddot{q}$. Substituting $\ddot{q}$ into the second equation we obtain the LCP for the sliding contact [10,14,13]

$$ 
\begin{aligned}
\ddot{g}_N &= w_N^T M^{-1} (w_N - W_T e \mu) \lambda_N + (w_N^T M^{-1} h + \bar{w}_N) \\
\ddot{g}_N \ge 0, \quad \lambda_N \ge 0, \quad \ddot{g}_N \lambda_N = 0.
\end{aligned} 
\quad (22) 
$$

Obviously, friction may cause the LCP-matrix (for one contact being a scalar) to become indefinite.

## 7. Friction pyramids for rolling contacts

Rolling contacts are characterized by vanishing tangential velocities $\dot{g}_T = 0$. In contrast to the previous section, the tangential contact law is now described by the second condition in (20) which holds for stiction as well as for the transition to sliding. Before discussing approximations of this friction law by polytopes we will rewrite (20) for $\dot{g}_T = 0$ in a more suitable form which will also be used in the following sections. Let

$$ D_T := \{ \lambda_T \mid \|\lambda_T\| \le \mu_0 \lambda_N \} \quad (23) $$

be the (convex) set of the admissible tangential contact forces. The tangential contact law is then given by

$$ -\dot{g}_T \in N_{D_T}(\lambda_T), \quad (24) $$

see e.g. [5,9], where $N_{D_T}(\lambda_T)$ denotes the normal cone to $D_T$ at $\lambda_T$,

$$ N_{D_T}(\lambda_T) = \begin{cases} \kappa e & \text{for } 0 < |\lambda_T| = \mu_0 \lambda_N \\ 0 & \text{for } 0 \le |\lambda_T| < \mu_0 \lambda_N \\ \mathbb{R}^2 & \text{for } 0 = |\lambda_T| = \mu_0 \lambda_N \end{cases} \quad (25) $$

with $e = \lambda_T / |\lambda_T|$ and $\kappa \ge 0$ as depicted in Fig. 2.

**PROPOSITION 1.** The contact laws (20) for $\dot{g}_T = 0$ and (24) are equivalent.

*PROOF.* By inversion, (24) may be also written as $\lambda_T \in \partial I_{D_T}^*(-\ddot{g}_T)$, where $\partial$ denotes the subdifferential from convex analysis and $I_{D_T}^*(-\ddot{g}_T)$ is the support function of $D_T$ defined by $I_{D_T}^*(\dot{\ddot{g}}_T) = \mu_0 \lambda_N |\dot{\ddot{g}}_T|$. Solving this inclusion for $\ddot{g}_T = 0$ we obtain $\lambda_T \in \partial I_{D_T}^*(0) = B(0, \mu_0 \lambda_N)$, where $B(0, r)$ is the closed disk with center 0 and radius $r$. Hence, $|\lambda_T| \le \mu_0 \lambda_N$ which is the first condition in (20). For $\ddot{g}_T \neq 0$ the support function is
<figure><img src="image_3.png" /><figcaption>Fig. 2. Admissible tangential forces for the friction cone.</figcaption></figure>

differentiable. Thus, $\lambda_T = -\mu_0 \lambda_N \dot{g}_T / |\ddot{g}_T|$, i.e. $|\lambda_T| = \mu_0 \lambda_N$ and $\lambda_T = -\beta \dot{g}_T$ with $\beta = \mu_0 \lambda_N / |\ddot{g}_T| \ge 0$ since $\mu_0 \ge 0$, $\lambda_N \ge 0$. This is the second condition in (20).

Following [5] we may approximate the set $D_T$ by a polytope generated by the intersection of $n$ closed half spaces $C_{Ti}$, each of which being represented by an affine inequality $\sigma_i(\lambda_T)$,

$$ C_T := \bigcap_{i=1}^{n} C_{Ti} \quad \text{with } C_{Ti} = \{\lambda_T \mid \sigma_i(\lambda_T) \ge 0\} \\ \text{and } \sigma_i(\lambda_T) = \mu_0 \lambda_N - e_i^T \lambda_T \qquad (26) $$

where $e_i$ denote the outward unit vectors to $D_T$ at some chosen points on the boundary $\partial D_T$ corresponding to angles $\alpha_i$, see Fig. 3. We now assume the contact law to be of the form

$$ -\ddot{g}_T \in N_{C_T}(\lambda_T) \qquad (27) $$

This is the friction cone approximation presented in [5]. For inner approximations of the friction cone we refer to [9]. Note that the original friction law (24) and its approximation (27) differ only in the definition of the sets $D_T$ and $C_T$. We call the $\sigma_i$ the friction saturation with respect to $C_{Ti}$ because it determines the distance between a

<figure><img src="image_4.png" /><figcaption>Fig. 3. Approximation of the friction disk by polytopes.</figcaption></figure>
given tangential force $\lambda_T$ and the boundary of $C_{T_i}$. It therefore describes the amount of forces which might additionally act on the contact in the direction $e_i$ until sliding starts.

In the following we will derive a system of complementarity conditions that describe the friction law (27). The half-spaces $C_{T_i}$ defined in (26) have at least the point 0 in common since $\sigma_i(0) \ge 0$ because $\lambda_N \ge 0$, $\mu_0 \ge 0$. Thus, the normal cone to $C_T = \cap_{i=1}^n C_{T_i}$ is given by

$$N_{C_T}(\lambda_T) = \sum_{i=1}^{n} N_{C_{T_i}}(\lambda_{T_i}), \qquad (28)$$

cf. [12], and the contact law (27) may be expressed by the normal cone conditions [12]

$$-\ddot{g}_T \in \sum_{i=1}^{n} -\kappa_i \partial \sigma_i(\lambda_T); \quad \sigma_i \ge 0, \quad \kappa_i \ge 0, \quad \sigma_i \kappa_i = 0 \qquad (29)$$

Since $\sigma_i(\lambda_T)$ is differentiable this inclusion may be further reduced to the equation

$$-\ddot{g}_T = \sum_{i=1}^{n} e_i \kappa_i; \quad \sigma_i \ge 0, \quad \kappa_i \ge 0, \quad \sigma_i \kappa_i = 0 \qquad (30)$$

providing already the desired complementarity conditions. Note that the tangential relative acceleration $\ddot{g}_T$ is expressed by a positive linear combination of given directions $e_i$. The scalars $\kappa_i$ are accelerations that are complementary to the friction saturations $\sigma_i$. For a given direction $e_i$, either the relative acceleration $\kappa_i$ or the friction saturation $\sigma_i$ vanishes.

We may now summarize the equations that describe the rolling contact problem: The Newton/Euler equations (16), the normal contact law for $g_N = 0$, $\dot{g}_N = 0$ from the third line in (18), the approximation of the tangent contact law consisting of the definition of the friction saturations $\sigma_i$ in (26), the acceleration equation in (30), and the complementarities in (30). Finally, we need from kinematics the two equations in the second line of (13) since both, $\ddot{g}_N$ and $\dot{g}_T$ enter the contact laws. Altogether we have the following complete description of the rolling contact problem:

$$\begin{align}
& M\ddot{q} - h - w_N\lambda_N - W_T\lambda_T = 0 \nonumber \\
& \ddot{g}_N = w_N^T\ddot{q} + \bar{w}_N \nonumber \\
& \ddot{g}_T = W_T^T\ddot{q} + \bar{w}_T \nonumber \\
& \lambda_N \ge 0, \quad \ddot{g}_N \ge 0, \quad \lambda_N \ddot{g}_N = 0 \tag{31} \\
& \sigma_i = \mu_0 \lambda_N - e_i^T \lambda_T, \quad i=1, \dots, n \nonumber \\
& -\ddot{g}_T = \sum_{i=1}^n e_i \kappa_i \nonumber \\
& \sigma_i \ge 0, \quad \kappa_i \ge 0, \quad \sigma_i \kappa_i = 0, \quad i=1, \dots, n \nonumber
\end{align}$$

The remaining part of the section is devoted to the formulation of (31) as a linear complementarity problem. There are many approaches spread through the literature how to perform this step, more or less suited for contact problems in multibody dynamics. In our opinion a good formulation should fulfill at least the following three demands: In the first place, any inversion of an expression $W^T M^{-1} W$ has to be avoided. For multi-contact problems the matrix $W$ consists of a selection of constraint vectors $w_{Ni}$, $w_{Ui}$, $w_{Vi}$ depending on the actual contact state which changes during time. In multibody dynamics there is nearly no way to avoid contact situations where these vectors become linearly dependent. Even for the simple case of a planar rectangular block resting on a surface this situation occurs when the unilateral frictional constraints are modeled by two contact points at the block's lower corners (see e.g. [3] and [11]). Thus, formulations as presented in [4] and [9] usually fail for dynamic multibody problems. This situation has been overcome in [3,11] by introducing additional variables and thus enlarging the dimension of the LCP which leads to our second demand: Since the LCP has to be solved at every time step during numerical integration its dimension must be kept as small as possible. In our case an optimal LCP would have dimension $n+1$, corresponding to the number of inequalities introduced for the approximation of the friction disk, plus the condition for the normal contact law. Thirdly, a standard LCP formulation is desirable.
The formulation of the LCP for the rolling contact problem is not as straightforward as for the sliding case in Section 6 where the Newton/Euler equations had just to be substituted into the kinematics equation in the normal direction. The reason for that is the missing complementarity between $\ddot{g}_T$ and $\lambda_T$. Thus, one has to eliminate $\lambda_T$ from the Newton/Euler equations and to reformulate the tangential acceleration $\ddot{g}_T$ such that complementarity appears in the same manner as for the normal contact law. Thereby, the tricky part of this process consists of avoiding the aforementioned inversion of $W^T M^{-1} W$. Following [13] we choose an indexing $i$ in the last three equations in (31) in a manner that the directions $e_i$ for $i=1, 2$ are linearly independent. Introducing the abbreviations

$$ \sigma_p = \begin{pmatrix} \sigma_1 \\ \sigma_2 \end{pmatrix}, \quad \kappa_p = \begin{pmatrix} \kappa_1 \\ \kappa_2 \end{pmatrix}, \quad \mu_p = \begin{pmatrix} \mu_0 \\ \mu_0 \end{pmatrix}, \quad I_p = (e_1 \ e_2) \qquad (32) $$

$$ \sigma_R = \begin{pmatrix} \sigma_3 \\ \vdots \\ \sigma_n \end{pmatrix}, \quad \kappa_R = \begin{pmatrix} \kappa_3 \\ \vdots \\ \kappa_n \end{pmatrix}, \quad \mu_R = \begin{pmatrix} \mu_0 \\ \vdots \\ \mu_n \end{pmatrix}, \quad I_R = (e_3 \cdots e_n) $$

we may then rewrite the fifth line in (31) in matrix notation,

$$ \sigma_p = \mu_p \lambda_N - I_p^T \lambda_T, \quad \sigma_R = \mu_R \lambda_N - I_R^T \lambda_T \qquad (33) $$

Since $I_p = (e_1 e_2)$ is invertible we may solve the left equation in (33) for $\lambda_T$,

$$ \lambda_T = I_p^{-T} \mu_p \lambda_N - I_p^{-T} \sigma_p \qquad (34) $$

and then substitute the result into the right equation,

$$ \sigma_R = (\mu_R - I_R^T I_p^{-T} \mu_p) \lambda_N + I_R^T I_p^{-T} \sigma_p \qquad (35) $$

The accelerations from the sixth line in (31) are solved for $\kappa_p$ which yields

$$ \kappa_p = -I_p^{-1} \ddot{g}_T - I_p^{-1} I_R \kappa_R \qquad (36) $$

Finally, we put $\lambda_T$ from (34) into the first equation in (31) and substitute $\ddot{g}_T$ from the third line into (36). This gives, together with (35), the new system of equations

$$ \begin{aligned} & M \ddot{q} - h - (w_N + W_T I_p^{-T} \mu_p) \lambda_N + W_T I_p^{-T} \sigma_p = 0 \\ & \ddot{g}_N = w_N^T \ddot{q} + \bar{w}_N \\ & \kappa_p = -I_p^{-1} W_T^T \ddot{q} - I_p^{-1} \bar{w}_T - I_p^{-1} I_R \kappa_R \\ & \sigma_R = (\mu_R - I_R^T I_p^{-T} \mu_p) \lambda_N + I_R^T I_p^{-T} \sigma_p \\ & \ddot{g}_N \ge 0, \quad \lambda_N \ge 0, \quad \ddot{g}_N \lambda_N = 0 \\ & \kappa_p \ge 0, \quad \sigma_p \ge 0, \quad \kappa_p^T \sigma_p = 0 \\ & \sigma_R \ge 0, \quad \kappa_R \ge 0, \quad \sigma_R^T \kappa_R = 0 \end{aligned} \qquad (37) $$

that describes the rolling contact problem. In a last step we solve the first equation in (37) for $\ddot{q}$ and substitute the result into the second and the third equation. With $W_p := W_T I_p^{-T}$ and $\bar{w}_p := I_p^{-1} \bar{w}_T$ we finally obtain

$$ \begin{pmatrix} \ddot{g}_N \\ \kappa_p \\ \sigma_R \end{pmatrix} = \begin{pmatrix} w_N^T M^{-1} (w_N + W_p \mu_p) & -w_N^T M^{-1} W_p & 0 \\ -W_p^T M^{-1} (w_N + W_p \mu_p) & W_p^T M^{-1} W_p & -I_p^{-1} I_R \\ \mu_R & I_R^T I_p^{-T} \mu_p & 0 \end{pmatrix} \begin{pmatrix} \lambda_N \\ \sigma_p \\ \kappa_R \end{pmatrix} + \begin{pmatrix} w_N^T M^{-1} h + \bar{w}_N \\ -W_p^T M^{-1} h - \bar{w}_p \\ 0 \end{pmatrix} \qquad (38) $$

which is, together with the inequalities and complementarities from the last three lines in (37), the desired linear complementarity problem of dimension $n+1$ in standard form [13].

If the friction disc $D_T$ is approximated, for example, by a square, the directions $e_i$, $i=1, \dots, 4$ are simply $e_1^T = (1 \ 0)$, $e_2^T = (0 \ 1)$, $e_3^T = (-1 \ 0)$, $e_4^T = (0 \ -1)$. In this case the matrices $I_p$ and $I_R$ as defined in (32) become $I_p = -I_R = E$, where $E$ is the $2 \times 2$ identity matrix, and the dimension of the LCP is 5. Thus, $W_p = W_T$ and $\bar{w}_p = \bar{w}_T$, and the entries in the last row and column of the LCP matrix in (38) become
$$ \mu_R - I_R^T I_P^{-T} \mu_P = 2\mu_T, \quad I_R^T I_P^{-T} = -E, \quad -I_P^{-1} I_R = E \quad (39) $$

with $\mu_T := \mu_P \equiv \mu_R = (\mu_0 \ \mu_0)^T$. Exactly this formulation but expanded to multi-contact configurations is used in [10].

In the case of one-dimensional friction the directions $e_1^T = (1 \ 0)$ and $e_3^T = (-1 \ 0)$ reduce to +1 and -1, respectively, whereas $e_2$ and $e_4$ are no longer needed. The LCP (38) then writes

$$ \begin{pmatrix} \ddot{g}_N \\ \kappa_1 \\ \sigma_3 \end{pmatrix} = \begin{pmatrix} \mathbf{w}_N^T \mathbf{M}^{-1} (\mathbf{w}_N + \mathbf{w}_U \mu_0) & -\mathbf{w}_N^T \mathbf{M}^{-1} \mathbf{w}_U & 0 \\ -\mathbf{w}_U^T \mathbf{M}^{-1} (\mathbf{w}_N + \mathbf{w}_U \mu_0) & \mathbf{w}_U^T \mathbf{M}^{-1} \mathbf{w}_U & 1 \\ 2\mu_0 & -1 & 0 \end{pmatrix} \begin{pmatrix} \lambda_N \\ \sigma_1 \\ \kappa_3 \end{pmatrix} + \begin{pmatrix} \mathbf{w}_N^T \mathbf{M}^{-1} \mathbf{h} + \bar{\mathbf{w}}_N \\ -\mathbf{w}_U^T \mathbf{M}^{-1} \mathbf{h} - \bar{\mathbf{w}}_U \\ 0 \end{pmatrix} \quad (40) $$

and is of dimension 3. The index $U$ denotes the first of the former two tangential directions and has to be understood as introduced in Section 3, Eq. (10). In this case the tangential relative acceleration (30) becomes

$$ \ddot{g}_U = \kappa_3 - \kappa_1, \quad (41) $$

and the friction saturations (26) are

$$ \sigma_1 = \mu_0 \lambda_N - \lambda_U, \quad \sigma_3 = \mu_0 \lambda_N + \lambda_U \quad (42) $$

with $\lambda_U$ being the tangential contact force. Due to the complementarities $\kappa_i \ge 0$, $\sigma_i \ge 0$, $\kappa_i \sigma_i = 0$ we finally note that $\kappa_3$ and $\kappa_1$ correspond to the positive and negative parts of $\ddot{g}_U$ for $\mu_0 \lambda_N > 0$, respectively, i.e. $\ddot{g}_U = \ddot{g}_U^+ - \ddot{g}_U^-$.

The presented method of formulating linear complementarity problems for dynamic multibody systems is general. Every polytope obtained from the intersection of a finite number of half-spaces may be treated in this way. The dimension of the resulting LCP is always equal to the number of half-spaced used. For example, the most general force laws in multibody dynamics are six-dimensional, corresponding to the number of degrees of freedom of an unconstrained body. In this case, a force vector $\lambda \in \mathbb{R}^6$ would enter the Newton/Euler equations, and the elimination of $\lambda$ would require six linearly independent directions $e_i$ to be arranged in $I_P$ from (32).

## 8. Friction cones and NCP-formulations

With the results from the preceding section one is able to treat spatial friction at least by approximations. However, due to the pyramid model, the collinearity of the negative friction forces $\lambda_T$ and the tangential relative accelerations $\ddot{g}_T$ is lost. This lack of consistency with respect to the physical model of friction in mind (and also with respect to the sliding friction problem from Section 6) may be reduced by increasing the number of vertices of the polytope which, however, leads to LCPs of higher dimensions. We therefore try to find a formulation which, on the one hand, incorporates the conic friction law (23), (24) and, on the other hand, requires as less equations as possible. Since this problem is obviously beyond linearity we try to use the next available standard formulation, i.e. the representation as a NCP.

Again, we consider the rolling contact problem which is described by the normal contact law from the third line in (18) and by the tangential contact law from the second condition in (20) or, equivalently, by (23) and (24). Introducing the friction saturation $\sigma_D$ of the friction disk we may rewrite (23) in the form

$$ D_T = \{ \lambda_T \mid \sigma_D(\lambda_T) \ge 0 \} \\ \text{with} \\ \sigma_D(\lambda_T) = \mu_0 \lambda_N - |\lambda_T| \quad (43) $$

From now, the derivation of the NCP follows exactly the procedure presented in Section 7. A crucial step necessary to avoid inversions of terms $W^T M^{-1} W$ had been performed in Eqs. (34) and (37) where the tangential forces $\lambda_T$ had been eliminated from the Newton/Euler equations. This, however, requires at least two equations which can be solved for $\lambda_T$, preferably being linear. We therefore introduce two additional linear inequalities on the tangential forces, not changing $D_T$ when intersection is performed: Let
$$C_{Ti} = \{ \lambda_T | \sigma_i(\lambda_T) \ge 0 \}$$

with

(44)

$$\sigma_i(\lambda_T) = \mu_0 \lambda_N - e_i^T \lambda_T, \quad i = 1, 2$$

where $e_1, e_2$ are two unit vectors as defined in Fig. 4. The set of the admissible tangential forces may now be written as

$$C_T := \bigcap_{i=1}^{2} C_{Ti} \cap D_T \equiv D_T$$

(45)

and the corresponding normal cone becomes

$$N_{C_T}(\lambda_T) = \sum_{i=1}^{2} N_{C_{Ti}}(\lambda_T) + N_{D_T}(\lambda_T) \equiv N_{D_T}(\lambda_T)$$

(46)

Since $C_T = D_T$ we also have $N_{C_T}(\lambda_T) = N_{D_T}(\lambda_T)$. Thus, the tangential contact law (24) may be equivalently expressed by

$$-\dot{g}_T \in N_{C_j}(\lambda_T)$$

(47)

Writing down the normal cone conditions for (47) we obtain

$$-\dot{g}_T \in \sum_{i=1}^{2} e_i \kappa_i - \partial \sigma_D(\lambda_T) \kappa_D$$

with

(48)

$$\kappa_j \ge 0, \quad \sigma_j \ge 0, \quad \kappa_j \sigma_j = 0, \quad j = 1, 2, D$$

where, in contrast to (30), no further simplification is possible because $\sigma_D$ is not differentiable at $\lambda_T = 0$, i.e. $\partial \sigma_D(\lambda_T = 0) = B(0, 1)$ is set-valued.

Instead of (48) we introduce a contact law which is expressed by an equality. We set

$$-\ddot{g}_T = \sum_{i=1}^{2} e_i \kappa_i + e_D(\lambda_T) \kappa_D$$

with

(49)

$$\kappa_j \ge 0, \quad \sigma_j \ge 0, \quad \kappa_j \sigma_j = 0, \quad j = 1, 2, D$$

<figure><img src="image_5.png" /><figcaption>Fig. 4. Friction disk with additional constraints.</figcaption></figure>
and

$$e_D(\lambda_T) = \begin{cases} \lambda_T / |\lambda_T| & \text{if } \lambda_T \neq 0 \\ e_3 & \text{if } \lambda_T = 0 \end{cases}$$

where $e_3$ is the unit vector as defined in Fig. 4.

**PROPOSITION 2.** *The contact laws (48) and (49) are equivalent.*

*PROOF.* For $\lambda_T \neq 0$ the friction saturation $\sigma_D$ as defined in (43) is differentiable. Thus, $\partial\sigma_D(\lambda_T) = \nabla\sigma_D(\lambda_T) = -\lambda_T/|\lambda_T|$ in (48) which is also (49). For $\lambda_T = 0$, Eqs. (48) and (49) become

$$-\ddot{g}_T \in \sum_{i=1}^{2} e_i \kappa_i + B(0, 1) \kappa_D, \quad -\ddot{g}_T = \sum_{i=1}^{2} e_i \kappa_i + e_3 \kappa_D \quad (50)$$

with $\sigma_i = \mu_0 \lambda_N \ge 0$, $\kappa_i \ge 0$, $\sigma_i \kappa_i = 0$ for $i = 1, 2, D$. Let $\mu_0 \lambda_N > 0$. In this case we have $\sigma_i > 0$, $\kappa_i = 0$ for $i = 1, 2, D$, and both equations in (50) reduce to $\dot{g}_T = 0$. Now let $\mu_0 \lambda_N = 0$. We obtain $\sigma_i = \mu_0 \lambda_N = 0$ and $\kappa_i \ge 0$ without any further restrictions. The left equation in (50) becomes $-\dot{g}_T \in \mathbb{R}^2$ which is also true for the right equation because $(e_1, e_2, e_3)$ have been chosen such that the union of the convex cones generated by $(e_i, e_j)$, $i, j = 1, 2, 3$ is $\mathbb{R}^2$.

As in (31) we may now state the rolling contact problem. The first four expressions consisting of the Newton/Euler equations, the kinematic expressions, and the normal contact law may be adopted unchanged from (31), whereas the tangential contact law is now described by Eq. (49) together with the definitions of $\sigma_D$ and $\sigma_i$ from (43) and (44), respectively. Altogether we have

$$\begin{align*}
& \mathbf{M}\ddot{\mathbf{q}} - \mathbf{h} - \mathbf{w}_N\lambda_N - \mathbf{W}_T\lambda_T = 0 \\
& \ddot{g}_N = \mathbf{w}_N^T\ddot{\mathbf{q}} + \bar{\mathbf{w}}_N \\
& \ddot{g}_T = \mathbf{W}_T^T\ddot{\mathbf{q}} + \bar{\mathbf{w}}_T \\
& \lambda_N \ge 0, \quad \ddot{g}_N \ge 0, \quad \lambda_N \ddot{g}_N = 0 \tag{51} \\
& \sigma_i = \mu_0 \lambda_N - \mathbf{e}_i^T\lambda_T, \quad i=1,2 \\
& \sigma_D = \mu_0 \lambda_N - |\lambda_T| \\
& -\ddot{g}_T = \sum_{i=1}^2 e_i \kappa_i + e_D(\lambda_T)\kappa_D \\
& \sigma_j \ge 0, \quad \kappa_j \ge 0, \quad \sigma_j \kappa_j = 0, \quad i=1,2,D
\end{align*}$$

Since the directions $e_1$ and $e_2$ are linearly independent we may use the same notation as in (32). We set

$$\boldsymbol{\sigma}_p = \begin{pmatrix} \sigma_1 \\ \sigma_2 \end{pmatrix}, \quad \boldsymbol{\kappa}_p = \begin{pmatrix} \kappa_1 \\ \kappa_2 \end{pmatrix}, \quad \boldsymbol{\mu}_p = \begin{pmatrix} \mu_0 \\ \mu_0 \end{pmatrix}, \quad \mathbf{I}_p = (\mathbf{e}_1 \ \mathbf{e}_2) \quad (52)$$

in order to solve the fifth equation in (51) for $\lambda_T$,

$$\lambda_T = \mathbf{I}_p^{-T} \boldsymbol{\mu}_p \lambda_N - \mathbf{I}_p^{-T} \boldsymbol{\sigma}_p \quad (53)$$

Putting (53) into the sixth and seventh equation in (51) yields

$$\sigma_D = \mu_0 \lambda_N - |\mathbf{I}_p^{-T}(\boldsymbol{\mu}_p \lambda_N - \boldsymbol{\sigma}_p)| \quad (54)$$

$$\boldsymbol{\kappa}_p = -\mathbf{I}_p^{-1} \ddot{g}_T - \mathbf{I}_p^{-1} \boldsymbol{e}_D(\lambda_N, \boldsymbol{\sigma}_p) \boldsymbol{\kappa}_D \quad (55)$$

where the accelerations $\ddot{g}_T$ have additionally been solved for $\boldsymbol{\kappa}_p$. Note that the direction $\boldsymbol{e}_D$ in (55) as defined in (49) now depends on $\lambda_N$ and $\boldsymbol{\sigma}_p$. As in (37) we eliminate $\lambda_T$ with the help of (53) from the Newton/Euler equations in (51), and we substitute $\ddot{g}_T$ from the third line in (51) into (55). Altogether we have now
$$ M\ddot{q} - h - (w_N + W_T I_p^{-T} \mu_p)\lambda_N + W_T I_p^{-T} \sigma_p = 0 $$

$$ \ddot{g}_N = w_N^T \ddot{q} + \bar{w}_N $$

$$ \kappa_p = -I_p^{-1} W_T^{-1} \ddot{q} - I_p^{-1} \bar{w}_T - I_p^{-1} e_D(\lambda_N, \sigma_p) \kappa_D $$

$$ \sigma_D = \mu_0 \lambda_N - |I_p^{-1}(\mu_p \lambda_N - \sigma_p)| \quad (56) $$

$$ \ddot{g}_N \ge 0, \quad \lambda_N \ge 0, \quad \ddot{g}_N \lambda_N = 0 $$

$$ \kappa_p \ge 0, \quad \sigma_p \ge 0, \quad \kappa_p^T \sigma_p = 0 $$

$$ \sigma_D \ge 0, \quad \kappa_D \ge 0, \quad \sigma_D \kappa_D = 0 $$

where in a last step the accelerations $\ddot{q}$ have to be eliminated from the second and the third line by using the first equation. Similarly to (38) the remaining three equations may be written in the form

$$ \begin{pmatrix} \ddot{g}_N \\ \kappa_p \\ \sigma_D \end{pmatrix} = \begin{pmatrix} w_N^T M^{-1} (w_N + W_p \mu_p) & -w_N^T M^{-1} W_p & 0 \\ -W_p^T M^{-1} (w_N + W_p \mu_p) & W_p^T M^{-1} W_p & 0 \\ \mu_0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \lambda_N \\ \sigma_p \\ \kappa_D \end{pmatrix} \\ + \begin{pmatrix} 0 \\ -I_p^{-1} e_D(\lambda_N, \sigma_p) \kappa_D \\ -|I_p^{-1}(\mu_p \lambda_N - \sigma_p)| \end{pmatrix} + \begin{pmatrix} w_N^T M^{-1} h + \bar{w}_N \\ -W_p^T M^{-1} h - \bar{w}_p \\ 0 \end{pmatrix} \quad (57) $$

where the abbreviations $W_p = W_T I_p^{-T}$ and $\bar{w}_p = I_p^{-1} \bar{w}_T$ have been used. Eq. (57) together with the inequalities and complementarities from the last three lines in (56) constitutes a nonlinear complementarity problem in standard form, i.e. a problem of the form $y = f(x), y \ge 0, x \ge 0, y^T x = 0$, with $f(x)$ being here $f(x) = Ax + g(x) + c$.

The dimension of the NCP in (57) is 4. We suppose that this is also the minimal number of equations necessary to describe spatial friction problems in NCP-standard form: In the case $\mu_0 \lambda_N = 0$ the friction law (47) reduces to $\ddot{g}_T \in \mathbb{R}^2$. In order to express any $\ddot{g}_T \in \mathbb{R}^2$ by a positive linear combination $\sum e_i \kappa_i$ at least three directions $e_i$ are necessary. The scalars $\kappa_i \ge 0$ enter the NCP as independent variables which makes, together with the normal contact law, an overall dimension of at least 4.

## 9. A differentiable NCP-formulation for rolling contacts

There is a serious lack in the NCP-formulation of the last section: The nonlinear function $g(x)$ contained in the right-hand side $f(x)$ of (57) is neither differentiable nor continuous. The direction $e_D$ defined in (49) may jump at points $\lambda_T = 0$, i.e. if $\mu_p \lambda_N = \sigma_p$. In order to overcome this problem we will finally present a NCP formulation with a continuous and differentiable right-hand side but with an overall dimension increased to 5.

For this purpose we characterize the friction disc $D_T$ by a parabolic function $\sigma_D$ which is differentiable everywhere

$$ D_T = \{ \lambda_T | \sigma_D(\lambda_T) \ge 0 \} \quad (58) $$

$$ \sigma_D(\lambda_T) = \mu_0^2 \lambda_N^2 - |\lambda_T|^2 $$

but we note that the gradient of $\sigma_D$ vanishes for $\lambda_T = 0$. The subdifferential $\partial \sigma_D(\lambda_T) = \nabla \sigma_D(\lambda_T)$ is therefore not capable of describing the normal cone to $D_T$ at $\lambda_T = 0$. We therefore introduce additional inequalities as in (44) but now representing the whole triangle that is depicted in Fig. 4, i.e. the index $i$ is now running from 1 to 3. Obviously the normal cones to $D_T$ and $\cap_{i=1}^3 C_{T,i}$ at $\lambda_T = 0$ are identical and may be represented by using the three directions $e_i$, $i = 1, 2, 3$ that are defined in Fig. 4. The set of the admissible tangential forces $C_T$ and the corresponding normal cone $N_{C_T}(\lambda_T)$ may be expressed in the same manner as in (45) and (46), and the tangential contact law (24) becomes

$$ -\ddot{g}_T \in N_{C_T}(\lambda_T) \quad (59) $$
since $C_T \equiv D_T$ and thus $N_{C_T}(\lambda_T) \equiv N_{D_T}(\lambda_T)$. Since the gradient of $\sigma_D$ vanishes for $\lambda_T = 0$, the normal cone conditions as in (48) are no longer valid in the general case. However, we formally write

$$ -\ddot{\mathbf{g}}_T = \sum_{i=1}^{3} \mathbf{e}_i \kappa_i + 2\mathbf{\lambda}_T \kappa_D \quad (60) $$
with
$$ \kappa_j \ge 0, \quad \sigma_j \ge 0, \quad \kappa_j \sigma_j = 0, \quad j = 1, 2, 3, D $$

and prove equivalence.

**PROPOSITION 3.** The friction laws (60) and (24) are equivalent.

**PROOF.** Let $0 < |\lambda_T| \le \mu_0 \lambda_N$. In this case the (relative) interior of $C_T$ is not empty, and the gradients of $\sigma_j$, $j = 1, 2, 3, D$ do not vanish. Therefore, (60) are the normal cone conditions for $N_{C_T}(\lambda_T)$. Since $N_{C_T}(\lambda_T) \equiv N_{D_T}(\lambda_T)$ equivalence of (60) and (25) holds. Now let $0 = |\lambda_T| \le \mu_0 \lambda_N$. In this case (60) reduces to $-\ddot{\mathbf{g}}_T = \sum_{i=1}^3 \mathbf{e}_i \kappa_i$ with $\sigma_i = \mu_0 \lambda_N \ge 0$, $\kappa_i \ge 0$, $\sigma_i \kappa_i = 0$ and $i = 1, 2, 3$. Assume that $\mu_0 \lambda_N > 0$. We obtain $\sigma_i > 0$, $\kappa_i = 0$ for $i = 1, 2, 3$. Thus, $\ddot{\mathbf{g}}_T = 0$, which is the second condition in (25) for $\lambda_T = 0$. Assume that $\mu_0 \lambda_N = 0$. We obtain $\sigma_i = 0$, $\kappa_i \ge 0$ for $i = 1, 2, 3$. Since the convex cone generated by $(\mathbf{e}_1, \mathbf{e}_2, \mathbf{e}_3)$ is $\mathbb{R}^2$ we have $\ddot{\mathbf{g}}_T \in \mathbb{R}^2$ which is the third condition in (25).

Using the tangential contact law from (60) together with the definitions of $\sigma_D$ and $\sigma_i$ from (58) and (44) we may now state the rolling contact problem in the form

$$ \begin{align} & \mathbf{M}\ddot{\mathbf{q}} - \mathbf{h} - \mathbf{w}_N\lambda_N - \mathbf{W}_T\lambda_T = 0 \notag \\ & \ddot{\mathbf{g}}_N = \mathbf{w}_N^\mathrm{T}\ddot{\mathbf{q}} + \bar{\mathbf{w}}_N \notag \\ & \ddot{\mathbf{g}}_T = \mathbf{W}_T^\mathrm{T}\ddot{\mathbf{q}} + \bar{\mathbf{w}}_T \notag \\ & \lambda_N \ge 0, \quad \ddot{\mathbf{g}}_N \ge 0, \quad \lambda_N \ddot{\mathbf{g}}_N = 0 \tag{61} \\ & \sigma_i = \mu_0 \lambda_N - \mathbf{e}_i^\mathrm{T}\lambda_T, \quad i = 1, 2, 3 \notag \\ & \sigma_D = \mu_0^2 \lambda_N^2 - |\lambda_T|^2 \notag \\ & -\ddot{\mathbf{g}}_T = \sum_{i=1}^3 \mathbf{e}_i \kappa_i + 2\mathbf{\lambda}_T \kappa_D \notag \\ & \sigma_j \ge 0, \quad \kappa_j \ge 0, \quad \sigma_j \kappa_j = 0, \quad i = 1, 2, 3, D \notag \end{align} $$

The derivation of the NCP is straightforward and will be performed in the same manner as in Section 8. Using the abbreviations (52) the fifth equation leads to the tangential forces (53) which are then substituted in the Newton/Euler equations and in the remaining friction saturations $\sigma_3$ and $\sigma_D$,

$$ \begin{align} \sigma_3 &= (\mu_0 - \mathbf{e}_3^\mathrm{T} \mathbf{I}_P^{-\mathrm{T}} \mu_P) \lambda_N + \mathbf{e}_3^\mathrm{T} \mathbf{I}_P^{-\mathrm{T}} \boldsymbol{\sigma}_P \notag \\ \sigma_D &= \mu_0^2 \lambda_N^2 - |\mathbf{I}_P^{-\mathrm{T}}(\mu_P \lambda_N - \boldsymbol{\sigma}_P)|^2 \tag{62} \end{align} $$

Solving the accelerations of the seventh equation in (61) for $\boldsymbol{\kappa}_p$,

$$ \boldsymbol{\kappa}_p = -\mathbf{I}_P^{-1} \ddot{\mathbf{g}}_T - \mathbf{I}_P^{-1} \mathbf{e}_3 \boldsymbol{\kappa}_3 - 2 \mathbf{I}_P^{-1} \mathbf{I}_P^{-\mathrm{T}} (\boldsymbol{\mu}_p \lambda_N - \boldsymbol{\sigma}_P) \boldsymbol{\kappa}_D, \quad (63) $$

we may rewrite the system of equations (61) to be

$$ \begin{align*} & \mathbf{M}\ddot{\mathbf{q}} - \mathbf{h} - (\mathbf{w}_N + \mathbf{W}_T\mathbf{I}_P^{-\mathrm{T}}\boldsymbol{\mu}_P)\lambda_N + \mathbf{W}_T\mathbf{I}_P^{-\mathrm{T}}\boldsymbol{\sigma}_P = 0 \\ & \ddot{\mathbf{g}}_N = \mathbf{w}_N^\mathrm{T}\ddot{\mathbf{q}} + \bar{\mathbf{w}}_N \\ & \boldsymbol{\kappa}_p = -\mathbf{I}_P^{-1}\mathbf{W}_T^{-\mathrm{T}}\ddot{\mathbf{q}} - \mathbf{I}_P^{-1}\bar{\mathbf{w}}_T - \mathbf{I}_P^{-1}\mathbf{e}_3\boldsymbol{\kappa}_3 - 2\mathbf{I}_P^{-1}\mathbf{I}_P^{-\mathrm{T}}(\boldsymbol{\mu}_p\lambda_N - \boldsymbol{\sigma}_P)\boldsymbol{\kappa}_D \\ & \sigma_3 = (\mu_0 - \mathbf{e}_3^\mathrm{T}\mathbf{I}_P^{-\mathrm{T}}\boldsymbol{\mu}_P)\lambda_N + \mathbf{e}_3^\mathrm{T}\mathbf{I}_P^{-\mathrm{T}}\boldsymbol{\sigma}_P \end{align*} $$
$$
\begin{align*}
\sigma_D &= \mu_0^2 \lambda_N^2 - |\boldsymbol{I}_P^{-\mathrm{T}}(\boldsymbol{\mu}_P \lambda_N - \boldsymbol{\sigma}_P)|^2 \tag{64} \\
\ddot{\lambda}_N &\ge 0, \quad \lambda_N \ge 0, \quad \ddot{\lambda}_N \lambda_N = 0 \\
\boldsymbol{\kappa}_P &\ge 0, \quad \boldsymbol{\sigma}_P \ge 0, \quad \boldsymbol{\kappa}_P^\mathrm{T} \boldsymbol{\sigma}_P = 0 \\
\sigma_3 &\ge 0, \quad \kappa_3 \ge 0, \quad \sigma_3 \kappa_3 = 0 \\
\sigma_D &\ge 0, \quad \kappa_D \ge 0, \quad \sigma_D \kappa_D = 0
\end{align*}
$$

Finally, the accelerations $\ddot{q}$ are eliminated which gives, together with the inequalities and complementarities
from the last four lines in (64), a NCP in standard form of dimension 5,

$$
\begin{pmatrix} \ddot{\lambda}_N \\ \boldsymbol{\kappa}_P \\ \sigma_3 \\ \sigma_D \end{pmatrix} = \begin{pmatrix} \boldsymbol{w}_N^\mathrm{T} \boldsymbol{M}^{-1} (\boldsymbol{w}_N + \boldsymbol{W}_P \boldsymbol{\mu}_P) & -\boldsymbol{w}_N^\mathrm{T} \boldsymbol{M}^{-1} \boldsymbol{W}_P & 0 & 0 \\ -\boldsymbol{W}_P^\mathrm{T} \boldsymbol{M}^{-1} (\boldsymbol{w}_N + \boldsymbol{W}_P \boldsymbol{\mu}_P) & \boldsymbol{W}_P^\mathrm{T} \boldsymbol{M}^{-1} \boldsymbol{W}_P & -\boldsymbol{I}_P^{-1} \boldsymbol{e}_3 & 0 \\ \boldsymbol{\mu}_0 - \boldsymbol{e}_3^\mathrm{T} \boldsymbol{I}_P^{-1} \boldsymbol{\mu}_P & \boldsymbol{e}_3^\mathrm{T} \boldsymbol{I}_P^{-1} & 0 & 0 \\ 0 & 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} \lambda_N \\ \boldsymbol{\sigma}_P \\ \boldsymbol{\kappa}_3 \\ \boldsymbol{\kappa}_D \end{pmatrix} \\
+ \begin{pmatrix} 0 \\ -2\boldsymbol{I}(\boldsymbol{\mu}_P \lambda_N - \boldsymbol{\sigma}_P) \boldsymbol{\kappa}_D \\ 0 \\ \boldsymbol{\mu}_0^2 \lambda_N^2 - \|\boldsymbol{\mu}_P \lambda_N - \boldsymbol{\sigma}_P\|_I^2 \end{pmatrix} + \begin{pmatrix} \boldsymbol{w}_N^\mathrm{T} \boldsymbol{M}^{-1} \boldsymbol{h} + \bar{\boldsymbol{w}}_N \\ -\boldsymbol{W}_P^\mathrm{T} \boldsymbol{M}^{-1} \boldsymbol{h} - \bar{\boldsymbol{w}}_P \\ 0 \\ 0 \end{pmatrix} \tag{65}
$$

Hereby, the abbreviations $W_P = W_T I_P^{-\mathrm{T}}$, $\bar{w}_P = I_P^{-1} \bar{w}_T$, and $I = I_P^{-1} I_P^{-\mathrm{T}}$ have been used, where $I$ is a symmetric and positive definite $2 \times 2$ matrix. The NCP (65) is of the form

$$
y = f(x), \quad y \ge 0, \quad x \ge 0, \quad y^T x = 0 \tag{66}
$$

where $f(x)$ in (65) is represented by a linear part, a nonlinear part, and a constant,

$$
f(x) = Ax + g(x) + c. \tag{67}
$$

Since $f(x)$ is continuous and differentiable we may finally derive its Jacobian. With

$$
D_x f = A + D_x g \tag{68}
$$

and $D_x g$ explicitly depending on $\lambda_N$, $\sigma_P$ and $\kappa_D$, we obtain

$$
D_x g = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 - 2I\boldsymbol{\mu}_P\boldsymbol{\kappa}_D & 2I\boldsymbol{\kappa}_D & 0 & -2I\boldsymbol{\lambda}_P \\ 0 & 0 & 0 & 0 \\ 2\boldsymbol{\mu}_0^2\boldsymbol{\lambda}_N - 2\boldsymbol{\lambda}_P^\mathrm{T}I^\mathrm{T}\boldsymbol{\mu}_P & 2\boldsymbol{\lambda}_P^\mathrm{T}I^\mathrm{T} & 0 & 0 \end{pmatrix} \quad (69)
$$

with

$$
\lambda_P(\lambda_N, \boldsymbol{\sigma}_P) = \boldsymbol{\mu}_P \lambda_N - \boldsymbol{\sigma}_P \quad \text{and} \quad I = I^T
$$

Finally, we draw attention to the similar structure of the LCP matrix from (38) and the Jacobian (68). In both
cases the second column of the corresponding matrix, multiplied by the friction coefficients μ<sub>P</sub>, appears in the
first column. This indicates that the tangential forces λ<sub>T</sub> have been eliminated by the use of certain friction
saturation<sub>s</sub>s.

10. Conclusion

In this paper, we have presented different formulations of the unilateral contact problem with Coulomb
friction in dynamic multibody systems. By using the same methodical approach we have formulated pyramidal
and conic friction laws in terms of linear and nonlinear complementarity problems, respectively. As an
important result, these formulations have been obtained in standard form and, moreover, in a way that also
dependent contact constraints can be treated. We have observed similar structures arising in the Jacobians of the
linear and the nonlinear complementarity problem, due to the fact that the formulation processes follow the same rules.

As we have already pointed out, the treatment of the Coulomb friction law as presented in this paper should be understood as an example. The way in which standard LCP formulations may be obtained when approximating arbitrary set-valued force interactions by polyhedral convex sets is general. The NCP formulation was simply based on the fact that the intersection of convex sets remains convex, and that additional inequalities may be added if convenient. This method works if the normal cone can be generated by a finite number of given directions.

We have not discussed any numerical aspects. The friction pyramid model as presented in Section 7 has been used for about two years and has been applied, for example, to the friction problem at the legs and in the gears of a tube crawling robot as well as to the investigation of the wheel/rail contact of the undercarriage of a roller coaster. The NCP formulation, however, is new and has not yet been tested.

## References

[1] P. Alart and A. Curnier, A mixed formulation for frictional contact problems prone to Newton like solution methods, Comput. Methods Appl. Mech. Engrg. 92(3) (1991) 353-357.

[2] Ch. Glocker, The principles of d'Alembert, Jourdain, and Gauss in nonsmooth dynamics. Part I: Scleronomic multibody systems, ZAMM 78(1) (1998) 21-37.

[3] Ch. Glocker, Dynamik von Starrkörpersystemen mit Reibung und Stößen, Fortschrittberichte VDI, Reihe 18, Nr. 182 (VDI-Verlag, Düsseldorf, 1995).

[4] A. Klarbring. A mathematical programming approach to three-dimensional contact problems with friction, Comput. Methods Appl. Mech. Engrg. 58 (1986) 175-200.

[5] A. Klarbring and G. Björkman, A mathematical programming approach to contact problems with friction and varying contact surface, Comput. Struct. 30(5) (1988) 1185-1198.

[6] P. Lötstedt, Numerical simulation of time-dependent contact and friction problems in rigid body mechanics, Technical Report TRITA-NA-8214, Dept. of Numerical Analysis and Computing Science (The Royal Institute of Technology, Sweden, 1982).

[7] Th. Meitinger, Dynamik automatisierter Montageprozesse, Fortschrittberichte VDI, Reihe 2, Nr. 476 (VDI-Verlag, Düsseldorf, 1998).

[8] J.J. Moreau, Unilateral contact and dry friction in finite freedom dynamics, in: J.J. Moreau and P.D. Panagiotopoulos, eds., CISM Courses and Lectures No. 302, Nonsmooth Mechanics and Applications (Springer-Verlag, Wien, New York, 1988) 1-82.

[9] P.D. Panagiotopoulos, Hemivariational Inequalities—Applications in Mechanics and Engineering (Springer-Verlag, Berlin, Heidelberg, New York, 1993).

[10] J.S. Pang and J.C. Trinkle, Complementarity formulations and existence of solutions of dynamic multi-rigid-body contact problems with Coulomb friction. Math. Prog. 30 (1996) 199-226.

[11] F. Pfeiffer and Ch. Glocker, Multibody Danamics with Unilateral Contacts (John Wiley and Sons, New York, 1996).

[12] R.T. Rockafellar, Convex Analysis (Princeton University Press, Princeton, New Jersey, 1970).

[13] Th. Rossmann and Ch. Glocker, Efficient algorithms for non-smooth dynamics, in: Proc. of the ASME International Mechanical Engineering Congress and Exposition, Dallas, Texas, Nov. 16-21, 1997, The American Society of Mechanical Engineers, 1997.

[14] M. Wösle, Dynamik von räumlichen strukturvarianten Starrkörpersystemen, Fortschrittberichte VDI, Reihe 18, Nr. 213 (VDI-Verlag, Düsseldorf, 1997).