Vincent Acary
email: vincent.acary@inria.fr

Keywords: Computational contact mechanics, flexible multibody dynamics, unilateral contact, Signorini's condition, impact, Coulomb's friction, geometric numerical integration, Gear-Gupta-Leimkuhler(GGL) technique, event-capturing schemes Computational contact mechanics, flexible multibody dynamics, unilateral contact, Signorini's condition, impact, Coulomb's friction, geometric numerical integration, Gear-Gupta-Leimkuhler(GGL) technique, event-capturing schemes

This work addresses the problem of the numerical time-integration of nonsmooth mechanical systems subjected to unilateral contacts, impacts and Coulomb's friction. The considered systems are the space-discretized continuous systems obtained by using a Finite Element Method (FEM) approach or the multi-body systems, or a mix of them as in flexible multibody dynamics. Up to now, two main numerical schemes are available for this purpose: the Moreau-Jean scheme which solves the constraints at the velocity level together with a Newton impact law and the Schatzman-Paoli scheme which directly considers the constraints at the position level. In both schemes, the position and velocity constraints are not both satisfied in discrete time. A first attempt is made by directly using the Gear-Gupta-Leimkuhler(GGL) approach for Differential Algebraic Equations (DAE), that solves, in discrete time, the constraints on both position and velocity levels. The obtained direct projection scheme succeeds in solving in discrete time both position and velocity constraints, but introduces some chattering at contact after a finite accumulation of impacts. The proposed scheme improves the direct projected scheme by combining several steps of activation and projection to avoid the chattering effect. The stability and the local order of the scheme will be discussed. The usefulness of the scheme is demonstrated on several academic examples and is illustrated on an industrial application : the modeling and simulation of an electrical circuit breaker.

Schémas à capture d'événements pour les systèmes mécaniques non réguliers en présence de contact et de frottement.

Résumé : Dans ce travail, on propose une nouvelle stratégie pour l'intégration numérique en temps des systèmes mécaniques non réguliers. Les systèmes considérés sont les systèmes mécaniques multicorps flexibles ou les systèmes mécaniques continus discretisés par éléments finis. Jusqu'à maintenant, deux schémas principaux existaient : le schéma de Moreau-Jean qui résoud les contraintes au niveau des vitesses avec une loi d'impact de type Newton et le schéma de Schatzman-Paoli qui considére directement les contraintes en position. Dans les deux schémas, les contraintes en vitesse et en position ne sont pas satisfaites simultanément. Une première solution est proposée utilisant directement la technique de Gear-Gupta-Leimkuhler(GGL) pour les équations différentielles algébriques, qui permet de résoudre en temps discret les contraintes à la fois en position et en vitesse. Le schéma proposé résoud le problème de respect des contraintes mais introduit des oscillations artificielles aux contacts stabilisés après une accumulation finie d'impacts.

Le second schéma proposé améliore le schéma direct en combinant des étapes successives de projection et d'activation des contraintes. Ce schéma élimine les oscillations au contact. La stabilité et l'ordre de la méthode sont discutés. L'utilité de ces schémas est démontrée sur différents exemples académiques ainsi que sur une application industrielle : le prototypage virtuel des disjoncteurs électriques.

Mots-clés : Mécanique numérique du contact, dynamique des systèmes mécaniques flexibles, contact unilatéral, condition de Signorini, forttement de Coulomb, intégration numérique géométrique, technique de Gear-Gupta-Leimkuhler(GGL), schéma à capture d'événements.

Introduction and motivations

This work addresses the problem of the numerical time-integration of nonsmooth mechanical systems subjected to unilateral contacts, impacts and Coulomb's friction. The targeted systems are the multi-body systems where interconnected rigid or flexible bodies interact trough perfect joints and ideal unilateral frictional interfaces. As flexibility may have an important role in the global dynamic behavior, we are also interested in considering discrete systems which result of a space-discretization of a solid, for instance, by finite element techniques. The simulation, and especially time-integration of nonsmooth mechanical systems with unilateral contacts is an active research domain due to the complexity of performing an efficient, accurate and robust simulation. The main issue is the inherent nonsmoothness of the time evolution as a result of the nonsmoothness of the models based on the unilateral contact (Signorini's condition) and Coulomb's friction. It is well-known that the presence of unilateral contact may imply the occurrence of impacts (velocity jumps and/or reaction impulses) and Coulomb's friction may also generate velocity jumps as in the well-known Painlevé example [START_REF] Génot | New results on Painlevé paradoxes[END_REF]. This demands for specific time-integration techniques which are usually classified into two categories: the event-tracking time-stepping schemes (also commonly and shortly called event-driven schemes) and the eventcapturing time-stepping schemes (shortly time-stepping schemes). The first family of schemes is based on an accurate detection of events (closing and opening contacts, changes in the direction of sliding, transition from sliding to sticking or vice-versa, . . . ). Such schemes are mainly dedicated to systems with a small number of events and mainly in the bi-dimensional configuration. For more details on such schemes, we refer to [START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF] and [4,Chapter 8].

