A proximal-like algorithm for vibro-impact problems with a non-smooth set of constraints

We consider a discrete mechanical system subjected to perfect unilateral contraints characterized by some geometrical inequalities $f_\alpha(q) >=0$, $\alpha \in {1,...,v}$, with $v >=1$. We assume that the transmission of the velocities 
at impacts is governed by a Newton’s impact law with a restitution coefficient 
$e \in [0, 1]$, allowing for conservation of kinetic energy if $e=1$, or loss of kinetic 
energy if $e \in [0, 1)$, when the constraints are saturated. Starting from a formulation of the dynamics as a first order measure-differential inclusion for the 
unknown velocities, time-stepping schemes inspired by the proximal methods 
can be proposed. Convergence results in the single-constraint case $(v = 1)$ 
are recalled and extended to the multi-constraint case $(v > 1)$, leading to new 
existence results for this kind of problems.

1. Description of the dynamics. We consider a mechanical system with a finite number of degrees of freedom. We denote by q ∈ R d the representative point of the system in generalized coordinates and by M (q) the inertia operator. The unconstrained dynamics is described by a second order Ordinary Differential Equation M (q)q = g(t, q,q). (1) We assume that the system is subjected to unilateral constraints described by some geometrical inequalities f α q(t) ≥ 0 ∀t, ∀α ∈ {1, . . . , ν}, ν ≥ 1 (2) with smooth (at least of class C 1 ) functions f α . The set of admissible configurations is then defined by . . , ν} and we denote by J(q) the set of active constraints at q given by J(q) = α ∈ {1, . . . , ν}; f α (q) ≤ 0 , ∀q ∈ R d .
When the constraints are saturated, i.e. when at least one of the inequalities (2) is an equality, a reaction force appears and should be added to the right hand side of equation (1): M (q)q = g(t, q,q) + R, Supp(R) ⊂ t; J q(t) = ∅ .
We assume moreover that the constraints are perfect i.e.

