Mihai Anitescu
email: anitescu@mcs.anl.

OPTIMIZATION-BASED SIMULATION OF NONSMOOTH RIGID MULTIBODY DYNAMICS

Keywords: Multibody Dynamics, Contact, Coulomb Friction, Convex Relaxation AMS subject classifications. 65K10, 90C33

We present a time-stepping method to simulate rigid multibody dynamics with inelastic collision, contact, and friction. The method progresses with fixed time step without backtracking for collision and solves at every step a strictly convex quadratic program. We prove that a solution sequence of the method converges to the solution of a measure differential inclusion. We present numerical results for a few examples, and we illustrate the difference between the results from our scheme and previous, linear-complementarity-based time-stepping schemes.

1. Introduction. Nonsmooth rigid multibody dynamics (NRMD) is an important paradigm of modern computational science. It consists of predicting the position and velocity evolution of a group of rigid particles that are subject to noninterpenetration, collision, adhesion, and dry friction constraints and to possibly global forces (such as electrostatic and gravitational forces). The dynamics of such a group of particles is nonsmooth because of the intermittent nature of noninterpenetration, collision, and adhesion constraints and because of the nonsmooth nature of the dry friction constraints at stick-slip transitions. NRMD has been successfully applied to a vast group of diverse applications such as granular [36] and rock dynamics [START_REF] Camborde | Numerical study of rock and concrete behaviour by discrete element modeling[END_REF], masonry stability analysis [START_REF] Thavalingham | Computational framework for discontinuous modeling of masonry arch bridges[END_REF], simulation of concrete obstacle response to explosion [START_REF] Camborde | Numerical study of rock and concrete behaviour by discrete element modeling[END_REF], tumbling mill design (for grinding ore in mineral-processing industries) [START_REF] Mishra | A review of computer simulation of tumbling mills by the discrete element method part ii-practical applications[END_REF], interactive virtual reality [START_REF]Karma user's guide[END_REF], and robot simulation and design [START_REF] Anitescu | Constraint stabilization for time-stepping ap[END_REF]. In civil engineering applications NRMD is often used under the name of the discrete element method; in the physics literature it is regarded as a particular instance of a molecular dynamics approach.

In addition to these applications, exciting new prospects for NRMD have been generated by recent investigation in dynamic self-assembly in electrostatically driven granular media [START_REF] Sapozhnikov | Dynamic selfassembly and patterns in electrostatically driven granular media[END_REF] and pedestrian and evacuation dynamics (PED) [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF][START_REF] Helbing | Simulation of pedestrian crowds in normal and evacuation situations[END_REF]. In PED, NRMD has provided a rational and systematic way to improve the design of evacuation paths and other areas of high pedestrian density [START_REF] Helbing | Simulating dynamical features of escape panic[END_REF]. For example, a stunning result obtained through NRMD simulations is that emergency evacuations complete faster and with fewer injuries if a column is placed in front of (and at some distance upflow from) the emergency exit [START_REF] Helbing | Simulation of pedestrian crowds in normal and evacuation situations[END_REF]. Clearly NRMD has a huge impact potential for significant applications with complex dynamics.

The NRMD approach presents several difficulties. It does not always have a classical solution if the dry friction coefficient is not zero [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF]. The nonsmoothness of its description and its constrained nature make direct application of numerical methods for ordinary differential equations and differential algebraic equations impossible. In response to these difficulties the vast majority of numerical approaches create a smooth and stiff approximation of the problem that is integrated explicitly with a time step inside the stability region. In turn, this creates other difficulties. First, the more accurate the approximation, the stiffer the problem, which may lead to prohibitively small time steps to achieve stability. Second, the approximation may introduce nonphysical artifacts. For example, it may systematically underestimate the particle velocities [START_REF] Hryciw | Video tracking for experimental validation of discrete element simulations of large discontinuous deformation[END_REF]. Thirdly, the smoothing parameters are notoriously difficult to tune to achieve the desired compromise between accuracy and computational efficiency.

An alternative to the smoothing approach is to incorporate contact and friction as hard constraints, either in acceleration-force approaches [START_REF] Featherstone | Robot Dynamics Algorithms[END_REF][START_REF] Glocker | An lcp-approach for multibody systems with planar friction[END_REF]9,[START_REF] Trinkle | On dynamic multi-rigid-body contact prob[END_REF] or in impulsevelocity approaches [START_REF] Stewart | An implicit time-stepping scheme for rigid-body dynamics with inelastic collisions and Coulomb friction[END_REF]6,[START_REF] Anitescu | A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF]. The latter approach, to which we will exclusively refer to in this paper, is related to the contact dynamics method [START_REF] Radjai | Force distributions in dense twodimensional granular systems[END_REF] used for granular dynamics simulation. It consists of using a backward Euler-based time-stepping scheme where at every time step all nonsmooth constraints are enforced by complementarity and inequality constraints and where the fundamental variables become the velocities and impulses, as opposed to the accelerations and forces. The method presents several advantages over classical approaches.

1. The time-stepping scheme has a solution for any choice of parameters, as opposed to its acceleration formulation, and produces an uniformly bounded velocity sequence over any finite time interval [6]. This is not the case for acceleration-force approaches, which would necessarily fail when a configuration that does not have a classical solution is encountered [9,[START_REF]Rigid-body dynamics with friction and impact[END_REF]. 2. The discrete solution converges to a weak solution of the problem, at least when adhesion is absent and the restitution coefficient is 0 [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF]. A classical solution does not always exist [START_REF]Rigid-body dynamics with friction and impact[END_REF]. 3. The time-stepping scheme progresses stably over any finite time interval without the need to tune any additional parameters [START_REF] Anitescu | A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF]6]. [START_REF]Solving nonconvex problems of multibody dynamics with joints, contact and small friction by sequential convex relaxation[END_REF]. The scheme can be modified to accommodate stiff external forces [7] and to ensure constraint stabilization while solving only one linearly constrained complementarity subproblem per time step [START_REF] Anitescu | A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF][START_REF] Anitescu | Constraint stabilization for time-stepping ap[END_REF]. To progress with a lower bounded time step, the algorithm is not stopped at collisions, but potentially active noninterpenetration constraints are added to the complementarity constrained subproblem, similar to the contact dynamics method. This has led to an event-driven single-processor version of the time-stepping scheme, which includes the author's contribution of using a linearized version of the constraints and a modification of the mass matrix that accommodates stiff external forces. That version was implemented by a commercial provider of physical content for small scale interactive virtual reality simulations [START_REF]Karma user's guide[END_REF]. For this class of applications the variety and heterogeneity of the scenarios required by the users make it impossible to appropriately tune a smoothing model.

The hard constraint techniques we use here result in algorithms that are stable for larger timesteps when compared to the penalty molecular dynamics approaches. This, however, comes at the cost of solving a linear complementarity subproblem that may have a nonconvex solution set and thus may be expensive to solve [START_REF]Solving nonconvex problems of multibody dynamics with joints, contact and small friction by sequential convex relaxation[END_REF]. This difficulty is related to the fact that the Coulomb friction law is not associative: the linear complementarity subproblem does not represent the optimality conditions for some optimization problem, as is the case for frictionless dynamics.

Many applications of interest consist of thousands of fairly tightly packed rigid bodies. An important source of very large scale problems with contact and friction is the direct simulation approach of granular matter, which was used for the simulation of 10 million rigid bodies [10]. In order to extend the benefits of the hard constraint approach to such large-scale systems, alternative approaches are needed that solve a simpler subproblem while maintaining the stability properties of the simulation. The need for such a requirement comes from the fact that, to date, there is no reasonable very large-scale approach for solving non-convex linear complementarity problems of the type that is needed in rigid multibody dynamics simulations [START_REF]A fixed-point iteration approach for multibody dynamics with contact and friction[END_REF]. Ideally, such approaches should preserve as much of the quality of the solution as the linear complementarity approaches, in the limit of the time step going to 0.