When a large number of events are expected in three-dimensional configuration, only eventcapturing time-stepping schemes are sufficiently efficient and robust. This is the case for structural dynamics, and multibody dynamics where the density of events with respect to time prevents the use of an accurate detection of the instants of events. Two main numerical schemes available to integrate the nonsmooth dynamics with impacts: the Moreau-Jean scheme [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid bodies collections[END_REF][START_REF] Moreau | Numerical aspects of the sweeping process[END_REF][START_REF] Jean | The non smooth contact dynamics method[END_REF] which solves the constraints at the velocity level together with a Newton impact law and the Schatzman-Paoli scheme [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF][START_REF] Paoli | A numerical scheme for impact problems I: The one-dimensional case[END_REF][START_REF] Paoli | A numerical scheme for impact problems II: The multidimensional case[END_REF] which directly considers the constraints at the position level. For these schemes, rigorous mathematical analysis have been carried out [START_REF] Monteiro Marques | Differential Inclusions in Nonsmooth Mechanical Problems. Shocks and Dry Friction[END_REF][START_REF] Stewart | Convergence of a time-stepping scheme for rigid-body dynamics and resolution of Painlevé's problem[END_REF][START_REF] Dzonou | A convergence result for a vibro-impact problem with a general inertia operator[END_REF] and numerous large scale applications have proven their own interests [START_REF] Radjaï | Discrete-element modeling of granular materials[END_REF]. These schemes have also numerous variants that have been presented in the literature (see [4] for details), but in any of these schemes, the position and velocity constraints are both satisfied in discrete time.

In computational contact mechanics of solids and structures, the Newmark family of schemes (HHT, α-scheme, . . . ) are generally used for the time integration of space-discretized structures with Signorini's condition and Coulomb's friction [START_REF] Laursen | Computational Contact and Impact Mechanics -Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis[END_REF][START_REF] Wriggers | Nonsmooth mechanical systems with unilateral contact 6[END_REF][START_REF] Deuflhard | A contact-stabilized Newmark method for dynamical contact problems[END_REF][START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF]. In these proposed approaches, it is implicitly assumed that the solutions (position/displacements, velocities, contact forces) are sufficiently smooth such that the Newmark family of schemes of order 2 can be applied without problems. In our case, the contact activation between finite-freedom mechanical models induces the nonsmoothness of the solutions. Therefore, the direct application of higher order schemes in this context may be hazardous. Attempts have recently been made to improve the global order of accuracy of time-stepping with nonsmooth events [START_REF] Studer | Step size adjustment and extrapolation for time stepping schemes in non-smooth dynamics[END_REF]3,[START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous galerkin methods: Definition and outlook[END_REF].

The motivations to build a scheme satisfying constraints both at the position level and at the velocity level in discrete time are :

• The study of multi-body systems with clearances in joints. If the joints with clearances are modelled with unilateral contact, we need to keep the drift of the constraints as smallest as possible with respect to the charateristic lengths of the clearances.

• For multibody systems with perfect ideal joints (bilateral constraints), we want to able to solve the well-known drift issue of the constraints if they are treated at the acceleration level or at the velocity level.

• In quasi-static applications, mostly when finite element method is involved, we want to avoid penetration between bodies, so we want to enforce the constraint at the position level, but a smooth evolution of the local relative velocities at contact. Spurious oscillations at contact of the local velocities are an extensively studied issue in the literature [START_REF] Laursen | Improved implicit integrators for transient impact problemsgeometric admissibility within the conserving framework[END_REF][START_REF] Deuflhard | A contact-stabilized Newmark method for dynamical contact problems[END_REF]. Mimicking the plastic impact law at the velocity level (quasi-collision) allows one to stabilize the velocity and then the stresses at contact.

• Finally, the last motivation is to maintain the consistency of the geometrical model required by the computational geometry system of Computer Aided Design (CAD) tools. For most of collision detection algorithms, avoiding penetrations of bodies is a requirement to guarantee the efficiency and robustness of the results of the contact detection process.

The aim of this work is to propose a new strategy based on the Moreau-Jean time-stepping scheme which enforces in discrete time the constraints on both position and velocity levels. The quest for such a scheme is connected to the approaches of geometric numerical integration theory of differential systems where the discrete approximation of the solution preserves some geometrical properties of the flow [START_REF] Hairer | Geometric numerical integration. Structure-preserving algorithms for ordinary differential equations[END_REF]. In our case, we want that the solution preserves the constraints and the impact law in discrete time.

A first solution is proposed in Section 4 which makes a direct use of the Gear-Gupta--Leimkuhler(GGL) approach. The GGL technique for Differential Algebraic Equations (DAE) [START_REF] Gear | Automatic integration of Euler-Lagrange equations with constraints[END_REF] solves, in discrete time, the constraints on both position and velocity levels. For this purpose, the authors introduced in the continuous time formulation, new kind of multipliers to enforce additional and redundant constraints and reducing by the way the index of the DAE. These multipliers can be understood as the multipliers associated with the projection onto the constraints. The direct application of this idea which was already suggested in [2,[START_REF] Studer | Numerics of Unilateral Contacts and Friction. -Modeling and Numerical Time Integration in Non-Smooth Dynamics[END_REF] results in a scheme able to satisfy the constraints requirements in a very efficient way in most of the configurations. Results on the local order of the scheme are given and qualitative properties are discussed on several academic and industrial examples. Its main drawback lies in the introduction of numerical spurious oscillations when a contact is kept closed after a finite accumulation of impacts. This spurious oscillations, termed as chattering in this paper, is mainly due to an harmful interaction between the unilateral condition and the increase of energy due to the projection onto constraints.

In Section 5, an improved solution is proposed to circumvent the chattering problem. This solution improves the direct projected algorithm by using a special combination of the activation of the constraints at the velocity level and the projection onto these activated constraints. The new combined scheme is mainly based on the direct projected scheme and shares the same favorable properties (respect of the constraints, order, straightforward implementation). Furthermore, it cancels the chattering at contact and avoids increasing energy due to the projection.

The outline of the article is as follows. In Section 2, basic concepts and equations that are used to model nonsmooth multi-body systems are introduced. In Section 3, the basic time-stepping schemes i.e., the Moreau-Jean scheme and the Schatzman-Paoli scheme are briefly reviewed. The direct application of the GGL idea is developed in Section 4. After the formulation of the scheme, some results on the local order of consistency are proven. The main drawback, the socalled chattering effect, is also exhibited on academic examples. A combined activation/projection procedure is investigated in Section 5, which gives the correct answer to our problem. Finally, applications are developed in Section 6 and we show on the slider-crank mechanism and a model of a circuit breaker that the approach is a promising solution in industrial prototying process of mechanisms with clearances. Section 7 concludes the article.

Notation

The following notation is used throughout the paper. The uniform norm for a function f is denoted by f ∞ and for a vector x ∈ IR n by x . A function f is said to be of class C p if it is continuously differentiable up to the order p. Let I denote a real time interval of any sort. The set of functions f : I → IR n of bounded variations (BV) is denoted by BV (I, IR n ). For f ∈ BV (I, IR n ), we denote the right-limit function by f + (t) = lim s→t,s>t f (s), and respectively the left-limit by f -(t) = lim s→t,s<t f (s). The value var(f, I) denotes its total variation on I. We denote by 0 = t 0 < t 1 < . . . < t k < . . . < t N = T a finite partition (or a subdivision) of the time interval [0, T ] (T > 0). For simplicity's sake, the length of a time step is considered to be constant and is denoted by h = t k+1 -t k . The value of a real function x(t) at the time t k , is approximated by x k . In the same way, the notation

x k+θ = (1 -θ)x k + θx k+1 is used for θ ∈ [0, 1].
The notation O(h) is to be understood as h → 0. The notation dt defines the Lebesgue measure on IR. The notation N C (x) is used for the normal cone in the Convex Analysis sense to a convex set C at the point x [START_REF] Rockafellar | Convex Analysis[END_REF]. The function prox M (C, x) returns the closest element of C to x in the metric defined by a definite positive matrix M . For any matrix A ∈ IR n×n and a set of indices α ⊂ {1, . . . , n} ⊂ IN , A αα denotes the submatrix composed of the rows and the columns indexed by the indices in α. The matrix A α• (respectively A •α ) stands for the matrix that collects all the rows (respectively the columns) indexed by α. The matrix I n×n denotes the identity matrix of IR n×n .

Nonsmooth mechanical systems with unilateral contact

In this section, we give the basic ingredients for the modeling of multibody systems with unilateral constraints. For more details on the modeling of multibody systems with unilateral constraints, we refer to [4,[START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF][START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF] and for the mathematical analysis, we refer to [START_REF] Schatzman | A class of nonlinear differential equations of second order in time[END_REF][START_REF] Monteiro Marques | Differential Inclusions in Nonsmooth Mechanical Problems. Shocks and Dry Friction[END_REF][START_REF] Stewart | Rigid body dynamics with friction and impact[END_REF][START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF].

The frictionless case in a pure Lagrangian setting

Let us first consider a pure Lagrangian setting. The equations of motion of multibody systems with unilateral constraints are

                 q(t 0 ) = q 0 , v(t 0 ) = v 0 , (1a) q
(t) = v(t), (1b)
M (q(t)) v(t) + F (t, q(t), v(t)) = G(t, q)λ(t), (1c)
g α (t, q(t)) = 0, α ∈ E, (1d)
g α (t, q(t)) 0, λ α 0, λ α g α (t, q) = 0 α ∈ I, (1e)
where

• q(t) ∈ IR n is the generalized coordinates and v(t) = q(t) the associated generalized velocities,

• the initial conditions are q 0 ∈ IR n and v 0 ∈ IR n ,

• M (q(t)) ∈ IR n×n is the inertia, F (t, q(t), v(t)) ∈ IR n the forces,

• the function g(t, q(t)) ∈ IR m defines the constraints in the dynamical system, and G (t, q(t)) = ∇ q g(t, q(t)) is the Jacobian matrix of g with respect to q,

• λ ∈ IR m is the Lagrange multiplier vector associated with the constraints, and

• the sets E ⊂ IN and I ⊂ IN respectively describe the set of bilateral constraints (joints) and unilateral constraints (contacts).

Remark 1 In the Newton/Euler formalism [START_REF] Goldstein | Classical Mechanics[END_REF][START_REF] Géradin | Flexible Multibody Dynamics: A finite element Approach[END_REF][START_REF] Bauchau | Review of contemporary approaches for constraint enforcement in multibody systems[END_REF], the vector of parameters q usually contains the position of the center of mass x and a parametrization of the finite rotation θ which models the orientation of the body with respect to a spatial frame. The velocity is usually composed the velocity of the center of mass ẋ and an angular velocity Ω expressed for instance in the inertial frame. Therefore, the velocity is not the time-derivative of the parameters q, but generally related to q thanks to an operator T (q) such that

q(t) = T (q(t))v(t). ( 2
)
The equations of motion (1) can be extended to the Newton/Euler formalism by considering (2) rather than (1b) and by defining G as G(t, q(t)) = ∇ q g(t, q(t)) T (q((t)).

(3)

In the remaining of the article, we will consider only the Lagrangian setting to enlighten the notation.

Remark 2 After a space-discretization of continuum solids by a finite element approach, the generalized coordinates q usually contains the nodal displacements, and possibly the nodal rotations if any. Nevertheless, the generalized velocity is most of the time derivative of the coordinates q.

For the sake of simplicity, we also restrict our presentation to holonomic perfect unilateral constraints, that is, we will consider in this paper that E = ∅ and that the constraints are scleronomic constraints, i.e. g(t, q(t)) = g(q(t)). Applications in Section 6 will however show more general cases. The constitutive law for the perfect unilateral constraints is given by the Signorini condition 0 g(q(t)) ⊥ λ(t) 0,

where the inequalities involving vectors are understood to hold component-wise and the x ⊥ y symbol means that y x = 0. Let us define the following variables relative to the constraints, called local variables: the local velocity U (t) and the (local) Lagrange multiplier λ(t) which is associated with the generalized reaction forces r(t) such that

U (t) = G (q) v(t), r(t) = G(q)λ(t). (5)
For finite-freedom mechanical systems, an impact law must be added to close the system of equations. The most simple impact law will be considered in this work given by Newton's impact law

U + (t) = -e U -(t), if g(q(t)) = 0, ( 6
)
where e is the coefficient of restitution. Throughout the paper, several academic test examples are chosen to highlight the properties of the considered numerical integration schemes.

Example 1 (The bouncing ball) This is the standard bouncing ball under gravity depicted in Figure 1(a). The dynamics is constant with a forcing term equal to f together with a unilateral contact on the ground,

v(t) = f (t) + λ(t), q(t) = v(t), 0 q(t) ⊥ λ(t) 0, v + (t) = -ρv -(t), if q(t) = 0, (7)
The interesting feature of the bouncing ball example is the presence of a finite accumulation of impact when 0 < e < 1 and f < 0. The analytical solution of this example can be found in [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control[END_REF].

A more pleasant analytical solution due to Ballard [START_REF] Ballard | Bounded variations and measure theory on the line[END_REF] for f = -2 and ρ = 1/2 is detailed in [2]. It will be used as a benchmark in the further sections.

Example 2 (The linear oscillator) The dynamics of this one-degree-of-freedom system depicted in Figure 1(b) example is similar to the dynamics [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF] but with a linear spring-damper internal force, that is

m v(t) + cv(t) + kq(t) = λ(t). (8)
The explicit analytical solution can be found in [START_REF] Janin | Comparison of several numerical methods for mechanical systems with impacts[END_REF]. The previous trivial free dynamics [START_REF] Ballard | The dynamics of discrete mechanical systems with perfect unilateral constraints[END_REF] with a null or a constant forcing term are exactly integrated with any first order scheme. With the linear, but nontrivial, dynamical term in [START_REF] Ballard | Bounded variations and measure theory on the line[END_REF], the order of accuracy of higher order schemes can be exhibited. at the corner B. The block is parametrized by the coordinates of the center of mass [x, y] and the angle with respect to the ground θ, that is q = [x, y, θ] . The unilateral constraints read as

f A (q) = y -l 2 cos θ + L 2 sin θ 0, for the contact point A, f B (q) = y -l 2 cos θ -L 2 sin θ 0, for the contact point B. ( 9
)
The equations of motion in the frictionless case are

     mẍ = 0 mÿ = -mg + λ A + λ B I θ = λ A [ l 2 sin θ + L 2 cos θ] + λ B [ l 2 sin θ -L 2 cos θ] ( 10
)
where m is the mass of the block and I = m 12 (l 2 + L 2 ) the inertia. Despite the fact that the Newton impact law might not be not the most appropriate law for reproducing the rocking behaviour of the block, we have chosen this example for the strong coupling between the contact points and the nonlinear constraints. Especially, the projection onto the constraints of one of the contact points can lead to a violation of the constraint for the other contact point if it has not been taken into account in a proper way.

Coulomb's friction

Let us consider now Coulomb's friction. In such a case when more complex contact laws are considered, the pure Lagrangian modelling of constraints is not sufficient. Indeed, the use of the Jacobian matrix of the constraints G (t, q(t)) in order to define the normal to the constraints is not necessarily convenient to introduce richer mechanical behaviors at the interface. Hence, we introduce for each contact α a local orthonormal frame at contact point C α composed of a normal vector n α and two tangent vectors t α and s α . In this frame, the local velocity at contact U α and the reaction force λ α are decomposed as

U α = U α N n α + U α T , U α N ∈ IR, U α T ∈ IR 2 , λ α = λ α N n α + λ α T , λ α N ∈ IR, λ α T ∈ IR 2 . ( 11
)
Note that the operator G(q) in (5) that links variables expressed in the local frame to generalized variables is no longer the gradient of some constraints. Coulomb's friction is expressed in a disjunctive form as

   if U T = 0 then λ ∈ C if U T = 0 then ||λ T || = µ|λ N |
and there exists a scalar a 0 such that λ T = -aU T [START_REF] Chen | A Newmark-Type Integrator for Flexible Systems Considering Nonsmooth Unilateral Constraints[END_REF] where C = {λ, ||λ T || µ|λ N | } is the Coulomb friction cone. Let us introduce the modified velocity U [START_REF] Saxcé | New inequality and functional for contact with friction : The implicit standard material approach[END_REF] defined by

U = U + µ ||U T ||n. ( 13
)
With the Signorini condition at the velocity level, this notation provides us with a synthetic form of the Coulomb friction as

-U ∈ N C (λ), (14)
where N C is the normal cone to C [START_REF] Rockafellar | Convex Analysis[END_REF], or equivalently,

C * U ⊥ λ ∈ C. ( 15
)
where

C * = {v ∈ IR n | r T v 0, ∀r ∈ C} is the dual cone of C.
For more details on this formulation and its theoretical interest, we refer to [5].

In this form, the numerical time integration of systems with Coulomb's friction is similar to case with only Signorini's condition written in terms of complementarity at the velocity level. The standard schemes and the new approaches developed in the sequel will directly apply to the case with Coulomb's friction. To improve the readability, only the Signorini condition case will be detailed in the sequel.

Time-integration methods for nonsmooth dynamics

Leaving aside the time-integration methods based on an accurate event detection procedure (event-tracking schemes or event-driven schemes [4, Chap. 8]), two main numerical schemes are available to date for integrating nonsmooth mechanical systems which are sound from the mathematical analysis point of view and which take advantage of a strong practical experience: the Moreau-Jean scheme and the Schatzman-Paoli scheme.

3.1 Moreau-Jean's scheme [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] The Moreau-Jean scheme [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] is based on the Moreau sweeping process which enables to write the unilateral constraints at the velocity level including Newton's impact law,

         M (q(t))dv = F (t, q(t), v + (t))dt + G(q(t))di, q(t) = v + (t), U (t) = G (q(t)) v(t) if g α (q(t)) 0, then 0 U α,+ (t) + eU α,-(t) ⊥ di 0. ( 16
)
where dt is the Lebesgue measure, dv is a differential measure associated with v and di is an impulse reaction measure. When the evolution is smooth, r(t) is considered as the density of di with respect to the Lebesgue measure, that is

r(t) = di dt (t). (17)
The associated local reaction measure is defined by dI = G(q)di.

The numerical time integration of the measure differential inclusion (MDI) ( 16) is performed on an interval (t k , t k+1 ] of length h as follows (θ ∈ [0, 1]):

                 M (q k+θ )(v k+1 -v k ) -hF (t k+θ , q k+θ , v k+θ ) = p k+1 = G(q k+1 )P k+1 , (18a)
q k+1 = q k + hv k+θ , (18b)
U k+1 = G (q k+1 ) v k+1 (18c) if ḡα k+1 0 then 0 U α k+1 + eU α k ⊥ P α k+1 0, otherwise P α k+1 = 0. , α ∈ I (18d)
where the following approximations are considered

v k+1 ≈ v + (t k+1 ); U k+1 ≈ U + (t k+1 ); p k+1 ≈ di(]t k , t k+1 ]), P k+1 ≈ dI(]t k , t k+1 ]). ( 19
)
The value ḡk+1 is a prediction of the constraint that manages the activation at the velocity level. Several formulas for this forecast will be discussed in Section 4.4.

The numerical scheme which solves [START_REF] Dzonou | A convergence result for a vibro-impact problem with a general inertia operator[END_REF] enforces in discrete time the Newton impact law at each time step. On the contrary, the constraints in position g(q(t)) 0 are not strictly satisfied. A violation of the constraints can occur at the activation of the contact and a drift of the constraints is generally observed if the constraints g(q) is non linear.

3.2 Schatzman-Paoli's scheme [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF][START_REF] Paoli | A numerical scheme for impact problems I: The one-dimensional case[END_REF][START_REF] Paoli | A numerical scheme for impact problems II: The multidimensional case[END_REF] The Schatzman-Paoli scheme [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF][START_REF] Paoli | A numerical scheme for impact problems I: The one-dimensional case[END_REF][START_REF] Paoli | A numerical scheme for impact problems II: The multidimensional case[END_REF] deals directly with the unilateral constraint g(q(t)) 0 in discrete time and incorporates the Newton impact law such that the law is satisfied over two or three time-steps. In this scheme for e = 0, the position constraints is satisfied in discrete time but not the impact law.

For a non trivial mass matrix, in the multi-contact case and with a θ-method, the following scheme can be viewed as an extension of the original Schatzman-Paoli scheme

           M (q k+1 )(q k+1 -2q k + q k-1 ) -h 2 F (t k+θ , q k+θ , v k+θ ) = p k+1 , (20a)
v k+1 = q k+1 -q k-1 2h , (20b)
-p k+1 ∈ N K q k+1 + eq k-1 1 + e . ( 20c
)
For an admissible set defined by a finite set of unilateral constraints,

K = {q ∈ IR n , y = g(q) 0}, (21)
the inclusion into the normal can be recast under some constraints qualification conditions as a nonlinear complementarity problem of the form

           g k+1 = g q k+1 + eq k-1
1 + e , p k+1 = G q k+1 + eq k-1

1 + e P k+1 ,

0 g k+1 ⊥ P k+1 0. ( 22
)
The convergence of Schatzman-Paoli's scheme is studied in [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF][START_REF] Paoli | A numerical scheme for impact problems I: The one-dimensional case[END_REF][START_REF] Paoli | A numerical scheme for impact problems II: The multidimensional case[END_REF][START_REF] Paoli | An existence result for non-smooth vibro-impact problems[END_REF] under various assumptions. When the impacts are perfectly inelastic (e = 0), we observe that the constraint g(q k+1 ) 0 is satisfied in discrete time.

Qualitative comparison of the schemes

For the sake of simplicity, let us consider a nonsmooth multi-body system subjected to simple linear constraints q 0. Providing that M is symmetric positive definite in order to define an associated metric and using some basics in Convex Analysis [START_REF] Brogliato | Nonsmooth Mechanics: Models, Dynamics and Control[END_REF], we can write:

M (x -y) -b ∈ -λN K (x), λ > 0 x = argmin z∈K 1 2 (z -y) M (z -y) -(z -y) b x = prox M K; y + M -1 b (23)
Moreau-Jean's time-stepping scheme can be written in terms of the proximal operator as

v k+1 + ev k = prox M (q k+1 ) T IR n + (q k+1 ); (1 + e)v k + hM -1 (q k + 1)F (t k+θ , q k+θ , v k+θ ) (24)
and Schatzman-Paoli's scheme as

q k+1 + eq k-1 = prox M (q k+1 ) IR n + ; 2q k -(1 -e)q k-1 + h 2 M -1 (q k + 1)F (t k+θ , q k+θ , v k+θ ) (25)
On the qualitative point of view, the main difference between these two schemes is the mechanical nature of the projected variable. In the Moreau-Jean scheme, the variable which is projected is homogeneous to a velocity. One of the interesting remarks is that the Newton impact law is respected for the discrete velocity in a very natural way by noting that

U k+1 = -eU k if P k+1 > 0. ( 26
)
This fact leads to a straightforward interpretation of the discrete multiplier as a mechanical impulse. However, the projection of the velocity onto the tangent cone of IR n + yields a slight violation of the constraints which occurs at the impact.

In the Schatzman-Paoli scheme, the generalized coordinates is directly projected onto the admissible set. The result is that there is no violation of the discrete constraints when e = 0. On the contrary, the discrete velocity does not satisfy the Newton impact law. Furthermore, the multiplier involved in the projection of the coordinates has no direct mechanical meaning. The Newton impact law is satisfied after several steps. On the other hand, on the simple linear oscillator example, the scheme does not generate artificial rebound in presence of flexibility.

A first solution : A direct projected scheme

In this section, we first propose a scheme which both satisfies in discrete time the position constraints and the velocity constraints, i.e., the impact law. This scheme is an adaption of the Moreau-Jean scheme based on the direct use of the Gear-Gupta-Leimkuhler (GGL) method [START_REF] Gear | Automatic integration of Euler-Lagrange equations with constraints[END_REF]. Since this scheme will serve as the basis for an improved version in Section 5, we detail its local order of accuracy, its implementation, the choice of the activation rule and its main drawback : the chattering at contact.

General presentation of the direct projected scheme

Let us start by considering the following "augmented" system

               M (q(t))dv = F (t, q(t), v + (t))dt + G(q(t))di, q(t) = v + (t) + G(q(t))µ(t), U (t) = G (q(t))v(t) if g α (q(t)) 0, then 0 U α,+ (t) + eU α,-(t) ⊥ di 0, 0 g(q(t)) ⊥ µ(t) 0. (27)
where µ(t) is a new multiplier which corresponds to the redundant constraints g(q(t)) 0. Thanks to Moreau's viability lemma [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF], we expect that the multiplier is identically zero and that the solution of ( 27) is equivalent to the solution of ( 16). The proposed time-stepping scheme, called the direct projected scheme reads as

                   M (q k+θ )(v k+1 -v k ) -hF k+θ = G(q k+1 )P k+1 , q k+1 = q k + hv k+θ + G(q k+1 )τ k+1 , U k+1 = G (q k+1 ) v k+1 , g k+1 = g(q k+1 ), if ḡα k+1 0, then 0 U α k+1 + eU α k ⊥ P α k+1 0. 0 g k+1 ⊥ τ k+1 0. ( 28
)
The discrete multiplier

τ k+1 ≈ t k+1 t k µ(t) dt (29)
ensures the constraints at the position level in discrete time that is g k+1 0.

Empirical convergence and order analysis

The global order of the scheme is shown on the two very simple systems described in Example 1 and Example 2. For the sake of simplicity we choose for the prediction of the constraints, ḡk+1 the fully explicit forecast ḡk+1 = g k . Other choices will be discussed in § 4.4. Figure 2 shows the global order of convergence of the direct projected scheme follows the order of the Moreau-Jean scheme and the Paoli-Schatzman scheme. We can notice that the projection does not improve the global quality of the solution. Let us now study the local order of accuracy of the scheme and let us start with the discrete multiplier τ . On Figure 3, it is shown that the discrete multiplier τ k+1 is of order O(h). This result is proven in Proposition 1 under the following assumptions: Assumption 1 (Existence and uniqueness) A unique global solution over [0, T ] for augmented Moreau's sweeping process ( 27) is assumed such that q( • ) is absolutely continuous and admits a right velocity v + ( • ) at every instant t of [0, T ] and such that the function

v + ∈ LBV ([0, T ], IR n ).
Assumption 2 (Smoothness of data) The following smoothness on the data will be assumed: a) the inertia operator M (q) is assumed to be of class C 0 and definite positive, b) the force mapping F (t, q, v) is assumed to be of class C 0 , c) the constraint functions g(q) are assumed to be of class C 1 and d) the Jacobian matrix G (q) = ∇ q g(q) is assumed to have full-row rank.