1
• there is no adhesion (R, v) ≥ 0 ∀v ∈ T K (q) • the contact is frictionless where (·, ·) denotes the Euclidean inner product in R d and T K (q) is the set of kinematically admissible right velocities at q given by Using Farkas's lemma we infer that If J q(t) = ∅ the inequalities (2) implẏ Thus the velocities may be discontinuous at impacts. It follows that R is a measure and (3)-(4) lead to a second order Measure Differential Inclusion We infer that the jumps of velocities at impacts satisfy but relations (5) and (7) do not define uniquelyq + (t) and we have to complete the description of the dynamics by a constitutive impact law. We will assume that the tangential part of the left velocity is conserved while its normal part is reversed and multiplied by a coefficient e ∈ [0, 1], i.e.
are defined respectively as the projection ofq − (t) on the convex cones T K q(t) and N * K q(t) = M −1 q(t) N K q(t) relatively to the kinetic metric at q(t). More precisely we define the kinetic metric at q by which is a Newton's law with a restitution coefficient e ∈ [0, 1]. When e = 0 we simply getq and we recognize the description of standard inelastic shocks introduced by J.J. Moreau ([9]). When e = 0, we can use the lemma of the two cones (see [8]) to rewrite (8) aṡ We can observe that this model is mechanically consistent. Indeed the kinetic and we have conservation of energy at impacts if e = 1 (elastic shocks).
Following J.J.Moreau's approach (see [9] or [10] for instance) we can describe the dynamics at the velocity level by the following first order Measure Differential Inclusion where ψ TK (q) denotes the indicator function of T K (q) and ∂ψ TK (q) its subdifferential given by For a more detailed study of the equivalence of the system (6)-(8) with the MDI (9) the reader is referred to [13].
More precisely, for any given admissible initial data (q 0 , u 0 ) ∈ K × T K (q 0 ) we will consider the Cauchy problem (P): (P4) there exists a non-negative measure µ such that the Stieltjes measure du and the Lebesgue's measure dt admit densities relatively to dµ, denoted respectively u µ and t µ , and dµ a.e. on (0, τ ), 2. Time-stepping scheme. Let h > 0 be a given time-step and (q 0 , u 0 ) be any given admissible initial data. Starting from (9) the discrete positions and velocities are defined by and and u h,i as the right and left velocities at t h,i+1 = (i + 1)h, this inclusion is a very natural discretization of the MDI (9). Moreover, using the definition of ∂ψ TK (q) , we can rewrite it as which leads to an approximation of the impact law (8).
This scheme has been introduced by J.J.Moreau in the 80's (see [10] for instance). It is inspired by sweeping process techniques and can be interpreted as a proximallike method for the MDI (9). In the single constraint case (i.e. ν = 1) its convergence has been established first in the case of a trivial mass matrix and standard inelastic shocks (i.e. M (q) ≡ Id R d and e = 0) by M.Monteiro Marques ([6], [7]). This result has been extended to the case of partially or totally elastic shocks (i.e. e ∈ [0, 1]) but still a trivial mass matrix by M. Mabrouk ([4]).
For a non-trivial mass matrix and a vanishing restitution coefficient, the convergence of the scheme has been established by B.Maury for a constant inertia operator M and by R.Dzonou and M.Monteiro Marques for a position-dependent inertia operator M = M (q), both in 2006 ( [5] and [2]). Finally, the general single constraint case (i.e. M = M (q) ≡ Id R d and e ∈ [0, 1]) has been considered in [3].
Motivated by applications to systems of rigid bodies (see for instance [11] or [12] for examples of implementation with granular materials), we can wonder if the previous convergence results can be extended to the multi-constraint case. Unfortunately we meet a new difficulty: continuity on data does not hold in general when ν > 1.
Indeed, let us consider the model problem of a material point moving in an angular domain of R 2 : it can be seen easily that continuity on initial data does not hold if the edge angle is obtuse (see figure 1 and [14] for detailed computations).
In such a case, even if we can prove a theoretical convergence result for the algorithm (10)-(11), round-up errors may lead to a kind of unpredictibility. So we have to introduce some geometrical assumptions on the active constraints to ensure continuity on data and avoid this difficulty. Keeping in mind the previous model problem, it has been proved that continuity on data holds if for all q ∈ ∂K for all (α, β) ∈ J(q) 2 such that α = β (see [1] and [14]). These relations mean that the active constraints create right or acute angles with respect to the momentum metric defined by M −1 (q) if e = 0 or right angles if e = 0 and, in this framework, we will establish the convergence of the approximate trajectories. with an equality when the constraints are never saturated and finite time explosion may occur. Nevertheless, we can establish that Proposition 3.1. ( [15], [16]) Let C > u 0 M(q0) . Then, there exists τ (C) ∈ (0, T ] such that, for any solution (q, u) of problem (P) defined on [0, τ ], we have q(t) − q 0 ≤ C ∀t ∈ 0, min τ (C), τ , u(t) M(q(t)) ≤ C dt a.e. on 0, min τ (C), τ .

and (q, u) is a solution of problem (P).
Let us emphasize that this convergence result provides an existence result for problem (P) under weaker regularity assumptions on the data than in [1]. Since the Cauchy problem may have several solutions (see [17] or [1] for counter-examples to uniqueness) we can not expect the convergence of the whole sequence of approximate solutions (q h , u h ) h>0 .
Sketch of the proof The proof is divided in four steps.
Step 1 We begin with an estimate of the discrete velocities.
With the lemma of the two cones we can check easily that