To that end, we present an optimization-based method for the simulation of nonsmooth rigid body dynamics: rigid body dynamics with collision, contact and friction. At every step, the method solves one convex quadratic program, and progresses with a fixed time step. Therefore, the number of quadratic programs that is needed to be solved is determined ahead of time once the time step is chosen. This represents a major advantage over integrate-detect-restart strategies [9, 1] for which there is no upper bound in the number of subproblems that need to be solved. The method has been introduced in [START_REF] Anitescu | A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF] following a fixed-point iteration approach introduced in [START_REF]Solving nonconvex problems of multibody dynamics with joints, contact and small friction by sequential convex relaxation[END_REF][START_REF]A fixed-point iteration approach for multibody dynamics with contact and friction[END_REF], but its convergence as h → 0 has not been studied prior to this work.

In this work we complete its analysis, and we show that, in the limit of the time step going to 0, the solution of the numerical scheme will converge to the solution of a measure differential inclusion [START_REF] Moreau | Standard inelastic shocks and the dynamics of unilateral constraints[END_REF][START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF][START_REF]Rigid-body dynamics with friction and impact[END_REF]. Given the recent success in solving very large scale quadratic programs, this research may open a new avenue in the hard constraint simulation of very large scale NRMD.

We also discuss the physical meaning of our approximation, and we present numerical solutions to a few widely discussed examples in the literature, to which we have also applied for comparison the equivalent LCP-based method [6,[START_REF] Stewart | An implicit time-stepping scheme for rigid-body dynamics with inelastic collisions and Coulomb friction[END_REF], as well as to a medium-scale granular matter simulation.

The Contact Model.

In the following, we present our contact model, and we compare it to previous approaches. The object of study is a system of rigid bodies, described by state variables and contact and frictional constraints.

2.1. System representation. At a time t, the position of the system is described by generalized coordinates q(t) (which may include rotational coordinates that cannot be defined over a subspace homeomorphic to R n , for some n), and generalized velocities v(t). In classical mechanics, v(t) is continuous, and we can write dq dt = v. Since we are attempting to accommodate impact, v(t) is not, in general, continuous. Therefore we will require that the position of editor and the velocity vector satisfied the following, weaker, condition

q(t) -q(0) = t 0 v(τ )dτ.
In two dimensions, each rigid body can be described by the position (x, y) of the center of mass and the angle φ of rotation around the center of mass. The result is a total of three coordinates per body. So a system of n rigid bodies can be described by the set of coordinates x 1 , y 1 , φ 1 , x 2 , y 2 , φ 2 , . . . , x n , y n , φ n , or 3n coordinates.

In three dimensions, the position of a rigid body is described by the position x, y, z of the center of mass and a 3 × 3 orthogonal matrix Q that represents the rotation of a frame attached to the body with respect to a fixed-world frame. Here we assume that Q can be represented smoothly by three parameters. This parameterization is valid only locally, but this problem can be easily remedied by an appropriate reparameterization, which does not affect the dynamics [START_REF] Haug | Computer Aided Kinematics and Dynamics of Mechanical Systems[END_REF]. To simplify our approach, we assume that the position of a three-dimensional rigid body can be represented by six parameters: x, y, z, φ, ζ, θ. Therefore a system with n bodies in three dimensions is represented by 6n coordinates.

Nonpenetration constraints.

Two rigid bodies should not penetrate, and, if they are in contact, there should be friction acting at the interface. To enforce the nonpenetration constraint, we assume that there exists a function Φ(q) that satisfies Φ(q) = > 0 if the bodies are separated = 0 if the bodies touch each other < 0 if the bodies are interpenetrating (2.1) For such a function, the nonpenetration constraint becomes Φ(q) ≥ 0.

An example for such a mapping is the signed distance function [START_REF] Kim | Deep: Dual-space expansion for estimating penetration depth between convex polytopes[END_REF], which is differentiable when the bodies are smooth and convex, at least up to some values of the interpenetration [START_REF] Anitescu | Formulating 3d contact dynamics problems[END_REF]. For most cases, even simple ones involving the relative position of two spheres, a differentiable signed distance function cannot be defined for all values of q. The fact that Φ(q) can be differentiably defined only on a neighborhood of the set Φ(q) ≥ 0 can be accommodated at the cost of making the analysis substantially more involved [START_REF] Anitescu | A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF]. To simplify our discussion, we make the following assumption.

Assumption A1 (Differentiability of geometrical constraint data). Any contact is described by a signed distance function Φ(q) that is everywhere twice continuously differentiable.

Frictional constraints.

In this work we describe the frictional constraints by hard constraints, which in turn can be reduced to complementarity models [6, 33].

2.3.1.

The Coulomb friction model. The model we represent and approximate is the Coulomb friction model. If a position q is feasible and the contact is active, that is, Φ(q) = 0, then at the contact we have a normal force and a tangential force.

Let n be the normal at the contact pointing from the second body to the first body, and let t 1 and t 2 be the tangents at the contact. Here n, t 1 , t 2 are mutually orthogonal vectors in three dimensions of length one. The vectors n, t 1 , and t 2 are a function of the position q, but we ignore this fact until the end of this section.

The reaction force is impressed on the system by means of multipliers c n ≥ 0, β 1 and β 2 . The normal component of the force is F N = c n n, and the tangential component of the force is

F T = β 1 t 1 + β 2 t 2 .
The Coulomb model consists of the following constraints:

c n ≥ 0, Φ(q) ≥ 0, Φ(q)c n = 0, μc n ≥ β 2 1 + β 2 2 , ||v T || μc n -β 2 1 + β 2 2 = 0, F T , v T = -||F T || ||v T || (2.2)
where v T is the relative tangential velocity at contact. The parameter that defines the effect of the friction over the dynamical system is μ ¿ 0, the friction coefficient.

The friction coefficient typically has values between 0 and 1.

The first constraint can be restated as

F = F N + F T = c n n + β 1 t 1 + β 2 t 2 ∈ K,
where K is a cone in three dimensions, whose slope is arctan(μ).

The constraint F T , v T = -||F T || ||v T || requires that the tangential force be opposite to the tangential velocity. This results in the reaction force being dissipative. In fact, an equivalent convenient way of expressing this constraint is by using the maximum dissipation principle [START_REF] Stewart | An implicit time-stepping scheme for rigid-body dynamics with inelastic collisions and Coulomb friction[END_REF][START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF][START_REF]Rigid-body dynamics with friction and impact[END_REF] (

β 1 , β 2 ) = argmin √ β 2 1 +β 2 2 ≤μc n β 1 t 1 + β 2 t 2 T v T .
These constraints are represented by mapping the vectors n, t 1 , t 2 from contact coordinates to generalized coordinates [START_REF] Anitescu | Formulating 3d contact dynamics problems[END_REF].

For example, if we have a two body system, then the generalized coordinates in the three-dimensional space are embedded in a twelve-dimensional space by using the coordinates:

x 1 , y 1 , z 1 , φ 1 , θ 1 , ζ 1 , x 2 , y 2 , z 2 , φ 2 , θ 2 , ζ 2 .
For a three-dimensional vector v, the mapping to generalized coordinates is

v ⇒ ⎛ ⎝ v r 1 × v -v -r 2 × v ⎞ ⎠ ,
where r 1 and r 2 are the relative positions of the contact point with respect to the centers of mass of the two bodies [START_REF] Anitescu | Formulating 3d contact dynamics problems[END_REF]. Using this mapping, we denote the generalized vector version of n, t 1 , t 2 by n, t 1 , t 2 . One unfortunate side effect of generalized coordinates mapping is that, in the new coordinates, n, t 1 , t 2 cease to be mutually orthogonal.

If v is the generalized velocity, the tangential velocity can be expressed by using the quantities in generalized coordinates as

v T = v T t 1 t 1 + v T t 2 t 2 .
In generalized coordinates, the Coulomb model thus becomes

F N = c n n, F T = β 1 t 1 + β 2 t 2 , (2.3)
c n ≥ 0, Φ(q) ≥ 0, c n Φ(q) = 0, (2.4
)

μc n ≥ β 2 1 + β 2 2 , (2.5)
v T = v T t 1 t 1 + v T t 2 t 2 , F T , v T = -||F T || ||v T || . (2.6)
The maximum dissipation principle can now be invoked in generalized coordinates to read

(β 1 , β 2 ) = argmin √ β 2 1 +β 2 2 ≤μc n (β 1 t 1 + β 2 t 2 ) T v,
where c n and v are considered fixed. With respect to the regularity of n, t 1 , t 2 as a function of the position variables q, we make the following assumption.