Proposition 1 Under Assumptions 1 and 2, the multiplier τ k+1 is of order h that is

τ k+1 = O(h). ( 30
)
proof Let us consider the first order Taylor expansion of the constraints, g(q) at q k

g k+1 = g(q k+1 ) = g(q k ) + ∇ q g(q k )(q k+1 -q k ) + O( q k+1 -q k 2 ) = g(q k ) + G (q k ) [hv k+θ + G(q k+1 )τ k+1 ] + O( q k+1 -q k 2 ) = g(q k ) + hG (q k )v k+θ + G (q k )G(q k+1 )τ k+1 + O( q k+1 -q k 2 ) ( 31
)
Let us denote the index sets of active constraints by β = {i | τ i k+1 > 0}. Since g β k+1 = 0 and G has full row rank, the solution of the LCP can be written as

τ β k+1 = -[G (q k )G(q k+1 )] -1 ββ g β (q k ) + hG β• (q k )v k+θ + O( q k+1 -q k 2 ) (32)
Since g β k+1 = 0, we also obtain

g β (q k ) = -∇ q g(q k )(q k+1 -q k ) + O( q k+1 -q k 2 ) = O( q k+1 -q k ) (33)
From ( 32) and ( 33), we get the following estimate

τ β k+1 = -[G (q k )G(q k+1 )] -1 ββ hG β• (q k )v k+θ + O( q k+1 -q k ) (34)
and more generally, since τ i k+1 = 0 for i / ∈ β, we get

τ k+1 = O( q k+1 -q k ) + O(h). (35)
since v k+θ is assumed to be bounded. Inserting this estimate in

q k+1 -q k = hv k+θ + G(q k+1 )τ k+1 (36)
we get that

O( q k+1 -q k ) = O(h) (37)
and we conclude that τ k+1 = O(h).

The following result is a straightforward extension of the Proposition 1 in [3].

Proposition 2 Under Assumptions 1 and 2, the local order of consistency of the Moreau-Jean time-stepping scheme with projection for the generalized coordinates is e q = O(h) and at least for the velocities e v = O(1).

