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

The time-evolution of a collection of nb rigid bodies interacting through friction and contact is described herein using Cartesian coordinates. The array of generalized coordinates $q=[r1T,ε1T,…,rnbT,εnbT]T∈ℝ7nb$, and its time derivative $q˙=[r˙1T,ε˙1T,…,r˙nbT,ε˙nbT]T∈ℝ7nb$, is used to represent the state of the system, where for body j, $1≤j≤nb$, rj and $εj$ are the absolute position of the center of mass and the body orientation Euler parameters, respectively. The time derivative of the Euler parameters $ε˙$ can be replaced with a different set of unknowns, i.e., the angular velocity in local coordinates $ω¯$. The unknown velocity $v=[r˙1T,ω¯1T,…,r˙nbT,ω¯nbT]T∈ℝ6nb$ is tied to $q˙$ via a linear transformation [1]
$q˙=L(q)v$
(1)

With the ground assigned by convention index 0, assume two bodies of index A and B, $0≤A 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 $Φi$, and a set of three orthonormal vectors: $ni, ui$, and $wi$. By convention, the normal vector $ni$ is oriented from the body of lower index to the body of higher index.

Any two bodies that are closer than a prescribed $δK≥0$ are considered to produce an active contact event. The gap function $Φi$ 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.

In each configuration $q(t)$, the collection of $NK$ contact events is denoted by $A(q(t),δK)$. The rotation matrices associated with bodies A and B are $AA=AA(ε1(t))$ and $AB=AB(ε1(t))$, respectively. The force acting on body B at point P is $Fi,B=γ̂i,nni+γ̂i,uui+γ̂i,wwi≡Aiγ̂i$. The location of point P on body B is $rBP=rB+ABs¯i,B$ and its virtual displacement is $δrBP=δrB−ABs¯̃i,Bδπ¯B$, where rB is the location of the center of mass of body B; a three-dimensional vector quantity with an over-bar, such as $s¯i,B$, indicates a representation of a geometric vector in the local (body-attached) centroidal and principal reference frame; the tilde operator produces the skew symmetric matrix associated with the vector it is used in conjunction with; and, the vector $δπ¯B$ is the virtual rotation associated with body B. The virtual work associated with the frictional contact force $Fi,B$ is
$δWi,B=[δrBP]T Fi,B=δrBTAiγ̂i+δπ¯BT s¯̃i,BABTAiγ̂i$
Similarly, the virtual work for body A is
$δWi,A=[δrAP]T (−Fi,B)=−δrATAiγ̂i−δπ¯ATs¯̃i,AAATAiγ̂i$
The virtual work that the presence of the frictional contact force $Fi,B$ imparts is
$δWi=Wi,A+Wi,B=δrTDiγ̂i$
where
$δr=[δr1δπ¯1⋮δrnbδπ¯nb]∈ℝ6nb Di=[03×3⋮03×3−Ai−s¯̃i,AAATAi03×3⋮03×3Ais¯̃i,BABTAi03×3⋮03×3]∈ℝ6nb×3$

and therefore the generalized force associated with the frictional contact force is $Diγi$.

The equations of motion assume the form [1]
$Mv˙=F(q,v,t)+Gλ̂B+Dγ̂K$
(2)

where $M=diag{m1I3×3,J¯1,…,mnbI3×3,J¯nb}$ is the constant mass matrix, $F(q,v,t)$ is the generalized applied and Coriolis forces, $Gλ̂B$ is the constraint reaction force associated with bilateral constraints, and $Dγ̂K$ is the frictional contact force associated with the presence of $NK$ contact events. In terms of notation $D≡[D1…DNK]∈ℝ6nb×3NK, G≡[G1…GNB]∈ℝ6nb×NB, γ̂K≡[γ̂1T,…,γ̂NKT]T∈ℝ3NK$, and $λ̂B≡[λ̂1T,…,λ̂NBT]T∈ℝNB$.

The bilateral constraint reaction forces are associated with the presence of bilateral constraints. These can be holonomic or nonholonomic; without loss of generality, they are assumed here to be holonomic and expressed as
$gj(q,t)=0$
(3a)
where $1≤j≤NB$. Their time derivative yields
$g˙j(q,v,t)≡GjTv+∂gj∂t=0$
(3b)
Finally
$DiTv=−AiTr˙A+AiTr˙B+AiTAAs¯̃i,Aω¯A−AiTABs¯̃i,Bω¯B=AiT(r˙B+ABω¯̃B s¯̃i,A−r˙A−AAω¯̃A s¯̃i,A)=AiT(r˙i,B−r˙i,A)≡[vi,n, vi,u, vi,w]T$

