Mikaël Barboteu
email: barboteu@univ-perp.fr

Serge Dumont

A primal-dual active set method for solving multi-rigidbody dynamic contact problems

Keywords: Granular media, unilateral constraint, rigid body, discrete element method, nonsmooth contact dynamics, semismooth Newton method, primal-dual active set, bipotential, augmented Lagrangian, multi-body contact, numerical simulations

In this work, an active set type method is considered in order to solve a mathematical problem that describes the frictionless dynamic contact of a multi-body rigid system, the so-called nonsmooth contact dynamics (NSCD) problem. Our aim, here, is to present the local treatment of contact conditions by an active set type method dedicated to NSCD and to carry out a comparison with the various well-known methods based on the bipotential theory and the augmented Lagrangian theory. After presenting the mechanical formulation of the NSCD and the resolution of the global problem concerning the equations of motion, we focus on the local level devoted to the resolution of the contact law. Then we detail the numerical treatment of the contact conditions within the framework of the primal-dual active set strategy. Finally, numerical experiments are presented to establish the efficiency of the proposed method by considering the comparison with the other numerical methods.

Introduction

The simulation of dynamic multi-body contact problems is involved in many engineering applications, such as granular media, masonry, tensegrities, geomaterials, robotic manipulation and mechanical deformable systems. In the case of both deformable bodies and rigid bodies, the numerical solution of multi-contact problems remains a difficult and nontrivial task because contact law is strongly nonlinear, nonsmooth and multivalued. Moreover, the numerical methods used to solve multi-body contact problems depend, in particular, on the mechanical behaviour (deformable or rigid) of the bodies involved.

In most cases, the bodies considered are assumed to be rigid. In a rigid-body model, the contact dynamics method, also called the nonsmooth contact dynamics (NSCD) method, is based on a discrete element method initially developed for the simulation of granular materials. For a survey on the subject, we can refer to Radjai and Richefeu [START_REF] Radjai | Contact dynamics method[END_REF]. This very well-known method emerged from a mathematical formulation of nonsmooth dynamics and subsequent algorithmic developments realized by Moreau and Jean [START_REF] Moreau | Application of convex analysis to some problems of dry friction[END_REF][START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Moreau | Numerical aspect of sweeping process[END_REF][START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF][START_REF] Jean | The non smooth contact dynamics method[END_REF]. In this context, it is important to note that the contact dynamics method deals with the following two major issues: the contact law based on complementarity relations between the contact (or impulse) forces and velocities, and nonsmooth motion involving velocity jumps with impulse forces. To give a consistent description of the dynamics, the discrete element method is based on two main computational level tasks: an implicit time-stepping scheme is used to solve the equations of motion at the global level, while an explicit local treatment deals with the evolution of the contact network between rigid bodies. The most commonly used discrete element method is based on the nonlinear Gauss-Seidel (NLGS) method developed by Jean and Moreau [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Moreau | Numerical aspect of sweeping process[END_REF][START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF] and it consists of an iterative process pertaining to the Gauss-Seidel method by considering, successively, each contact until convergence is achieved. More sophisticated methods, such as conjugate gradient-type methods, have also been developed (see, e.g., Renouf and Alart [START_REF] Renouf | Conjugate gradient type algorithms for frictional multicontact problems: applications to granular materials[END_REF]) to solve nonsmooth contact dynamics problems. For the local treatment of the contact conditions, several approaches can be considered, such as, for instance, methods based on the bipotential theory and the augmented Lagrangian theory [START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Jean | The non smooth contact dynamics method[END_REF][START_REF] Jourdan | A Gauss-Seidel like algorithm to solve frictional contact problems[END_REF][START_REF] Fortin | Simulation numérique de la dynamique des systèmes multicorps appliquée aux milieux granulaires[END_REF][START_REF] Fortin | An improved discrete element method based on a variational formulation of the contact law[END_REF][START_REF] Fortin | Numerical simulation of granular materials by an improved discrete element method[END_REF][START_REF] Dumont | On enhanced descend algorithms for solving frictional multi-contact problems: applications to the discrete element method[END_REF]. Furthermore, note that when dynamic contact systems are considered without friction, the resulting problem can be formulated as a convex quadratic problem and several other efficient methods based on quadratic programming can be used to solve these nonsmooth dynamical systems [START_REF] Nocedal | Numerical optimization[END_REF][START_REF] Acary | Numerical methods for nonsmooth dynamical systems: applications in mechanics and electronics[END_REF][START_REF] Cottle | The linear complementarity problem[END_REF].

When deformable bodies are considered, a time-stepping method and a finite-element method are generally used for the full discretization of the contact problem. A numerical treatment of the contact conditions can then be realized by several numerical approaches. A wide bibliography exists on this subject [START_REF] Wriggers | Computational contact mechanics[END_REF][START_REF] Laursen | Computational contact and impact mechanics[END_REF][START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]. For instance, Alart and Curnier [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF] have introduced an augmented Lagrangian formulation combined with a generalized Newtonian method to solve nondifferentiable but continuous equations arising from frictional contact problems. Several other methods have emerged during the last few years. Amongst them, active set strategies are very successful and are widely used because of their efficiency and simplicity of implementation. In the standard active set approach, which can be found in well-known books of optimization [START_REF] Nocedal | Numerical optimization[END_REF][START_REF] Cottle | The linear complementarity problem[END_REF][START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF], the methods are divided into primal active set methods and dual active set methods. When finding a feasible starting point, primal active set methods generate a sequence of primal feasible iterations until dual feasibility is achieved; hence, an optimal solution is obtained. Dual active set methods for convex quadratic problems generate a sequence of dual feasible iterations until primal feasibility is achieved; hence, an optimal solution is obtained. More recently, primal-dual active set methods have been considered to solve variational problems with unilateral constraints [START_REF] Hintermuller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF][START_REF] Hintermuller | The primal-dual active set strategy as a semismooth Newton method[END_REF][START_REF] Hueber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF][START_REF] Kunisch | Generalized Newton methods for the 2D-Signorini contact problem with friction in function space[END_REF][START_REF] Hintermuller | Constrained optimization for interface cracks in composite materials subject to non-penetration conditions[END_REF][START_REF] Hintermuller | Obstacle problems with cohesion: a hemi-variational inequality approach and its efficient numerical solution[END_REF][START_REF] Abide | Analysis of two active set type methods to solve unilateral contact problems[END_REF]. These approaches are characterized by the fact that the active set is defined by a relation described by both the primal and the dual feasibilities, which are enforced together during each iteration. Furthermore, within the framework of semismooth Newtonian techniques [START_REF] Hintermuller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF][START_REF] Hintermuller | The primal-dual active set strategy as a semismooth Newton method[END_REF][START_REF] Kunisch | Generalized Newton methods for the 2D-Signorini contact problem with friction in function space[END_REF], it can be proven that the primal-dual active set type methods lead to direct enforcement of the exact contact conditions on nodes related to the active and inactive contact sets.