proof The estimate e v on the velocity is trivial if we recall that M -1 k+θ (F k+θ + p k+1 ) is bounded on [t k , t k+1 ]. The BV function v( • ) is also bounded on [t k , t k+1 ] then we have that e v is also bounded. Therefore, we obtain e v = O(1). Using Lemma 1 in [3] for v + ∈ BV (I, IR n ) , we get

t k+1 t k v(s) ds -h(θv + (t k+1 ) + (1 -θ)v + (t k )) C(θ)h var(v + , I), (38)
with

C(θ) = θ if θ 1/2 and C(θ) = 1-θ otherwise. Since v k = v + (t k ) and v k+1 = v(t k+1 )+O(1)
, we obtain for (38)

t k+1 t k v(s) ds -h(θv k+1 + (1 -θ)v k ) -θO(h)) = e q +G(q k+θ )τ k+1 -θO(h) C(θ)h var(v + , I).
(39) Using the result of Proposition 1, τ k+1 = O(h), we get

e q + O(h) C(θ)h var(v + , I), (40)
which completes the proof.

Implementation

In this section, several possible implementations and variants of the direct projected scheme [START_REF] Jourdan | A Gauss Seidel like algorithm to solve frictional contact problems[END_REF] are detailed. Let first us discuss the strict implementation of [START_REF] Jourdan | A Gauss Seidel like algorithm to solve frictional contact problems[END_REF]. The nonlinear residual function is defined as

R(v, q) = M (θq + (1 -θ)q k )(v -v k ) -hF (t k+θ , θq + (1 -θ)q k , θv + (1 -θ)v k ) -G(q)P k+1 q -q k -h(θv + (1 -θ)v k ) -G(q)τ k+1 .
(41) For given values of P k+1 and τ k+1 , the unknowns v k+1 and q k+1 which solve the first two lines of (28) satisfy R(v k+1 , q k+1 ) = 0. This set of nonlinear equations can be solved by standard solvers [START_REF] Ortega | Iterative Solution of Nonlinear Equations in Several Variables[END_REF]. For Newton's method, the solution is sought as a limit of the sequence {v n , q n } n∈IN such that

v 0 = v k , q 0 = q k , R L (v n+1 , q n+1 ) ∆ = R(v n , q n ) + ∇ v R(v n , q n )(v n+1 -v n ) + ∇ q R(v n , q n )(q n+1 -q n ) = 0.
(42) In order to be self-contained but without enter into deepest details, we will describe a semi-Newton procedure in which the inertia matrix and the Jacobian of the constraints are only updated in a fixed point way and then evaluated at (v n , q n ). It amounts to neglecting the Jacobian of M and G in the Newton loop. We solve

R(v n , q n )+ M (q n θ ) -hθ∇ v F (t k+θ , q n θ , v n θ ) -hθI n×n (v n+1 -v n )+ -hθ∇ q F (t k+θ , q n θ , v n θ ) I n×n (q n+1 -q n ) = 0. ( 43
) where q n θ (resp. v n θ ) denotes θq n + (1 -θ)q k (resp. θv n + (1 -θ)v k ) and M n = M (q n θ ).
Let us denote the values of P k+1 and τ k+1 at the Newton iteration n by P n+1 and τ n+1 . After simple algebraic manipulations, we obtain

     [M n + hθC n ] (v n+1 -v n ) = -hθK n (q n+1 -q n ) -M (q n θ )(v n -v k ), +hF (t k+θ , q n θ , v n θ ) + G(q n )P n+1 , q n+1 = q k + hv n+1 θ + G(q n )τ n+1 , (44)
where C n = -∇ v F (t k+θ , q n θ , v n θ ) denotes the tangent damping matrix and K n = -∇ q F (t k+θ , q n θ , v n θ ) the tangent stiffness matrix. A substitution in [START_REF] Radjaï | Discrete-element modeling of granular materials[END_REF] of the second equation into the first one yields

         M n + hθC n + h 2 θ 2 K n (v n+1 -v n ) = -M (q n θ )(v n -v k ) + hF (t k+θ , q n θ , v n θ ) + hθK n (q k -q n + h(1 -θ)v k ) + hθK n G(q n )τ n+1 + G(q n )P n+1 , q n+1 = q k + hv n+1 θ + G(q n )τ n+1 . (45)
In condensed matrix notation we obtain

M 0 h I v n+1 -v n q n+1 -q n = f -hv n + G(q n ) hθK n G(q n ) 0 G(q n ) P n+1 τ n+1 , (46)
with the iteration matrix M defined by

M = M n + hθC n + h 2 θ 2 K n , (47)
and

f = -M n (v n -v k ) + hF (t k+θ , q n θ , v n θ ) + hθK n (q k -q n + h(1 -θ)v k ). ( 48
)
The nonlinear constraints g(q k+1 ) 0 are also linearized by considering the following nonlinear residual function

R g (y, q) = y -g(q), (49)
and its first order expansion,

R g,L (y n+1 , q n+1 ) ∆ = y n -g(q n ) + G (q n ) (q n+1 -q n ) + y n+1 -y n . (50)
By denoting y n+1 = g n+1 , we get the following system of linearized constraints

g n+1 = g(q n ) + G (q n ) (q n+1 -q n ) 0. ( 51
)
Let us note the set of active constraints in velocity by

I v = {α | ḡα k+1 0}
The unknown local velocity vector at step n for this set of constraints is shortly denoted as written as U n+1 = [U ,n+1,α , α ∈ I v ] . The following Mixed Linear Complementarity Problem (MLCP) must be solved at each Newton's loop

    M 0 0 0 h I 0 0 -G (q n ) 0 I 0 0 -G (q n ) 0 I         v n+1 -v n q n+1 -q n U n+1 -U n g n+1     =     f -hv n 0 g(q n )     +     G(q n ) hθK n G(q n ) 0 G(q n ) 0 0 0 0     P n+1 τ n+1 , 0 U n+1 ⊥ P n+1 0, 0 g n+1 ⊥ τ n+1 0. ( 52
) for v n+1 , U n+1 , P n+1 , g n+1 , τ n+1 . Since the first matrix in ( 52) is lower block triangular and invertible if M is invertible, we can build a condensed Linear Complementarity Problem (LCP) as follows

           U n+1 g n+1 = W P n+1 τ n+1 + a, 0 U n+1 ⊥ P n+1 0, 0 g n+1 ⊥ τ n+1 0, (53)
with

W = G (q n ) M -1 G(q n ) h 2 θ 2 G (q n ) M -1 K n hθG (q n ) M -1 G(q n ) G (q n )G(q n ) + h 2 θ 2 G (q n ) M -1 K n , (54)
and

a = U n + G (q n ) M -1 f g(q n ) + G (q n ) q k -q n + hθ(v n + M -1 f ) . ( 55
)

A decoupled implementation

In this section, we propose a simplified implementation of the method by updating higher order terms in h in the Newton iterations in a fixed point way. Let us consider the following modified version of (45):

             M (q n θ ) + hθC n + h 2 θ 2 K n (v n+1 -v n ) = -M (q n θ )(v n -v k ) + hF (t k+θ , q n θ , v n θ ) + hθK n (q k -q n + h(1 -θ)v k ) +hθK n G(q n )τ n + G(q n )P n+1 (56a) q n+1 = q k + hv n+1 θ + G(q n )τ n+1 . ( 56b
)
Note that the value of τ k+1 in the right hand side of the equation is taken at step n that is τ n . Doing that way, the following LCP can be first solved

     U n+1 = G (q n ) M -1 G(q n ) P n+1 + U n + G (q n ) M -1 b + hθG (q n )K n G(q n )τ n , 0 U n+1 ⊥ P n+1 0. ( 57
) Knowing the value P n+1 , hence the value of v n+1 given by (56a), a second LCP is solved as follows

   q n+1 = G (q n )G(q n ) τ n+1 + hv n+1 θ , 0 q n+1 ⊥ τ n+1 0. ( 58
)
The main interest of this decoupled implementation lies in the formulation of two smaller LCPs [START_REF]CASCADE Technology[END_REF] and [START_REF] Wriggers | Nonsmooth mechanical systems with unilateral contact 6[END_REF] rather than the larger one [START_REF] Stewart | Rigid body dynamics with friction and impact[END_REF]. Furthermore, the matrix involved in the first LCP at the velocity level ( 57) is identical to the LCP matrix that is used in the standard Moreau-Jean scheme. This decoupled implementation needs only two slight modifications of the standard Moreau-Jean scheme adding a new term in the right hand side of (56a) and performing the decoupled projection by solving the LCP [START_REF] Wriggers | Nonsmooth mechanical systems with unilateral contact 6[END_REF]. The algorithm is detailed in Algorithm 1.

Remark 3 It can also be interested to consider the following rule for the correction of the position

q(t) = v + (t) + N (q(t))G(q(t))µ(t), (59)
where N (q) ∈ IR n×n is a positive definite matrix. It amounts to choosing a special metric or the projection onto the constraints. Since the projection of the velocity is based on the kinetic metric, it may convenient to select the same kind of projection for the position. Applying the proposed discretization and in the decoupled implementation, the following LCP equivalent to ( 58)

is obtained    q n+1 = G (q n )N (q n )G(q n ) τ n+1 + hv n+1 θ , 0 q n+1 ⊥ τ n+1 0. ( 60
)
Choosing the N (q n ) = M -1 allows one to form the same matrix for the first LCP in velocity and the second one in position.

In the numerical practise the decoupled implementation performs very efficiently. This is manly due to the fact that the neglected coupling terms are of order h 2 and hence have a weak influence in the behavior of the scheme. Therefore, the decoupling strategy does only slightly change the average number of Newton's iterations and fixed point iterations (see Table 3 for an example).

Algorithm 1 Direct Projected Algorithm for one time-step (decoupled implementation)

Require: h time-step, I = [t k , t k+1 ], tol ∈ IR a user tolerance Require: q k , v k initial conditions of the step. Ensure: q k+1 , v k+1 , P k+1 , τ k+1 // update the index Set

I v ← {α | ḡα k+1
0} with one the rule Eq. ( 61) or (63).

v n ← v k ; q n ← q k , P n ← 0; τ n ← 0 // Start Newton's loop while R(v n+1 , q n+1 ) > tol or R g (g n+1 , q n+1 ) > tol do
// Solve the first LCP (Eq. ( 57)) for P n+1 .

     U n,u+1 = G (q n ) M -1 G(q n ) P n+1 + U n + G (q n ) M -1 b + hθG (q n )K n G(q n )τ n 0 U n+1 ⊥ P n+1 0
// Update the velocity v n+1 (Eq. (56a))

M (q n θ ) + hθC n + h 2 θ 2 K n (v n+1 -v n ) = -M (q n θ )(v n -v k ) + hF (t k+θ , q n θ , v n θ ) + hθK n (q k -q n + h(1 -θ)v k ) +hθK n G(q n )τ n + G(q n )P n+1
// Solve the second LCP (Eq. ( 58)) for q n+1 and τ n+1 .

   q n+1 = G (q n )G(q n ) τ n+1 + hv n+1 θ 0 q n+1 ⊥ τ n+1 0 n ← n + 1 end while v k+1 ← v n , P k+1 ← P n q k+1 ← q n , τ k+1 ← τ n

Convergence criteria

