Formulatin of spatial contact situations in rigid multibody systems

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.


I. 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 occuring 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 0045-7825/99/$ -see front matter © 1999 Elsevier Science BY. All rights reserved. PII: S0045-7825(98)0038 1-8 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.

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 Vpi, ~Qi (i = 1, 2) where Vp~ denotes the velocity of the body-fixed point P of body i and ~ 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 rp,~ = rp,=~(q, ~, 71i), see Fig. 1 [7].
Note that rp,. only depends on the surface parameters ~: and ?7 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 := 0¢ rp, and t := O,Trp, -in order to obtain ~rpE i V~P 11  (1) where we agree that (s ~, 7/) are ordered such that n becomes the outward normal. The (absolute) changes in time of the contour vector fez and the trihedral (t/, tJ, t/) are i e z = 1 2 × r p z + R ~ with R = O~r e , = = ( s , t ) where ~" = (~, ~/)T and (k, K) is any of (u, U), (v, V), (n, N). The velocity v z of a point moving on the surface in terms of the rigid body motion (re, 12) is thus given by v , = Vp + ip, which yields Here, v c denotes the rigid body velocity portion of v,.. Analogously, we may find the acceleration a c = v c by differentiation of the second equation in (3) Finally, we express v c and a c by the generalized magnitudes q and ~, Vc = J c 4 + i~, a c = J c q + ic (5) where tc is given by i c = J c q + tc, and J c is the Jacobian of the contact point as it will turn out later.

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 rt) : = ro,_ 2 -ro~ ~ , (6) cf. [3], where ro,=i are vectors pointing from a common inertially fixed point O to the surfaces under consideration.
Next, the two tangent planes spanned by (u i, vi) are adjusted parallel to each other such that the distance vector becomes parallel to each of the normals. The surface parameters (~'*, ~'*) corresponding with that situation solve the four nonlinear equations and are called the contact parameters; the corresponding points at the surfaces given by rp,i(~'*) are called the contact points. When the bodies are moving the conditions (7) for the contact points should not change [7,3], thus must hold which constitutes a system of four linear equations in ~* when substituting (t/i, t/i, t~i) from (2) and After the solution (~'*, ~'*) of (7) has been found the distance gN of the contact points can be calculated as T gN = n j rl) (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 (u~, v~, n~) are then given by where it can be easily verified that differentiation of (9) with respect to time results in (10) for (K, k~ ) : (X, n~). The temporal changes of the relative velocities gK in (10) are given by where ac~ is known from (4), k" 1 may be taken from (2), and ~* involved in ac~ and k'l is the solution of (8).
Substituting Vc~ and acz from (5) into (10), (11) yields a representation of the relative velocities and accelerations in the form T -T . .
gK : W K q q-I~K ' g K : W K q q-WK (12) = -J c~) k~, WK = (ic2 I£CI )Tkl, and ~ = w Kq "Jr-fie x. Finally we, arrange the two tangential relative velocities in the vector gT, i.e. gT : = (g;v'gv) which gives us the sliding direction of the contact points in the common tangent plane. With the abbreviations W r := (w v, W v ), w ~ -rr = ( g'v, v~v ), g' TT = ( g'v, V~'v ) we rewrite (12) in the form g r = W r q + f i r which will serve as the contact kinematic equations in standard form in the sequel.

Kinetics
The dynamics of a multibody system with n bodies and one contact is described by the Newton-Euler equations in the form -J c~F c l -J c 2 F c 2 = O, i~l see e.g. [3,1 1]. The terms Pi and L i denote the momentum and the moment of momentum of body i, and F i and Mi are the applied forces and moments, respectively. The contact forces acting at the contact points C~ and C 2 are denoted by Fcl and Fc2. In order to perform a projection into the configuration space all terms have to be premultiplied by the corresponding Jacobians. In particular, JCl and Jc2 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 j, v~, n~) admits, together with the cutting principle, a unique representation in the form Fc2 = -Fc~ : -: n l A N + u j z~ U + Vj/~ v Putting (15) into (14) and using the abbreviations w K below we arrive with : ( J c2 -Jc,)Tkl and W r = (Wv, Wv) from Eq. (12) and (16) where the vector ~t r of the tangential contact forces is defined to be ATr = (A U, Av). The sum in (14) is represented by the term (M~/-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.

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 gN being the distance of the contact points and A N the corresponding normal force one might express impenetrability of rigid bodies by the Signorini-Fichera condition [1, 3,5,6,9,10] This condition may be transferred to the velocity level and, in a second step, to the acceleration level [3,11,10] gN>O : 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 gN = 0, gN < 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 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 q(t) = u+(t) = u-(t) since q(t) is differentiable. This also applies to gN(t) through (10). Similarly, the accelerations 4/(t) (and ~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] = 0 laTI 0AN (19) where /x denotes the coefficient of friction which depends smoothly on gr, and /z o = # ( g r = 0). Similarly we obtain from (19) the friction law on the acceleration level in the form ]10,14] with e being the unit vector of the sliding direction, i.e. e = gr/IgTI. 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 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 gT = 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 q. In the case of closed contacts this is usually done by two steps [3,10]: First, the accelerations 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 g and the contact forces A. When this solution has been found, A is re-substituted into the Newton/Euler equations (16) which permits finally the calculation of 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.

Sliding contacts
In this section we derive the linear complementarity problem for a closed contact (gN = 0, gN : 0) which is sliding (gr ~ 0). This situation is described by the Newton/Euler equations (16), the normal contact law for gN = O, g N = 0 from the third line in (18), and the tangential contact law for sliding contacts gT ~0 from the first equation in (20). Additionally as kinematical relation the third equation in (13) is needed, since g u is involved in the normal contact law. Altogether these equations write Using the fourth equation A r is substituted into the first one which is then solved for q'. Substituting ~/ into the second equation we obtain the LCP for the sliding contact [10,14,13] .
Obviously, friction may cause the LCP-matrix (for one contact being a scalar) to become indefinite.

Friction pyramids for rolling contacts
Roiling contacts are characterized by vanishing tangential velocities g r = 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 gT = 0 in a more suitable form which will also be used in the following sections. Let DT :--{at I laTI maN} (23) be the (convex) set of the admissible tangential contact forces. The tangential contact law is then given by see e.g. [5,9], where NDT(AT) denotes the normal cone to D T at A T, Following [5] we may approximate the set D r by a polytope generated by the intersection of n closed half spaces Cri, each of which being represented by an affine inequality o;,(Ar), 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 r and C r. We call the ~ the friction saturation with respect to Cr~ because it determines the distance between a Fig. 3. Approximation of the friction disk by polytopes.
given tangential force '#~T and the boundary of CTi. It therefore describes the amount of forces which might additionally act on the contact in the direction ei until sliding starts. In the following we will derive a system of complementarity conditions that describe the friction law (27). The half-spaces Cri defined in (26) have at least the point 0 in common since ~(0) i> 0 because A N t> 0, /%/> 0. Thus, the normal cone to C r = f-I'/_~ Cr~ is given by i-I cf. [12], and the contact law (27)  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 -JW has to be avoided. For multi-contact problems the matrix W consists of a selection of constraint vectors w m, wu~, Wv~ 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 comers (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. 8 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 gT and A r. Thus, one has to eliminate A T from the Newton/Euler equations and to reformulate the tangential acceleration gT 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 WTM -~W. Following ]13] we choose an indexing i in the last three equations in (31)  with/x r :=/zp ~-/z R = (#o /z0) T. Exactly this formulation but expanded to multi-contact configurations is used in [10].
In the case of one-dimensional friction the directions el T = (1 0) and e T = (-1 0) reduce to + 1 and -1, respectively, whereas e 2 and e 4 are no longer needed. The LCP (38) then writes

gN) i / w~M-'(wN+wu~,,)
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 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 A E ~6 would enter the Newton/Euler equations, and the elimination of A would require six linearly independent directions e~ to be arranged in Ip from (32).

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 A r and the tangential relative accelerations gT is lost. This lack of consistence 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 NCE 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) (37) where the tangential forces A r had been eliminated from the Newton/Euler equations. This, however, requires at least two equations which can be solved for A r, preferably being linear. We therefore introduce two additional linear inequalities on the tangential forces, not changing D y when intersection is performed: Let Instead of (48) we introduce a contact law which is expressed by an equality. We set 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 0-t~ and ~ from (43)  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/zoA m = 0 the friction law (47) reduces to g r E ~2. In order to express any gT E ~2 by a positive linear combination S, e;r, at least three directions e; are necessary. The scalars x~/> 0 enter the NCP as independent variables which makes, together with the normal contact law, an overall dimension of at least 4.

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 , defined in (49) may jump at points h~ r = 0, i.e. if [d~p/~ u ~-O'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 r by a parabolic function or, which is differentiable everywhere but we note that the gradient of or;, vanishes for A r = 0. The subdifferential O¢r,(/tT-) = ~7(T;:,(~[7.) is therefore not capable of describing the normal cone to D r at ,A 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 (q)_~ C7.; at A r = 0 are identical and may be represented by using the three directions e;, i = 1,2, 3 that are defined in Fig. 4. The set of the admissible tangential forces C r and the corresponding normal cone NeT(l[ r) may be expressed in the same manner as in (45) and (46), and the tangential contact law (24) becomes -g r E Nc, ( Ar) 13 since C r -~D r and thus Nc~(A r) =--N.,(Ar). Since the gradient of o-o vanishes for A r = 0, the normal cone conditions as in (48)  Ei= j e~Ki with o], =/~oAu ~>0, t~ ~0 , ~K z = 0 and i = 1,2,3. Assume that /xoAN > 0 . We obtain ~> 0 , K~=0 for i = 1, 2, 3. Thus, gT = 0, which is the second condition in (25) for A T = 0. Assume that #oaN = 0. We obtain = 0, K~ ~> 0 for i = 1,2, 3. Since the convex cone generated by (e~, e 2, e3) is ~2 we have g r ~ I~-which is the third condition in (25).
Using the tangential contact law from (60)  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 lap, appears in the first column. This indicates that the tangential forces A r have been eliminated by the use of certain friction saturations trp.

I0. 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 pyramidic 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 n o n l i n e a r c o m p l e m e n t a r i t y problem, due to the fact that the formulation processes follow the same rules.
As we have already pointed out, the treatment of the C o u l o m b friction law as presented in this paper should be u n d e r s t o o d as an example. T h e way in which standard L C P formulations m a y be obtained w h e n a p p r o x i m a t i n g arbitrary set-valued force interactions by polyhedral c o n v e x sets is general. The N C P formulation was simply based on the fact that the intersection of c o n v e x sets remains convex, and that additional inequalities m a y be added if convenient. This m e t h o d works if the normal cone can be generated by a finite n u m b e r of given directions.
We have not discussed any numerical aspects. T h e friction p y r a m i d model as presented in Section 7 has been used for about two years and has been applied, for example, to the friction p r o b l e m at the legs and in the gears of a tube crawling robot as well as to the investigation of the w h e e l / r a i l contact of the undercarriage of a roller coaster. The N C P formulation, however, is new and has not yet been tested.