The contribution of this paper is to propose a primal-dual active set method for the numerical treatment of the contact conditions within the resolution of dynamic multi-rigid-body contact problems. To our knowledge, this is the first time that a primal-dual active set method has been considered for solving such types of contact dynamics problem. Very few works have been devoted to this topic. We can cite, for instance, the work of Sharaf [START_REF] Sharaf | An active set algorithm for a class of linear complementarity problems arising from rigid body dynamics[END_REF], in which a very particular class of rigid-body dynamics problems is considered, and the work of Koziara and Bicanic [START_REF] Koziara | Semismooth Newton method for frictional contact between pseudo-rigid bodies[END_REF], for which a semismooth Newtonian method is proposed to solve problems dealing with pseudo-rigid bodies. Our goal in this work is to show the performance and efficiency of the primal-dual active set method for NSCD problems, compared with other effective methods based on the bipotential and augmented Lagrangian theories that have been recently developed [START_REF] Dumont | On enhanced descend algorithms for solving frictional multi-contact problems: applications to the discrete element method[END_REF]. For this purpose, several numerical examples of rigid-body contact problems are considered. Furthermore, a comparison of performances of the whole methods can be realized by the use of only one stopping criterion [START_REF] Dumont | On enhanced descend algorithms for solving frictional multi-contact problems: applications to the discrete element method[END_REF].

The remainder of the article is organized as follows. In Section 2, we recall the usual contact conditions commonly used in contact dynamics, and the main traits of the primal-dual active set type method within the framework of a deformable body system. In Section 3, the discrete element method context for solving a dynamic multi-rigid-body contact problem is presented by considering first the equations of motion and then the general algorithm for NSCD problems. After presenting methods based on the bipotential and the augmented Lagrangian theories, Section 4 is devoted to a numerical treatment of the dynamic contact conditions by a primal-dual active set method within the framework of dynamic multi-rigid-body contact problems. After that, in Section 5, we present several numerical simulations to illustrate the performances of the active set method compared with other effective methods. Finally, in Section 6, we conclude by discussing some prospects and future works in continuation of the present study.

Contact conditions and active set type methods

The aim of this section is to recall the usual contact conditions commonly used in contact dynamics. First, a brief presentation of the contact laws in the case of a deformable body is given; then the dynamic contact conditions for rigid-body systems are provided. Finally, the main traits of the primal-dual active set type method are also briefly recalled to solve dynamic contact problems in the framework of deformable bodies.

Contact conditions

First of all, to simplify the writing, let us consider the contact without friction of a deformable body with a rigid obstacle, the so-called foundation. Note that, in what follows, the contact between two or several bodies can be generalized quite easily by considering an effort on the formalism. For mathematical convenience, we assume that the body occupies a bounded domain ⊂ R d (d = 1, 2, 3), with a Lipschitz continuous boundary Ŵ. We denote by n the unit inner normal vector on Ŵ and t the associated tangent vector. We use the notation u and σ for the local displacement and the local stress tensor at one point of the domain , respectively. Also, we denote by u n and u t the normal and tangential components of u on Ŵ given by u n = u • n, u t = uu n n. Finally, σ n and σ t will represent the normal and the tangential stress on Ŵ, defined by σ n = (σ n) • n and σ t = σ nσ n n. Furthermore, a dot superscript represents the time derivative with respect to the time variable t, e.g., u = ∂u/∂t. Now, let us consider a potential contact between one point of Ŵ and the foundation. A potential contact point has the following dynamic content. As long as the normal distance u n between the body and the obstacle remains positive (corresponding to a gap), no force is activated and the normal force σ n is identically zero. But when u n = 0, a nonnegative (repulsive) normal force σ n is mobilized at the contact point and can take indefinitely large values, depending on the forces acting on the body. These conditions define a complementary relation, called Signorini's conditions [START_REF] Signorini | Sopra alcune questioni di elastostatica[END_REF], between u n and σ n . These frictionless contact conditions can be written following the Karush-Kuhn-Tucker conditions, as follows:

   u n ≥ 0, σ n ≥ 0, u n σ n = 0, (1)
σ t = 0. (2)
In addition, when contact dynamics are considered, the persistency condition has to be added for energy conservation purposes [START_REF] Laursen | Computational contact and impact mechanics[END_REF]. Indeed, to eliminate the work of the normal contact reaction at time t (W cont = Ŵ σ n un ), the following persistency condition has to be considered [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF][START_REF] Armero | Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems[END_REF]:

un σ n = 0. (3)
This condition means that the normal contact reaction can only appear during persistent contact. One can easily prove [START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Moreau | Numerical aspect of sweeping process[END_REF][START_REF] Laursen | Computational contact and impact mechanics[END_REF] that the addition of the persistency condition (equation ( 3)) to the unilateral contact law (equation ( 1)) gives the following contact dynamic conditions:

if u n > 0 , σ n = 0 (4) if u n = 0 ,    un ≥ 0, σ n ≥ 0, un σ n = 0, (5)
In the following, a couple (u n , σ n ) verifying this set of conditions is denoted by contact_law(u n , σ n ) = .true.

The previous alternatives in equations ( 4) and ( 5) lead to the condition of a complete contact law formulated by Moreau [START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Moreau | Numerical aspect of sweeping process[END_REF]. In the context of nonsmooth motion, the time derivative u is not unique and then the left-limit velocity and the right-limit velocity must be distinguished. Since the actual velocity at time t is immaterial, the question of velocity jump is crucial and is problematic when the step of time discretization is addressed. As a consequence, in a time-stepping formulation, the discrete forms of the complementarity condition formulated in displacement (u n σ n = 0) and in velocity (u n σ n = 0) are incompatible and cannot be enforced at the same discrete instant.

At this stage, we can distinguish two cases related to the nature of the dynamic contact system considered. In the case where deformable bodies are concerned, the contact conditions formulated in displacement (equation ( 1)) are usually preferred, with the use of an additional numerical treatment to take into account the persistency condition and thus the energy conservation properties [START_REF] Laursen | Computational contact and impact mechanics[END_REF][START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF][START_REF] Armero | Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems[END_REF][START_REF] Simo | The discrete energy-momentum method. Part I: Conserving algorithms for nonlinear elastodynamics[END_REF][START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in non linear elasticity[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF][START_REF] Ayyad | Formulation and analysis of two energy-consistent methods for nonlinear elastodynamic frictional contact problems[END_REF][START_REF] Acary | Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact[END_REF]. In the case of multi-rigid-body systems, the question of nondissipation of energy is essential to ensure reliable dynamics of the system. Then the contact conditions formulated in velocity (equation ( 5)) are always used and relate the impulse forces to the velocities. Here, the issue is to predict the velocities of the bodies and the impulse forces acting on the simultaneous multicontacts. To resume, the choice of the contact model and, more precisely, the choice of the complementarity relation depend on the nature of the physics of the problem (a system of deformable bodies or a multi-rigidbody system) involved in the context of nonsmooth contact dynamics. In the following, we recall the main traits of the primal-dual active set type method in the case of the unilateral contact conditions formulated in displacement to solve a system of deformable bodies. The description of the primal-dual active set method for the solution of multi-rigid-body system will be dealt with in Section 4.3.

Primal-dual active set type method

As mentioned by Hintermuller and colleagues [START_REF] Hintermuller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF][START_REF] Hintermuller | The primal-dual active set strategy as a semismooth Newton method[END_REF], the primal-dual active set type method can be seen as a semismooth Newtonian method, for which the contact conditions can be reformulated in terms of a fixed point problem [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and timemarching schemes[END_REF]. Let us recall this usual result in the case of the local unilateral contact conditions formulated in displacement.

Proposition 1. Let γ > 0, the contact conditions in displacement (equation ( 1)) are equivalent to:

σ n = [σ n -γ u n ] + on Ŵ c . ( 6
)
The proof of Proposition 1 can be found in Abide et al. [START_REF] Abide | Analysis of two active set type methods to solve unilateral contact problems[END_REF] or Chouly et al. [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and timemarching schemes[END_REF]. Here, Ŵ c denotes the boundary part of in contact.

The main trait of the active set type method is to consider separately the solution of the equation of motion (or the equation of equilibrium) R(u, σ ) = 0 and the solution of the fixed point σ n = [σ nγ u n ] + . Here, R(u, σ ) is the generalized nonlinear operator, which includes all the operators related to the acceleration, the constitutive behaviour of the material and the external solicitations. Let us denote by S the set of contact nodes of Ŵ c and by c a node of S. The solution of equation ( 6) leads us directly to enforce u c n = 0 if the nodes c ∈ S are in a status of contact by checking an active set condition, that is σ c nγ u c n ≥ 0 for all c ∈ S. Furthermore, the calculus of the fixed point of equation ( 6) leads also to the condition σ c n = 0 in the case of noncontact [START_REF] Hintermuller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF][START_REF] Hintermuller | The primal-dual active set strategy as a semismooth Newton method[END_REF]. Thereby, let us consider the active subset A of contact nodes defined by A = {c ∈ S : σ c nγ u c n ≥ 0}, and the associated inactive subset: I = S \ A. The general form of the iterative active set algorithm of index i can be as follows:

(i) Choose (u 0 , σ 0 ), set i = 0. (ii) Set:

A i+1 = {c ∈ S : σ c,i n -γ u c,i n ≥ 0}, (7)
I i+1 = S \ A i+1 . ( 8
) (iii) Find (u i+1 , σ i+1 ), such that R(u i+1 , σ i+1 ) = 0 u c,i+1 n = 0 for all c ∈ A i+1 , (9)
σ c,i+1 n = 0 for all c ∈ I i+1 . ( 10
) (iv) If A i+1 = A i stop, else go to (ii).
The goal of the active set type strategy is to find the correct subset A of all nodes that are currently in contact by considering the status of contact that is derived directly through computation of the fixed point of equation [START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF]. The main trait of this method is to consider exactly the contact conditions at the contact nodes without the need for using additional nodes to determine the normal contact stress, as is the case for the augmented Lagrangian approach. In the case of a nonlinear constitutive behaviour of the material, a standard Newtonian method can be coupled to the active set method to solve the equation R(u i+1 , σ i+1 ) = 0. Therefore, at each Newton iteration, stopping criteria related to the pair (u, σ ) and the residue R(u, σ ) must be added.

To solve a dynamic multi-rigid-body system with contact, the unilateral contact conditions must be expressed in terms of velocity. Then Proposition 1 can be naturally generalized in the case of a local unilateral contact law expressed in velocity. Proposition 2. Let γ > 0; the contact conditions in velocity (equation ( 5)) are equivalent to:

r n = [r n -γ un ] + ( 11
)
where r n is the reaction impulse force between two particles in contact.

Proof. We use here similar arguments already employed for the proof of Proposition 1 [START_REF] Abide | Analysis of two active set type methods to solve unilateral contact problems[END_REF].

First, let us assume that the following conditions hold: un ≥ 0, r n ≥ 0 and un r n = 0. We consider successively the cases un > 0 and un = 0. If un > 0, the condition un r n = 0 implies that r n = 0. Thus, [r n -

γ un ] + = [-γ un ] + = 0 = r n , since γ > 0. If un = 0 and r n ≥ 0 then [r n -γ un ] + = r n .
Conversely, we assume now that equation [START_REF] Fortin | An improved discrete element method based on a variational formulation of the contact law[END_REF] holds. This implies that r n ≥ 0. Next if r n = 0, we have [-γ un ] + = 0 and this leads to un ≥ 0, since γ > 0. Finally, if r n > 0 then r n = r nγ un > 0; so γ un = 0 , which implies un = 0, since γ > 0.

The goal of the next section is to propose a general algorithm for NSCD, for which the dynamic contact condition (equation ( 11)) is treated by an active set type method.

Nonsmooth contact dynamics (NSCD)

After presenting the equations of motion of a dynamic multi-rigid-body contact problem, the general algorithm for solving the full problem is described. This general algorithm is based on two main computational level tasks: the global level, with the solution of the equations of motion, and the local level, devoted to the contact resolution. At this stage, we briefly recall several methods dedicated to the numerical contact resolution. These methods are based on the bipotential theory and the augmented Lagrangian theory.

Equations of motion

Classically [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Moreau | Some numerical methods in multibody dynamics: application to granular materials[END_REF][START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections[END_REF][START_REF] Jean | The non smooth contact dynamics method[END_REF], the motion of a multi-contact system is described using a global generalized coordinate q describing the centre position and the rotation of each particle (for N p particles, q ∈ R d×N p , where d = 6 for a three-dimensional problem and d = 3 for a two-dimensional problem). Owing to the possible shocks between particles, it is necessary to introduce the generalized velocity denoted by q as a function of bounded variations and its associated differential measure d q. Then the equations of motion must be formulated in terms of a differential measure equation:

Md q + F int (t, q, q)dt = F ext (t, q, q)dt + dR, ( 12
)
where • M represents the generalized mass matrix; • F int and F ext represent the internal and external forces, respectively; • dR is a nonnegative real measure, representing the reaction forces and impulses between particles in contact.

For the sake of simplicity and without lost of generality, only the external forces are considered in the following. The internal forces are neglected because the general case can be easily derived through a linearizing procedure.

Then, for the numerics in a time-stepping approach, we consider that the time interval of interest [0, T] is discretized by introducing uniform time instants t k defined by t k+1 = t k + t for k = 0, . . . , N T -1, where t = T/N T is the time step and N T is the number of time steps. Then equation ( 12) is integrated over each time interval [t k , t k+1 ], and approximated using a θ-method with θ ∈] 1 2 , 1] for stability reasons [START_REF] Moreau | Numerical aspect of sweeping process[END_REF][START_REF] Renouf | Conjugate gradient type algorithms for frictional multicontact problems: applications to granular materials[END_REF]. Therefore, the classical approximation of equation [START_REF] Fortin | Numerical simulation of granular materials by an improved discrete element method[END_REF] yields

M( qk+1 -qk ) = t(θF k+1 + (1 -θ)F k ) + P k+1 q k+1 = q k + tθ qk+1 + t(1 -θ) qk (13)
where P k+1 represents the value of the total impulsion over the time step, which contains the contribution of smooth contact and the local percussion densities exerted during shocks, and F k (or F k+1 ) is the external force computed at time t k (or t k+1 ). We will denote qfree k = qk + M -1 t(θF k+1 + (1θ)F k ) as the free velocity (the velocity when the contact impulses vanish). Then, the first equation in equation ( 13) becomes

qk+1 = qfree k + M -1 P k+1 . ( 14
)
To write the contact law, for a contact node x c (note that c, 1 ≤ c ≤ N c , is used as a label for contact nodes, where N c is the total number of contact nodes), it is necessary to define the local-global mapping

v c = H * (q, c) q P = H(q, c)p c ( 15
)
where v c is the local relative velocity between the two bodies in contact and p c is the contact impulse (v c , p c ∈ R d where d is the dimension of the problem) and H(q, c) is the matrix of the local-global mapping that permit to compute the variables v c and p c in the local frame at a contact node x c from the global variable q and P. For duality reasons, H * is the transpose of matrix H. We recall that p c can be decomposed into the sum of a normal component p c n and a tangential component p c t , as p c = p c n n + p c t . We also denote by H(q) the total-global mapping, for v and p in R d×N c (vectors composed of all relative velocities and contact impulses, respectively):

v = H * (q) q P = H(q)p (16)
Let us notice that even if H(q, c) and H * (q, c) have good theoretical properties, this is not necessarily the case for H and H * . In the discretization, a prediction of q is computed to estimate the mapping H(q) (see equations [START_REF] Laursen | Computational contact and impact mechanics[END_REF] and ( 19)).

Using equations ( 13) and ( 16), the discretization of the motion of a multi-contact system, with contact between particles can be written:

ṽk+1 = ṽfree k + Wp k+1 contact_law(ṽ c k+1 , p c k+1 ) = .true ∀c ∈ {1, 2, . . . , N c } (17)
where W = H * M -1 H is the Delassus operator, and ṽfree k = P * qfree k is the relative free velocity. Notice that a Newtonian impact law is also considered in the first part of equation ( 17) (see equation ( 21)) [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], which modifies v k and v free k as ṽk and ṽfree k , respectively. The second part of equation ( 17) is the implicit frictionless contact law that is in our case the classical Signorini condition with p c t = 0.

Resolution of the global problem: the nonlinear Gauss-Seidel method (NLGS)

In this section, we describe the algorithm used at the global level to solve the problem of equation [START_REF] Wriggers | Computational contact mechanics[END_REF]. Following the ideas of Jean and Moreau [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Jean | The non smooth contact dynamics method[END_REF][START_REF] Jourdan | A Gauss-Seidel like algorithm to solve frictional contact problems[END_REF], we use the NLGS algorithm, which is the most commonly used. It involves considering each contact successively until convergence is achieved. The numerical criterion used to state the convergence will be studied later in this paper. This method is intrinsically sequential but it is possible to use a simple multi-threading technique, which consists of splitting the contact loop onto several threads. This method has been studied by Reouf et al. [START_REF] Renouf | A parallel version of the non smooth contact dynamics algorithm applied to the simulation of granular media[END_REF] for the case of a local algorithm based on the augmented Lagrangian method. Notice that it is also possible at this stage to consider more sophisticated methods, such as conjugate gradient-type methods [START_REF] Renouf | Conjugate gradient type algorithms for frictional multicontact problems: applications to granular materials[END_REF]. However, these techniques do not significantly speed up the convergence. This is why they are not considered in the rest of this paper. The time-stepping method combined with the NLGS algorithm takes the following form:

•

Loop on the step time k:

-Prediction of a position (for computation of the local-global mapping):

q k+ 1 2 = q k + t 2 qk ; (18)
v i = H * (q k+ 1 2
, c) qk ; velocity at the beginning of the time step ( 19)

v c,j,f = H * (q k+ 1 2 , c) qj k+1 ;
predicted velocity at the end of the time step (20)

• Newton shock law (using formal Moreau velocity)

ṽc,j+1 n = v c,j,f n + e n v i n 1 + e n ; (21)
• Computation of the contact law contact_law(ṽ c,j+1 n , p c,j+1 ) = .true;

• Actualization of the generalized displacement

qj+1 k+1 = qj k+1 + M -1 H(q k+ 1 2 , c)p c,j+1
. * End of the loop on contacts c. -End of the loop on j of NLGS. When the convergence is reached, actualization of the velocity: qk+1 = qj+1 k+1 -Actualization of the generalized displacements: q k+1 = q k+ 1 2

+ t 2 qk+1 •
End of the loop on the step time k.

Numerical computation of the local contact impulses

We present in this section three methods to compute the local contact forces in the previous algorithm (equation [START_REF] Hintermuller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF]). The first two are based on the work of Alart and Curnier [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF], de Saxcé and Feng [START_REF] De Saxcé | New inequality and functional for contact with friction: the implicit standard material approach[END_REF] and Fortin et al. [START_REF] Fortin | Numerical simulation of granular materials by an improved discrete element method[END_REF] and are classical, and will be useful to evaluate the third method, which is based on the primal-dual active set method. Moreover, the first two methods have been developed in the more general problems of frictional contact problems, which are no longer convex. Let us notice that for all of the proposed methods, an assumption on the existence and the uniqueness of the one-contact problem is needed. This property is verified here by the positive sign of the reduced element of the Delassus matrix.

A Newtonian or augmented Lagrangian technique to solve the contact problem

The first method, developed by Alart and Curnier [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF], is based on the optimization of an augmented Lagrangian problem. In the case of a contact problem with friction, the method can be written as follows. For a given contact point c (1 ≤ c ≤ N c ), one can define the Cartesian product of an infinite half cylinder with section equal to a ball B(0, µp c ) of radius µp c by C(µp c ) = R + × B(0, µp c ), where µ is the friction coefficient. Then the granular type contact problem is given by p c ∈ argmin s∈C(µp c ) J c (s), [START_REF] Hintermuller | The primal-dual active set strategy as a semismooth Newton method[END_REF] where

J (s) = 1 2 s • W cc s +   v c,free -ṽc + N c α=1,α =c W cα p α   • s
and W cα is the reduced element of the Delassus matrix, related to the contacts c and α. Then the projected gradient method can be applied to solve the minimization problem (equation ( 23)). Therefore, for each iteration j of the NLGS algorithm, this approach leads to the following form of the contact reaction as a saddle point of an augmented Lagrangian problem:

p c,j+1 = proj   p c,j -ρ   v c,free -ṽc + N c α=1,α =c W cα p α + W cc p c,j   , C(µp c )   , ( 24
)
where the function proj(v, C) is the orthogonal projection of v on the convex set C. Equation ( 24) can be formulated equivalently as follows:

p j+1 = proj(τ j+1 , C(µp c )) where τ j+1 = p j -ρv j and v j = v c,free -ṽc + N c α=1,α =c W cα p α + W cc p c,j . (25)
Here,τ j+1 is the augmented contact reaction and ρ is the arbitrary positive parameter. This method, which is based on an augmented Lagrangian technique with a projected gradient solver, will be referred to as the standard augmented Lagrangian (SAL) method.

It is also possible to use an iterative Newtonian method to find the minimum of J by seeking the solution as a zero of the function f (χ ) where, for a contact c, χ = (p c , ṽc ) and

fc (χ ) =      ṽc -ṽc,free - N c α=1 W cα p α Zc     
, in which the vector Zc is the error on the prediction of the reaction

Zc (p c , ṽc ) = p c -proj(τ c , C(µp c )), (26)
This method, which is based on an augmented Lagrangian technique and a quasi-Newtonian solver, will be referred to as the enhanced augmented Lagrangian (EAL) method. The idea of the enhancement is to consider, for the approximation of the gradient of fc , not only the part of fc that is always differentiable, as in the standard case, but all the part of fc that is differentiable, especially when τ c is in the interior of C(µp c ). This idea is also used in the next section with the bipotential, where it is explained in more detail [START_REF] Dumont | On enhanced descend algorithms for solving frictional multi-contact problems: applications to the discrete element method[END_REF].

A bipotential technique to solve the contact problem

It is also possible to use the notion of a bipotential, introduced by Fortin, de Sacxé and colleagues [START_REF] Fortin | Simulation numérique de la dynamique des systèmes multicorps appliquée aux milieux granulaires[END_REF][START_REF] Fortin | An improved discrete element method based on a variational formulation of the contact law[END_REF][START_REF] Fortin | Numerical simulation of granular materials by an improved discrete element method[END_REF][START_REF] De Saxcé | New inequality and functional for contact with friction: the implicit standard material approach[END_REF]. Using this framework, a couple (v, p) satisfies the Signorini's contact conditions if and only if

b c (w, s) + w • s ≥ b c (v, p) + v • p = 0 ∀w, s (27)
where b c is the bipotential

b c (-v, p) = R + (v n ) + K µ (p) + µp n v t ( 28
)
C stands for the indicator function of the set

C: C (x) = 0 if x ∈ C, C (x) = +∞ if x / ∈ C, and K µ is the set defined by K µ = {p = p n n + p t : p t -µp n ≤ 0}. ( 29
)
The minimization of equation ( 27) is classically realized using a Uzawa method without considering the singular term R + (ṽ c,i n ). This minimization can also be viewed as the proximal point of the augmented force pρ ṽ, with respect to the function p → ρb c (-ṽ, p) [START_REF] Fortin | Simulation numérique de la dynamique des systèmes multicorps appliquée aux milieux granulaires[END_REF][START_REF] Fortin | Numerical simulation of granular materials by an improved discrete element method[END_REF][START_REF] De Saxcé | New inequality and functional for contact with friction: the implicit standard material approach[END_REF]:

p = prox(p -ρ ṽ, ρb c (-ṽ, p))
The bipotential method with an Uzawa solver will be referred to as the standard bipotential (SBP) method.

Then, following the ideas developed by Joli and Feng [START_REF] Joli | Uzawa and Newton algorithms to solve frictional contact problems within the bi-potential framework[END_REF] for the case of continuum mechanics and adapted by Dumont [START_REF] Dumont | On enhanced descend algorithms for solving frictional multi-contact problems: applications to the discrete element method[END_REF] for the case of granular materials, to solve the local problem with a Newton-like method, it is necessary to reformulate the problem. Therefore, for each contact c, the local problem to be solved can be written as

     ṽc k+1 = ṽc,free k + N c α=1 W cα p α p c = proj(τ c , K µ ) ∀c = 1, . . . , N c (30)
where τ c = p cρ ṽ is the augmented reaction used in equation ( 24). This problem can be written equivalently

as      ṽc k+1 -ṽc,free k - N c α=1 W cα p α = 0 p c -proj(τ c , K µ ) = 0 ∀c = 1, . . . , N c (31)
Remembering now that we want to use a Newton algorithm to solve theses equations inside the NLGS loop on the variable j, we define now, for each contact c = 1, . . . , N c , the function

f j c (χ ) =      ṽc,j -ṽc,free k - N c α=1 W cα p α,j Z c,j     
where:

• the vector Z c is the error in the prediction of the reaction

Z c,j (p c,j , ṽc,j ) = p c,j -proj(τ c,j , K µ ). (32)
•

χ c = (p c,j , ṽc,j ) t . • χ = (χ 1 , χ 2 , . . . , χ N c ) t .
The first equality in the relation f (χ ) = 0 is the equation of motion for the bodies in contact, and the second relation is related to the friction condition between the bodies in contact, written within the bipotential framework [START_REF] Dumont | On enhanced descend algorithms for solving frictional multi-contact problems: applications to the discrete element method[END_REF][START_REF] Feng | The bi-potential method applied to the modeling of dynamics problems with friction[END_REF]. Then a Newtonian algorithm of index l is used to solve the problem f (χ ) = 0. This algorithm can be written, for a contact c, as follows:

• Initialization: χ 0 c = p 0 = p c,j , ṽ0 = ṽc,j t , ℓ = 0 • Loop on ℓ, until convergence: - τ c ℓ = p ℓ -ρ ṽℓ - Resolution: ∂f c ∂χ c (χ ℓ ) χ c = -f c (χ ℓ ) (33)
-

Actualization: χ ℓ+1 c = χ ℓ c + χ c •
End of the loop on ℓ until convergence, ṽc,j+1 = ṽℓ and p c,j+1 = p ℓ .

For a two-dimensional problem, the matrix

∂f c ∂χ c (χ ℓ ) is equal to ∂f c ∂χ c (χ ) = -W Id 2×2 A c B c ( 34
)
where

A c = ∂Z c ∂p n ; ∂Z c ∂r t ( 35
)
B c = ∂Z c ∂v n ; ∂Z c ∂v t ( 36
)
In the case with no friction (µ = 0), the matrices A c and B c take different forms according to the contact status:

• First case: contact A c = [0; t] B c = [ρn; 0]
• Second case: no contact. In this case

A c = Id 2×2 B c = 0 2×2
The previous bipotential technique with the use of a Newtonian solver will be referred to as the enhanced bipotential (EBP) method.

A primal-dual active set method for NSCD

This section is devoted to the numerical treatment of the contact conditions by a primal-dual active set method within the framework of multi-rigid-body dynamic contact problems.

For this purpose, it is necessary to define the active set defined in equations ( 7) and ( 8), and to compute the contact conditions on each set only in terms of contact impulses, using the local general equations of motion (equation ( 17)) in the form of equation [START_REF] Fortin | An improved discrete element method based on a variational formulation of the contact law[END_REF]. We recall that S denotes the set of potential contact particles and c ∈ S represents a potential contact between two particles. Moreover, let us consider the active subset A of contact defined by A = {c ∈ S : p c nγ ṽc n ≥ 0}, and the associated inactive subset: I = S \ A. With these notations, the numerical computation of the local contact step inside the NLGS iteration loop of index j leads to the following primal-dual active set algorithm:

(i) Compute: τ c,j+1 n = p c,j n -γ ṽc,j n for each c ∈ S. (ii) Set: A j+1 = {c ∈ S : τ c,j+1 n ≥ 0}, (37)
I j+1 = S \ A j+1 . ( 38
(iv) Convergence obtained for the active sets when : A j+1 = A j and I j+1 = I j .

We can remark that the term τ c,j+1 n computed in the first step is similar to the augmented expression (equation ( 25)) used for the methods presented in the two previous sections. The main trait of the active set type method is to enforce directly and exactly the contact conditions that are found when the problems considered in equations ( 24) and ( 31) are solved with a Newtonian type solver. Unlike the augmented Lagrangian approach, another trait of the active set type method is the nonuse of additional nodes for determination of the normal contact stress. Indeed, the computation of p c,j+1 n for all c ∈ A j+1 is obtained a posteriori, see equation [START_REF] Renouf | A parallel version of the non smooth contact dynamics algorithm applied to the simulation of granular media[END_REF]. The method described in this section, which is based on a primal-dual active set approach, will be referred to as the primaldual active set (PDAS) method.

Let us notice that only one iteration of the active set method is realized for one contact at each global iteration of the NLGS method, as in the case of the SBP method, and unlike the other methods. A numerical study showed that considering several iterations of the active set method at each NLGS step did not improve the method.

Stopping criterion

Since each numerical solution of the presented method ends the NLGS iteration with different properties (some of the methods are exact for the dynamic, some of the methods are exact for the contact law, and the others are not exact for either), it is necessary to use an identical stopping criterion for all the methods, which is able to take into account all types of error. Even if such very robust stopping criteria exist in the literature [START_REF] Facchinei | Finite-dimensional variational inequalities and complementarity problems[END_REF], we propose here a criterion that permits each type of error to be exhibited separately. This criterion, developed from that proposed by Fortin et al. [START_REF] Fortin | An improved discrete element method based on a variational formulation of the contact law[END_REF] has been extended in the case of the Newtonian and bipotential (EBP) method, where some terms are naturally vanishing in the original Uzawa and bipotential (SBP) method. This stopping criterion ε glob has been proposed in such a way that both the error in the equation of motion ε c motion , the error in the Signorini contact law ε b c and the error of penetration ε c pen must be sufficiently small. Therefore, the numerical solution provides good properties related to both the equation of motion and the Signorini contact law.

This criterion can be stated as

ε glob = 1 N c N c c=1 ε c motion + ε b c + ε c pen ( 42
)
where:

• ε c motion = ṽc -ṽc m , where ṽc m = ṽc,i + N c α=1 W cα p α , so ε motion measures the error on the equation of motion (see equation [START_REF] Koziara | Semismooth Newton method for frictional contact between pseudo-rigid bodies[END_REF]; this term vanishes for the SBP and SAL methods);

• ε b c = |ṽ c • p c |
is the absolute value of the bipotential, which must vanish if and only if the couple (ṽ c , p c ) verifies the Signorini contact law; • ε c pen = -min(0, ṽc n ) is the value of the penetration. One can notice that it is absolutely necessary to verify in the criterion that there is no penetration, because nothing in the presented algorithm ensures that this condition is satisfied at the end of the loop. Moreover, if this condition is not satisfied, the other part of the bipotential can be negative or vanish, even if the couple (ṽ, p) is not a solution.

Numerical simulations

In this section, several numerical examples are solved to illustrate the performances of the PDAS method, compared with various other methods (SBP, SAL, EBP and EAL). Three numerical examples of increasing complexity are considered: the sliding of one ball on a plane obstacle, the sedimentation of a collection of balls and a poly-disperse collection of balls falling over an inclined plane.

Sliding of one ball on plane obstacle

In this first example, a ball posed on a horizontal plane is considered (see Figure 1) with a nonvanishing initial horizontal velocity (v 0 = 1.5 m• s -1 ). There is only one point contact between the ball and the plane. Since there is no friction, the ball slides on the horizontal plane, without rolling. The final time is equal to T = 0.2 s and the number of time steps is fixed to 2000 with t = 10 -4 s. The results for this example with sliding are presented in Table 1. In this example, one can notice (see Table 1) that all the methods need a similar total CPU time to compute the solution. However, the time necessary to compute the contact forces is smaller for the three last methods. Moreover, the PDAS method needs only one iteration of NLGS, meaning that good status for the contact is found directly. For the EBP and EAL methods, two iterations of NLGS are necessary, for the standard methods, since there is no convergence for the contact law within each NLGS iteration, 18 iterations are necessary to obtain convergence.

Let us notice that this value also strongly depends on the penalization parameter, unlike the other methods, which are relatively independent of the parameters [START_REF] Dumont | On enhanced descend algorithms for solving frictional multi-contact problems: applications to the discrete element method[END_REF].

Sedimentation of a collection of balls

In this second example, the sedimentation of a collection of 500 rigid balls, with radii ranging from 0.25 mm to 0.5 mm, in a box is considered (see Figure 2). The stopping criterion is equal to ε glob = 10 -7 . The time step is equal to t = 5.10 -5 s and the number of time steps is equal to 1000.

This example is representative of multi-rigid-body contact problems, owing to the large number of rigid bodies considered and the large proportion of computing time necessary for the local treatment of the contact compared with the total CPU time. According to Table 2, one can observe in this example that the active set method provides the best results in terms of computing time, both for the total time and the time devoted to the treatment of the contact. The EBP and the EAL methods need similar computing time and the SBP and SAL methods also provided comparable results but with a larger computing time, essentially resulting from the computation of the contact forces.

Example of a poly-disperse collection of balls falling over an inclined plane

The third example concerns the falling of a poly-dispersed collection of 500 balls over an inclined plane, with a slope equal to 10 • . The radii of the balls range from 0.6 mm to 1 mm (see Figure 3). The balls in contact with the plane are sticking. The stopping criterion is equal to ε glob = 10 -7 . The time step is equal to t = 10 -5 s and the number of time steps is equal to 10000. The interest of this example lies in the fact that the balls can have very different velocities, which makes simulation using the SBP and SAL methods difficult, since the convergence and the optimal parameters depend strongly on the cinematics. Indeed, the slope of the plane has been chosen in such a way that the balls close to the bottom have small velocity, by contrast with the balls on the top.

As expected, one can observe in Table 3 that the two first methods are not efficient, since they are 75% more computationally expensive than the others. The EBP and EAL methods provide very similar results in terms of time consumption. Once again, the PDAS method is slightly faster than the other methods, with a gain of 17.3% compared with the EAL method.

To study the convergence of the NLGS method at one time step, Figure 4 shows the convergence of the NLGS method for three different treatments of the contact, namely the SBP, the EAL and the PDAS methods.

One can observe in Figure 4 that convergence of the NLGS algorithm associated with the SBP method is very slow, unlike that associated with the two other methods, which provide very close results. Note that the convergence is faster with the PDAS method.

Conclusion

In this article, we have proposed an adaptation of the active set method to solve the contact dynamics problem with several rigid particles. Numerical tests show that this method is efficient and fast compared with standard methods. These first encouraging results leads us to envisage taking into account friction, using Coulomb's law within the framework of the active set method. We refer to Hintermuller et al. [START_REF] Hintermuller | Constrained optimization for interface cracks in composite materials subject to non-penetration conditions[END_REF][START_REF] Hintermuller | Obstacle problems with cohesion: a hemi-variational inequality approach and its efficient numerical solution[END_REF] and Kunisch and Stadler [START_REF] Kunisch | Generalized Newton methods for the 2D-Signorini contact problem with friction in function space[END_REF] for the treatment of frictional and cohesive contact.

-

  Initialization of the motion: q0 k+1 = qfree k (initialization of the contact forces with P = 0); -Loop on j ≥ 0 (NLGS), until convergence: * Loop on the contacts c: • Computation of the local-global mapping

  ) (iii) Find (v j+1 , p j+1 ) such that if c ∈ A j+1 then ṽc,j+1 n = 0 and p c,j+1 t = 0, (contact) (39) Computation of the local contact force: p c,j+1 n = p c,j n -ṽc,j n /W nn . (40) if c ∈ I j+1 then p c,j+1 n = 0, and p c,j+1 t = 0 (no contact).

Figure 1 .

 1 Figure 1. Sliding of one ball on plane obstacle.

Figure 2 .

 2 Figure 2. Sedimentation of a collection of balls (left: initial state; right: final state).

Figure 3 .

 3 Figure 3. Balls falling over an inclined plane (initial state).

Figure 4 .

 4 Figure 4. Convergence of the NLGS method for various treatments of the contact law at the 1000th time step.

Table 1 .

 1 Results for example 1. For each method are provided the number of nonlinear Gauss-Seidel iterations for each time step (second column), the total CPU time (third column) and the CPU time devoted to the computation of the contact forces (fourth column).

		Iterations	CPU time (s)	Contact forces CPU time (s)
	Standard bipotential	18	4.10	1.29×10 -1
	Standard augmented Lagrangian	17	4.07	1.11×10 -1
	Enhanced bipotential	2	4.00	8.40×10 -2
	Enhanced augmented Lagrangian	2	4.02	9.5×10 -2
	Primal-dual active set	1	3.91	7.7×10 -2

Table 2 .

 2 Results for example 2. For each method is provided the total CPU time (second column) and the total CPU time necessary to compute the contact forces (third column).

	CPU time (s)	Contact forces CPU time (s)

Table 3 .

 3 Results for example 3. For each method is provided the total CPU time (second column) and the total CPU time necessary to compute the contact forces (third column).

	Standard bipotential	6790	6724
	Standard augmented Lagrangian	6704	6660
	Enhanced bipotential	3967	3956
	Enhanced augmented Lagrangian	3699	3676
	Primal-dual active set	3058	2977

CPU time (s)

Contact forces CPU time (s)

Funding

This research received no specific grant from any funding agency in the public, commercial or not-for-profit sectors.