The convergence of the scheme with the coupled or decoupled implementation is ensured by checking against the prescribed tolerance the norm of the residual term (41) R(v n+1 , q n+1 ) and the residual term (49) R g (g n+1 , q n+1 ). Furthermore, we check that the complementarity condition are satisfied inside the LCP solver.

4.4 How to choose a good prediction scheme for the activation of constraints ?

In this section, we discuss the prediction of the constraints, i.e., the computation of the value of ḡk+1 . The choice of the forecast ḡk+1 will be discussed in details as it plays a leading part in the global behavior of the scheme. The following implementation will be discussed:

• The fully explicit forecast consists in evaluating the constraints with the values of the previous time-step,

ḡk+1 = g k + γhU k , ( 61
)
where γ is usually chosen in [0, 2]. The main interest of the method lies in its simplicity. For γ = 0, the constraints on velocity are activated when the constraints in position are violated. Without projection, this rule is very robust, but may yield negative violations of order O(h).

For γ > 0, this scheme uses an extrapolation of the trajectory to guess if a constraint will be violated within the time-step. Let us recall that if the trajectory is absolutely continuous, the velocity is not. Therefore, the extrapolation can not be better that O(h). One of the interest of this approach is an activation of the constraints can be set before the violation of the constraints. Unfortunately, as the constraints are treated at the velocity level, a reaction force can be imposed even if g k+1 > 0. The only thing that can be said is that g k+1 = O(h).

• The free-position forecast is based on the evaluation of the position without any reaction forces due to unilateral constraints, that is

ḡk+1 = g k + h(θU free + (1 -θ)U k ) ( 62
)
where U free is the relative free velocity at contact and the parameter θ is chosen such that θ ∈ (0, 1].

Remark 4

The question of a semi-implicit forecast based on the update of the position inside the Newton loop without the projection onto the constraints, i.e.

g k+1 = g(q k + hv α k+θ ) ( 63
)
or a fully implicit forecast based on the implicit evaluation of the position without the projection onto the constraints, that is

g k+1 = g(q k + hv k+θ ) ( 64
)
is a difficult question. Since the activation of the constraints will modify the velocity at the end of the step, the previous choices in (63) and (63) are most of time inconsistent and yields the cycling in the activation of constraints without any convergence. However, we will see in Section 5 how such a similar idea can be used without cycling by augmenting in a unique way the set of active constraints.

Artificial oscillations, chattering and energy balance

Let us study the behavior of the projected scheme [START_REF] Jourdan | A Gauss Seidel like algorithm to solve frictional contact problems[END_REF] together with the explicit forecast (61) on the bouncing ball example 1. In Figure 4, a cycling behavior is observed after the finite accumulation when the ball would come to rest on the ground. This behavior was already observed in [START_REF] Studer | Numerics of Unilateral Contacts and Friction. -Modeling and Numerical Time Integration in Non-Smooth Dynamics[END_REF]. One of the explanation is that the energy brings to the system by the projection exactly compensates the energy that is dissipated during the impact. One of the consequences of this cycle is that the constraint is never stabilized and the reaction forces does not converge towards the constant value which counteracts the weight of the ball. Note that the chattering is observed in Figure 4 whatever the choice of the prediction parameter γ in (61). However, we note in Figure 5 that the chattering is not observed with the free position forecast (62). Indeed, when this latter forecast rule is used on the bouncing ball example, the constraint is never violated, i.e g k+1 > 0 and therefore τ k+1 is identically zero. In this particular example, the projected Moreau scheme ( 28) is equivalent to the original one.

In the case of the bouncing ball, it clearly appears that the projection brings some energy to the system. In Figure 6, the discrete kinetic energy, potential energy and the total mechanical energy are plotted versus time. Each projection onto the constraints affects the energy balance. At the instants of the projection, an increase of the potential energy is shown. In Figure 7, the coefficient of restitution is chosen equal to one such that the continuous system is conservative. We observe an increase of the total amount of energy of the system. Although the convergence is not rigorously proved in the paper, it does seem to call into question since the amplitude of the oscillatory artifact goes to zero as the time-step vanishes. The main weakness generated by the chattering is the quality of the approximation for a finite time-step. With a fixed time-step, the bouncing ball example never show a still equilibrium and the chattering never stops. Let us give another example in this section where the oscillations prevent to reach a static equilibrium. Let us consider the mechanism described in Figure 8. The system is composed of a roller of radius R submitted to an external applied torque Γ and a slider of mass m which hits the roller. The contact of the roller and the slider is modeled by the Signorini condition, the Newton impact law with a coefficient e < 1 and the Coulomb's friction of coefficient µ. The whole system is submitted to gravity. If the applied torque is less than µmg, a static equilibrium must be reached.

If the simulation is achieved with the standard Moreau-Jean scheme, the equilibrium is correctly approximated and observed with a finite time-step. With the direct projected scheme, the slider never stops to bounce and each time the contact is lost, we observe a slip of the roller under the slider. Therefore, the roller never reaches its static equilibrium.

Conclusion on the direct projected scheme

As often with a time-discretization method which is assumed to be convergent, the approximate solution in discrete-time does not keep all the properties of the continuous time solution. This is the goal of the geometric time-integration methods to ensure the conservation of properties in discrete time [START_REF] Hairer | Geometric numerical integration. Structure-preserving algorithms for ordinary differential equations[END_REF]. With the direct projected scheme, the constraints in velocity and position are satisfied in discrete time. The algorithm keeps the order of the standard Moreau-Jean scheme and the multiplier associated with the projection vanishes at the order O(h). With the decoupled approach, the implementation is straightforward and requires only slight modifications of the standard Moreau-Jean scheme.

The main drawback is the occurrence of the chattering at contact in some special configurations. The chattering can have major consequences as we can lose in several situations the existence of equilibria. This is not satisfactory for our purpose and the goal of the next section is to remedy to this drawback, retaining the favorable properties of the direct projected scheme.

A combined projection/activation algorithm

In this section, we present a scheme to circumvent a part of the problem listed in the previous section. The main goal of the scheme is to activate consistently the constraints at the velocity level with respect to the set of constraints which will be projected in the current time-step. Especially, we want to avoid the projection onto the constraints if the constraint at the velocity level is not activated. We have seen in Section 4 that these phenomena causes chattering in the direct projection scheme.

Presentation of the combined scheme

The combined scheme is based on the iterations denoted by ν of the following two steps :

1. The projection step is based on the solution of the following system

                       M (q k+θ )(v k+1 -v k ) -hF k+θ = G(q k+1 )P k+1 , q k+1 = q k + hv k+θ + G(q k+1 )τ k+1 , U k+1 = G (q k+1 ) v k+1 , g k+1 = g(q k+1 ), for all α ∈ I ν      0 U α k+1 + eU α k ⊥ P α k+1 0, g α k+1 = 0, τ α k+1 , if P α k+1 > 0, 0 g α k+1 ⊥ τ α k+1 0 otherwise . ( 65
)
for a given index set I ν of active constraints.

2. The activation step computes the index set I ν of active constraints by checking for a given value of g k+1 if the constraint is satisfied or not. Starting form I 0 = ∅, at each iteration ν, the activation performs the following operation

I ν+1 = I ν ∪ α | g α k+1 0 ( 66
)
The iterates of the solution q k+1 , v k+1 depends on the iteration number ν. In order to avoid useless complexity in the notation, we skip the superscript ν when there is no ambiguity. With this convention, the algorithm is described in Algorithm 2.

Algorithm 2 Activation/Projection Algorithm for one time-step

Require: h time-step, I = [t k , t k+1 ] Require: q k , v k initial conditions of the step. Ensure: q k+1 , v k+1 , P k+1 , τ k+1 // Initialization ν ← 0 I 0 ← ∅, I -1 ← {-1} while I ν = I ν-1 do // Solve the projection step (Eq. ( 65)).

                       M (q k+θ )(v k+1 -v k ) -hF k+θ = G(q k+1 )P k+1 , q k+1 = q k + hv k+θ + G(q k+1 )τ k+1 , U k+1 = G (q k+1 ) v k+1 , g k+1 = g(q k+1 ), for all α ∈ I ν      0 U α k+1 + eU α k ⊥ P α k+1 0, g α k+1 = 0, τ α k+1 , if P α k+1 > 0, 0 g α k+1 ⊥ τ α k+1 0 otherwise .
// Update the index set (Eq. ( 66))

I ν+1 ← I ν ∪ α | g α k+1 0 ν ← ν + 1 end while

Comments

Let us first note that the results in Proposition 1 and 2 are still valid for this combined projection/activation scheme.

The first step (65) is very similar to the scheme presented in Section 4. Note that only the index set onto we project is modified. This a priori minor modification in the implementation is nevertheless crucial for the qualitative behavior of the scheme.

Let us give now some insight on the behavior of the scheme. The first iteration is performed with I 0 = ∅ amounts to computing the free velocity and the free position of the system. The goal is to perform the first activation of the constraints. In other terms we performed first a step disregarding the constraints and we check what are the constraints which are not satisfied. The indices of these constraints that are violated compose the set I 1 , that is the first set of forecast activated constraints.

Let us introduce a new index set I ν c , subset of I ν , i.e. I ν c ⊂ I ν such that

I ν c = {α ∈ I ν | P α k+1 > 0}. ( 67
)
In the projection step, the following rule is used for the projection

for all α ∈ I ν      0 U α k+1 + eU α k ⊥ P α k+1 0, g α k+1 = 0, τ α k+1 , if α ∈ I ν c , 0 g α k+1 ⊥ τ α k+1 0 otherwise . ( 68
)
We note that the projection is only performed for the active constraints whose index belong to I ν c (positive contact impulse P α k+1 > 0). For the constraints α ∈ I ν ⊂ I ν c , we project on the manifold defined by {g α k+1 = 0, α ∈ I ν c }. Otherwise, we only require to have non penetration. Even if there is no rigorous mathematical proof, this correction rule shown the best compromise in practice. Indeed, imposing g α k+1 = 0 for all constraints in I ν may lead to unfeasible problem. At the end of the time step, we ensure that there are no violated constraints and no projected constraints without satisfying the jump rule 0 U α k+1 + eU α k ⊥ P α k+1 0.

RR n°8029

Concerning the implementation, the projection step is very similar to the direct projection scheme presented in Section 4. Its implementation follows the same line as in Section 4.3. It can be decoupled or not. As for the evaluation of I ν , we need to precise the rule for the computation of I ν c in order to obtain a proper LCP without switched-off constraints triggered by a conditional statement. The rule that we have chosen is

I ν c = {α ∈ I ν | P n,α k+1 > 0} (69)
for a coupled scheme and

I ν c = {α ∈ I ν | P n+1,α k+1 > 0} (70)
for the decoupled scheme since we known the value of P n+1 from the solution of the first LCP.

Rocking Block example

The efficiency of Algorithm 2 is firstly demonstrated on the block (Example 3). In terms of spurious oscillations, the others academic examples are simpler to deal with. Since the rocking block example has two strongly coupled nonlinear constraints, we focus our attention in this section on it. More complex examples will be treated in Section 6.