Assumption A2 (Differentiability of contact data). The mappings n(q), t 1 (q) and t 2 (q) are continuously differentiable.

2.4. The overall dynamical model. The other dynamical data needed to fulfill the model are the mass matrix M (q), the external force k(t, v, q), and the inertial force f c (q, v). The latter contains the centrifugal and Coriolis force. The mapping f c (q, v) is continuously differentiable and satisfies

[7] v T f c (q, v) = 0 ∀q, v.
With respect to these quantities, we invoke the following assumptions. Assumption A3 (Constant mass matrix). The mass matrix M (q) is positive definite and constant. This assumption is satisfied in two dimensions and three dimensions if we use the Newton-Euler formulation in body coordinates [START_REF] Murray | A mathematical introduction to robotic manipulation[END_REF]. Assumption A4 (At most linearly increasing external force). The external force is continuous and increases at most linearly with the position and the velocity and is uniformly bounded in time.

Hence,

k(t, v, q) = k 0 (t, v, q) + f c (v, q) + k 1 (v) + k 2 (q), (2.7)
and there exists c K ≥ 0 such that

||k 0 (t, v, q)|| ≤ c K , ||k 1 (v)|| ≤ c K ||v|| , ||k 2 (q)|| ≤ c K ||q|| . (2.8)
Assume now that we have p potential contact constraints, which are enforced by the nonpenetration constraints Φ (j) (q) ≥ 0, j = 1, 2, . . . , p.

In the following, we denote by the superscript (j) the data associated to the potential contact (j). The continuous model is the following differential variational inequality [START_REF] Stewart | Differential variational inequalities, Mathematical Programming, Series A[END_REF] M dv dt = j=1,2,...,p c (j)

n n (j) + β (j) 1 t (j) 1 + β (j) 2 t (j) 2 + f c (q, v) + k(q, v) q = v c (j) n ≥ 0 ⊥ Φ (j) (q) ≥ 0, j = 1, 2, . . . , p β (j) 1 , β (j) 2 = argmin μ (j) c (j) n ≥ q β (j) 1 +β (j) 2 , j = 1, 2, . . . , p v T t (j) 1 t (j) 1 + v T t (j) 2 t (j) 2 T β 1 t (j) 1 + β 2 t (j) 2
.