By definition of the algorithm
Thus, It follows that there exists τ ∈ (0, T ] such that the sequence (u h ) h>0 is uniformly bounded in L ∞ (0, τ ; R d ) and (q h ) h>0 is uniformly Lipschitz continuous on [0, τ ]. Using Ascoli's theorem, and possibly extracting a subsequence denoted (q hn , u hn ) n∈N , we get Furthermore, observing that the constraints are satisfied at each time step at the velocity level by the average discrete velocity u h,i+1 + eu h,i 1 + e (which belongs to T K (q h,i+1 )), we can prove that Thus there exists a compact set B such that q(t) ∈ K ∩ B for all t ∈ [0, τ ]. We infer that there exists a compact neighbourhood K of K ∩ B such that (H3) holds on K and q hn (t) ∈ K for all t ∈ [0, τ ] and for all n large enough. As a consequence we can prove Lemma 3.3. For all i ∈ 0, . . . , τ /h n , there exist non-positive real numbers µ α hn,i+1 ) α∈J(q hn ,i+1 ) such that and there exists a constant C 1 (independent of n and i) such that |µ α hn,i+1 | ≤ C 1 .
Step 2 Next we prove an a priori estimate for the discrete accelerations. More precisely Indeed, with lemma 3.3, we only need to estimate τ /hn i=1 α∈J(q hn ,i ) |µ α hn,i |. But, for all i ∈ 1, . . . , τ /h n , the family ∇f α (q hn,i+1 ) α∈J(q hn ,i+1 ) is linearly independent, so we can complete it as a basis of R d denoted by v j (q hn,i+1 ) 1≤j≤d and we define by w j (q hn,i+1 ) 1≤j≤d its dual basis. Then we introduce µ β hn,i = 0 if β ∈ J(q hn,i ) and we have We rewrite it in order to get a telescopic sum which allows us to conclude.
We infer that the sequence (u hn ) n∈N is uniformly bounded in BV (0, τ ; R d ). Thus we can pass to the limit with Helly's theorem: possibly extracting another subsequence, still denoted (q hn , u hn ) n∈N and possibly modifying v on a negligible subset of [0, τ ] we obtain and v ∈ BV (0, τ ; R d ). We define u ∈ BV (0, τ ; R d ) by Clearly u satisfies properties (P1)-(P2) and we have Furthermore, using the same techniques as in [3] we prove that (P4) holds at the continuity points of u and that u + (0) = u 0 .
Step 3 We study now the jumps of the velocities. First we observe that Next we consider t ∈ (0, τ ) such that u is discontinuous at With the previous results, we already know thaṫ It follows that J q(t) = ∅ and there exist non-positive real numbers (µ α ) α∈J(q(t)) such that µ α ∇f α q(t) .

Hence (13) is satisfied if and only if
v and for all α ∈ J q(t) . Recalling lemma 3.3, we have with µ α hn,i+1 ≤ 0 for all i ∈ 0, . . . , τ /h n − 1 , for all n ∈ N. So, if µ α = 0, we can prove with a contradiction argument that, in any neighbourhood V of the impact instant t, there exists at least one discrete impact t hn,i+1 such that α ∈ J(q hn,i+1 ) and µ α hn,i+1 < 0 for all n large enough. With the definition of the scheme we get also u hn,i+1 + eu hn,i 1 + e = Proj M(q hn ,i+1 ) T K (q hn,i+1 ), u hn,i + h n 1 + e g hn,i+1 .
We distinguish now two cases.
So, if µ α = 0, for any neighbourhood V of the impact instant t, we may consider the last discrete impact in V: we have ∇f α (q hn,i+1 ), u hn,i+1 = 0.
Using assumption (12) we obtain where i + = max i ∈ N; ih n ∈ V . So we can pass to the limit as n tends to +∞ first, then as |V| tends to zero, and we finally obtain ∇f α q(t) , v + (t) ≤ 0. But v + (t) ∈ T K q(t) , thus ∇f α q(t) , v + (t) ≥ 0 and we may conclude.
Once again, for any given neighbourhood V of t, we know that there exists at least one discrete impact t hn,i+1 ∈ V for all n large enough. It follows that ∇f α (q hn,i+1 ), u hn,i+1 + eu hn,i = 0.
Step 4 Let C > u 0 M(q0) . With the previous steps of the proof we already know that the convergence holds on a non-trivial time interval [0, τ ], with τ ∈ (0, T ] and it remains to extend it to 0, τ (C) . With proposition 3.1 we get lim n→+∞ u hn (t) M(q hn (t)) = u(t) M(q(t) ≤ C dt a.e. on 0, min τ (C), τ .
Then we use a contradiction argument to conclude.