In Figure 9, the results of the scheme based on a direct projection ( 28) is compared with Algorithm 2 based on the combination of projection and activation steps. The simulation parameters are as follows : l = 1.5 m, L = 1 m, x 0 = 0.0 m, y 0 = 1.0 m, θ 0 = 0.2 rad, ẋ0 = 0.0 m.s -1 , y 0 = 0.0 m.s -1 , θ 0 = 0.2 rad.s -1 , m = 1.0 kg, e = 0.5, t 0 = 0.0 s, T = 2.0 s, h = 10 -02 s, θ = 1/2 . In Figure 9(c), the spurious oscillations are observed when the rocking block reaches its equilibrium after a finite accumulation on impact. In Figure 9(d), we remark that the spurious oscillations are suppressed and the block reaches its equilibrium without any troubles. On the energetic point of view (Figures 9(e) and 9(f)), the combined approach dissipates more energy than the standard direct approach. With the requirement that the time-step length does not depend on events, dissipation of energy is the price to pay to avoid chattering with a projection onto the constraints.

6 Demonstrative applications 6.1 Software aspects Algorithms 1 and 2 are implemented in the open-source Siconos software [6,[START_REF] Siconos | A software for modeling and simulation of nonsmooth dynamical systems[END_REF]. This software provides a general framework for implementing numerical time integration schemes of nonsmooth dynamical systems. The solver for the discrete frictional contact problem that is used in this paper is a projected Gauss-Seidel solver [START_REF] Jourdan | A Gauss Seidel like algorithm to solve frictional contact problems[END_REF] developed in the Siconos/Numerics library of solvers. In the most of the following examples, the Siconos/Multibody library is used to model and simulate multibody systems with tri-dimensional contact, impacts and Coulomb's friction. This library allows the user to instantiate a Newton/Euler model linked to a geometrical representation in a industrial CAD library. In our examples, we use the open-source Software CAD library OpenCascade [START_REF]CASCADE Technology[END_REF] and its python wrapper PythonOCC [START_REF]PythonOCC. 3D CAD/CAE/PLM development framework for the Python programming language[END_REF].

The impacting elastic bar

The example which is considered in this section is the problem of a one-dimensional bar which hits a rigid wall with a constant initial velocity v 0 . We assume that the assumption of small displacements holds and the constitutive law is linear elastic. The problem depicted in Figure 10 consists of a linear elastic bar of cross section area S, with a Young modulus E and density ρ and of length L. This problem has been widely used in the literature (see for instance [START_REF] Hughes | A finite element method for a class of contact-impact problems[END_REF][START_REF] Hughes | The Finite Element Method, Linear Static and Dynamic Finite Element Analysis[END_REF][START_REF] Wriggers | Nonsmooth mechanical systems with unilateral contact 6[END_REF]) because of its interest from the mathematical and computational point of view.

From the mathematical point of view, this example is one of the rare example in elastodynamics with unilateral contact for which the mathematical properties are known in terms of existence and uniqueness. The problem of an elastic bar are indeed discussed from the mathematical point of view in [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF] with an associated numerical scheme. It shown that the problem can be uniquely solved without requiring to an additional energetic condition. More general cases are complex.

The vibration of a string on a point-shaped obstacle enjoys also existence and uniqueness. It has been proven in [START_REF] Schatzman | Un problème hyperbolique du 2ème ordre avec contrainte unilatérale: La corde vibrante avec obstacle ponctuel[END_REF] without asking for an additional energy condition. For a concave continuous obstacle, an energy condition has to be added to retrieve uniqueness (see [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF]). This means that we are as in the rigid body case where an impact law is added to solve the problem uniquely.

For more complex geometries, the problem remains open. For the elastic impacting bar without external forces, a simple solution can be constructed as it has been done in [START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF]. Let us denote by c 0 = E/ρ the wave speed. The contact time denoted by T corresponds to twice the time of rhe travelling of the elastic wave in the bar, that is T = 2L/c 0 . Within the contact time, the contact force is constant and equal to r = ESv 0 /c 0 .

From the computational point of view, the problem of an impacting elastic bar has also another special interest: it enables us to exhibit possible spurious oscillations of the standard numerical schemes (Newmark, HHT, α-schemes) of the contact velocity and forces. This spurious oscillations which are very different in nature with the chattering observed with the direct projected scheme can explain at least by two causes. The first one is the nonsmoothness of solutions. When the tip of the bar reaches the wall, a jump in the velocity of the tip and in the contact forces occurs. If the scheme approximates these unknown with a second-order approximation, oscillatory artifacts can be observed [START_REF] Laursen | Improved implicit integrators for transient impact problems -Dynamical frcitional dissipation within an admissible conserving framework[END_REF][START_REF] Deuflhard | A contact-stabilized Newmark method for dynamical contact problems[END_REF][START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF]. Most of the time, this first cause is circumvented by using a first-order fully implicit treatment of the contact forces [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF][START_REF] Laursen | Improved implicit integrators for transient impact problemsgeometric admissibility within the conserving framework[END_REF][START_REF] Laursen | Improved implicit integrators for transient impact problems -Dynamical frcitional dissipation within an admissible conserving framework[END_REF], or a direct use of the impulse [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid bodies collections[END_REF]. The second cause of oscillations is the index of the DAE resulting after a time and space discretization. If the unilateral constraint is written at the position level, the index of the DAE for a closed contact is 3. It is well known that the direct time-integration of DAEs of index 3 generates spurious oscillations in the time derivative of the constraints. To remove this oscillation, some attempts have been based on a) a writing of the constraints at the velocity level [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF], that is to perform an index reduction similarly as in the Moreau-Jean scheme, b) a mass redistribution which consists in removing the mass from the contact boundaries [START_REF] Khenous | On the discretization of contact problems in elastodynamics[END_REF][START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF] or c) a contact stabilization of the relative velocities at contact [START_REF] Deuflhard | A contact-stabilized Newmark method for dynamical contact problems[END_REF].

In the article, we propose to solve the impacting elastic bar problem by using an index reduction technique with a stabilization of the constraint based on a projection as in the seminal work of Gear, Gupta and Leimkuhler [START_REF] Gear | Automatic integration of Euler-Lagrange equations with constraints[END_REF]. The constraint is written at the velocity level and the dynamics is time integrated solved with the help of the Moreau-Jean scheme. Using an index reduced formulation at velocity level avoids the spurious oscillations at contact of the forces and the relative velocity. By the way, since the structure is space-discretized, it appears as a finite-freedom mechanical system with unilateral contact for which we know that we have to provide an energy condition under the form of an impact law. We choose a perfect inelastic impact law to mimic the continuous time solution where the contact stays closed for a finite time interval. In more complex situations, this question is also open. The price to pay in using a velocity based formulation for the constraint is the drift the constraint at the position level. This drift, or violation of the constraints is fixed by the projection onto the constraints and the additional multiplier.

In the results that follow, the bar is space discretized by N linear rod finite elements. The elementary mass and stiffness matrices are

M e = 1 6 ρSl e 2 1 1 2 , K e = 1 l e ES 1 -1 -1 1 , ( 71
)
where l e = L/N is the length of an element. The material parameters are summarized in Table 1.

In Figure 11, the spurious oscillations of the contact velocity and the contact force are depicted for the trapezoidal rule with a contact condition at the position level and an implicit treatment of the contact force. The scheme is similar to the Moreau-Jean scheme with θ = 1/2, but with a position level constraint. Although the constraint at the position level is perfectly satisfied on the time of contact equal to T = 3.8545 × 10 -04 s, we note that the contact velocity oscillates between two extreme values at each time step. The oscillations of the contact forces are also observed around the solution value r = 1271.472 N. Note that spurious oscillations are also observed for the contact velocity after the contact time, around the solution value -v 0 . This is mainly due to the jump of the contact force that excites the artificial high frequency modes of the bar induced by the space discretization. Using the pure trapezoidal rule do not enable us to introduce a small amount of numerical damping which would allow us to damp these latter oscillations. With the Moreau-Jean scheme, the spurious oscillations within the contact time are not observed in Figure 12. The constraint at the velocity level with a coefficient of restitution e = 0 yields a perfect stabilization of the velocity. The post contact oscillations due to the high-frequencies modes are still observed due the lack of numerical damping. A very small oscillation of the contact force occurs in the first step after the bar reaches the obstacle can be also observed. Is is mainly due to the fact that we deal with a finite-freedom system and the flexible structure is subjected to an impact. In order to understand a little bit further this phenomenon, Figures 13 and14 provide us with an analysis of the contact force with respect to the element size and the time-step. In Figure 13, it appears that for a decreasing number of elements, an increasing peak appears in the contact force. This peak reveals the occurrence of an impact when the bar reaches the obstacle. The peak increases since the finite mass of the last element involved in the contact increases as well. In Figure 14, we observe that the peak increases with a decreasing time-step for a constant number of elements. In the limit, we may expect that the value of this peak goes to infinity which is another expression of the occurrence of an impact. For a vanishing time-step, we converge to a finite-freedom mechanical systems with impact and the contact forces goes to infinity. The right unknown is then the impulse. As we said earlier, the presence of an impact and an impulse calls for the introduction of an impact law. In this simple case, we known a priori that the bar should in contact for a finite-time interval. This is the reason why we choose a Newton impact law with a coefficient of restitution equal to 0. In Figure 15, we report the results of the same scheme with e = 0.95. Since the discretized model is a finite-freedom mechanical system, the choice of a coefficient of restitution in [0, 1]| yields a well-posed problem. However, in the limit when the mesh size vanishes, we cannot expect to retrieve the elastic bar problem as the oscillations of the contact velocity shows.

This section is completed with the application of the combined projected scheme to the elastic bar example. In the previous example, the violation of the constraints ranges from 2.5 × 10 -05 m to 2.5 × 10 -08 m for a time-step from h = 5 × 10 -05 to h = 5 × 10 -07 . When the projected scheme is applied, the violation is equal to zero at the machine accuracy. In Figures 16 and17 scheme. For the contact force, spurious oscillations are observed which are due to the jump in position that excites the high-frequency modes of the structure. Nevertheless, and contrary to the real impact at the velocity level, these oscillations vanish as the time-step decreases, since the projection multiplier also vanishes. To conclude, we are able, for a flexible structure, to satisfy at the same time the constraints in position and in velocity. The high-frequency mode excitation of the structure calls for the introduction of numerical damping. These can be achieved by using 1 θ 1/2 as it is illustrated in Figure 18. Nevertheless, the use of numerical damping with the θ-method implies a dissipation over the whole range of frequencies. As we can observed, the response of the structure is largely damped. In [START_REF] Chen | A Newmark-Type Integrator for Flexible Systems Considering Nonsmooth Unilateral Constraints[END_REF], a consistent adaptation of the HHT scheme and the α-schemes to the impact mechanics is proposed to this end.

Pendulum in a ring.

This very simple example is a pendulum hung to a ring that is depicted as a schematic view (Figure 19(a)). The ring of radius R = 1 m and thickness r = 2 × 10 -3 m can only rotate about the y-axis (its axis of revolution) due to a perfect revolute joint. The mass of the ring is 1 kg and its moment of inertia about the y-axis is 0.8 kg • m -2 . The mechanical model of the pendulum is a rod of mass 5kg and the moment of inertia about the y-axis of 1kg • m -2 . The contact of the pendulum in the interior of the ring is punctual. The model of contact is the Signorini condition, the Coulomb friction with a coefficient µ = 0.3 and the perfectly inelastic newton impact with e = 0. The system is submitted to the gravity load g = 9.81m • s -2 .