(2.9) The Coulomb model, that we use in this work, is the predominant model that is used in the engineering literature to describe dry friction. Unfortunately, the model may be inconsistent: there exist configurations for which the model does not have a solution [9,[START_REF]Rigid-body dynamics with friction and impact[END_REF]. This situation has led to the need to explore weaker formlations of the model for dynamics.

We will consider that all collisions that appear during the simulation are of the inelastic type. Therefore, they are naturally treated by the time stepping scheme through a change of active set without the need to modify the algebraic expression of the scheme.

3. Measure Differential Inclusions. In the following, we explore weaker formulations of (2.9) based on the concept of measure differential inclusions. Such an investigation is necessary because the continuous model (2.9) does not necessarily have a classical solution, that is a solution q(t) that is twice continuously differentiable [9,[START_REF]Rigid-body dynamics with friction and impact[END_REF]. In effect, q(t) sometimes is not even once continuously differentiable.

The connection between measure differential inclusions and time stepping schemes like the ones in this work comes from the fact that it was recently proven that the linear complementarity-based time-stepping scheme converges to a solution of the measure differential equation [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF].

To introduce the concept of measure differential inclusion (for which we use [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF]), we describe the friction cone associated to the potential contact (j), where j ∈ {1, 2, . . . , p} by

F C (j) (q) = c (j) n n (j) + β (j) 1 t (j) 1 + β (j) 2 t (j) 2 c (j) n ≥ 0, β (j) 1 2 + β (j) 2 2 ≤ μ (j) c (j) n .
The total friction cone becomes

F C(q) = j=1,2,...,p c (j) n n (j) + β (j) 1 t (j) 1 + β (j) 2 t (j) 2 β (j) 1 2 + β (j) 2 2 ≤ μ (j) c (j) n , j = 1, 2, . . . , p c (j) n ≥ 0 ⊥ Φ (j) (q) = 0, j = 1, 2, . . . , p .
It is immediate that

F C(q) = j=1,2,...,p, Φ (j) (q)=0 F C (j) (q).
Hence a solution of (2.9) must satisfy

M dv dt = f C (q, v) + k(q, v) + ρ dq dt = v. ρ (j) (t) ∈ F C (j) (q(t)), j = 1, 2, . . . , p Φ (j) (q) ≥ 0, j = 1, 2, . . . , p ρ (j) Φ (j) (q) = 0, j = 1, 2, . . . , p, ρ = p j=1 ρ (j) (t). (3.1)
Clearly, from (3.1) it follows that ρ(t) ∈ F C(q). We are interested in v(t) and q(t) that satisfy (3.1). Since v(t) must be discontinuous to allow for collisions, dv dt is meaningless in a classical sense. Since the velocity is discontinuous, the force multipliers c n , β 1 , β 2 also cannot exist in a classical sense. We therefore allow them to be distributions or vector measures.

3.1. Vector measures and measure differential inclusions. In the following formalism we use the setup and some of the results from [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF].

A vector measure ν is defined in terms of its action on a continuous function φ:

ν, φ = φdν = φ(t)ν(dt).
Equation (3.1) can now be interpreted in terms of vector measures where v(t) is a function of bounded variation by means of the Riemann-Stjeltjes integral: φ(t)dν(t). For continuous φ, [a,b] φdν can be approximated by finite Riemann sums:

N -1 i=0 φ(τ i )[v(t i+1 ) -v(t i )],
where

a = t 0 < τ u < t 1 < . . . τ N -1 < t n = b.
For the Riemann-Stieltjes integral to be well defined, the function v(•) must have finite variation T 0 . Here T 0 is the supremum of the sums

N -1 i=0 ||v(t i+1 ) -v(t i )|| over all finite sequences a = t 0 < t 1 < . . . < t N -1 < t N = b over all N .
We now extend the following form of Newton's law,

M dv dt -f c (q(t), v(t)) -k(t, q(t), v(t)) ∈ F C(q(t)). (3.2)
which is valid for the classical case, to the case where v(t) is a function of bounded variation.

Definition (Measure Differential Inclusion [START_REF]Rigid-body dynamics with friction and impact[END_REF]) If ν is a measure and K(•) is a convex-set valued mapping, we say that dv dt ∈ μ K(t), if, for all continuous φ ≥ 0 with compact support, not identically 0, we have that

φ(t)ν(dt) φ(t)dt ∈ τ :φ(τ ) =0 K(τ ).
Definition: (Weak solution of (3.2)) We therefore say that v(t), q(t) is a weak solution of (3.2) if 1. v(t) is a function of bounded variation. 2. q(•) is a continuous, locally Lipschitz function that satisfies

q(t) = q(0) + t 0 v(τ )dτ.
3. The measure dv(t), which exists as a result of v being a bounded variation function, must satisfy,

d(M v) dt -k(t, v) -f c (q, v) ∈ μ F C(q(t)).
4. Φ (j) (q) ≥ 0, ∀j = 1, 2, . . . , p. Note that, from our definition of F C(q), the statement ∈ μ F C(q(t)) also contains the implication about the active set, because ρ(t) ∈ F C(q(t)) means that ρ(t) = p j=1, Φ (j) (q)=0 ρ (j) (t).

4. Regularity assumptions and polyhedral friction cone approximation. The regularity assumptions refer to the quality of the friction cone. We work with several types of friction cones.

One of these types is the -active friction cone, where ≥ 0:

F C (q) = Φ (j) (q)≤ F C (j) (q).
Clearly, F C 0 (q) = F C(q), and, because of the continuity of Φ (j) (q), we also have that

F C (q) →0 =⇒ F C(q).
4.1. Polyhedral approximations of the friction cone. Many of the existing approaches [6, [START_REF] Stewart | An implicit time-stepping scheme for rigid-body dynamics with inelastic collisions and Coulomb friction[END_REF][START_REF] Trinkle | On dynamic multi-rigid-body contact prob[END_REF] use polyhedral approximations of the friction cone. The reason is that the polyhedral description results in linear complementarity problems, for which reliable software, such as PATH [START_REF] Munson | Algorithms and Environments for Complementarity[END_REF][START_REF] Dirkse | The PATH solver: A non-monotone stabilization scheme for mixed complementarity problems[END_REF], exists.

We define a polyhderal approximation to F C (j) (q) as follows. We define a new set of tangent vectors

d (j) k = cos 2πk m (j) t (j) 1 + sin 2πk m (j) t (j) 2 , k = 1, 2, . . . , m (j) .
The friction cone approximation that we use, for a fixed (j), is

F C (j),m (j) (q) = n (j) c (j) n + m (j) k=1 β k d (j) k c (j) n ≥ 0; β (j) k ≥ 0, k = 1, 2, . . . , m (j) ; μ (j) c (j) n ≥ m (j) k=1 β (j) k .
Using the inequality

m (j) i=1 (β i cos(θ i )) 2 + m (j) i=1 (β i sin(θ i )) 2 ≤ ( β i ) 2 , which holds for β 1 , β 2 , β 3 , . . . , β N ≥ 0, it follows that F C (j),m (j) ⊂ F C (j) (q).
It is also clear that, for fixed (j), we must have

F C(q) (j),m (j) m (j) →∞ ---------→ F C (j) (q). (4.1)
We denote m = m

(j) 1 , m (j) 2 , . . . , m (j) p
. We define the 0 polyhedral approximation to F C(q) as

F C m (q) = Φ (j) (q)=0 F C(q) (j),m (j)
and the polyhedral approximation cone as

F C m (q) = Φ (j) (q)≤ F C(q) (j),m (j)
.

Dual cones.

In this work, an important role is played by the dual of the friction cone. For one contact, we have that

F C (j) (q) = t = c (j) n n (j) + β (j) 1 t 1 + β (j) 2 t (j) 2 μc (j) n ≥ β (j) 2 1 + β (j) 2 2 , c (j) n ≥ 0 .
The dual cone is

F C (j) * (q) = v ∈ R I n |n (j) T v ≥ μ (j) (d (j) T 1 v) 2 + (d (j) T 2 v) 2 .
For the polyhedral approximation with m (j) facets for one contact, we have that

F C (j) (q) m (j) = ⎧ ⎨ ⎩ m (j) k=1 d (j) k β (j) k + n (j) c (j) n |c (j) n ≥ 0, β (j) ≥ 0, m (j) k=1 β (j) k ≤ μ (j) c (j) n ⎫ ⎬ ⎭ .
The dual cone of the polyhedral cone is

F C (j) (q) m (j) * = v ∈ R I n |n (j) T v ≥ μ (j) max k=1,m (j) d (j) T k v .
Similar relations follow for the total cone and its dual.

4.3. The pointed friction cone assumption. In the following, we assume that we have a uniformly pointedness assumption that holds for both the primal and the dual cones.

Assumption A5 (the Uniformly pointed friction cone assumption). ∃ K , K * , and t(q, ) ∈ F C (q) and v(q, ) ∈ F C * (q), such that, ∀q ∈ R I n , and ∀ ∈ [0, ], we have that

• t(q, ) T w ≥ K ||t(q, )|| ||w||, ∀w ∈ F C (q).
• n (j) T v(q, ) ≥ μ t

(j) T 1 v(q, ) + t (j) T 2
v(q, ) + K * ||v(q, )||, for j = 1, 2, . . . , p. With this assumption, we can state the following result. Corollary 4.1. If Assumption A5 holds, then the following are true:

• t(q, ) T w ≥ K ||t(q, )|| ||w||, ∀w ∈ F C m (q, ).

• n (j) T v(q, ) ≥ μ (j) max k=1,2,...,m (j) {d (j) T k v(q, )}, for j = 1, 2, . . . , p. Proof Since F C m (q) ⊂ F C (q), and

d (j) k (q) T v(q, ) ≤ t (j)T 1 v(q, ) 2 + t (j)T 2 v(q, ) 2 , k = 1, 2, . . . , m (j) ,
for j = 1, 2, . . . , p, the conclusion follows.

The pointed friction cone assumption has one important consequence.

Lemma 4.2. Assume Assumption A5 holds. Let τ = Φ (j) (q)≤ c

(j)

n n + β (j) 1 t (j) 1 + β (j) 2 t (j) 2 . Then ∃K c such that c (1) n , β (1) 1 , β (1)
2 , c (2) n , β

≤ K c ||τ || , where the impulses corresponding to contact constraints that are not active are replaced with 0.

Proof Let v(q, ) be the vector that defines the pointedness property. We multiply with it through the definition of τ to get that

v(q, ) T τ = p j=1 c (j) n n (j) T v(q, ) + β (j) 1 t (j) T 1 v(q, ) + β (j) 2 t (j) T 2 v(q, ) ≥ Φ (j) ≤ μ (j) c (j) n t (j)T 1 v(q, ) 2 + t (j)T 2 v(q, ) 2 - β (j) 1 , β (j) 2 2 t (j)T 1 v(q, ) 2 + t (j)T 2 v(q, ) 2 + Φ (j) ≤ K * ||v(q, )|| c (j) n ≥ K * ||v(q, )|| Φ (j) (q)≤ c (j) n .
On the other hand, from the Cauchy-Schwarz inequality, we obtain that v(q, ) T τ ≤ ||v(q, )|| ||τ ||, which leads to

p j=1 c (j) n ≤ 1 K * ||τ || . (4.2) Since β (j) 2 1 + β (j) 2 2 ≤ μc (j)
n , j = 1, 2, . . . , p, the conclusion follows, after choosing

K C = 1 K * 1 + max j=1,2,...,p μ (j) 2 .
5. The time-stepping scheme. We now present the time-stepping scheme. This scheme was defined in [START_REF] Anitescu | A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF] and was shown to result in constraint stabilization.

We are going to define it first for the polyhedral approximation to the friction cone, and then for the full circular friction cones. Rewriting the classical continuous model coresponding to the polyhedral friction cone, we obtain [START_REF] Stewart | An implicit time-stepping scheme for rigid-body dynamics with inelastic collisions and Coulomb friction[END_REF][START_REF]Rigid-body dynamics with friction and impact[END_REF]

M dv dt = p j=1 n (j) ĉ(j) n + m (j) j=1 β(j) k d (j) k + f c (q, v) + k(q, v). q = v ĉ(j) n ≥ 0 ⊥ Φ (j) (q) ≥ 0, j = 1, 2, . . . , p. β(j) = argmin || β(j) || 1 ≤μc n , β(j) ≥0 m (j) k=1 v T d (j) k β(j) k , where β(j) = β(j) 1 , β(j) 2 , . . . , β (j)
m . In this work we do not go through the process of deriving the scheme as a convex approximation of a nonconvex LCP-based scheme, which was detailed in [START_REF]A fixed-point iteration approach for multibody dynamics with contact and friction[END_REF]. We simply define the scheme and then explain our choice through a particular microscopical interpretation of Coulomb friction.

We define the scheme for the polyhedral approximation. We start at the time t (l) , position q (l) , and velocity v (l) with time step h. The scheme is expressed by the following LCP:

M (v (l+1) -v l ) = j∈A(q (l) , ) m (j) k=1 β (j) k (n (j) + μ (j) d (j) k ) + hf c (q (l) , v (l) ) + hk(q (l) , v (l) ) 0 ≤ 1 h Φ (j) (q (l) ) + ∇Φ (j) T v (l+1) + μ (j) d (j) T k v (l+1) ⊥ β (j) k ≥ 0, j ∈ A(q (l) , ), k = 1, 2, . . . m (j) c (j) n = m (j) j=1 β (j) k q (l+1) -q (l) = hv (l+1) (5.1)
where A(q, ) = j j ∈ {1, 2, . . . , p} , Φ (j) (q) ≤ .

(5.2)

Note that β

(j)
k have meaning of impulses now and that the value of the normal impulse is determined as a function of the values of the tangential impulses. To simplify our notation, we will sometimes use the aggregate normal impulse vector and tangential impulses vector:

cn = c (1) n , c (2) n , . . . , c (p) n T , β = β (1) T , β (2) T , . . . , β (p) T T .
Under the pointed friction cone assumption, it has been shown that the scheme will always produce a solution [START_REF] Anitescu | A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF]. The scheme (5.1) produces the same solution as the original unrelaxed scheme when there is no slip at the contacts [START_REF] Anitescu | A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF].

