Multibody Dynamics with Unilateral Constraints: Computational Modelling of Soft Contact and Dry Friction

. We consider a system of rigid bodies subjected to unilateral constraints with soft contact and dry friction. When the constraints are saturated, velocity jumps may occur and the dynamics is described in generalized coordinates by a second-order measure diﬀerential inclusion for the unknown conﬁgurations. Observing that the right velocity obeys a minimization principle, a time-stepping algorithm is proposed. It allows to construct a sequence of approximate solutions satisfying at each time-step a discrete contact law which mimics the behaviour of the sys-tem in case of collision. In case of tangential contact, dry friction may lead to indeterminacies such as the famous Painlev´e’s paradoxes. By a precise study of the asymptotic properties of the scheme, it is shown that the limit of the approximate trajectories exhibits the same kind of indeterminacies.


Description of the Problem
We consider a discrete mechanical system with d degrees of freedom.We denote by q ∈ R d its representative point in generalized coordinates and by M (q) its inertia operator.We assume that the system is subjected to unilateral constraints characterized by the geometrical inequality g(q) ≤ 0 (non penetration condition) with a smooth (at least C 1 ) function g such that ∇g does not vanish in a neighborhood of q ∈ R d ; g(q) = 0 .Let us denote by •, • the Euclidean inner product in R d .If for some instant t > 0 the constraints are saturated, i.e. g q(t) = 0, then q− (t), ∇g q(t) ≥ 0, q+ (t), ∇g q(t) ≤ 0 and the velocity may be discontinuous.It follows that u = q is a function of Bounded Variations and the dynamics is described by the Measure Differential Equation M (q)du = f (t, q, u)dt + r where du is the Stieljes measure associated to u and r is the reaction due to the constraints.Of course a reaction is applied to the system only when a contact occurs and we have a complementarity condition g q(t) < 0 =⇒ r = 0.
Furthermore, we assume that the contact is non-adhesive which yields r, ∇g q(t) ≤ 0 if g q(t) = 0.
In the frictionless case we get r ∈ −N q(t) if g q(t) = 0 where N q(t) is the normal cone to the set of admissible configurations at q(t) and in the frictional case where C q(t) is the so-called friction cone at q(t).Hence we may rewrite (1) as a Measure Differential Inclusion ( [17], [10]) where and n(q) = − ∇g(q) ∇g(q) and D 1 (q) is the disc of center 0 and radius µ in Rn(q) ⊥ with µ = 0 (frictionless constraints) or µ > 0 (Coulomb's friction).
Let us assume also soft contact i.e.
We infer that In the frictionless case we may decompose u ± (t) as u ± (t), n q(t) n q(t) , M −1 q(t) n q(t) .
With ( 3)-( 4) we get where Proj q(t) denotes the projection relatively to the kinetic metric at q(t).
In the frictional case, when M (q) ≡ mId R d , m > 0, Coulomb's law ( [3], [11], [12]) yields where r N and r T are respectively the projection relatively to the Euclidean metric of r on Rn q(t) and T q(t) and ∂ψ r N D1(q(t)) is the indicatrix function of the disc of radius µr N and center 0 in T q(t) .Reminding that R q(t) = R + n q(t) + D 1 q(t) , ( 5) is equivalent to which can be rewritten as Let us emphasize that (4) and (6) imply that u + (t) = u − (t) only if g q(t) = 0 and u − (t), n q(t) < 0, i.e.only in case of collision.Moreover u + (t) is defined as the Argmin of the kinematically admissible right velocities.The same property holds when M (q) ≡ mId R d and µ > 0 and we still have ( [12], [6]) if g q(t) = 0 and u − (t), n q(t) < 0. On the contrary, when g q(t) = 0 and u − (t), n q(t) = 0, (4) yields and velocity jumps without collision may occur if M −1 q(t) R q(t) ∩T q(t) = {0}.Such phenomena can easily be observed when we consider the model problem of a slender rod in contact at one edge with an horizontal obstacle, leading to the famous Painlevé's paradoxes ( [13], [14]): there exists a subset A q(t), u − (t) , containing u − (t) but not reduced to this single point, such that any value of u + (t) ∈ A q(t), u − (t) solves the problem (see for instance [8], [4], [1] or more recently [12], [2], [7]).Hence, in case of tangential contact with dry friction and non-trivial inertia operator, the dynamics exhibits indeterminacies.

Computational Modelling: the Contact Dynamics Approach
In order to solve numerically the problem, the Contact Dynamics approach has been introduced by J.J. Moreau in the mid 80's ( [10], [11], [12]).The core idea is to avoid any regularization of the unilateral constraints and to build a timestepping scheme by combining an Euler discretization of the measure differential inclusion (2) on each interval [t i , t i+1 ] with an impulsional form of the contact law at t i+1 .More precisley the approximate position is updated as and a "free" left velocity at t i+1 is defined by Then u i+1 is the right velocity at t i+1 given by and where S is a discrete analogous of the contact law.Starting from the definition of R(q i+1 ), we get immediately If g(q i+1 ) ≥ 0 and v i+1 , n(q i+1 ) > 0, the left velocity v i+1 points inward and If g(q i+1 ) ≥ 0 and v i+1 , n(q i+1 ) < 0, we may interpret t i+1 as a collision instant and with (6) we get Finally, if g(q i+1 ) ≥ 0 and v i+1 , n(q i+1 ) = 0, we get Hence u i+1 is uniquely defined if M −1 (q i+1 )R(q i+1 )∩T (q i+1 ) = {0} and we have u i+1 = S(q i+1 , v i+1 ) = v i+1 .On the contrary, if M −1 (q i+1 )R(q i+1 ) ∩ T (q i+1 ) = {0}, t i+1 may be interpreted as a discrete tangential contact with possible indeterminacies and there is not any natural choice for u i+1 .
Otherwise, v i+1 ∈ V (q i+1 ) and v i+1 , n(q i+1 ) < 0. By definition of S we get Observing that the inequality u i+1 , r i+1 ≤ 0 corresponds to a dissipativity property, we infer that the scheme reproduces at the discrete level the main features of the dynamics except the indeterminacies of Coulomb's law.Indeed, S(q, v) = v if g(q) = 0, v, n(q) = 0 and the discrete contact law does not lead to any velocity jump in case of tangential contact.
Starting from Proposition 1, the convergence of the scheme has been established by M. Monteiro-Marques in [9] when M (q) ≡ mId R d (m > 0) and µ ≥ 0 and by R. Dzonou and M. Monteiro-Marques in [5] when M (q) ≡ mId R d and µ = 0.In both cases we have M −1 (q)R(q) ∩ T (q) = {0} for all q ∈ R d such that g(q) = 0 and the difficulty due to indeterminacies is avoided.Nevertheless the convergence has also been proved recently when M (q) ≡ mId R d and µ > 0 ( [15]), so a natural question arises: is it possible to recover with such a scheme velocity jumps without collisions at the limit when h tends to zero?A first answer has been given by J.J. Moreau in [12]: numerical simulations show that the approximated trajectories exhibit the plurality of solutions given by Coulomb's law.

Asymptotic Properties of the Discrete Contact Law
Let us assume from now on that M (q) ≡ mId R d (m > 0) and µ > 0. Then the limit trajectory will satisfy in case of convergence the following property We may recover the indeterminacies of Coulomb's law if, for all (q, u − ) such that g(q) = 0, u − , n(q) = 0 and M −1 (q)R(q) ∩ T (q) = {0}, we have with Let us recall that the Hausdorff distance between two subsets A and B of R d is defined as Hence ( 8) can be decomposed as which can be interpreted as an asymptotic consistency property of the discrete contact law S and lim ε→0 sup dist v, A ε (q, u − ) ; v ∈ A(q, u − ) = 0 (10) which can be interpreted as an asymptotic indeterminacy of the scheme.
In the one-dimensional friction case, i.e. when Dim Span D 1 (q) = 1, then (see [12]).Then we can prove that ( 9) is satisfied while (10) is not always true and depends on the evolution of the mappings q ε → R(q ε ) and q ε → n(q ε ) in a neighborhood of the contact point.More precisely let us assume that (H1) the mapping M is of class C 1 from R d to the set of symmetric positive definite d × d matrices; (H2) the function g belongs to C 1 (R d ), ∇g is locally Lipschitz continuous and does not vanish in a neighbourhood of {q ∈ R d ; g(q) = 0}; (H3) for all q ∈ R d , D 1 (q) is a closed, bounded, convex subset of R d such that 0 ∈ D 1 (q) and the multivalued mapping q → D 1 (q) is Hausdorff continuous.Furthermore, ∇g(q) ∈ Span D 1 (q) for all q ∈ R d such that ∇g(q) = 0. We denote by K the set of admissible configurations i.e.
With assumption (H2) there exists r q > 0 such that the mapping n : is well defined and Lipschitz continuous.Let us assume moreover that, possibly reducing r q , we have (H4) dim Span D 1 (q ) = 1 for all q ∈ B(q, r q ).
We infer that the discrete contact law always satisfies the asymptotic consistency property while the asymptotic indeterminacy of the scheme holds only if γ(q) < 0 or γ(q) = 0 and for all ε ∈ (0, r q ) there exists q ε ∈ B(q, ε)\ Int(K)∪{q} such that γ(q ε ) > 0 if ũ = u − .In the latter case, any v ∈ [u − , ũ] \ {u − , ũ} is the limit of a sequence of post-collision velocities (S(q εn , u − εn ) = ũεn ) n≥1 with (ε n ) n≥1 decreasing to zero and (q εn , u − εn ) ∈ B(q, ε n ) × B(u − , ε n ) for all n ≥ 1.Hence, for all n ≥ 1, S(q εn , u − εn ) = ũεn is defined as the Argmin of the kinetic norm of the admissible right velocities at (q εn , u − εn ) but v is not the Argmin of the kinetic norm of the admissible right velocities at (q, u − ).It means that the minimization property (7) defining post-collision velocities is not continuous at (q, u − ) and it appears that this mathematical propery is deeply related to the indeterminacies of Coulomb's law.
Finally let us emphasize that (H3) allows to take into account both isotropic and anisotropic friction.Moreover the conclusions of Theorem 1 are still valid when (H4) is replaced by (H'4) D 1 (q ) is strictly convex for any q ∈ B(q, r q ) i.e. for any w 1 and w 2 belonging to D 1 (q ) such that w 1 = w 2 , and for any γ ∈ (0, 1), γw 1 + (1 − γ)w 2 belongs to the relative interior of D 1 (q ), i.e. there exists a open subset O of R d such that which is always true for the classical isotropic Coulomb's friction characterized by a friction coefficient µ > 0.
B), e(B, A) where e(A, B) = sup a∈A dist(a, B) = sup a∈A inf b∈B a − b (the excess of A from B), e(B, A) = sup b∈B dist(b, A) = sup b∈B inf a∈A b − a (the excess of B from A).