At the initial time, the pendulum is left in a vertical position with an initial gap equal to 10 -2 m and the ring is launched with an initial angular velocity about the y-axis equal to 10rad • s -1 . In a first phase, the pendulum freely falls under gravity. When it touches the moving ring the ring describes some large oscillatory motions due to Coulomb's friction. We have an alternance of slip and stick phases with reversal sliding motions. The most interesting phenomenon for our purpose is the violation of the constraint. There are two causes of constraints violation with the standard Moreau-Jean scheme: 0.0e+00 The first cause is due to a first penetration when the pendulum hits the ring. We use in this example the explicit rule (61) which γ = 0 to activate the constraints at the velocity level. This first penetration is not recovered by the scheme. If the constraint is linear, this first penetration do not become more pronounced as the time passing because the relative velocity is set to zero by the impact law.

2. In our case, the contact constraint is nonlinear. The fact that we treat the constraints at the velocity level causes some drift from the constraints at the position level.

The result of the constraints violation is illustrated in Figure 20 where the contact distance is plotted with respect to time. We can observe the first penetration that is not too large due to the chosen reasonable time-step (h = 5 × 10 -4 s). From Proposition 1, we can infer that the penetration is of the same order of the time-step, that is O(h). We also see easily the result of the drift from the constraint when the pendulum slides on the ring. This drift increases as the contact point oscillates in the ring up to losing the contact when the pendulum reaches the external boundary of the ring. On the right side in Figure 20, we have depicted the gap as it is computed by the CAD library, which minimizes the distance between the contact point of the rod and the external or the internal boundary of the ring. The gap is positive when we are inside the ring and negative when the pendulum enter into the thickness of the ring. In this zone, the computed gap is decreasing in a first period and then increasing again since the contact is computed on the external boundary. Such type of troubles frequently occurs with thin structures if the penetrations are too large and not corrected.

In the case of the combined projection/activation algorithm, the penetration is kept less than the machine precision due to the projection. There is no drift in the constraints at the position level and therefore no loss of contacts.

Classical slider-crank mechanism

In this section, we are interested in the application of the proposed schemes to a classical mechanism : the slider-crank. The slider crank mechanism depicted in Figure 21 is composed of 3 mobile bodies: the crank ( 1 in light dray), the connecting rod ( 2 in dark gray) and the slider ( 3 in white). Concerning the clearances in joints, two configurations will be studied. In the first one, there is only play in the transitional joint between the slider and the cylinder. In the second configuration, a play is also introduced in the revolute joint between the crank and the connecting rod. In each configuration, the clearance is modeled by a unilateral contact with Coulomb's friction and Newton's impact law. The modeling of clearances is detailed in Figure 22. The size of the play between the crank and the connecting rod is denoted as r.

In the sequel, we will consider that the first configuration is obtained with a zero play, i.e., r = 0. This configuration is identical to those developed in [START_REF] Flores | Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach[END_REF] where the slider-crank is studied in the nonsmooth dynamics framework with a model based on unilateral contact, Coulomb's friction and Newton's impact law. The time integration in [START_REF] Flores | Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach[END_REF] is performed with the Moreau-Jean timestepping scheme. In a first step ( § 6.4.1), the equations of motion and the unilateral constraints are written in a pure Lagrangian setting with minimal coordinates chosen as in [START_REF] Flores | Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach[END_REF] : the generalized coordinates is defined by the crank angle, the connecting rod angle and the slider angle with respect to the x-axis. In the second step ( § 6.4.2), we use the full Newton/Euler formalism with maximal coordinates to facilitate the introduction of clearances in every joint of the mechanism.

In order to validate our approach and our model, the geometrical and mechanical properties for the system are rigorously the same as in [START_REF] Flores | Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach[END_REF]. They are listed in Table 2. The main discrepancy Geometrical properties l 1 = 0.1530m Table 2: Geometrical, mechanical and numerical parameters used in the slider-crank model with [START_REF] Flores | Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach[END_REF] is the adaptation of the geometry for the full three-dimensional case in § 6.4.2: the crank, the connecting rod and the slider are considered to be slender rods rather than planar laminates (see Figure 27 for an illustrtion). In their initial positions, the cylinders are aligned with the x-axis and the moments of inertia in Table 2 are given along their principal axis of inertia which coincides with the (x, y, z) frame at the initial time.

l 2 = 0.3060m a = 0.0500m b = 0.0250m c = 0.0010m d = 0.0520m Inertial properties m 1 = 0.0380kg m 2 = 0.0380kg m 2 = 0.0760kg J y,1 = J z,1 = 7.4 × 10 -5 kg • m 2 J x,1 = 1 × 10 -5 kg • m 2 J y,2 = J z,2 = 8.9 × 10 -4 kg • m 2 J x,2 = 1 × 10 -5 kg • m 2 J y,3 = J z,3 = 2.7 × 10 -6 kg • m 2 J x,
The system is under the action of the gravity acceleration equal to g = 9.81s • m -2 . The initial conditions are also chosen as in [START_REF] Flores | Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach[END_REF] that is the slider-crank is aligned with the x-axis and initial angular velocities are imposed to the crank ω 1,z = 150rad/s and the connecting rod ω 2,z = -75rad/s.

Lagrangian setting with no clearance between the crank and the connecting rod

As we recalled earlier in this section, we use the pure Lagrangian setting with minimal coordinates and without any play in the revolute joint. In Figure 23, the results of the simulation with the classical Moreau-Jean scheme are reproduced for the time-step h = 2 × 10 -6 over two crank revolutions. The results corroborates the results in [START_REF] Flores | Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach[END_REF]. In Figure 24, the motion of the corners of the slider are depicted. Note that the motion of the four corners are identical, due to the fact that the reaction forces results in a vanishing torque at the center mass. The motion of the slider is therefore a pure translation. We note also some violation of the constraints and drift of the constraints which are proportional to the time-step. This remark is confirmed by the numerical values reported in Table 3 where the maximum of violation of the constraints is of the same order as the time-step. It also corroborates the result of Proposition 1. In order to minimize the violation of constraints, the discrete frictional contact problem is solved at the velocity with a numerical tolerance of 10 -12 . With a smaller accuracy (a larger tolerance), the symmetry of the problem may be lost and a drift in the slider angle is observed leading to larger drifts in the constraints.

In Figure 26, the slider portrait (y position with respect to the x position of the slider's center of mass) is depicted for the various schemes : the Moreau-Jean scheme (18d) in Figure 25(a), the direct projected scheme (Algorithm 1) in Figure 25(b) and the combined scheme (Algorithm 2) in Figure 25(c). A quite large time-step h = 10 -04 s is used to exhibit large violations for the Moreau-Jean scheme. These violations can be observed in the kinks of the curve in Figure 25( 3 that the maximum of violations drops from 1.324 × 10 -04 to 8.099 × 10 -11 by using Algorithm 1. In Figure 26, we give some details on the relative velocity at the corner 1 of the slider for values of the crank revolution comprised in [0, 0.25]. We can observed in Figures 26(a) and 26(b) that the stabilization on the constraints with the Moreau-Jean scheme is smooth. This smoothness is lost with the direct projected scheme as we can observed in Figure 26(d) some oscillatory artifacts on the relative velocity when the contact should be closed. This is exactly the same phenomenon as we observed in Figure 4 for the bouncing ball example. In Figure 26(f), the effect of the combined scheme of Algorithm 2 is to retrieve a smooth stabilization of the constraints by keeping the satisfaction of the constraints at the position level.

In Table 3, we give some details on the computational efficiency of the projected algorithms for various time-steps. The first conclusion that can be drawn on this example is that the number of projection steps is quite low and it does not perturb the convergence of the Newton algorithm. On the contrary, we can observe that the convergence of the Newton algorihtm is improved by the projection onto the constraints. We can also observe that the average number of projection iterations with the combined scheme is lower than the average number of iterations for the projected scheme for the time-steps h = 10 -05 s and h = 10 -06 s. This is mainly due to the modification of the procedure for activating the constraints. In the combined scheme, the constraints are most the time activating before the contact is closed.

Newton/Euler setting with clearance between the crank and the connecting rod

In this section, the Newton/Euler equations are used for each body of the mechanism using a maximal set of coordinates (3 translations and a unit quaternion to parametrize the finite rotations). The formulation allows us to introduce some clearances in every joint of the mechanism without redefining a new choice of coordinates. In Figure 27, a view of the CAD is model is given where the arrow represents the reaction force in the revolute joint. The clearance is geometrically modelled by two cylinders that allows out-of-plane motions. With no play in the perfect revolute joint, i.e, the joint is modelled by an ideal revolute joint, the results, which are not reproduced here for a sake of space, perfectly corroborates the curves obtained in [START_REF] Flores | Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach[END_REF] and in the previous section (see Figures 23 and 24). In Figure 28, we report the result of the simulation with the standard Moreau-Jean timestepping scheme with a time-step h = 10 -05 s. A play of size r = 5 × 10 -04 m is introduced in the revolute joint with a perfectly plastic impact law. We can first observe in Figure 29(a), 29(b) and 29(c) that in a first phase, the angular speeds of the crank and the connecting are constant and equal to the prescribed values at the initial time. This reveals that the contact between the crank and the connecting rod is not active up to a first perfectly plastic impact occurs. We can also observe that the presence of clearances in the revolute joint introduces a higher frequency motion around the periodic motion of the slider. Most importantly, we observe a complete change of the periodic motion just before the end of the second crank revolution. This reveals that the contact in the revolute joint is lost during the simulation mainly due to too large constraints violations. During the motion, the CAD library is not able to follow correctly the contact point and the contact detection failed because the geometries interpenetrate. This problem can be fixed with a smaller time-step, for instance h = 10 -07 s that limits the violation of constraints. In Figure 29, the problem is solved with the same time-step h = 10 -05 s by using the combined projection algorithm described in Algorithm 2. In Figure 30, the simulation is performed with a one-order larger clearance in the revolute joint r = 5 × 10 -03 m and shows larger oscillations around the nominal motion.

Electrical circuit breaker's mechanisms

In this final section, we give an insight of the usefulness of the proposed projected schemes for the virtual prototyping of electrical circuit breakers designed by the company Schneider Electric. The model that we considered in a C60 circuit breaker, which is a domestic low voltage circuit breaker depicted on Figure 31. The mechanism is only composed of 7 moving bodies, but 12 contacts come into play when the breaker switches off. Furthermore, when the mechanism is in closed position, the equilibrium is guaranteed by means of Coulomb's friction in the contact of the two bodies described in Figure 31(c). When the breaker opens the circuit, a lot of events (impacts, stick-slip transitions, . . . ) are observed in experimental setups. A rather complete description of its behaviour and its nonsmooth modeling in 2D can be found in [START_REF] Abadie | Dynamic simulation of rigid bodies: Modelling of frictional contact[END_REF]. The study in 3D with clearances that is performed with the scheme described in Algorithm 2 allows us to accurately study the effect of the clearances in joints on the out-of-plane motion of the breaker. Furthermore, it helps to state on the stability and the robustness of the fundamental properties of the circuit breaker with respect to the size of the clearances. Such studies are not possible with standard event-driven schemes which have a lot of difficulties to deal with 3D frictional contacts and a bunch of events. Indeed, the presence of clearances in joints generates a lot of finite accumulation of impacts and numerous stick-slip transitions. Such studies are also difficult with standard eventcapturing schemes for which the violation of constraints come into play with the characteristic lengths of the clearances.

Concluding remarks

In this paper, an efficient strategy is proposed to perform the time-discretization of nonsmooth multibody systems that satisfies in discrete time both the constraints at the position level and at the velocity level. This strategy consists in

• a first direct project scheme which both satisfies in discrete time the position constraints and the velocity constraints, i.e., the impact law. This scheme based on the Moreau-Jean timestepping scheme is a direct extension of Gear-Gupta-Leimkuhler (GGL) method [START_REF] Gear | Automatic integration of Euler-Lagrange equations with constraints[END_REF] to unilateral constraints and impacts. The algorithm keeps the order of the standard Moreau-Jean scheme and the multiplier associated with the projection vanishes at the order O(h).

With the decoupled approach, the implementation is straightforward and requires only slight modifications of the standard Moreau-Jean scheme. This direct projected scheme efficiently performs on most applications. Nevertheless, in the special case of the stabilization on the constraints after a finite accumulation of impacts, the direct application of the GGL technique yields chattering at contact.

• an improved combined projection/activation scheme has been proposed to circumvent this problem and to make robust the simulation with the respect to the activation strategy of constraints. We end up with an event-capturing time-stepping scheme which retains the most favorable properties of the direct scheme (order, simple implementation and efficiency, respect of position and velocity constraints) avoiding the artificial oscillations at contact of the relative velocity. Further works must be done on the filtering of artificial high frequencies modes that occurs when a jump of velocity travels into a flexible structure discretized by finite elements. A first step has been done in [START_REF] Chen | A Newmark-Type Integrator for Flexible Systems Considering Nonsmooth Unilateral Constraints[END_REF] by adapting standard schemes for elastodynamics (HHT, α-schemes) that provides a targeted numerical damping. The question of shock wave capturing methods in such applications has also to be studied in the light of recent work on time-integration techniques for shock wave propagation (see [START_REF] Park | A method for computation of discontinuous wave propagation in heterogeneous solids: basic algorithm description and application to one-dimensional problems[END_REF] and references therein).

Example 3 (

 3 The rocking block) The rocking block of length L and thickness l is depicted in Figure1(c). Let us consider that the contact with the rigid ground can occur at the corner A and

Figure 1 :

 1 Figure 1: Simple archetypal test examples

Figure 2 :Figure 3 :

 23 Figure2: Empirical order of convergence of time-stepping schemes. Paoli-Schatzman scheme: I) uniform norm, IV) l 1 norm, VII) l 2 norm. Moreau-Jean's scheme: II) uniform norm, V) l 1 norm, VIII )l 2 norm. Direct projected scheme[START_REF] Jourdan | A Gauss Seidel like algorithm to solve frictional contact problems[END_REF]. III) uniform norm, VII) l 1 norm, IX) l 2 norm.