An important observation, from a computational point of view, is that the solution at step (l) of (5.1) is the primal-dual solution of the following quadratic program:

min 1 2 v T M v + k (l) T v subject to 1 h Φ (j) (q (l) ) + ∇Φ (j) T v (l+1) + μ (j) d (j) T k v (l+1) ≥ 0, j ∈ A(q (l) , ), k = 1, 2, . . . , m (j) , (5.3)
where the linear term is defined by k (l) = -M v (l)hf c (q (l) , v (l) )hk(q (l) , v (l) ).

(5.4)

One advantage of this formulation is that it also suggests the corresponding formulation when the full circular cone is considered:

min 1 2 v T M v + k (l) T v subject to ∇Φ (j) T v (l+1) -μ (j) t (j) T 1 v (l+1) 2 + t (j) T 2 v (l+1) 2 + 1 h (l) Φ (j) (q (l)
) ≥ 0 j ∈ A(q (l) , ), k = 1, 2, . . . , m (j) .

(5.5) Lemma 5.1. Let v (l+1),m be the solution of (5.3) and v (l+1) be the solution of (5.5). Then lim m→∞ v (l+1),m = v (l+1) . Here by m (1) , m (2)

, • • • , m (p) = m → ∞ we mean that each component goes to ∞.
Proof It is an immediate consequence of the pointed friction cone assumption and Lemma 4.2.

5.1. Physical meaning of the scheme. The equivalence between (5.3) and (5.1) is immediate on the basis of the optimality conditions for (5.3). Hence, we focus our attention on (5.1). The essential difference between the time stepping scheme from this work (5.1) and previous LCP time stepping schemes [START_REF] Stewart | An implicit time-stepping scheme for rigid-body dynamics with inelastic collisions and Coulomb friction[END_REF]6] resides in the constraint

1 h Φ (j) (q (l) ) + ∇Φ (j) T v (l+1) + μ (j) d (j) T k v (l+1) ≥ 0, j ∈ A(q (l) , ), k = 1, 2, . . . , m (j) .
This constraint is assembled from two effects: the constraint linearization portion and friction treatment potion.

1. The constraint linearization portion addresses the replacement of Φ (j) (q) ≥ 0 by its linearization, based on the update q (l+1) = q (l) + hv (l+1) . Here we use the approximation Φ (j) (q (l+1) ) ≈ Φ (j) (q (l) ) + h∇Φ (j) T (q (l) )v (l+1) ≥ 0. This is the basis for the constraint stabilization effect. Under the assumptions from this work, we have shown that the total energy of the system stays bounded on every finite time interval and that the geometrical constraints are satisfied with one order higher than the order of the method [START_REF] Anitescu | A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF], which is one way to describe constraint stabilization [8]. 2. For the friction portion, we consider a two-dimensional configuration with one contact. It is well known that the friction coefficient models the average behavior of asperities at the contact surface between two bodies. In an idealized approach presented in Figure 5.1, we assume that the bodies in contact are very rigid and that the friction coefficient appears from a see-saw-like profile of the two surfaces. Then the constraints of motion are precisely

∇Φ (j) T v + μ (j) d (j) T 1 v ≥ 0, ∇Φ (j) T v + μ (j) d (j) T 2 v ≥ 0,
where d 1 = -d 2 are the tangent vectors (in generalized coordinates). This approach creates a vertical motion that is inexistent in the original friction and contact constraints. This motion is in effect existent in the microscale but will be highly exaggerated by our approach, since the timesteps taken will be orders of magnitude larger the ones for which the real vertical motion occurs. Nevertheless, we arctg( μ) show that this effect disappears for many interesting examples as the timestep goes to 0 and that the solution of the relaxed time-stepping scheme satisfies the measure differential inclusion.

In some sense, our model is excessively tuned toward static friction. This is also visible from our microscopic interpretation, where we assume that the asperities have no backward/forward motion with respect to the bodies on which they appear. That is not true for real contact: clearly, as the objects separate at the scale of the asperity, there is a pressure that would cause the asperity to bend toward the direction of motion and thus lowering the friction coefficient. This is a well-known effect in dynamics: the dynamic coefficient of friction is lower than the static coefficient of friction.

In effect, a behavior where the friction coefficient is dynamically dependent can be accommodated by solving at every time step the optimization problem min

1 2 v T M v + k (l) T v subject to 1 h Φ (j) (q (l) ) + ∇Φ (j) T v (l+1) + μ (j) ( (t T 1 v) 2 + (t T 2 v) 2 ) ≥ 0, j ∈ A(q (l) , ),
where μ(•) is a continuous concave mapping. Such an approach would dramatically reduce the amount of vertical motion due to our approach.

Unfortunately, the feasible set now becomes nonconvex, which creates a new difficulty. Because of the lack of appropriate software (since the problem is not convex and has conic constraints), we have not used this model for our simulations.

Boundedness of the velocity solution.

The key toward proving convergence of the numerical solution of our scheme to the solution of the appropriate measure differential inclusion is to show that the numerical solution is uniformly bounded.

We introduce a measure of the geometrical infeasibility. Of main concern is the behavior of the infeasibility of the noninterpenetration constraint, since it is divided by h in our model. We define the following measure of constraint infeasibility:

I(q) = max 1≤j≤p Φ (j) -(q) . (6.1)
Here for a function g(•), we define by g -(•) its negative part, that is, g -(•) = max{0, -g(•)}. The boundedness of the velocity sequence of the numerical solution is provided by the following result. Theorem 6.1. Consider the time-stepping algorithm defined above with the choice of active set defined by (5.2). The algorithm is applied over a finite time interval [0, T ], with a time step 0 < h. It is assumed that the system satisfies Assumptions A1-A5 and that the system is initially feasible, that is, I(q (0) ) = 0.

Then, there exist H > 0, V > 0, and C c > 0 such that, whenever in addition to the requirements above we have that h < H, ∀ l, 0 ≤ l ≤ N -1, we also have that [START_REF] Anitescu | A constraint-stabilized time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF]Theorem 4.1]. It can immediately be seen that Assumptions A1-A5 imply Assumptions (A1) and (D1-D3) from the same reference.

1. v (l) ≤ V , ∀1 ≤ l ≤ N and 2. I(q (l) ) ≤ C c v (l) 2 h 2 l-1 , ∀1 ≤ l ≤ N. Proof This result follows from
7. Closed cone-valued mapping. In the following section, we use Assumption A2. We have the following result.

Theorem 7.1. Under Assumptions A2 and A5, the mapping (q, ) → F C (q) and (q, ) → F C m (q) are closed.

Proof Choose q (n) , (n) , and let

t (n) ⊂ F C (n) (q (n)
). Recall that we defined the active set at q as A(q, ) = j ∈ {1, 2, . . . , m} |Φ (j) (q) ≤ .

Therefore, ∃ c

(j),n n , j ∈ A(q (n) , ), β (j),n 1 , β (j),n 2
, j ∈ A(q (n) , (n) ) and k = 1, 2, . . . , m (j) such that β

(j),n 1 2 + β (j,n) 2 2 ≤ μ (j) c (j),n n and t n = Φ (j) (q (n) )≤ (n) nc (j),n n (q n ) + β (j),n 1 t (j),n 1
(q n ) + β (j,n) 2 t (j) 2 (q n ).
Assume now that the sequence q n , n , t n has an accumulation point in R I n , R I , R I n , which we denote by q * , * , t * . Since the mappings Φ (j) (q) are continuous, we will have that j|Φ (j) (q n ) ≤ n ⊂ j|Φ (j) (q) ≤ * for n sufficiently large. By eventually padding the set with 0 we can write that

t n = Φ (j) (q (n) )≤ * nc (j),n n (q n ) + β (j),n 1 t (j),n 1 (q n ) + β (j,n) 2 t (j) 2 (q n ).
Since F C(q) is pointed, it follows from Lemma 4.2 that the sequence

(c (j),n n ), β (j),n 1 , β (j),n 2
is bounded, and thus it admits a convergent subsequence to c

