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]T7nb, and its time derivative q˙=[r˙1T,ε˙1T,,r˙nbT,ε˙nbT]T7nb, is used to represent the state of the system, where for body j, 1jnb, 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]T6nb is tied to q˙ via a linear transformation [1]
(1)

With the ground assigned by convention index 0, assume two bodies of index A and B, 0A<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 Φ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.

Fig. 1
Bodies A and B in contact; a local reference frame {ni, ui, wi} 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 s¯i,A and s¯i,B vectors.
Fig. 1
Bodies A and B in contact; a local reference frame {ni, ui, wi} 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 s¯i,A and s¯i,B vectors.
Close modal

Any two bodies that are closer than a prescribed δK0 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,wwiAiγ̂i. The location of point P on body B is rBP=rB+ABs¯i,B and its virtual displacement is δrBP=δrBABs¯̃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
Similarly, the virtual work for body A is
The virtual work that the presence of the frictional contact force Fi,B imparts is
where

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

The equations of motion assume the form [1]
(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[D1DNK]6nb×3NK,G[G1GNB]6nb×NB,γ̂K[γ̂1T,,γ̂NKT]T3NK, and λ̂B[λ̂1T,,λ̂NBT]TNB.

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
(3a)
where 1jNB. Their time derivative yields
(3b)
Finally

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:
(4a)
(4b)
(4c)

where μi0 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:
(5a)
(5b)
(5c)
and, iAδK(l),
(5d)

In Eq. (5), f(l)ΔtF(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)+Δtvi,n(l+1),iAδ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

(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]T3:0xμixy2+z20}. Similarly, the polar cone Ki° associated with the friction cone Ki is defined as Ki°{[a,b,c]T3:ax+by+cz0[x,,y,z]TKi}. Based on Eq. (5), γi(l+1)[γi,n(l+1),γi,u(l+1),γi,w(l+1)]TKi. Define next di[1ΔtΦi(l)+vi,n(l+1),vi,u(l+1),vi,w(l+1)]T3. Then, using Eq. (5)

and therefore γi(l+1)di.

Next, we show that diKi°, i.e., that diT·p0,p=[xyz]TKi. If x = 0, then y=z=0 and diT·p0. 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
We thus conclude that we can equivalently express the conditions in Eq. (5) as the following CCP:
(7)

The “Bilateral Constraints” Component.

Let 0jNB. For each bilateral kinematic constraint equation j, define
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·y0xBj}. Note that this polar cone set has only one element: Bj°={0}. Therefore, we can reformulate the condition in Eq. (5c) as
(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
(9)

Let di,0[(1/Δt)Φi(l),0,0]T3,di,1di,0+Di(l),T(v(l)+M1f(l))3, and bj,0(1/Δt)gj(l)+gj(l)/t+Gj(l),T(v(l)+M1f(l)). Therefore, di=di,0+Di(l),Tv(l+1)=di,1+Di(l),TM1G(l)λB,(l+1)+Di(l),TM1D(l)γK,(l+1), and bj=bj,0+Gj(l)M1G(l)λB,(l+1)+Gj(l)M1D(l)γK,(l+1).

Next, define P[D(l)G(l)]6nb×(3NK+NB),ν(l+1)[γK,(l+1),T,λB,(l+1),T]T3NK+NB, and p[d1,1T,,dNK,1T,b1,0,,bNB,0]T3NK+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),TM1Pν(l+1) and bj=bj,0+Gj(l)M1Pν(l+1). Therefore, we have a collection of CCPs that can be generically represented as
(10)

where NPTM1P,CK1KNKB1BNB, and C°K1°KNK°B1°BNB°.

The Quadratic Problem Angle

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
(11)
To that end, formulate the Karush–Kuhn–Tucker (KKT) optimality conditions via the Lagrangian [9]
where ψ and ϕ are dummy Lagrange multipliers. The first-order optimality conditions assume the form
(12)
The first condition earlier leads to two sets of equalities. First, for 1iNK, the gradient with respect to γi yields
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)
Next, let a vector α[abc]TKi. We want to show that αTdi0, or equivalently, αTdi0. First, note that if a = 0, then b=c=0, and therefore αTdi=0. Otherwise, a > 0, and then

As μiab2+c2, using the Cauchy–Schwartz inequality, μiaγi,u2+γi,w2+bγi,u+cγi,wb2+c2γi,u2+γi,w2+bγi,u+cγi,w0, which proves that as far as γi is concerned, KiγidiKi°. 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 λj0, necessarily bj = 0. In other words, we have that BjλjbjBj°, 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 δK0 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, K1KNKB1BNB

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)]T3

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)

εi4 =

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

f(l) =

notation for impulse ΔtF(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 PTM1P

P =

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

p =

coefficient of the linear term in the CCP

q7nb =

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

ri3 =

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]T6nb

ω¯i3 =

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

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
.