represents the relative velocity at the contact point between the two bodies expressed in the local reference frame {$ni, ui, wi$}.

## Frictional Contact Model

Herein, all impacts are considered inelastic; i.e., the restitution coefficient is zero. A contact event i is captured as a nonpenetration condition posed as a complementarity equation. Friction factors in via the Coulomb dry friction model. Accordingly, for contact event i, the frictional contact model requires that the following three conditions hold simultaneously:
$0≤Φi⊥γ̂i,n≥0$
(4a)

$0≤vi,u2+vi,w2⊥(μiγ̂i,n−γ̂i,u2+γ̂i,w2)≥0$
(4b)

$∃ αi≥0:{vi,u=−αiγ̂i,uvi,w=−αiγ̂i,w ,$
(4c)

where $μi≥0$ is the friction coefficient.

## The Discretized Equations of Motion

The discretization scheme adopted is a half-implicit symplectic Euler method; see Ref. [2]. It is used to discretize the kinematic differential equations in Eq. (1), the Newton–Euler equations of motion in Eq. (2), and the Coulomb friction model stated in Eq. (4). This yields the following nonlinear complementarity problem:
$q(l+1)︷generalized positions=q(l)+Δt︷step sizeL(q(l))︸velocity transformation matrixv(l+1)$
(5a)

$M(v(l+1)︷gen. speeds−v(l))=f(l)+G(l)λB,(l+1)︷reaction impulse+D(l)γK,(l+1)︸frict. contact impulse$
(5b)

$0=1Δt g(l)︸stabilization term+G(l),Tv(l+1)+gt(l)$
(5c)
and, $∀ i∈AδK(l)$,
${0≤γi,n(l+1) ⊥ ( 1ΔtΦi(l)︷stabilization term+vi,n(l+1)−μivi,T(l+1)︷relaxation term )≥00≤vi,T(l+1) ⊥ ( μiγi,n(l+1)−γi,F(l+1) )≥0∃ αi≥0:{vi,u(l+1)=−αiγi,u(l+1)vi,w(l+1)=−αiγi,w(l+1)$
(5d)

In Eq. (5), $f(l)≡Δt F(t(l),q(l),v(l))$; $γK,(l+1)≡Δt γ̂K,(l+1)$; and, $λB,(l+1)≡Δt λ̂B,(l+1)$. Moreover, $G(l)≡G(q(l),t(l))$, and $D(l)≡D(q(l),t(l))$. In Eq. (5c), $g(l)≡g(q(l),t(l))$ and $gt(l)≡(∂g(q(l),t(l))/∂t)$. Finally, in Eq. (5d), $Φi(l)≡Φi(q(l)), vi,T(l+1)≡(vi,u(l+1))2+(vi,w(l+1))2$ is the magnitude of the tangential velocity at the point of contact, $γi,F(l+1)≡(γi,u(l+1))2+(γi,w(l+1))2$ is the magnitude of the friction impulse associated with the contact event i, and $AδK(l)$ represents the set of active contact events $A(q(t),δK)$ at time $t(l)$ and in configuration $q(l)$.

There are two notable aspects tied to the discretization of the differential variational inclusion problem above.

1. (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 [36]. Second, since the method is half-implicit, we chose to evaluate the partial time derivate gt in the configuration $(q(l),t(l))$.

2. (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 $t(l+1)$. That is, $Φi(l+1)≈Φi(l)+Δt vi,n(l+1), ∀ i∈AδK(l)$. Given that $Δt>0$, the condition that the gap between two bodies at the next time step can be at most zero translates into the inequality

$Φi(l)Δt+vi,n(l+1)≥0$
(6)

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 $−μivi,T(l+1)$. 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 $−(Φi(l)/Δt)$. Note that a negative gap represents penetration and the numerical method will produce a solution that seeks to enforce a non-negative gap $Φi(l)$. 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 $Δt$ 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 $Φi(l)$ to a non-negative gap, in one time step, i.e., from $t(l)$ to $t(l+1)$. 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 $Δt$, which divides the gap $Φi(l)$, 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 $AδK(l)$ 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 $t(l)$, i.e., choosing a strictly positive $δK>0$.

The discussion covering the penetration cases (i) through (iii) above was built around the $Φi(l)$ 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 $Φ̃i≡Φi(l)−μiΔtvi,T(l+1)$. Qualitatively, little changes when replacing $Φi(l)$ by $Φ̃i$. Quantitatively though, the use of $Φ̃i$ might lead to a subtle numerical artifact. If relative to the size of $Φi(l)$, the term $μiΔtvi,T(l+1)$ 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 $Δt$. The quantity $vi,T(l+1)$ represents the tangential relative velocity at the contact point. If it is large enough, then the net effect of $μiΔtvi,T(l+1)$ 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 $−μiΔtvi,T(l+1)$. 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 $−μiΔtvi,T(l+1)$. Once the value of the term $−μiΔtvi,T(l+1)$ 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=0$; 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 $AδK(l)$ dictates that two bodies that are δK-close to each other lead to a contact event. As such, the complementarity condition $(Φ̃i/Δt)+vi,n(l+1)≥0$ will enforce the condition that the numerical solution produce a velocity that leads to no penetration at $t(l+1)$. 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.

The friction cone $Ki$ associated with contact event i is defined as $Ki≡{[x, y, z]T∈ℝ3: 0≤x ∧ μix−y2+z2≥0}$. Similarly, the polar cone $Ki°$ associated with the friction cone $Ki$ is defined as $Ki°≡{[a, b, c]T∈ℝ3:ax+by+cz≤0 ∀[x, ,y, z]T∈Ki}$. Based on Eq. (5), $γi(l+1)≡[γi,n(l+1), γi,u(l+1), γi,w(l+1)]T∈Ki$. Define next $di≡[1ΔtΦi(l)+vi,n(l+1), vi,u(l+1), vi,w(l+1)]T∈ℝ3$. Then, using Eq. (5)
$diT·γi(l+1)=γi(l+1)(1ΔtΦi(l)+vi,n(l+1))+γi,u(l+1) vi,u(l+1)+γi,w(l+1) vi,w(l+1)=μiγi,n(l+1)vi,T(l+1)+vi,u(l+1)γi,u(l+1)+vi,w(l+1)γi,w(l+1)=γi,F(l+1) vi,T(l+1)+vi,u(l+1)γi,u(l+1)+vi,w(l+1)γi,w(l+1)=αi γi,F(l+1)−αi γi,F(l+1)=0$

and therefore $γi(l+1)⊥di$.

Next, we show that $−di∈Ki°$, i.e., that $diT·p≥0, ∀ p=[x y z]T∈Ki$. If x = 0, then $y=z=0$ and $diT·p≥0$. If x > 0, then we can scale $p$ by a constant $β>0$ such that $x=γi,n(l+1)$. Note that this scaling does not change the sign of the dot product $diT·p$. We assume $γi,n(l+1)>0$ since the case $γi,n(l+1)=0$ is trivial. Then, using the Cauchy–Schwartz inequality and that $μiγi,n(l+1)≥b2+c2$
$diT·p=γi,n(l+1) (1ΔtΦi(l)+vi,n(l+1))+vi,u(l+1) b+vi,w(l+1) c=μiγi,n(l+1)vi,T(l+1)+vi,u(l+1) b+vi,w(l+1) c≥b2+c2 vi,T(l+1)+vi,u(l+1) b+vi,w(l+1) c≥0$
We thus conclude that we can equivalently express the conditions in Eq. (5) as the following CCP:
$Ki∋γi(l+1)⊥−di∈Ki°$
(7)

#### The “Bilateral Constraints” Component.

Let $0≤j≤NB$. For each bilateral kinematic constraint equation j, define
$bj≡1Δt gj(l)+Gj(l),Tv(l+1)+∂gj(l)∂t$
In the light of Eq. (5c), one has $0=bj⊥λj(l+1)∈ℝ$. Following in the steps of the argument made for the unilateral constraints, one can define the cone $Bj≡ℝ$ and the polar cone $Bj°≡{y: x·y≤0 ∀ x∈Bj}$. Note that this polar cone set has only one element: $Bj°={0}$. Therefore, we can reformulate the condition in Eq. (5c) as
$Bj∋λj(l+1)⊥−bj∈Bj°$
(8)

which is the “bilateral constraint” CCP analog of the condition in Eq. (7).

### Reformulating the CCP.

The next goal is to eliminate any dependency on the unknown velocity $v(l+1)$ in the cone complementarity problems of Eqs. (7) and (8). To this end, using the force balance condition stated in Eq. (5b), one has that
$v(l+1)=v(l)+M−1 f(l)+M−1 G(l)λB,(l+1)+M−1 D(l)γK,(l+1)$
(9)

Let $di,0≡[(1/Δt)Φi(l), 0, 0]T∈ℝ3, di,1≡di,0+Di(l),T(v(l)+M−1 f(l))∈ℝ3$, and $bj,0≡(1/Δt) gj(l)+∂gj(l)/∂t+Gj(l),T(v(l)+M−1 f(l))∈ℝ$. Therefore, $di=di,0+Di(l),Tv(l+1)=di,1+Di(l),TM−1 G(l)λB,(l+1)+Di(l),TM−1 D(l)γK,(l+1)$, and $bj=bj,0+Gj(l)M−1 G(l)λB,(l+1)+Gj(l)M−1 D(l)γK,(l+1)$.

Next, define $P≡[D(l) G(l)]∈ℝ6nb×(3NK+NB), ν(l+1)≡[γK,(l+1),T, λB,(l+1),T]T∈ℝ3NK+NB$, and $p≡[d1,1T,…, dNK,1T, b1,0,…, bNB,0]T∈ℝ3NK+NB$. The terms entering the CCPs both for unilateral and bilateral constraints can then be expressed without any recourse to the velocity $v(l+1)$: $di=di,1+Di(l),TM−1Pν(l+1)$ and $bj=bj,0+Gj(l)M−1Pν(l+1)$. Therefore, we have a collection of CCPs that can be generically represented as
$Ck∋νk(l+1)⊥−(p+Nν(l+1))k∈Ck°$
(10)

where $N≡PTM−1P, C≡K1⊗…⊗KNK⊗B1⊗…⊗BNB$, and $C°≡K1°⊗…⊗KNK°⊗B1°⊗…⊗BNB°$.

We show next that the CCP stated in Eq. (10) represents the first-order optimality conditions [9] for the convex quadratic optimization problem with conic constraints
$ν(l+1)=minν12νTNν+pTνsubject to νk∈Ck$
(11)
To that end, formulate the Karush–Kuhn–Tucker (KKT) optimality conditions via the Lagrangian [9]
$L(ν,ψ,ϕ)=12νTNν+pTν+∑i=1NKψi(γi,u2+γi,w2−μiγi,n)+∑j=1NBϕjλi$
where $ψ$ and $ϕ$ are dummy Lagrange multipliers. The first-order optimality conditions assume the form
${∇νL=03NK+NB1≤i≤NK: 0≤ψi ⊥ μiγi,n−γi,u2+γi,w2≥0 1≤j≤NB: ϕj λj=0$
(12)
The first condition earlier leads to two sets of equalities. First, for $1≤i≤NK$, the gradient with respect to $γi$ yields
$Di(l),TM−1Pν+di,1+ψi[−μiγi,uγi,u2+γi,w2γi,wγi,u2+γi,w2]=03$
which leads to $diT=ψi[−μi (γi,u/γi,u2+γi,w2)(γi,w/γi,u2+γi,w2) ]$. Therefore, using the first set of complementarity conditions in Eq. (12)
$diTγi=−ψi(−μiγi,n+γi,u2+γi,w2γi,u2+γi,w2)=ψi(μiγi,n−γi,u2+γi,w2)=0$
Next, let a vector $α≡[a b c]T∈Ki$. We want to show that $−αTdi≤0$, or equivalently, $αTdi≥0$. First, note that if a = 0, then $b=c=0$, and therefore $αTdi=0$. Otherwise, a > 0, and then
$diT·[abc]≥0⇔aμi+bγi,u+cγi,wγi,u2+γi,w2≥0⇔μiaγi,u2+γi,w2+bγi,u+cγi,w≥0$

As $μia≥b2+c2$, using the Cauchy–Schwartz inequality, $μiaγi,u2+γi,w2+bγi,u+cγi,w≥b2+c2γi,u2+γi,w2+bγi,u+cγi,w≥0$, which proves that as far as $γi$ is concerned, $Ki∋γi⊥−di∈Ki°$. In other words, if $γi$ 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 $bj=−ϕj$. Using the last complementarity condition in Eq. (12), one has that whenever $λj≠0$, necessarily bj = 0. In other words, we have that $Bj∋λj⊥−bj∈Bj°$, which indicates that a $λj$ that satisfies the first KKT conditions in Eq. (12) is a solution of the CCP problem in Eq. (8).

Note that the dynamics step if essentially done once $νk(l+1)$ is computed. Indeed, the new velocity is evaluated using Eq. (9), while the new position is obtained via Eq. (5a).

## 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 $N$ 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.

Results of an experimental validation of the solution methodology discussed are reported in Refs. [1012]. This solution methodology is embedded in the simulation software Chrono [13].

## 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

•
• $A(q(t),δK)$ =

active contact set; two bodies contribute a contact event to this set if they are closer than $δK≥0$ from each other

•
• $AδK(l)$ =

the set of active contact events $A(q(t),δK)$ at time $t(l)$ and in configuration $q(l)$

•
• AA =

rotation matrix associated with body A

•
• $bj$ =

composite velocity associated with bilateral constraint j

•
• $Bj, Bj°$ =

cones associated with bilateral constraint j

•
• $C$ =

system level friction cone, $K1⊗…⊗KNK⊗B1⊗…⊗BNB$

•
• $C°$ =

system level polar cone, $K1°⊗…⊗KNK°⊗B1°⊗…⊗BNB°$

•
• di =

composite velocity; $di≡[(1/Δt)Φi(l)+vi,n(l+1), vi,u(l+1), vi,w(l+1)]T∈ℝ3$

•
• $D(l)$ =

notation for $D(q(l),t(l))$

•
• $D$ =

projection matrix used to generate the set of generalized forces $Dγ̂K$ induced by the contact events present in the system (unilateral constraints)

•
• $εi∈ℝ4$ =

set of four Euler parameters associated with body i. Provides orientation with respect to a global reference frame

•
• $f(l)$ =

notation for impulse $Δt F(t(l),q(l),v(l))$

•
• $F(q,v,t)$ =

the set of generalized applied and Coriolis forces

•
• $g(l)$ =

notation for $g(q(l),t(l))$

•
• $gt(l)$ =

notation for $(∂g(q(l),t(l))/∂t)$

•
• $gj(q,t)$ =

function expression that defines the jth bilateral constraint active in the system at time t, namely $gj(q,t)=0$

•
• $G$ =

projection matrix used to generate the set of generalized forces $Gλ̂B$ induced by the bilateral constraints

•
• $G(l)$ =

notation for $G(q(l),t(l))$

•
• $L(q)$ =

transformation matrix that links the time derivative of the Euler parameters to the angular velocity

•
• $M$ =

constant, system level, mass matrix

•
• nb =

number of rigid bodies in the system

•
• $NB$ =

number of bilateral constraints present in the system

•
• $ni$ =

unit normal vector at point of contact; oriented from the body of lower index to the body of higher index

•
• $ni,ui,wi$; combine to form the frictional contact vector $γ̂i$

•
• $NK$ =

number of contact events present in the system at a given time (explicit dependency on time dropped for convenience)

•
• $N$ =

quadratic term matrix, optimization problem; computed as $PTM−1P$

•
• $P$ =

system level projection matrix, defined as $[D(l) G(l)]$

•
• $p$ =

coefficient of the linear term in the CCP

•
• $q∈ℝ7nb$ =

set of generalized coordinates associated with the nb bodies in the system

•
• $ri∈ℝ3$ =

Cartesian space location of body i

•
• $s¯̃i,B$ =

vector that for contact event i provides the location of the contact point expressed in the local reference frame associated with body B

•
• $t(l)$ =

time associated with the integration time step l; also, $Δt=t(l+1)−t(l)$

•
• $ui, wi$ =

two unit normal vectors at point of contact; span the tangent contact plane; together with $ni$ form a orthogonal reference frame

•
• $vi,n, vi,u, vi,w$ =

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 {$ni, ui, wi$}

•
• $vi,T(l+1)$ =

magnitude of the tangential velocity at the point of contact; notation for $(vi,u(l+1))2+(vi,w(l+1))2$

•
• $v$ =

collection of all body velocities $[r˙1T,ω¯1T,…,r˙nbT,ω¯nbT]T∈ℝ6nb$

•
• $ω¯i∈ℝ3$ =

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

•
• $Φi$ =

signed gap function: negative if bodies penetrate, positive if they don't touch, zero if one point of contact

•
• $γ̂i,n, γ̂i,u, γ̂i,w$ =

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

•
• $γ̂K$ =

set of Lagrange multipliers associated with the unilateral constraints present in the system

•
• $γK,(l+1)$ =

notation for frictional contact impulse $Δt γ̂K,(l+1)$

•
• $γi,F(l+1)$ =

magnitude of the friction impulse associated with the contact event i; notation for $(γi,u(l+1))2+(γi,w(l+1))2$

•
• $Δt$ =

numerical integration time step

•
• $Ki$ =

friction cone associated with contact event i

•
• $Ki°$ =

polar cone associated with the friction cone $Ki$

•
• μi =

friction coefficient associated with contact event i

•
• $λ̂B$ =

set of Lagrange multipliers associated with the bilateral constraints

•
• $λB,(l+1)$ =

notation for bilateral constraint reaction impulse $Δt λ̂B,(l+1)$

•
• $Φi(l)$ =

notation for $Φi(q(l))$; i.e., gap for contact event i at time $t(l)$

## References

References
1.
Haug
,
E. J.
,
1989
,
Computer-Aided Kinematics and Dynamics of Mechanical Systems Volume-I
,
Prentice Hall
,
Englewood Cliffs, NJ
.
2.
Tasora
,
A.
, and
Anitescu
,
M.
,
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
.
3.
Baumgarte
,
J.
,
1972
, “
Stabilization of Constraints and Integrals of Motion in Dynamical Systems
,”
Comput. Methods Appl. Mech. Eng.
,
1
(
1
), pp.
1
16
.
4.
Ostermeyer
,
G.
,
1990
, “
On Baumgarte Stabilization for Differential Algebraic Equations
,”
Real-Time Integration Methods for Mechanical System Simulation
(Nato ASI Subseries F:),
Springer-Verlag
,
Berlin
, pp.
193
207
.
5.
Ascher
,
U. M.
,
Chin
,
H.
, and
Reich
,
S.
,
1994
, “
Stabilization of DAEs and Invariant Manifolds
,”
Numer. Math.
,
67
(
2
), pp.
131
149
.
6.
Kikuuwe
,
R.
, and
Brogliato
,
B.
,
2017
, “
A New Representation of Systems With Frictional Unilateral Constraints and Its Baumgarte-Like Relaxation
,”
Multibody Syst. Dyn.
,
39
(3), pp.
267
290
.
7.
Mazhar
,
H.
,
Heyn
,
T.
,
Tasora
,
A.
, and
Negrut
,
D.
,
2015
, “
Using Nesterov's Method to Accelerate Multibody Dynamics With Friction and Contact
,”
ACM Trans. Graphics
,
34
(
3
), pp. 1–14.
8.
Smith
,
B.
,
Kaufman
,
D. M.
,
Vouga
,
E.
,
Tamstorf
,
R.
, and
Grinspun
,
E.
,
2012
, “
Reflections on Simultaneous Impact
,”
ACM Trans. Graphics
,
31
(
4
), p. 106.
9.
Bertsekas
,
D. P.
,
1995
,
Nonlinear Programming
,
Athena Scientific
,
Belmont, MA
.
10.
Kwarta
,
M.
, and
Negrut
,
D.
,
2016
, “
Using the Complementarity and Penalty Methods for Solving Frictional Contact Problems in Chrono: Validation for the Cone Penetrometer Test
,” Simulation-Based Engineering Laboratory, University of Wisconsin-Madison, Madison, WI, Technical Report No.
TR-2016-16
.http://sbel.wisc.edu/documents/TR-2016-16.pdf
11.
Kwarta
,
M.
, and
Negrut
,
D.
,
2016
, “
Using the Complementarity and Penalty Methods for Solving Frictional Contact Problems in Chrono: Validation for the Triaxial Test
,” Simulation-Based Engineering Laboratory, University of Wisconsin-Madison, Madison, WI, Technical Report No. TR-2016-17.
12.
Kwarta
,
M.
, and
Negrut
,
D.
,
2016
, “
Using the Complementarity and Penalty Methods for Solving Frictional Contact Problems in Chrono: Validation for the Shear-Test With Particle Image Velocimetry
,” Simulation-Based Engineering Laboratory, University of Wisconsin-Madison, Madison, WI, Technical Report No. TR-2016-18.
13.
Tasora
,
A.
,
Serban
,
R.
,
Mazhar
,
H.
,
Pazouki
,
A.
,
Melanz
,
D.
,
Fleischmann
,
J.
,
Taylor
,
M.
,
Sugiyama
,
H.
, and
Negrut
,
D.
,
2016
, “
Chrono: An Open Source Multi-Physics Dynamics Engine
,”
High Performance Computing in Science and Engineering
(Lecture Notes in Computer Science),
T.
Kozubek
, ed.,
Springer
,
Cham, Switzerland
, pp.
19
49
.