(j), * n , β (j), * 1 , β (j), * 2 such that t * = Φ (j) (q * )≤ * nc (j), * n (q * ) + β (j), * 1 t (j),n 1 (q * ) + β (j), * 2 t (j) 2 (q * ) ⊂ F C(q * ).
This proves that (q, ) → F C (q) is a closed mapping. The fact that (q, ) → F C m (q, ) is a closed mapping for fixed m is proved much in the same way.

8. Measure differential inclusions. In this section we prove the convergence of our scheme. The approach used here follows, in general lines, [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF], with one important exception explained below.

8.1. The pointwise quantities whose convergence is investigated and the main theorem. Let q (l),h be the value of the position vector computed at step (l) of the method using step size h, and similarly for velocities v (l),h and similarly for the quantities c(l+1),h n and β(l+1),h . The value q (h) (t) is the linear interpolant of q (h) (lh) = q (l);h and q (h) ((l + 1)h) = q (l+1);h for t ∈ [lh, (l + 1)h]. For the velocities, set v (h) (t) = v (l+1);h for t ∈ (lh, (l + 1)h]. For the velocities, set v h (t) = v (l+1);h for t ∈ (lh, (l + 1)h]. Then

q (h) (t) = q (h) (0) + t 0 v (h) (τ )dτ
for all t > 0 and h > 0.

The functions c(h) n and β(h) are defined for all h > 0 as follows. c(h) n is a sum of Dirac δ functions with strength c(l+1);h n at the time t l = lh. Similarly, β(h) is a sum of Dirac δ functions with strength β(l+1),h at time t (l);h . Let

n (h) (t) = n(q (h) (t)), D (h) (t) = D(q (h) (t)), and k (h) (t) = k(q (h) (t (l) ), v (h) (t (l) )) for t l ≤ t ≤ t l+1 .
Since the non-penetration constraint is enforced differently from other time stepping schemes it does not immediately follow that the constraint

c (j) n ≥ 0 ⊥ Φ (j) ≥ 0, (8.1)
which is supposed to occur in the limit of h → 0, will be satisfied. The difficulty is that coefficient μ now enters the linearization of Φ (j) . However, we now show that, in the limit, the effect of the friction coefficient on the complementarity relation (8.1) drops to 0 and that the limit solution satisfies the same differential inclusion as the time stepping schemes from [START_REF] Stewart | An implicit time-stepping scheme for rigid-body dynamics with inelastic collisions and Coulomb friction[END_REF]6,[START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF]. Theorem 8.1. Let q h (t) and v h (t) be the pointwise position and velocities sequence that can be defined from the numerical sequence produced over a finite time interval [0, T ] by the time stepping scheme (5.1)-(5.2). Under Assumptions A1-A5, there exists a subsequence h k → 0 where

• q h k (•) → q(•) uniformly. • v h k (•) → v(•) pointwise a.e.
• dv h k (•) → dv(•) weak * as Borel measures in [0,T], and every such subsequence converges to a solution (q(•), v(•)) of the measure differential inclusion (3.2). Therefore, q(t), v(t) is a weak solution of our model.

Proof The proof follows the same outline as in [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF]. The only difference is that the friction cone considered here, which also includes the almost active constraints, is larger when h = 0.

Boundedness of the velocities. Following Theorem 6.1, we have that v (l) is uniformly bounded. This means that there exists V such that

v (l),h ≤ V ; ∀h < H.
Convergence of a subsequence of q h (•) that is feasible for the geometrical constraints. From dq h dt = v h , the boundedness of the velocities implies that the function q h (•) is uniformly Lipschitz in [0,T]. Combined with the condition q h (0) = q(0), this implies that the family q h (•) is equicontinuous and equibounded. By the Arzela-Ascoli theorem, there exists a uniformly convergent subsequence, which we also denote by q h (•), that converges q h (•) → q(•) uniformly in [0, T ]. In addition, from Theorem 6.1 we have that I(q(t)) = 0, that is, Φ (j) (q) ≥ 0, j = 1, 2, . . . , p.

Uniform boundedness of the variation of velocities. We now show that the numerical velocity function v h (•) has bounded variation uniformly in h > 0.

Lemma 8.2.

T 0 v h (•) is uniformly bounded as h → 0. Proof The proof is essentially the same as the one in [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF]Lemma 6]. The cones used are different, but the only information used by the cones is the closure of the mapping (q, ) → F C (q) and the pointed friction cone assumption.

Consider an arbitrary time instant t ∈ [0, T ]. Since F C(q(t)) is a pointed cone, so is F C m (q(t)) ⊂ F C(q(t)). Thus (since pointedness is an affine property) so are the cones M -1 F C(q(t)) and M (-1) F C m (q(t)). Thus there exists an unit vector n 0 (t) and a parameter (t) > 0 such that in any vector norm in R n ,

z ∈ M -1 F C(q(t)) =⇒ n 0 (t) T z ≥ (t) ||z|| .
By the closed graph property of F C(•), that follows from Theorem 7.1 there exists an η(t) > 0 and an h 0 > 0 such that h 0 > h > 0 and ||t -t|| ≤ (t) implies that

z ∈ M -1 F C(q h (t )) =⇒ n 0 (t) T z ≥ 1 2 (t) ||z|| .
Then, provided that both lh and (l + 1)h lie in [tη(t), t + η(t)], we obtain from the time stepping scheme (5.1)-(5.2) that

v (l+1);h -v (l);h ∈ hM -1 k (l);h + M -1 F C(q (l);h ). Write v (l+1);h -v (l);h ∈ hM -1 k (l);h + z (l+1);h (8.2)
for some z (l+1);h ∈ M -1 F C(q (l);h ). Then

n 0 (t) T v (l+1);h -v (l);h = hn 0 (t) T M -1 k (l);h + n 0 (t) T z (l+1);h ≥ hn 0 (t) T M -1 k (l);h + 1 2 (t) z (l+1);h . Set l min = (t -η(t))/h , l max = (t + η(t))/h . Then l=l max -1 l=l min n 0 (t) T v (l+1);h -v (l);h ≥ l=l max -1 l=l min hn 0 (t) T M -1 k (l);h + 1 2 (t) l=l max -1 l=l min z (l+1);h .
The left-hand sum telescopes to give the following inequality:

l=l max -1 l=l min n 0 (t) T v (l+1);h -v (l);h = n 0 (t) T v (l max );h -v (l min );h ≤ v (l+1);h + v (l);h .
By choosing η(t) > 0 so small that t + η(t) < T , we obtain that v h (•) ≤ C 2 uniformly as h ↓ 0 on [tη(t), t + η(t)]. Similarly, from Assumption 3, there is a bound

M -1 k l;h ≤ C 3 uniformly as h ↓ 0 on [t -η(t), t + η(t)]. Thus, 1 2 (t) l=l max -1 l=l min z (l+1);h ≤ 2C 2 + h(l max -l min )C 3 ≤ 2C 2 + 2η(t)C 3 .
This gives the bound, (uniformly as h ↓ 0),

l=l max -1 l=l min z (l+1);h ≤ 2 (t) (2C 2 + 2η(t)C 3 ) .
By (8.2), we obtain that

l=l max -1 l=l min v (l+1);h -v (l);h ≤ 2η(t)C 3 + 2 (t) (2C 2 + 2η(t)C 3 ) , which implies that t+η(t)/2 t-η(t)/2 v h (•)
is uniformly bounded as h ↓ 0. Here, we have used the usual notation for the total variation.

b a v h (•) = sup m; a=t 0 <t 1 <...t m =b m i=1 v h (t i ) -v h (t i-1 ) . Since t -η(t) 2 , t + η(t) 2 |t ∈ [0, T ] , is a covering of [0, T ], there exists a finite subcovering t i -η(t i ) 2 , t i + η(t i ) 2 |i = 0, 1, 2, .
. . , m . Summing the contributions over these subintervals give an uniform bound on T 0 v h (•) as h ↓ 0. Using the fact that v h (•) has bounded variation, then, by Helly's selection theorem, there exists a subsequence of v h k (•) of v h (•) that converges pointwise to v(•) and has bounded variation. The coresponding functions q (h k ) (•) converge to the indefinite integral of v(•) by the pointwise convergence theorem for Lebesgue integrals. We assume for simplicity that this is the entire sequence and therefore q h (•) → q(•) and

