This technical brief revisits the method outlined in Tasora and Anitescu 2011 [“A Matrix-Free Cone Complementarity Approach for Solving Large-Scale, Nonsmooth, Rigid Body Dynamics,” Comput. Methods Appl. Mech. Eng., 200(5–8), pp. 439–453], which was introduced to solve the rigid multibody dynamics problem in the presence of friction and contact. The discretized equations of motion obtained here are identical to the ones in Tasora and Anitescu 2011 [“A Matrix-Free Cone Complementarity Approach for Solving Large-Scale, Nonsmooth, Rigid Body Dynamics,” Comput. Methods Appl. Mech. Eng., 200(5–8), pp. 439–453]; what is different is the process of obtaining these equations. Instead of using maximum dissipation conditions as the basis for the Coulomb friction model, the approach detailed uses complementarity conditions that combine with contact unilateral constraints to augment the classical index-3 differential algebraic equations of multibody dynamics. The resulting set of differential, algebraic, and complementarity equations is relaxed after time discretization to a cone complementarity problem (CCP) whose solution represents the first-order optimality condition of a quadratic program with conic constraints. The method discussed herein has proven reliable in handling large frictional contact problems. Recently, it has been used with promising results in fluid–solid interaction applications. Alas, this solution is not perfect, and it is hoped that the detailed account provided herein will serve as a starting point for future improvements.
Notation: Problem setup
With the ground assigned by convention index 0, assume two bodies of index A and B, are in contact. As in Fig. 1, let i identify this contact event. A collision detection process produces the point of contact P, a signed distance function , and a set of three orthonormal vectors: , and . By convention, the normal vector is oriented from the body of lower index to the body of higher index.