1 - 2 -

 12 Zoom on t ∈[3, 5] for γ = 0 Zoom on t ∈[3, 5] for γ = Zoom on t ∈[3, 5] for γ = 2

Figure 4 : 2 - 2 - 2 Figure 5 :

 42225 Figure 4: The chattering behavior for the bouncing ball (Example 1) with the explicit position forecast (61). Position and velocity vs. time for h = 5 × 10 -2

  Zoom on t ∈ [3, 5] Energy diagram with θ = 1/2. Zoom on t ∈ [3, 5] Energy diagram with θ = 1/2. Free position forecast (62).

Figure 6 :

 6 Figure 6: Energy for the bouncing ball (Example 1). h = 5 × 10 -2

Figure 7 : 2 Frictional contact µ, e < 1 RFigure 8 :

 7218 Figure 7: Energy in the elastic case (e = 1) for the bouncing ball (Example 1). h = 5 × 10 -2

Figure 9 :

 9 Figure 9: The Rocking Block (Example 3). Comparison of the direct projection scheme (Algorithm 1) and the combined activation/projection approach (Algorithm 2).

v 0 LFigure 10 :

 010 Figure 10: The elastic impacting bar problem

  Geometrical properties L = 1.0 m S = 3.14 × 10 04 m 2 Material properties ρ = 7800 kg.m -3 E = 210 GPa Initial Conditions v 0 = -0.1 m.s -1

Figure 11 :

 11 Figure 11: Trapezoidal rule with a position-based constraints. Position q, velocity v and reaction force r at the contact point. h = 2 × 10 -06 s. Number of elements N = 1000.

Figure 12 :

 12 Figure 12: Moreau-Jean scheme. Position q, velocity v and reaction force r at the contact point. θ = 1/2, h = 5 × 10 -05 s. Number of elements N = 1000.

  , the bar tip position and velocity are depicted. They are very similar to those obtained with the Moreau-Jean Number of elements N = 10 13: Contact force r for Moreau-Jean's scheme. θ = 1/2, h = 2 × 10 -06 s. Effect of the element length. Time-step h = 2 × 10 -05 s. Time-step h = 2 × 10 -06 s. 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60 1.80 2.00 N time [ms] r (c) Time-step h = 2 × 10 -07 s. Time-step h = 2 × 10 -08 s.

Figure 14 :

 14 Figure 14: Contact force r for Moreau-Jean's scheme. θ = 1/2. Number of elements N = 100. Effect of the time-step.

Figure 15 :

 15 Figure 15: Moreau-Jean's scheme with a restitution coefficient e = 0.95. h = 2 × 10 -06 s. Number of elements N = 100.

Figure 18 :

 18 Figure 18: Combined scheme (Algorithm 2). Position q, velocity v and reaction force r at the contact point. θ = 1, h = 2 × 10 -07 s. Number of elements N = 1000.

Figure 19 :Figure 20 :

 1920 Figure 19: Pendulum in a ring

Figure 21 :

 21 Figure 21: Classical slider-crank mechanism with clearances . 1. the crank, 2. the connecting rod and 3. the slider. Clearances in the translational joint and in the revolute joint between the crank and the connecting rod

Figure 22 :

 22 Figure 22: Details of the clearance modeling

2 = 1 ×

 21 10 -5 kg • m 2 Contact parameters e = 0.4, µ = 0.01 (slider/cylinder contact) e = 0.0, µ = 0.01 (crank/connecting rod joint) Numerical parameters θ = 1/2, Newton tolerance 10 -10 , violation tolerance 10 -10

Figure 23 :

 23 Figure 23: Details of the motion of the slider crank with the standard Moreau-Jean time-stepping scheme for h = 10 -6 s

Figure 24 : 1 -Figure 25 : 1 -

 241251 Figure 24: Dimensionless motion of the slider corners with the standard Moreau-Jean timestepping scheme for different time-steps

Figure 26 :

 26 Figure 26: Details on the relative velocity vs. the crank revolution at corner 1 with different time-stepping schemes. h = 10 -4 s

Figure 27 :

 27 Figure 27: CAD view of the slider-crank mechanism using OpenCascade [57] and PythonOCC [43].

1 -

 1 (rad/s) connecting-rod angle (rad) (c) connecting-rod portrait-

Figure 28 :

 28 Figure 28: Details of the motion of the slider crank with the standard Moreau's time-stepping scheme for h = 10 -05 s and a play in the revolute joint r = 5 × 10 -04 m

Figure 29 :

 29 Figure 29: Details of the motion of the slider crank with the combined projection time-stepping scheme for h = 10 -05 s and a play in the revolute joint r = 5 × 10 -04 m

Figure 30 :

 30 Figure 30: Details of the motion of the slider crank with the combined projection time-stepping scheme for h = 10 -05 s and a play in the revolute joint r = 5 × 10 -03 m

Figure 31 :

 31 Figure 31: C60 electrical circuit breaker mechanism. Courtesy of Schneider Electric

  The efficiency and the robustness of the schemes have been shown on several academic examples and illustrates on an industrial applications. The schemes are freely available in the open-source platform Siconos with the examples developed in this article.

Table 1 :

 1 Geometrical and material parameters and initial conditions for the impacting bar example

	m	-2.0e-05 0.0e+00 2.0e-05 4.0e-05 1.2e-04 1.0e-04 8.0e-05 6.0e-05

  Figure 16: Combined scheme (Algorithm 2). Position q, velocity v and reaction force r at the contact point. = 1/2, h = 2 × 10 -06 s. Number of elements N = 1000.Figure 17: Combined scheme (Algorithm 2). Position q, velocity v and reaction force r at the contact point. θ = 1/2, h = 2 × 10 -07 s. Number of elements N = 1000.

	m m/s N m m/s N	2.0e-05 4.0e-05 6.0e-05 8.0e-05 1.0e-04 1.2e-04 -1.0e-01 -5.0e-02 0.0e+00 5.0e-02 1.0e-01 1.5e-01 2.0e-01 0.0e+00 5.0e+02 1.0e+03 1.5e+03 2.0e+03 2.5e+03 5.0e-05 8.0e-05 9.0e-05 1.0e-04 7.0e-05 6.0e-05 4.0e-05 3.0e-05 2.0e-05 1.0e-05 0.0e+00 -2.0e-02 0.0e+00 2.0e-02 6.0e-02 8.0e-02 1.0e-01 4.0e-02 -4.0e-02 -6.0e-02 -8.0e-02 -1.0e-01 0.0e+00 2.0e+02 4.0e+02 6.0e+02 8.0e+02 1.0e+03 1.2e+03 1.4e+03	q v 0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60 1.80 2.00 time [ms] time [ms] r q v 0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60 1.80 2.00 r
	m m/s N	1.2e-04 1.0e-04 8.0e-05 6.0e-05 4.0e-05 2.0e-05 0.0e+00 5.0e-02 2.0e-01 1.5e-01 1.0e-01 0.0e+00 -5.0e-02 -1.0e-01 0.0e+00 2.0e+02 4.0e+02 6.0e+02 8.0e+02 1.0e+03 1.2e+03 1.4e+03	q v 0.00 0.20 0.40 0.60 0.80 1.00 1.20 1.40 1.60 1.80 2.00 r
			time [ms]

Table 3 :

 3 Details on the computational efficiency of the projection algorithms in the Lagrangian setting

	-6 s
	Method	time-step	violation	Newton iteration projection iteration index sets iteration
		[s]	[m] (max.)	(avg./max.)	(avg./max.)	(avg./max.)
	Moreau-Jean (18d)	10 -04	1.324 × 10 -04	2.95/4	N/A	N/A
	Moreau-Jean (18d)	10 -05	1.234 × 10 -05	1.99/3	N/A	N/A
	Moreau-Jean (18d)	10 -06	1.119 × 10 -06	1.96/2	N/A	N/A
	Algorithm 1	10 -04	8.099 × 10 -11	1.96/2	0.09/2	N/A
	Algorithm 1	10 -05	4.833 × 10 -11	1.01/2	0.11/1	N/A
	Algorithm 1	10 -06	9.99 × 10 -11	1.001/2	0.04/1	N/A
	Algorithm 2	10 -04	8.410 × 10 -11	2.15/4	0.20/2	1.18/2
	Algorithm 2	10 -05	9.940 × 10 -11	1.127/2	0.019/1	1.127/2
	Algorithm 2	10 -06	8.650 × 10 -11	1.12/2	0.00036/1	1.12/2

RR n°8029

Acknowlegments

This work has been supported by the French National Research Agency (ANR) through COSINUS program (project SALADYN ANR-08-COSI-014)http://saladyn.inria.gforge.fr/. The author warmly acknowledges Olivier Bonnefon for the first implementation of the direct projected scheme in Siconos.