v h (•) → v(•).
Weak * convergence. Since T 0 v h (•) are uniformly bounded as h → 0 and v h (0) = v(0) and since v h (•) → v(•) pointwise, it follows that dv h → dv weakly *, that is

T 0 φ(t) T dv h (t) → T 0 φ(t) T dv(t)
for all continuous functions φ(t). Therefore, dv h (•) → dv(•) weak * as Borel measures.

Convergence to the measure differential inclusion. From the time-stepping scheme definition (5.1) we see that, from Assumption A2, it follows that

n T v (l+1) + μd T v (l+1) ≤ V K,
where K is some fixed constant.

Using the complementarity conditions from (5.1), we obtain that Φ (j) (q) ≥ hV K =⇒ β (j) = 0. Therefore, the time-stepping scheme (5.1) is equivalent to the following time-stepping scheme:

M (v (l+1) -v l ) = Φ (j) (q (l) )≤V Kh m (j) k=1 β (j) k (n (j) + μ (j) d (j) k ) + hf c (q (l) , v (l) ) + hk(q (l) , v (l) ) 0 ≤ 1 h Φ (j) (q (l) ) + ∇Φ (j) T v (l+1) + μ (j) d (j) T k v (l+1) ⊥ β (j)
k ≥ 0, j = 1, 2, . . . p, k = 1, 2, . . . , m (j) q (l+1)q (l) = hv (l+1) .

(8. [START_REF]A fixed-point iteration approach for multibody dynamics with contact and friction[END_REF] In the rest of the proof we will follow very closely the proof of [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF]Lemma 7], except at one crucial point, that will be outlined bellow.

Following (8.3), we obtain that the numerical solutions (q h (•), v h (•)) satisfy the measure differential inclusion

M dv h dt -f h c (q h (•)) -k h (t)dτ h (t) ∈ μ F C m V Kh (q h (•)) ⊂ F C V Kh (q h (•)) dq h dt = v h ,
where k h (t) = k(q h (t), v h (t)), and τ h (t) = (t/h)h . Let k(t) = k(q(t), v(t)). Letting dt denote the Lebesgue measure on R, we show that

k h (•)dτ h (•) → k(•)dt weak * Consider φ a C 1 function on R with compact support. Then φ(•)k h (•)τ h (•)
is a function of bounded variation, with associated measure

dφk h τ h = φ k h τ h dt + φτ h dk h + φk h dτ h ,
where τ h is continuous from the left and k h is continuous from the right. Thus

φk h dτ h = d(φk h τ h ) -φ k h τ h dt -φτ h dk h . Now k h (t) → k(t) for all t, with k h (•) L ∞ uniformly bounded as h ↓ 0. Further, τ h (•) → (•) uniformly as h ↓ 0. Thus φ (•)k h (•)τ h (•) -→ (•)
pointwise, and thus is in L 1 by the dominated convergence theorem. Also, φ (•)τ h (•) converges uniformly as h ↓ 0. Noting that

d(φk h τ h ) = 0 = d(φkτ h )
because φ has compact support, we obtain

lim h↓0 φk h dτ h = -lim h↓0 φ k h τ h dt -lim h↓0 φτ h dk h = -φ (t)k(t)tdt -φ(t)tdk(t) = φ(t)k(t)dt.
To show that lim h↓0 φk h dτ h = φkdt for all continuous φ with compact support, we take C 1 approximations φ where the support of φ is within distance of the support [t 1 , t 2 ] of φ and φ → φ uniformly as → 0. This can be done by using convolution with mollifiers, for example. Then

lim sup h↓0 φk h dτ h -φkdt ≤ lim sup h↓0 φ k h dτ h -φ kdt + lim sup h↓0 φ k h dτ h -φk h dt + lim sup h↓0 φ kdτ h -φkdt ≤ 2 ||φ -φ|| L ∞ lim sup h↓0 k h (•) ∞ L (t 2 -t 1 + 2 ).
The norm k h (•) L ∞ is taken over the interval [t 1 -, t 2 + ]. Now taking → 0 gives lim sup h↓0 φk h dτ hφkdt = 0 and thus k h dτ h → kdt as h → 0, as required. At this point, however, our proof departs from the proof of [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF]Lemma7], because in this work we consider a larger friction cone than the one in [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF]Lemma7].

Same as in [30, Lemma 7] we will invoke the following result, whose statement we included here.

[32, Theorem4] Suppose that q n (•) are continuous, v n (•) have uniformly bounded variation and k n (•) are uniformly bounded, all on [0, T ], and q n (•) → q(•) uniformly,

v n (•) → v(•) pointwise and k n (•) → k(•) pointwise. Suppose also that K : R n ⇒ C(R n ) has closed graph, min {||z|| |z ∈ K(w)} is uniformly bounded and K(w) is pointed for all w ∈ R n . Then if dv n dt (t) ∈ μ K(q n (t)) -k n (t)
for all n, then the limit satisfies dv dt (t) ∈ μ K(q(t))k(t).

In our case, the requirement that min {||z|| |z ∈ K(w)} is uniformly bounded is immediately satisfied since K(w) are cones that always contain the 0 element. Nevertheless, as opposed to [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF]Lemma 7], this result cannot be applied directly since the cone from (8.1) also depends on the timestep.

We introduce the following sequences:

q h (•) = q (h) (•) V Kh v h (•) = v (h) (•) 0 M = M 0 0 1 k h (•) = k (h) (•) 0 F C( q h (•)) = F C V Kh (q h (•)) × {0}
and the corresponding limits

q(•) = q(•) 0 v(•) = v(•) 0 k(•) = k (h) (•) 0
It is then immediate that P1 q h (•) converges uniformly to q(•).

P2 de q h (•) dt = v h (•), v h (•)
is uniformly bounded and has bounded variation, and it converges pointwise to v(•). P3 k h (•) converges pointwise k(•) and is uniformly bounded. P4 From Theorem 7.1, the cone-valued mapping F C( q) has a closed graph. P5 For any q, the cone F C( (q)) is pointed. P6 The quantities q h , v h , k h , satisfy the following differential inclusion.

M de v h dt -k h (t)dτ h (t) ∈ μ F C m ( q h (•)) ⊂ F C( q h (•)) e q h
dt = v h , Since the properties P1-P6 hold, we can use the above result, [START_REF]Reformulations of measure differential inclusions and their closed graph property[END_REF]Theorem4] to claim that the limits will satisfy the similar measure differential equation. Discarding the last component of q(•) and v(•), and using the relations

F C( q(t)) = F C 0 (q(t)) = F C(q(t)),
we obtain that the limits satisfy the measure differential inclusion.

M dv dt -k(t)dt ∈ μ F C m (q(t)) ⊂ F C(q(t)) (8.4)
To complete this subsection, note that for any T ≥ t 2 > t 1 ≥ 0, we have that

q h (t 2 ) -q h (t 1 ) = t 2 t 1 v h (τ )dτ.
Since this holds for all t 2 > t 1 in [0, T ], we obtain that