Bodies A and B in contact; a local reference frame {} is generated at the contact point based on contact detection information. The contact point is located in the centroidal and principal reference frames via the and vectors.
Any two bodies that are closer than a prescribed are considered to produce an active contact event. The gap function is negative if the two bodies share more than one point; it is zero, if they share one point; it is greater than zero, if they share no point. The geometry of the bodies is assumed to be convex in a neighborhood of the contact area.
and therefore the generalized force associated with the frictional contact force is .
where is the constant mass matrix, is the generalized applied and Coriolis forces, is the constraint reaction force associated with bilateral constraints, and is the frictional contact force associated with the presence of contact events. In terms of notation , and .
represents the relative velocity at the contact point between the two bodies expressed in the local reference frame {}.
Frictional Contact Model
where is the friction coefficient.
The Discretized Equations of Motion
In Eq. (5), ; ; and, . Moreover, , and . In Eq. (5c), and . Finally, in Eq. (5d), is the magnitude of the tangential velocity at the point of contact, is the magnitude of the friction impulse associated with the contact event i, and represents the set of active contact events at time and in configuration .
There are two notable aspects tied to the discretization of the differential variational inclusion problem above.
- (1)
The bilateral kinematic constraint equations, see Eq. (3a), are not used. Instead, we use the velocity-level set of kinematic constraints in Eq. (3b). However, the latter are modified in two respects. First, to account for violation in satisfying the kinematic constraints at position level, a “stabilization term” is considered in the discretized form of the equation [3–6]. Second, since the method is half-implicit, we chose to evaluate the partial time derivate gt in the configuration .
- (2)
When discretizing the expression of the signed gap function in Eq. (5d), there are two approximations involved in the process. First, the complementarity conditions are imposed using an approximation of the signed gap function at . That is, . Given that , the condition that the gap between two bodies at the next time step can be at most zero translates into the inequality
Second, in order to render the nonlinear complementarity problem obtained upon time discretization tractable, a relaxation of the approximation above is introduced via the term . As shown in Sec. 4, this change enables one to pose the problem in Eq. (5) as a cone complementarity problem (CCP).
The “stabilization term” in Eqs. (5c) and (5d) leads to a numerical penalty force that penalizes the constraint violation. This penalty force is devoid of physical meaning; for instance, it is not tied to the stiffness of the bodies participating in the contact event. To better understand the implications of the stabilization and relaxation performed, we focus next on the unilateral constraints; i.e., bullet 2 above. Note that Eq. (6) effectively places a lower bound on the new velocity. Indeed, to enforce nonpenetration, the velocity should be at least larger than . Note that a negative gap represents penetration and the numerical method will produce a solution that seeks to enforce a non-negative gap . Could then the gap ever become negative in this method? There are three situations when this can happen: (i) the choice of initial conditions leads to bodies penetrating; (ii) discretization errors lead to a negative gap—the approximation that led to Eq. (6) is a first-order Taylor expansion and large values of impact the quality of the approximation; and (iii) bodies that come at each other very fast and are on a penetration course. In all these three cases though, the solution described seeks to correct the wrong, i.e., move from a negative gap to a non-negative gap, in one time step, i.e., from to . This can be a problem under (i), since one can have large penetrations that will lead to large corrections forces needed to correct the penetration in one time step. This is particularly problematic when the step size , which divides the gap , is very small. Then, in one small time step, the method attempts to correct the penetration gap, which could lead to large normal forces. In practice, one might place an upper bound on the value of the resulting normal velocity, effectively capping the value of the normal force that attempts to address the gap violation in one step. For (ii), large negative gaps are unlikely. For (iii), a scenario that can lead to large penetration is one in which two bodies moving very fast toward each other and are on a collision course are not colliding yet but are about to. At the next time step, the penetration will be large, and as such a large normal force will subsequently correct the large penetration in one step. This situation can be mitigated by selecting the active set in Eq. (5d) to include not only bodies that are in contact at the current time step, but also bodies that are close to each other at , i.e., choosing a strictly positive .
The discussion covering the penetration cases (i) through (iii) above was built around the gap value. There is a clear intuition behind this quantity that facilitates a better grasp of the arguments put forth. However, Eq. (5d) makes it clear that the numerical method works with the modified gap . Qualitatively, little changes when replacing by . Quantitatively though, the use of might lead to a subtle numerical artifact. If relative to the size of , the term becomes large, then its size will play a role in the economy of the numerical solution. By and large, μi is small and so is . The quantity represents the tangential relative velocity at the contact point. If it is large enough, then the net effect of is to make the penetration more negative, i.e., deeper penetration. Indeed, if one bowling ball would slide fast without rolling on the floor and there would be no physical penetration, the numerical solution nevertheless registers a numerical gap equal in size to . As such, the numerical scheme will produce an excess normal force whose purpose is to compensate for this nonphysical gap. This artifact has been demonstrated/discussed in Ref. [7]; see Figs. 2 and 3 therein. Indeed, when the bowling ball slides fast on the floor, there is a vertical force caused by the non-negligible sinkage . Once the value of the term becomes small, so does the values of this force. Note that once the bowling ball rolls without slip, this “numerical artifact” force vanishes. Moreover, this force is zero if ; i.e., in the absence of friction.
Finally, note that the solution methodology outlined here leads to an inelastic treatment of impact. The use of dictates that two bodies that are δK-close to each other lead to a contact event. As such, the complementarity condition will enforce the condition that the numerical solution produce a velocity that leads to no penetration at . Bar small penetrations that are due to numerical approximations, the numerical solution will maintain this zero gap until the forces applied lead to a “lift-off” condition at which point the two bodies separate. Ways to change this behavior from inelastic to elastic are not discussed here—the focus of this contribution is exclusively on handling of frictional contact, which is ubiquitous. To the best of our knowledge, deriving a general purpose elastic impact scheme for rigid bodies of arbitrary shape that experience simultaneous impacts remains an open problem. One approach that is promising is described in Ref. [8], albeit therein the authors resort to a decoupling of the impact computation from contact computation, as well as the computation of the friction force from the normal force.
The Cone Complementarity Problem
Posing the Problem
The “Unilateral Constraints” Component.
and therefore .
The “Bilateral Constraints” Component.
which is the “bilateral constraint” CCP analog of the condition in Eq. (7).
Reformulating the CCP.
Let , and . Therefore, , and .
where , and .
The Quadratic Problem Angle
As , using the Cauchy–Schwartz inequality, , which proves that as far as is concerned, . In other words, if satisfies the KKT conditions in Eq. (12), it is also a solution of the CCP problem in Eq. (7).
A similar result can be obtained in the bilateral constraints case. Indeed, in this case . Using the last complementarity condition in Eq. (12), one has that whenever , necessarily bj = 0. In other words, we have that , which indicates that a that satisfies the first KKT conditions in Eq. (12) is a solution of the CCP problem in Eq. (8).
Conclusions and Future Work
A differential inclusion approach is used to formulate the rigid multibody dynamics problem in the presence of mutual contact and friction. Its salient attribute is reliance on complementarity conditions to pose both the contact non-penetration condition and the Coulomb dry friction model. Time discretization of the resulting differential, algebraic, and complementarity equations yields a nonlinear complementarity problem. In the discretization process, the unilateral and bilateral kinematic constraints are imposed at the velocity level. Drift in the position-level constraints is prevented in Eq. (5) via “stabilization terms.” Moreover, the nonpenetration unilateral constraint, formulated at the velocity level, is further modified via a “relaxation term” to morph what would otherwise be a nonlinear complementarity problem into a cone complementarity problem. The latter has a solution that is produced by solving of a convex quadratic optimization problem with conic constraints. The same numerical method is proposed in Ref. [1] but following a different set of intermediate steps. As such, this contribution does not focus on the numerical method, but rather on the steps to obtain this numerical method. It is hoped that this presentation is detailed enough to reveal the strengths and weaknesses of the solution methodology.
There are three aspects in which the method described can be improved. First, a better approach would enforce the unilateral and bilateral constraints at the position level to impose a tight and numerically robust control on constraint drift. Second, the relaxation in Eq. (5d) was shown to lead to “lift-off forces” at high sliding speed, an artifact that can introduce noise in the solution. Third, the rigid body assumption leads to scenarios in which, due to the presence of redundant constraints, the matrix in Eq. (11) is symmetric positive semi-definite. As such, a solution of the convex optimization problem with conic constraints, while global, is not unique. There are early indications that these three limitations can be addressed. A discussion of this issue falls outside the scope of this document.
Acknowledgment
This research was supported in part by National Science Foundation grant CMMI–1635004. We would like to thank Michał Kwarta and two anonymous reviewers for their suggestions, which improved the quality of this contribution.
Funding Data
National Science Foundation (Grant No. 1635004).
Nomenclature
The list below summarizes the meaning of the main symbols used in the manuscript. The list is not exhaustive; it includes only the more important symbols that were used beyond the immediate location where defined
- A, B =
dummy indexes for arbitrary bodies
- =
active contact set; two bodies contribute a contact event to this set if they are closer than from each other
- =
the set of active contact events at time and in configuration
- AA =
rotation matrix associated with body A
- =
composite velocity associated with bilateral constraint j
- =
cones associated with bilateral constraint j
- =
system level friction cone,
- =
system level polar cone,
- di =
composite velocity;
- =
notation for
- =
projection matrix used to generate the set of generalized forces induced by the contact events present in the system (unilateral constraints)
- =
set of four Euler parameters associated with body i. Provides orientation with respect to a global reference frame
- =
notation for impulse
- =
the set of generalized applied and Coriolis forces
- =
notation for
- =
notation for
- =
function expression that defines the jth bilateral constraint active in the system at time t, namely
- =
projection matrix used to generate the set of generalized forces induced by the bilateral constraints
- =
notation for
- =
transformation matrix that links the time derivative of the Euler parameters to the angular velocity
- =
constant, system level, mass matrix
- nb =
number of rigid bodies in the system
- =
number of bilateral constraints present in the system
- =
unit normal vector at point of contact; oriented from the body of lower index to the body of higher index
; combine to form the frictional contact vector
- =
number of contact events present in the system at a given time (explicit dependency on time dropped for convenience)
- =
quadratic term matrix, optimization problem; computed as
- =
system level projection matrix, defined as
- =
coefficient of the linear term in the CCP
- =
set of generalized coordinates associated with the nb bodies in the system
- =
Cartesian space location of body i
- =
vector that for contact event i provides the location of the contact point expressed in the local reference frame associated with body B
- =
time associated with the integration time step l; also,
- =
two unit normal vectors at point of contact; span the tangent contact plane; together with form a orthogonal reference frame
- =
components of the relative velocity at the contact point between the two bodies involved in contact event i; components associated with the local reference frame {}
- =
magnitude of the tangential velocity at the point of contact; notation for
- =
collection of all body velocities
- =
angular velocity of body i expressed in its centroidal and principal reference frame
- δK =
user defined threshold value that defines when two bodies yield an active contact event
- =
signed gap function: negative if bodies penetrate, positive if they don't touch, zero if one point of contact
- =
Lagrange multipliers associated with the three unknown components of the frictional contact force for contact event i. Components expressed in conjunction with the reference frame
- =
set of Lagrange multipliers associated with the unilateral constraints present in the system
- =
notation for frictional contact impulse
- =
magnitude of the friction impulse associated with the contact event i; notation for
- =
numerical integration time step
- =
friction cone associated with contact event i
- =
polar cone associated with the friction cone
- μi =
friction coefficient associated with contact event i
- =
set of Lagrange multipliers associated with the bilateral constraints
- =
notation for bilateral constraint reaction impulse
- =
notation for ; i.e., gap for contact event i at time