dq dt = v,
10. Numerical Results. In this section we compare the optimization-based method that we have analyzed in this work with the coresponding LCP-based method. We consider two types of examples. The first group of examples compares the optimization-based method developed in this work with the LCP-based method [6,[START_REF] Stewart | An implicit time-stepping scheme for rigid-body dynamics with inelastic collisions and Coulomb friction[END_REF].

The second group investigates the potential of this method for the simulation of granular materials that can easily become very large scale computations. 10.1. Comparison between the LCP-based method and the optimization based method. We consider the following two-dimensional examples. All the examples include the effect of gravity.

1. Particle of 0 diameter falling on a flat, horizontal line. The initial position of the particle is x = 0, y = 3. The initial velocity is ẋ = 3, ẏ = 0. There is friction between the particle and the line with friction coefficient μ = 0.3. The mass of the particle is 1. We have simulated this with both the LCP algorithm and the optimization algorithm with time steps shown in Table 10.1. In Table 10.1 we also give the scaled values of the 2-norm, which show that in the limit both algorithms produce the same solution sequence. A graph of the trajectory is presented in Figure 10.1. 2. Particle of 0 diameter starting in contact with a flat, horizontal line. The initial position of the particle is x = 0, y = 0. The initial velocity is ẋ = 3, ẏ = 0. There is friction between the particle and the line with friction coefficient μ = 0.3. The mass of the particle is 1. A comparison between the y-positions produced by the LCP algorithm and the QP algorithm for time step 0.001 is presented in Figure 10.2. Although the differences are small in absolute terms, there is a difference between the solutions that does not disappear as the time step goes to 0. The trajectory of the particle is presented in Figure 10.2. 3. Bar of length 2 dropped on a flat, horizontal line. At the collision point, the data obtained are the same as for a well-known example of the Painleve Paradox [9,[START_REF]Rigid-body dynamics with friction and impact[END_REF]. Much like in the dropped-ball case, we see that when there is no initial contact, the limit solution produced by the optimization-based scheme is the same as the solution of the LCP-based scheme, as can be seen from the Table 10.2. The trajectory of the bar with friction that starts without initial contact is presented in Figure 10.3. 4. Bar of length 2 dropped on a flat, horizontal line and starting in contact with the horizontal line. The data obtained are the same as for a well-known example of the Painleve Paradox [9, 31]. The optimization scheme and the LCP scheme produce close but different solutions, as can be seen from the Figure 10.5, where we present the prediction of the position of the tip of the bar with both methods. In either case, the figures show that the difference between the two methods is small in all cases and 0 when the contacts are not initially active. The fact that the two simulation methods produce different results is not surprising. It was well established that the measure differential inclusion will not necessarily have unique solutions [START_REF]Rigid-body dynamics with friction and impact[END_REF].

An intriguing fact, however, is the fact that the two simulations converge to the same solution when there is no initial contact. We plan to further investigate this observation to determine whether it can be stated a with reasonable generality. 10.2. Size-based segregation of granular materials. In the example of the Brazil nut effect [START_REF] Rosato | Why the brazil nuts are on top: Size segregation of particulate matter by shacking[END_REF], also communicated in [START_REF] Hart | A hard constraint time-stepping approach for rigid multibody dynamics with joints, contact and friction[END_REF], we have smaller particles shaken together with a large particle that, after a while, emerges on top. This behavior is characteristic of granular matter. This effect is seen in four frames of the simulation in Figure 10.6. That example contains 201 bodies and is simulated with the optimization method for 75 seconds with time step of 0.1 s and friction coefficient of 0.5 at all interactions. Our time step compares favorably with traditional molecular dynamics approaches that use a penalty method and that need time steps on the order of microseconds [START_REF] Hryciw | Video tracking for experimental validation of discrete element simulations of large discontinuous deformation[END_REF]. Currently, we do not have any molecular dynamics result to compare with, since most molecular dynamics simulations of size-based segregation use periodic boundary conditions [START_REF] Rosato | Why the brazil nuts are on top: Size segregation of particulate matter by shacking[END_REF], whereas we use hard walls. In any event, the number of shakes needed to get the larger body emerge to the top is comparable to the number of shakes needed by a molecular dynamics simulation that uses periodic boundary conditions and temperature to model the dynamics of the bodies (between 30 and 50) [START_REF] Rosato | Why the brazil nuts are on top: Size segregation of particulate matter by shacking[END_REF].

11. Conclusions and future work. We have shown that the solution of a scheme that progresses with fixed time step and solves only a convex quadratic program per step converges to the weak solution of the rigid multibody dynamics problem with contact and friction as defined in [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF][START_REF]Rigid-body dynamics with friction and impact[END_REF]. In addition, the scheme decreases the energy in the limit if the external force is a potential force. This shows that solution to the most general weak formulations of this problem can be produced relatively efficiently.

The only result from [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF] that needs to be extended to our scheme is the weak form of the equation ||v T || μc nβ 2 1 + β 2 2 = 0 from the Coulomb model (2.2). We note that even in [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF] the result was proven in the general case only for one inelastic contact, therefore proving such a result for an arbitrary number of contacts in our case would mean much more than simply expanding the results from [START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF] to this work. Two of our simulations, whose trajectory is plotted in Figures 10.2 and 10.5 are indicating that that result will not hold for our scheme. However, though it is desirable for this result to hold, not having proved it does not invalidate our conclusion of convergence to a measure differential inclusion, since the latter does not depend on this part of the Coulomb model. We nevertheless plan to provide a more refined investigation of the differences between the schemes in this work and the schemes in [START_REF] Stewart | An implicit time-stepping scheme for rigid-body dynamics with inelastic collisions and Coulomb friction[END_REF]6,[START_REF] Stewart | Convergence of a time-stepping scheme for rigid body dynamics and resolution of painleve's problems[END_REF].

Additional extensions that we are investigating are (1) expanding the results in this work for the case of partially elastic and totally elastic collisions and (2) accommodating piecewise differentiable constraints that are widely used in the computer graphics community. Also, we are interested in investigating whether our observation that the optimization-based scheme and LCP-based scheme produce the same initial Table 10.1 Differences between the LCP method and the optimization method for one particle with friction and no initial contact . results when all bodies start without contact can be stated with some generality. The submitted manuscript has been created by the University of Chicago as Operator of Argonne National Laboratory ("Argonne") under Contract No. W-31-109-ENG-38 with the U.S. Department of Energy. The U.S. Government retains for itself, and others acting on its behalf, a paid-up, nonexclusive, irrevocable worldwide license in said article to reproduce, prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the Government.

Fig. 5 . 1 .

 51 Fig. 5.1. Successive blowup of the microscopic interpretation of our model

Fig. 10. 1 .Fig. 10. 2 .

 12 Fig. 10.1. Trajectory for particle with friction example; no initial contact.

Fig. 10. 3 .

 3 Fig. 10.3. Trajectory for bar with friction example; no initial contact.

Fig. 10. 4 .

 4 Fig. 10.4. Trajectory for the Painleve example; bar with friction with initial contact and nonzero tangential velocity (computed by the LCP method).

Fig. 10. 5 .

 5 Fig. 10.5. Comparison between the y positions of the tip that is initially in contact for the Painleve example.

  lems with coulomb friction, Zeithschrift fur Angewandte Mathematik und Mechanik, 77 (1997), pp. 267-279. [36] D. Volfson, L. S. Tsimring, and I. S. Aranson, Partially fluidized shear granular flows: Continuum theory and MD simulations, Phys. Rev. E, 65 (2003).

Acknowledgments. We thank Michael Ferris and Todd Munson for providing and maintaining PATH [24,11], a package for solving general linear complementarity problems. The size segregation of granular materials example was developed in colaboration with Gary Hart from the University of Pittsburgh. This work was supported by the Mathematical, Information, and Computational Sciences Division subprogram of the Office of Advanced Scientific Computing Research, Office of Science, U.S. Department of Energy, under Contract W-31-109-ENG-38.

as required.

Since q h (•) → q(•) uniformly and v h (•) → v(•) pointwise on [0, T ] and dv h (•) → dv(•) weak *, and using the fact that the limit satisfies (8.4), we obtain that q(•), v(•) satisfy the measure differential inclusion M dvk(q, v)dt ∈ μ F C(q), dq dt = v.

This completes the proof of Theorem 8.1.

9. Exact dissipation of solutions. We now assume that the external force k c (q, v) originates in a potential, that is, k(q, v) = -∇V (q), where V is a smooth function of q. In this section we show that the numerical scheme satisfies the exact dissipativity in the limit.

Lemma 9.1. If q(•) and v(•) are the limit of a subsequence

Proof Consider a sequence for which l h k is the largest integer that satisfies l h k h k ≤ t. From the definition of the time stepping scheme (5.1)-(5.2), we have that v h (t) = v (l+1);h and

After rewriting the last sum, we obtain

We have

as a result of the smoothness of V (q) and the uniform boundedness of V (q). We replace this in (9.1) and obtain

In the last equation, we used assumption A4, with C V an upper bound on ||∇V (q)||.

Adding from s = 0 to s = l all equations like the preceding one, we obtain that

We take the limit as h → 0, and since we consider piecewise constant velocities v l,h that have bounded variation, we obtain that, at all points where v (h) (t) → v(t) and q (h) (t) → q(t), and since v(•) is continuously differentiable,

The proof is therefore complete.

The preceding result shows that the energy does not increase and that the scheme is indeed dissipative. Therefore, the scheme mimics an important feature of the physical phenomenon that we are simulating. Table 10.2 Differences between the LCP method and the optimization method for the bar with friction friction and no initial contact.