Keywords: timestepping scheme, discontinuous Galerkin method, nonsmooth dynamics, flexible multibody system, friction, impact, index reduction Schémas d'intégration en temps, méthodes de Galerkin discontinues, dynamique nonrégulière, frottement, impact, réduction d'index

. It is based on time-discontinuous Galerkin methods to carry over higher order trial functions of event-driven integration schemes to consistent timestepping schemes for nonsmooth dynamical systems with friction and impacts. Splitting separates the portion of impulsive contact forces from the portion of non-impulsive contact forces. Impacts are included within the discontinuity of the piecewise continuous trial functions, i.e., with first-order accuracy. Non-impulsive contact forces are integrated with respect to the local order of the trial functions. In order to satisfy the constraints, a set of nonsmooth equations has to be solved in each time step depending on the number of stages; the solution of the velocity jump together with the corresponding impulse yields another nonsmooth equation. All nonsmooth equations are treated separately by semi-smooth Newton methods. The integration scheme on acceleration level was first introduced in [51] labeled "forecasting trapezoidal rule". It was analyzed and applied to a decoupled bouncing ball example concerning principal suitability without taking friction into account. In this work, the approach is algorithmically specified, improved and applied to nonlinear multi-contact examples with friction. It is compared to other numerical schemes and it is shown that the newly proposed integration scheme yields a unified behavior for the description of contact mechanical problems.

Introduction

In this paper, we study numerical integration schemes for the simulation of nonsmooth mechanical systems. Rigid and flexible bodies, joints as well as contacts and impacts with dry friction constitute the mechanical models. Thereby, we formulate the contact conditions as constraints and do not allow any penetration, e.g. due to penalty techniques. As a result, velocity jumps occur during the transient simulation of semi-discrete models and we have to be cautious in the formulation of efficient and stable time-discretization schemes.

We distinguish two cases.

Inria 1. non-impulsive contact forces -For the contact between flexible bodies, the contact force is finite in continuum models, i.e., in not semi-discrete models, although velocity jumps may arise. Classic (implicit) time integration schemes for computational mechanical problems, i.e., members of the Newmark family [START_REF] Newmark | A method for computation for structural dynamics[END_REF][START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF][START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-a method[END_REF][START_REF] Simeon | Computational Flexible Multibody Dynamics[END_REF][START_REF] Bauchau | Flexible Multibody Dynamics[END_REF], have been adapted to these demands and extended with respect to contact/velocity updates (Laursen-Love scheme) [START_REF] Laursen | Computational contact and impact mechanics[END_REF][START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF]. Another strategy to preferably get a well-posed problem is the application of energy-momentum paradigms like in energy-momentum schemes [START_REF] Simo | The discrete energy-momentum method. conserving algorithms for nonlinear elastodynamics[END_REF][START_REF] Bauchau | Flexible Multibody Dynamics[END_REF], i.e., modifications of the midpoint rule [START_REF] Simeon | Computational Flexible Multibody Dynamics[END_REF], to impact problems (Laursen-Chawla scheme) [START_REF] Laursen | Computational contact and impact mechanics[END_REF]. A contact-stabilized Newmark scheme is proposed in [START_REF] Deuflhard | A contact-stabilized newmark method for dynamical contact problems[END_REF].

2. impulsive contact forces -For the contact between rigid bodies, the reaction forces are impulsive and the classic time integration schemes do not work anymore [START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF]. The application of mass redistribution techniques [START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF] is a procedure, which reminds of penalty approaches with the benefit of having a theoretical foundation [START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF]. However in [START_REF] Krause | Presentation and comparison of selected algorithms for dynamic contact based on the newmark scheme[END_REF] it is shown, that all these schemes suffer from oscillations in the relative contact velocities. Event-driven schemes and timestepping schemes are further concepts to simulate rigid multibody systems or semi-discrete systems consistently by applying impact laws. Thereby, event-driven schemes resolve impact events to a high precision.

In-between impact events, standard integration schemes are used. Classic timestepping schemes do not resolve impact events, but include their possible existence directly in the discretization. Thus, they have low accuracy in non-impulsive intervals, whereas event-driven schemes may get inefficient and inconsistent for many impact events [START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF].

The aim of the present paper is to improve the consistent and robust concept of timestepping schemes for semi-discrete mechanical systems. The main drawback is the lack of problem adaptive accuracy in non-impulsive intervals. Classic timestepping schemes can be embedded within the context of timediscontinuous Galerkin methods, when we choose piecewise constant trial functions for the velocity approximation [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook[END_REF]. In the aforementioned paper, two different families of timestepping schemes on acceleration level based on discontinuous Galerkin methods have been introduced and analyzed. They differ in the interpretation, if one assumes the velocity jump at the beginning, called D + timestepping schemes, or at the end, called D timestepping schemes, of each discretization interval. First, using higher-order but piecewise continuous trial functions for the velocity, and, second, splitting of impulsive and nonimpulsive contact reactions, offer the opportunity to both stay consistent and benefit from a higher-order integration of non-impulsive contact reactions. The "forecasting trapezoidal rule" is the D representative for piecewise linear velocity trial functions; concerning function evaluations and implementation complexity it is the easiest scheme. Its basic practical applicability apart from theoretical propositions has been shown with a decoupled bouncing ball example. It is the basis for extensions in the present paper. We algorithmically specify and improve it and apply it to nonlinear multi-contact examples with friction.

Thereby for an automatic switching between non-impulsive and impulsive intervals, constraint equations for non-impulsive reactions are also formulated on velocity level in an augmented Lagrangian setting [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF] like the constraint equations for contact impulse and jumping velocity. Therefore, we focus on halfexplicit methods and evaluate all other magnitudes but the constraints explicitly [START_REF] Brasey | Hem 5 user's guide[END_REF][START_REF] Brasey | Methodes demi-explicites pour equations differentielles algebraiques semi-explicites d'index 2[END_REF][START_REF] Murua | Phem 56 user's guide[END_REF][START_REF] Murua | Partitioned half-explicit runge-kutta methods for differential-algebraicsystems of index 2[END_REF][START_REF] Arnold | Non-stiff integrators for differential algebraic systems of index 2[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic-Problems[END_REF].

The number of nonsmooth equations for the constraints depends on the number of stages, i.e., the local order of the underlying trial functions. These three nonsmooth equations are solved separately by semi-smooth Newton methods [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF][START_REF] Christensen | Formulation and comparison of algorithms for frictional contact problems[END_REF][START_REF] Qi | A survey of some nonsmooth equations and smoothing newton methods[END_REF].

The recent proposal for higher-order integration of non-impulsive contact reactions on basis of the generalized-a scheme [START_REF] Brüls | Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized-a scheme[END_REF] is also footing on a specific splitting approach. Splitting is interpreted as separating the "flow" of applied forces, bilateral contact reactions and unilateral contact reactions in the classic "Strang-sense" [START_REF] Mclachlan | Splitting methods[END_REF]. In contrast, we split non-impulsive and impulsive reactions not depending on their source. In particular, closed, i.e., active unilateral contacts, are treated with first-order accuracy in [START_REF] Brüls | Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized-a scheme[END_REF] but benefit from a higher-order approximation in the present approach.

The paper introduces the new time-discretization scheme step-by-step. The proposed (intermedi-ate) schemes are directly evaluated using an impacting slider-crank mechanism which we introduce in Sect. 2.1. The paper is organized as follows. Section 2 introduces the equations of motion and common time-discretization schemes for nonsmooth mechanical systems. Half-explicit timestepping schemes on acceleration level are introduced in Sect. 3. They generalize the approach in [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook[END_REF] algorithmically but suffer from drift-off effects. In Sect. 4, a first attempt for a velocity-level approximation is discussed to reduce the drift-off effect. It follows the constraint evaluation ideas of the half-explicit timestepping schemes on acceleration level and yields an explicit but not half-explicit scheme on velocity level, with even increased drift-off effect. The final half-explicit scheme on velocity level is derived and described in Sect. 5. It is tested within the example of the slider-crank mechanism. In Sect. [START_REF] Arnold | Numerical methods for simulation in applied dynamics[END_REF], we analyze the performance of our proposed timestepping scheme with the example of the slider-crank mechanism and a bouncing ball example. In Sect. 7, we extend this analysis to flexible systems with the examples of an elastic bar and a rotor. Section 8 concludes the paper.

Equations of motion and time-discretization

This section introduces the continuous equations of motion, the impacting slider-crank mechanism as an accompanying example and basic time-discretization schemes. We consider an impacting mechanical system with dry friction. For non-impulsive motion in time, it is described by the following set of equations

q = v , (1)
and

Mv = h +W N l N +W T l T , (2)
0  g N ? l N 0 , (3)
( kl T k  µl N for ġT = 0 ^gN  0 l T = ġT kġ T k µl N for ġT 6 = 0 ^gN  0 . ( 4
)
For mechanical systems, these equations typically occur after semi-discretization in space of respective mathematical models, e.g. by the finite element method [START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF]. If one does not use any classic concept of regularization by introducing contact stiffnesses or contact potentials [START_REF] Wriggers | Inovallée 655 avenue de l'Europe Montbonnot 38334 Saint Ismier[END_REF][START_REF] Fetecau | Nonsmooth lagrangian mechanics and variational collision integrators[END_REF], an impact may occur at time t j 0 and we have to use

M (q (t j )) ⇣ v + j v j ⌘ = W N (q (t j )) L N, j +W T (q (t j )) L T, j , (5)
0  ġ+ N, j + e N ġ N, j ? L N, j 0 , (6) 8 <
:

L T, j  µL N, j for ġ+ T, j + e T ġ T, j = 0 L T, j = ġ+ T, j +e T ġ T, j ġ+ T, j +e T ġ T, j
µL N, j for ġ+ T, j + e T ġ T, j 6 = 0 [START_REF] Arnold | Non-stiff integrators for differential algebraic systems of index 2[END_REF] instead of ( 2)-(4). We want to calculate the generalized position and generalized velocity functions on the time interval I := [0, T ] with T > 0:

q : [0, T ] ! IR N d , t 7 ! q (t) , (8)
v : [0, T ] ! IR N d , t 7 ! v (t) . (9)
Inria With the initial conditions

q (0) = q 0 , (10)
v (0) = v 0 , (11)
the evolution of q and v is non-impulsive almost everywhere due to (2)-( 4). The generalized mass matrix satisfies

M : IR N d ! IR N d ,N d , q 7 ! M (q) (12)
and the nonlinear generalized force has the structure

h : IR N d ⇥ IR N d ! IR N d , (q, v) 7 ! h (q, v) . ( 13
)
The matrices of generalized force directions W N and W T are the derivatives of the nonlinear normal and tangential gap functions

g N : IR N d ! IR N c , q 7 ! g N (q) , (14)
g T : IR N d ! IR N c , q 7 ! g T (q) (15)
with respect to q. These describe normal distances and tangential displacements of contacting bodies in a mechanical system and are therefore constraining its motion. The Lagrange multipliers

l N : [0, T ] ! IR N c , t 7 ! l N (t) (16)
can be identified with normal contact forces in the equality of forces (2) and in the non-penetration condition [START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF]. The Lagrange multipliers

l T : [0, T ] ! IR N c , t 7 ! l T (t) (17)
can be identified with tangential contact forces in the equality of forces [START_REF] Acary | Higher order event capturing time-stepping schemes for nonsmooth multibody systems with unilateral constraints and impacts[END_REF] and in (4), which contains sticking and sliding conditions with the friction coefficient µ 0 and the time-derivative denoted by a dot. The index set I 0 contains all constraints on the system.

I q 1 = k 2 I 0 : g N k (q)  0 is the continuous index set of the closed constraints , (18)
I q,v± 2 = k 2 I q 1 : ġN k (q, v ± )  0 is the continuous index set of the closed constraints, which stay closed .

For countable time instances t j , the evolution of the mechanical system might get impulsive, e.g. if g N k ⇤ (q (t j )) = 0 for some component k ⇤ , but g N k ⇤ (q (t)) > 0 for t j d  t < t j with an appropriate d > 0. The impact equations (5), Newton's impact law [START_REF] Arnold | Numerical methods for simulation in applied dynamics[END_REF] with the coefficient of restitution e N 2 [0, 1] and Newton's impact law [START_REF] Arnold | Non-stiff integrators for differential algebraic systems of index 2[END_REF] with the coefficient of restitution e T 2 [0, 1] have to be solved instead of ( 2)-( 4). This possibly leads to jumps in the velocity variables. Their derivatives do not exist anymore in the classical sense. One has to define the left-hand and right-hand limits

ġ N j := lim t"t j ġN (t) , ġ T j := lim t"t j ġT (t) , (20)
ġ+ N j := lim t#t j ġN (t) , ġ+ T j := lim t#t j ġT (t) . (21)
Then, the Lagrange multipliers

L N j = lim d #0 Z t j t j d l N dt , L T j = lim d #0 Z t j t j d l T dt (22)
describe the finite impulsive interaction in the sense of distributions. Equations ( 1)-( 7) describe impacting mechanical systems in general. Their mechanical, mathematical and numerical behaviour have been detailed in the monographs [START_REF] Brogliato | Nonsmooth mechanics : models, dynamics and control[END_REF][START_REF] Glocker | Set-valued force laws in rigid body dynamics : dynamics of non-smooth systems[END_REF][START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF][START_REF] Remco | Stability and convergence of mechanical systems with unilateral constraints[END_REF][START_REF] Pfeiffer | Mechanical system dynamics[END_REF][START_REF] Johnson | Contact Mechanics[END_REF][START_REF] Stewart | Dynamics with inequalities[END_REF] even discussing more sophisticated impact laws of kinetic or energetic type.

Slider-crank mechanism

The slider-crank mechanism in Fig. 1 is a nonlinear benchmark example for an impacting mechanical system with dry friction. The generalized positions q consist of the angles q 1 , q 2 and q 3 as well as the C r a n k

S li d e r R o d y x q 1 q 2 q 3
Figure 1: Slider-crank mechanism with unilateral constraints and friction [22,[START_REF] Schoeder | Discussion of the Gear-Gupta-Leimkuhler method for impacting mechanical systems[END_REF].

generalized velocities v consist of the angular velocities w 1 , w 2 and w 3 .

The generalized mass matrix satisfies

M : IR 3 ! IR 3,3 , q 7 ! M (q) = 0 @ J 1 + l 2 1 m 1 4 + m 2 + m 3 l 1 l 2 cos (q 1 q 2 ) m 2 2 + m 3 0 l 1 l 2 cos (q 1 q 2 ) m 2 2 + m 3 J 2 + l 2 2 m 2 4 + m 3 0 0 0 J 3 1 A , (23)
and the nonlinear generalized force can be obtained by

h : IR 3 ⇥ IR 3 ! IR 3 , (q, v) 7 ! h (q, v) = 0 @ l 1 l 2 sin (q 1 q 2 ) m 2 2 + m 3 w 2 2 gl 1 cos q 1 m 1 2 + m 2 + m 3 l 1 l 2 sin (q 1 q 2 ) m 2 2 + m 3 w 2 1 gl 2 cos q 2 m 2 2 + m 3 0 1 A . (24)
Thereby, it is l 1 the length of the crank and l 2 the length of the connecting rod. The inertia of crank, connecting rod and slider consist of translational masses, i.e., m 1 , m 2 and m 3 , as well as of rotational inertia values, i.e., J 1 , J 2 and J 3 . The system is subject to gravitation g.

g N 3 g N 4 g N 1 g N 2 q 3 y x g T 1 g T 3 g T 4 g T 2 2b 2a d

Inria

Considering the geometry of the slider according to Fig. 2, the nonlinear normal and tangential gap functions split up for each corner:

g N 1 (q) = d 2 l 1 sin q 1 l 2 sin q 2 + a sin q 3 b cos q 3 , (25)
g N 2 (q) = d 2 l 1 sin q 1 l 2 sin q 2 a sin q 3 b cos q 3 , (26)

g N 3 (q) = d 2 + l 1 sin q 1 + l 2 sin q 2 a sin q 3 b cos q 3 , (27)
g N 4 (q) = d 2 + l 1 sin q 1 + l 2 sin q 2 + a sin q 3 b cos q 3 , (28)
g T 1 (q) = l 1 cos q 1 + l 2 cos q 2 a cos q 3 b sin q 3 , (29)
g T 2 (q) = l 1 cos q 1 + l 2 cos q 2 + a cos q 3 b sin q 3 , (30)

g T 3 (q) = l 1 cos q 1 + l 2 cos q 2 a cos q 3 + b sin q 3 , (31)
g T 4 (q) = l 1 cos q 1 + l 2 cos q 2 + a cos q 3 + b sin q 3 . ( 32
)
It is 2a the length and 2b the height of the slider. The height of the notch is given by d and the gap can be defined by c. The matrices of generalized force directions are given by

W N : IR 3 ! IR 3,4 , q 7 ! W N (q) = 0 B B @ l 1 cos q 1 l 2 cos q 2 a cos q 3 + b sin q 3 l 1 cos q 1 l 2 cos q 2 a cos q 3 + b sin q 3 l 1 cos q 1 l 2 cos q 2 a cos q 3 + b sin q 3 l 1 cos q 1 l 2 cos q 2 a cos q 3 + b sin q 3 1 C C A T , (33)
W T : IR 3 ! IR 3,4 , q 7 ! W T (q) = 0 B B @ l 1 sin q 1 l 2 sin q 2 a sin q 3 b cos q 3 l 1 sin q 1 l 2 sin q 2 a sin q 3 b cos q 3 l 1 sin q 1 l 2 sin q 2 a sin q 3 + b cos q 3 l 1 sin q 1 l 2 sin q 2 a sin q 3 + b cos q

3 1 C C A T . ( 34
)
Only the i-th column of the matrix 4,3 does not vanish. We summarize this information in the matrix

∂W T N ∂ q i 2 IR
Q N = 0 B B @ l 1 sin q 1
l 2 sin q 2 a sin q 3 + b cos q 3 l 1 sin q 1 l 2 sin q 2 a sin q 3 + b cos q 3 l 1 sin q 1 l 2 sin q 2 a sin q 3 + b cos q 3 l 1 sin q 1 l 2 sin q 2 a sin q 3 + b cos q

3 1 C C A (35)
and define

v ⌅ = v 2 1 v 2 2 v 2 3 T . ( 36
)
Fixed characteristics used in this work are given in Table 1.

We focus on numerical integration and therefore point out particular previous work just in this direction. Classically, we distinguish event-driven respectively event-tracking schemes and timestepping respectively event-capturing schemes. q 10 = 0.0 q 2 0 = 0.0 q 3 0 = 0.0 w 1 0 = 150.01/s w 2 0 = 75.01/s w 3 0 = 0.01/s Table 1: Characteristics of the slider-crank mechanism with unilateral constraints and friction [22,[START_REF] Schoeder | Discussion of the Gear-Gupta-Leimkuhler method for impacting mechanical systems[END_REF].

Event-driven schemes

or stick (ġ T = 0)-slip (ġ T 6 = 0) transitions. In such cases, the exact transition time is resolved by a root-finding algorithm with respect to a given tolerance because the structure, i.e., the actual degree of freedom of the mathematical model, changes and therefore also the numerical discretization has to be adapted. In the case of an open-close transition, the impact equations ( 5)- [START_REF] Arnold | Non-stiff integrators for differential algebraic systems of index 2[END_REF] have to be solved separately and subsequently the non-impulsive integration has to be restarted. LSODAR 1 for index 0 or DASKR 2for index 2 formulations are typical examples for the derivation of the numerical model by standard integration schemes. Thereby, the relation of the index of the non-impulsive equations ( 1)-(4) and of their numerical counterpart play a crucial part for the discretization [START_REF] Eich | Numerical methods in multibody dynamics[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic-Problems[END_REF]. We discuss this concept for timestepping schemes in the following subsection. As event-driven schemes resolve the exact transition times, they, first, cannot resolve Zeno phenomena, i.e., an infinite number of impacts in a finite time interval. Second, even the resolution of the transition time itself may be crucial as the used tolerance is depending on the specific mechanical problem, e.g. on penetration velocities. These are two main drawbacks of event-driven schemes for their application to the mathematical description of impacting mechanical systems. However for non-impulsive intervals, event-driven schemes are extremely efficient because they are based on sophisticated numerical integration schemes with e.g. high accuracy for the numerical representation of the non-impulsive equations.

Classic timestepping schemes

Classic timestepping schemes do not distinguish impulsive and non-impulsive equations, i.e., they do not resolve exact transition times. Hence, no additional tolerances are required. Classic timestepping schemes discretize the equations of motion including the constraints, impact laws and impact equations in a physically consistent and uniform way [START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF]Chap. 10]. This means that possible impacts are the determining factor for the order of classic timestepping schemes. In comparison to the maximum possible order in non-impulsive phases, the actual order may be low. However, convergence results are available for classic timestepping schemes and not for event-driven schemes. Two mainstream concepts for classic timestepping schemes exist: Schatzman-Paoli [START_REF] Paoli | A numerical scheme for impact problems i: the onedimensional case[END_REF][START_REF] Paoli | A numerical scheme for impact problems i: the multidimensional case[END_REF] and Moreau-Jean schemes [START_REF] Jacques | Numerical aspects of the sweeping process[END_REF][START_REF] Jean | The nonsmooth contact dynamics method[END_REF]. The Schatzman-Paoli scheme discretizes normal gap functions on the level of positions and satisfies impact laws after several time steps. That is why, we prefer to formulate the impact law on the level of velocities in discrete time, which is the core of Moreau-Jean schemes.

To exemplary formulate the classic explicit Moreau-Jean scheme, we construct the approximations at time t i q i ⇡ q (t i ) ,

v i ⇡ lim t"t i v (t) , v + i ⇡ lim t#t i v (t) , (37)
L N i ⇡ Â j:t j 2I i L N, j , L T i ⇡ Â j:t j 2I i L T, j (38)
and use them in the abbreviations

M 1 i := M 1 (q i ) , (40)
h ± i := h q i , v ± i , (41)
W N i := W N (q i ) ,W T i := W T (q i ) . (42)
Then, the classic explicit Moreau-Jean scheme on an arbitrary time interval I i := [t i 1 ,t i ) with the time step-size Dt i := t i t i 1 is given by

q i = q i 1 + Dt i v + i 1 , (43)
v + i = v + i 1 + Dt i M 1 i 1 h + i 1 + M 1 i [W N i L N i +W T i L T i ] (44)
together with the active normal impact equations on velocity level

L N i ,I i 1,M proj IR + 0 ✓ L N i ,I i 1,M r(ġ + N i ,I i 1,M + e N ġ+ N i 1 ,I i 1,M ) ◆ = 0 (45)
and the active tangential impact equations on velocity level

L T i ,I i 1,M proj C T (L N i ,I i 1,M ) ✓ L T i ,I i 1,M r(ġ + T i ,I i 1,M + e T ġ+ T i 1 ,I i 1,M ) ◆ = 0 . (46)
To simplify the nonlinearity, usually one estimates active contact equations in advance by checking if a prediction, e.g. the leapfrog rule, is negative:

I i 1,M := ⇢ k 2 I 0 : gi N k := g N k (q i 1 ) + Dt 2 ġNk (q i 1 , v + i 1 )  0 . ( 47
)
We have used the projection to a convex set C as a special case of proximations [START_REF] Tyrrell | Convex Analysis[END_REF] to equivalently transform complementarity formulations for normal and tangential contacts. It has to be interpreted rowby-row [START_REF] Schindler | Understanding the difference between prox and complementarity formulations for simulation of systems with contact[END_REF]:

proj C : IR ! IR , x 7 ! proj C (x) = arg min x ⇤ 2C kx x ⇤ k . ( 48
)
The convex set to describe dry friction is the friction disc with the local normal force as radius:

C T : IR ! P IR 2 : y 7 ! C T (y) = x 2 IR 2 | kxk  µ|y| . (49)
The parameter r > 0 is mathematically arbitrary. It physically transforms units and can also be interpreted as a stabilization in an augmented Lagrangian setting [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF]. Therefore, it can be used to improve convergence of numerical solution schemes, which are applied to solve ( 43)- [START_REF] Paoli | A numerical scheme for impact problems i: the onedimensional case[END_REF] for the unknowns

q i , v + i , L N i ,I i 1,M and L T i ,I i 1,M
. The classic explicit Moreau-Jean scheme does not adapt the time step-size and naturally combines non-impulsive and impulsive equations, contacts as well as impacts by calculating implicitly the finite mean impulses

L N i ,I i 1,M and L T i ,I i 1,M
within I i ; their trajectory within I i in this sense is an assumption. Hence, several extensions are of possible interest and are current research topics:

• Constraint stabilization [START_REF] Schoeder | Discussion of the Gear-Gupta-Leimkuhler method for impacting mechanical systems[END_REF][START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF][START_REF] Brüls | Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized-a scheme[END_REF] Moreau-Jean schemes in general are formulated on the level of velocities. Hence, a linear drift-off effect from the non-penetration invariant will occur [START_REF] Führer | Numerical solution of differential-algebraic equations for constrained mechanical motion[END_REF][START_REF] Arnold | Numerical methods for simulation in applied dynamics[END_REF]. The Gear-Gupta-Leimkuhler formulation [START_REF] William Gear | Automatic integration of Euler-Lagrange equations with constraints[END_REF] treats the non-penetration invariant as additional constraint and therefore inserts a derivative projection; applied to Moreau's timestepping interpretation of the midpoint rule [START_REF] Jacques | Numerical aspects of the sweeping process[END_REF], this can be understood as a symmetric projection [START_REF] Hairer | Symmetric projection methods for differential equations on manifolds[END_REF][START_REF] Simeon | Computational Flexible Multibody Dynamics[END_REF].

• High-frequency damping [START_REF] Chen | A nonsmooth generalizeda scheme for flexible multibody systems with unilateral constraints[END_REF][START_REF] Brüls | Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized-a scheme[END_REF] Artificial high-frequency oscillations occur due to finite element discretizations and e.g. excitation in contact problems. One possible remedy is applying the generalized-a method [START_REF] Chung | A time integration algorithm for structural dynamics with improved numerical dissipation: the generalized-a method[END_REF] in the sense of timestepping schemes.

• Variational integration [START_REF] Kaufman | Geometric numerical integration of inequality constrained, nonsmooth Hamiltonian systems[END_REF] Sophisticated integration schemes can often be derived from a discrete variational principle [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF]. This concept is tried to be extended to timestepping schemes for impacting mechanical systems. However, some kind of splitting seems to be necessary because a symplectic method will not stay symplectic for an impact event [START_REF] Stewart | Rigid-body dynamics with friction and impact[END_REF].

• Higher order timestepping and step-size adaptation There are two types of higher order timestepping schemes which can consistently deal with impacts.

-Augmented timestepping schemes [START_REF] Studer | Numerics of unilateral contacts and friction : modeling and numerical time integration in non-smooth dynamics[END_REF][START_REF] Huber | Higher order integration of non-smooth dynamical systems using parallel computed extrapolation methods based on time-stepping schemes[END_REF] Within an augmented timestepping scheme, one applies classic augmentation strategies like order extrapolation whenever one heuristically expects no impacts during an integration step. Extrapolation induces instabilities in closed contacts because of chattering in Aitken-Neville schemes or because of the absence of splitting of non-impulsive and impulsive force propagations. It is conceptional a serial process. That is why, usually one applies a fixed increased order when we anticipate a non-impulsive phase. Also from a practical point of view, i.e., not infinitely-differentiable input data, an extrapolation-based augmentation strategy like in [START_REF] Lubich | Extrapolation integrators for constrained multibody systems[END_REF] seems to be over-designed.

-Mixed timestepping schemes [START_REF] Acary | Higher order event capturing time-stepping schemes for nonsmooth multibody systems with unilateral constraints and impacts[END_REF][START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook[END_REF][START_REF] Esefeld | A coupling approach for the numerical integration of non-smooth multibody systems[END_REF] A mixed timestepping scheme couples standard integration methods for non-impulsive differential algebraic equations with timestepping schemes for impulsive phases usually by heuristics.

We extend [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook[END_REF] to multi-contact problems with friction. In [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook[END_REF], heuristics for switching between impulsive and non-impulsive phases can be avoided by embedding in discontinuous Galerkin methods. The synchronization of fixed increased order integration and consistent low order integration is automatic but as we will see perhaps not the most robust way for an implementation.

Inria

Step size adaptation can be applied for both higher order timestepping strategies. It extends the classic approach [START_REF] Deuflhard | Scientific computing with ordinary differential equations[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic-Problems[END_REF] by heuristics which meet the timestepping idea. We mention foreseeing gap-estimations, retrospective time-step bisection, time-step size switching maintaining the overall integration order and error estimation based on Richardson extrapolation. For the latter, one could exclude velocities, discuss appropriate norms, separate impacts interval-by-interval or include the penetration depth of closed contacts. However to the knowledge of the authors, these strategies have not been theoretically analyzed in the literature.

Half-explicit timestepping schemes on acceleration level

This section generalizes timestepping schemes based on discontinuous Galerkin methods [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook[END_REF] to multicontact examples and shows their advantages and disadvantages. The theoretical benefit of these schemes over mixed timestepping schemes is that they do not rely on heuristics. Thereby as in [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook[END_REF], we start by introducing contact forces on acceleration level and by splitting non-impulsive and impulsive force propagation. We try to get some experience with linear trial functions for the velocity discretization, i.e., half-explicit "forecasting" trapezoidal rules, not regarding friction:

q i 1,0 = q i 1 , (50)
q i 1,1 = q i 1 + Dt i 2 {v i 1,0 + v i 1,1 } , (51)
q i = q i 1 + Dt i 2 {v i 1,0 + v i 1,1 } , (52)
v i 1,0 = v + i 1 , (53)
v i 1,1 = v + i 1 + Dt i M 1 i 1 h h + i 1 +W N i 1 l + N i 1 i , (54)
v + i = v + i 1 + Dt i 2 M 1 i 1 h h + i 1 +W N i 1 l + N i 1 i + Dt i 2 M 1 i 1,1 ⇥ h i 1,1 +W N i 1,1 l N i ⇤ + M 1 i W N i L N i (55)
together with the active contact equations on acceleration level

l + N i 1 ,I (i 1),0 2 proj IR + 0  l + N i 1 ,I (i 1),0 2 r g+ N i 1 ,I (i 1),0 2 = 0 , (56)
l N i ,I (i 1),1 2 proj IR + 0  l N i ,I (i 1),1 2 r g N i ,I (i 1),1 2 = 0 (57)
and the active impact equations on velocity level similar to (45) but referring to another index set, which we have to define. It is

I i,l 1 = k 2 I 0 : g N k (q i,l )  0 the discrete index set of the closed constraints , (58)
I i,l± 2 = n k 2 I i,l 1 : ġN k (q i,l , v ± i,l )  0 o the discrete index set of the closed constraints, which stay closed , (59)
and

l N i ⇡ lim t"t i l N (t) , l + N i ⇡ lim t#t i l N (t) . (60)
Now, what is active for this half-explicit trapezoidal rule? We assume that the generalized position q i 1 and the right-hand limit of the generalized velocity v + i 1 are known at t i 1 as illustrated in the graphical

t t i t i 1 t i+1 t i+2 v
Figure 3: Interpretation of velocity jumps [START_REF] Schoeder | Discussion of the Gear-Gupta-Leimkuhler method for impacting mechanical systems[END_REF].

interpretation in Fig. 3. These values are set to the first stage of the generalized position ( 50) and of the generalized velocity [START_REF] Schindler | Understanding the difference between prox and complementarity formulations for simulation of systems with contact[END_REF], respectively. The generalized velocity is assumed to be continuous in the interior of the time interval. In the specific case of linear trial functions, there is a linear velocity propagation until its left-hand limit v i 1,1 at t i . Hence, the second stage (54) of the generalized velocity is calculated with a step of the explicit Euler method. Everything is known but for the right-hand limit of the contact force l + N i 1 . We postpone the question whether to incorporate (56) for its calculation. From a general point of view subsequently, the second stage of the generalized position at t i is calculated with the trapezoidal rule [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook[END_REF]. As the jump at t i occurs in such a way that the generalized velocity at t i is right-continuous, the left-hand limit of the generalized velocity v i is corrected by a velocity jump leading to the right-hand limit of the generalized velocity v + i . The respective trapezoidal rule (55) involves the unknown left-hand limit of the contact force l N i and the impulse L N i as well as ( 57) and ( 45) if necessary. The effects of all impulses that would technically occur in the interior of the time interval are summarized at the end of the time interval.

The procedure of incorporating ( 56)-( 57) and ( 45) can be implemented straightforward for a decoupled bouncing ball example (cf. Sect. 6.2 for gravitational acceleration in y-direction) because there is only one single contact possibility [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook[END_REF]:

• If on the one side g N i 1  0 at the beginning of the time interval, we will not consider the impulse L N i = 0 at the end of the time interval to take advantage of higher order integration possibilities. If also ġ+ 56) is active for the explicit evaluation of l + N i 1 with a root-finding algorithm. Thereby, we substitute the circumstantial local acceleration because of (2):

N i 1 > 0, it is l + N i 1 = 0. Otherwise, Equation (
g+ N i 1 := W T N i 1 v+ i 1 = W T N i 1 M 1 i 1 h + i 1 +W T N i 1 M 1 i 1 W N i 1 l + N i 1 . (61)
If additionally 57) is active for the explicit evaluation of l N i with a root-finding algorithm and the local acceleration

g N i 1,1 > 0 or ġN i 1,1 > 0, it is l N i = 0. Otherwise, Equation (
g N i := W T N i 1,1 v i = W T N i 1,1 M 1 i 1,1 h i 1,1 +W T N i 1,1 M 1 i 1,1 W N i 1,1 l N i . (62)
• If on the other side g N i 1 > 0 at the beginning of the time interval, clearly l + N i 1 = 0. However, we will not consider the contact force l N i = 0 at the end of the time interval because analytically it does not exist during an impact. Concerning the impulse, the additional condition g N i 1,1 > 0 yields L N i = 0. Otherwise, Equation ( 45) is active for the implicit solution of both v + i and L N i with a root-finding algorithm and the local velocity

ġ+ N i := W T N i v + i = ġ N i +W T N i M 1 i W N i L N i (63)
Inria according to [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF].

For the decoupled bouncing ball example, the proposed half-explicit trapezoidal rule works quite well and shows the expected characteristics in numerical experiments, i.e., first order accuracy in impulsive phases and second order accuracy in non-impulsive phases [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook[END_REF]. For multi-contact situations it is not clear how to decide which contacts shall be considered active concerning the three different constraint equations ( 56)-( 57) and [START_REF] Newmark | A method for computation for structural dynamics[END_REF]. The half-explicit trapezoidal rule needs further interpretation for a consistent extension to simultaneous contact occurrences and nonlinear dynamics.

Calculation of contact forces on acceleration level

It seems not to be a bad idea to use (3) together with (2) for the calculation of a contact force l N at a certain time also for more general nonlinear mechanical systems. We again interpret (3) on acceleration level

l N = ( 0 if g N > 0 _ ġN > 0 proj IR + 0 [l N r gN ] else (64)
using the projection-formulation row-by-row. Often, we have to deal with holonomic-scleronomic constraints:

ġN = W T N v , (65) gN
= Ẇ T N |{z} Â i ∂W T N ∂ q i v i v +W T N v . (66)
We can simplify (66) to

gN = Q N v ⌅ +W T N v (67)
with v ⌅ containing squared generalized velocity combinations. In comparison to the decoupled bouncing ball example, the nonlinearity in (67) and in the right-hand side h newly arises. As before, we eliminate v in (67) by using (2):

gN = Q N v ⌅ +W T N M 1 h + G N l N (68)
with the Delassus matrix [START_REF] Remco | Stability and convergence of mechanical systems with unilateral constraints[END_REF], i.e., the mass action matrix [START_REF] Pfeiffer | Mechanical system dynamics[END_REF],

G N = W T N M 1 W N . (69)
We focus on active contacts, I 2 on acceleration level, and transform (64) formally using row-by-row interpretation:

l N,I 2 = proj IR + 0 2 6 6 4 l N,I 2 r(Q N,I 2 v ⌅ +W T N,I 2 M 1 h + G N,I 2 l N,I 2 | {z } gN,I 2 ) 3 7 7 5 . (70)
In the multi-contact case, active contacts might be depending. Hence if we decide to use e.g. a semismooth Newton method as root-finding algorithm to solve (70), we have to switch to semi-smooth variants of the Gauss-Newton method. This means that we systematically have to choose an approximate root l N,I 2 of the function

f : IR |I 2 | ! IR |I 2 | , l N,I 2 7 ! f (l N,I 2 ) = l N,I 2 proj IR + 0 ⇥ l N,I 2 r gN,I 2 (l N,I 2 ) ⇤ (71)
with the Moore-Penrose pseudoinverse operator pinv. The Gauss-Newton algorithm reads

l N,I 2 = 0 , f = f ( l N,I 2 ) while f > TOL -f ( l N,I 2 ) = I Q ⇥ l N,I 2 r gN,I 2 ( lN,I 2 ) ⇤ (I rG N,I 2 ) l N new ,I 2 = l N,I 2 pinv(-f ( l N,I 2 )) f l N,I 2 = l N new ,I 2 f = f ( l N,I 2 ) end ( 72
)
The occurring Heaviside function

Q : IR ! IR , x 7 ! Q(x) = ( 0 if x < 0 1 else (73)
is interpreted row-by-row and we use a fixed parameter r [START_REF] Schindler | Understanding the difference between prox and complementarity formulations for simulation of systems with contact[END_REF].

Smooth position and velocity prediction

After having calculated the right-hand limit of the contact force l + N i 1 at t i 1 according to Sect. 3.1, we predict the left-hand limit of the generalized velocity at t i , i.e., the second stage of the generalized velocity:

v i 1,1 = v + i 1 + Dt i M 1 i 1 h h + i 1 +W N i 1 l + N i 1 i . ( 74
)
Subsequently, we calculate the second stage of the generalized position:

q i 1,1 = q i 1 + Dt i 2 {v i 1,0 + v i 1,1 } . ( 75
)

Assignment of contact forces and impulses

Now, we have a pair of generalized position and generalized velocity at the left-hand limit of t i that allows us to calculate the gap function g N i .

• If on the one hand there is no gap function which has been inactive (g N k ⇤ q i 1,0 > 0) at t i 1 and which now becomes active (g N k ⇤ q i 1,1  0) at t i , no impact has to be considered. We set L N i = 0 to fully profit of higher integration possibilities and calculate the left-hand limit of the contact force l N i at t i according to Sect. 3.1. Knowing l N i , we improve the predicted left-hand limit of the generalized velocity at t i :

v i := v + i 1 + Dt i 2 M 1 i 1 h h + i 1 +W N i 1 l + N i 1 i + Dt i 2 M 1 i 1,1 ⇥ h i 1,1 +W N i 1,1 l N i ⇤ . ( 76
)
Inria • If on the other hand there is at least one inactive gap function (g N k ⇤ q i 1,0 > 0) at t i 1 that now becomes active (g N k ⇤ q i 1,1  0), an impact occurs in the whole rigidly connected component of the impact source. Hence in this rigidly connected component, the calculation of l N i is not consistent and does not improve the accuracy of the scheme. For the ease of description, we set l N i = 0 and also l + N i 1 = 0 everywhere in the multibody system and calculate the left-hand limit of the generalized velocity at t i , i.e., the free velocity:

v i := v + i 1 + Dt i 2 n M 1 i 1 h + i 1 + M 1 i 1,1 h i 1,1 o . ( 77
)
This global consequence does not occur in the decoupled bouncing ball example. For the computation of v + i and L N i , we write ( 6) on velocity level row-by-row:

L N i = 8 > > > > > < > > > > > : 0 if g N q i 1,1 > 0 proj IR + 0 2 6 6 6 4
L N i r(ġ + N i + e N ġ N i | {z } ġ+ N i ) 3 7 7 7 5 else . (78)
Equation ( 5) can be used to eliminate ġ+ N i , row-by-row resulting in

L N i ,I i 1,1 1 = proj IR + 0 2 6 6 6 6 6 6 4
L N i ,I i 1,1 1 r(G N i ,I i 1,1 1 L N i ,I i 1,1 1
+ (I + e N,I i 1,1 1 )W T N i ,I i 1,1 1 v i | {z } ġ+ N i ,I i 1, 1 1 ) 3 7 7 7 7 7 7 5 (79)
with G N i := G N (q i ). We look for the roots of

f : IR |I i 1,1 1 | ! IR |I i 1,1 1 | , L N i ,I i 1,1 1 7 ! f (L N i ,I i 1,1 1 ) = L N i ,I i 1,1 1 proj IR + 0  L N i ,I i 1,1 1 r ġ+ N i ,I i 1, 1 1 (80)
with a semi-smooth Gauss-Newton method as in (72). The derivative of f at

L N i ,I i 1,1 1 is given by -f (L N i ,I i 1,1 1 ) = I Q  L N i ,I i 1,1 1 r ġ+ N i ,I i 1,1 1
(I rG N i ,I i 1, 1 1
) .

We have to use structures from classic timestepping schemes and event-driven schemes to define the half-explicit trapezoidal rule.

State increment

Finally, we update the generalized velocity:

v + i = v i + M 1 i W N i L N i . (82)
Hence in the impact-free case, it is v + i = v i as expected. For the generalized position, it is

q i = q i 1 + Dt i 2 {v i 1,0 + v i 1,1 } . ( 83
) RR n°8623
Specify characteristics, start time t = 0, end time T , time step size Dt

Initialize i = 1 For t < T Evaluate left-hand side Evaluate q i 1,0 = q i 1 (50) Evaluate v i 1,0 = v + i 1 (53) Evaluate l + N i 1 on acceleration level (I (i 1), 0 2
)

Predict right-hand side Compute v i 1,1 = v + i 1 + Dt i M 1 i 1 h h + i 1 +W N i 1 l + N i 1 i (54)
Compute

q i 1,1 = q i 1 + Dt i 2 {v i 1,0 + v i 1,1 } (51) Decide strategy if @k ⇤ : g N k ⇤ q i 1,0 > 0 ^gN k ⇤ q i 1,1  0 (14) Set L N i = 0 Evaluate l N i on acceleration level (I (i 1), 1 2
)

Compute v i = v + i 1 + Dt i 2 M 1 i 1 h h + i 1 +W N i 1 l + N i 1 i + Dt i 2 M 1 i 1,1 ⇥ h i 1,1 +W N i 1,1 l N i ⇤ (76) if 9k ⇤ : g N k ⇤ q i 1,0 > 0 ^gN k ⇤ q i 1,1  0 (14) Set l + N i 1 = 0 and l N i = 0 Compute v i = v + i 1 + Dt i 2 n M 1 i 1 h + i 1 + M 1 i 1,1 h i 1,1 o (77)
Compute

L N i on velocity level (I i 1, 1 1
)

State increment Compute v + i = v i + M 1 i W N i L N i (82) Compute q i = q i 1 + Dt i 2 {v i 1,0 + v i 1,1 } (

Overview

The overall algorithm can be summarized as shown in Fig. 4.

First analysis of simulation results

Simulation results for the slider-crank mechanism in Fig. 1 of Sect. 2.1 with e N = 0.4 and Dt = 10 4 s are depicted in Figs. 5 and6. The curves show a nearly perfect behavior. The benefits from the timediscontinuous Galerkin schemes presented in [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook[END_REF] concerning theoretical investigations and a decoupled bouncing ball example can also be carried over to multi-contact examples. Impulsive and non-impulsive periods are seperated; hence, higher-order trial functions can be used in non-impulsive periods whereas consistency is preserved for the impacts. Nevertheless, we improve the curves concerning the following observations. The penetration of the bordering wall occurs because timestepping schemes in general do not detect but only capture events, e.g. the time instant when a gap function closes. Contact activity means g N i,l  0 ^ġ ± N i,l  0 and not g N i,l = 0 ^ġ ± N i,l = 0. We will see that a velocity level discretization of the contact forces is a remedy for the drift-off effect and opens the path to include friction naturally. On acceleration level, the implementation of the friction transitions would also be rather complicated. On velocity level, the nonlinearity in the equations of motion, e.g. v ⌅ , will be reduced. Nonlinearity in general induces the not always decreasing energy trend of the trapezoidal rule which differs from the midpoint rule in the nonlinear regime [START_REF] Simeon | Computational Flexible Multibody Dynamics[END_REF]. On acceleration level, contacts and impacts are defined by different activity rules using I 2 and I 1 . As a compromise between acceleration and position level discretizations, the velocity level discretization will unify contact and impact activation rules.

Explicit timestepping schemes on velocity level

In this section, we reduce the drawbacks, i.e., the degree of nonlinearity and the drift-off effect, for halfexplicit timestepping schemes on velocity level from the preceding section by introducing a velocity level discretization: together with the active normal contact equations on velocity level

q i 1,0 = q i 1 , (84)
q i 1,1 = q i 1 + Dt i 2 {v i 1,0 + v i 1,1 } , (85)
q i = q i 1 + Dt i 2 {v i 1,0 + v i 1,1 } , ( 86
)
v i 1,0 = v + i 1 , (87)
v i 1,1 = v + i 1 + Dt i M 1 i 1 h h + i 1 +W N i 1 l + N i 1 +W T i 1 l + T i 1 i , (88)
v + i = v + i 1 + Dt i 2 M 1 i 1 h h + i 1 +W N i 1 l + N i 1 +W T i 1 l + T i 1 i + Dt i 2 M 1 i 1,1 ⇥ h i 1,1 +W N i 1,1 l N i +W T i 1,1 l T i ⇤ + M 1 i [W N i L N i +W T i L T i ] (89
l + N i 1 ,I i 1,0 1 proj IR + 0 ✓ l + N i 1 ,I i 1,0 1 r ġ+ N i 1 ,I i 1,0 1 ◆ = 0 , (90)
l N i ,I i 1,1 1 proj IR + 0 ✓ l N i ,I i 1,1 1 r ġN i 1,1 ,I i 1,1 1 ◆ = 0 (91)
and the active normal impact equations on velocity level (78). Then, we can easily include friction which naturally is stated on velocity level. For the active tangential contact equations on velocity level, we get

l + T i 1 ,I i 1,0 1 proj C T (l + N i 1 ,I i 1,0 1 ) ✓ l + T i 1 ,I i 1,0 1 r ġ+ T i 1 ,I i 1,0 1 ◆ = 0 , (92)
l T i ,I i 1,1 1 proj C T (l N i ,I i 1, 1 1
)

✓ l T i ,I i 1,1 1 r ġT i 1,1 ,I i 1,1 1 ◆ = 0 (93) with l T i ⇡ lim t"t i l T (t) , l + T i ⇡ lim t#t i l T (t) . (94)
The active tangential impact equations on velocity level are similar to [START_REF] Paoli | A numerical scheme for impact problems i: the onedimensional case[END_REF].

Calculation of contact forces on velocity level

As in Sect. 3.1, we use ( 3) and (4) together with (2) for the calculation of contact forces l N and l T . We interpret ( 3) and ( 4) on velocity level

l N = ( 0 if g N > 0 proj IR + 0 [l N r ġN ] else , (95)
l T = ( 0 if g N > 0 proj C T (l N ) [l T r ġT ] else (96)
Inria using the projection-formulation row-by-row. Thereby now in contrast to Sect. 3.1, the necessary local kinematics

ġN = W T N v , ġT = W T T v (97)
is already known and we cannot directly calculate the contact forces with (95), (96). Hence, we start calculating equivalent forces l N and l T such that the given local velocities are projected into their respective admissible space, e.g. during one time step -which is an assumption of an explicit representation. That is why, we substitute

ġN,proj = ġN + Dt W T N M 1 h + G N l N +W T N M 1 W T l T , (98)
ġT,proj

= ġT + Dt W T T M 1 h +W T T M 1 W N l N + G T l T ( 99
)
in ( 3) and (4) using

h : IR N d ⇥ IR N d ! IR N d , (q, v) 7 ! h (q, v) = ( h (q, v) option 1: more implicit trend 0 option 2: less implicit trend . ( 100
)
Thereby, it is

G N = W T N M 1 W N , G T = W T T M 1 W T . (101)
We again focus on active contacts and transform (95) and (96) formally using row-by-row interpretation:

l N,I 1 = proj IR + 0 ⇥ l N,I 1 r ġN,proj,I 1 (l N,I 1 , l T,I 1 ) ⇤ , (102)
l T,I 1 = proj C T (l N,I 1 ) ⇥ l T,I 1 r ġT,proj,I 1 (l N,I 1 , l T,I 1 ) ⇤ . (103)
As in Sect. 3.1, we exemplary choose a semi-smooth variant of the Gauss-Newton method, i.e., we systematically select an approximate root l N,I 1 , l T,I 1 of the function

f : IR |I 1 | ⇥ IR 2|I 1 | ! IR |I 1 | ⇥ IR 2|I 1 | , l N,I 1 , l T,I 1 7 ! f (l N,I 1 , l T,I 1 ) = l N,I 1 proj IR + 0 ⇥ l N,I 1 r ġN,proj,I 1 (l N,I 1 , l T,I 1 ) ⇤ l T,I 1 proj C T (l N,I 1 ) ⇥ l T,I 1 r ġT,proj,I 1 (l N,I 1 , l T,I 1 ) ⇤ ! (104)
with the Moore-Penrose pseudoinverse operator pinv. The Gauss-Newton algorithm reads

l N,I 1 , l T,I 1 = (0, 0) , f = f l N,I 1 , l T,I 1 while f > TOL -f l N,I 1 , l T,I 1 = ✓ I 0 0 I ◆ ✓ Q N (I rDtG N,I 1 ) Q N ( rDtW T N,I 1 M 1 W T,I 1 ) Q T ( rDt(W T T,I 1 M 1 W N,I 1 ) Q T (I rDtG T,I 1 ) ◆ ⇣ l N new,I 1 , l T new,I 1 ⌘ = l N,I 1 , l T,I 1 pinv(-f l N,I 1 , l T,I 1 ) f l N,I 1 , l T,I 1 = ⇣ l N new,I 1 , l T new,I 1 ⌘ f = f l N,I 1 , l T,I 1 end (105)
The occurring Heaviside functions

Q N : IR |I 1 | ⇥ IR 2|I 1 | ! diag |I 1 |,|I 1 | , Q N kk l N,I 1 , l T,I 1 = ( 0 if l N k r ġN k ,proj < 0 1 else , (106)
Q T : IR |I 1 | ⇥ IR 2|I 1 | ! diag |I 1 |,|I 1 | , Q T kk l N,I 1 , l T,I 1 = ( 0 if l T k r ġT k ,proj > µ l N k 1 else . ( 107
)
are interpreted row-by-row and we use a fixed parameter r.

Smooth position and velocity prediction

According to the acceleration level case, we predict the left-hand limit of the generalized velocity at t i , which now also includes the friction contribution:

v i 1,1 = v + i 1 + Dt i M 1 i 1 h h + i 1 +W N i 1 l + N i 1 +W T i 1 l + T i 1 i . (108)
Subsequently, we calculate the second stage of the generalized position with (75).

Assignment of contact forces and impulses

With the gap function g N i 1,1 , we can distinguish non-impulsive and impulsive phases.

• Non-impulsive:

v i = v + i 1 + Dt i 2 M 1 i 1 h h + i 1 +W N i 1 l + N i 1 +W T i 1 l + T i 1 i + Dt i 2 M 1 i 1,1 ⇥ h i 1,1 +W N i 1,1 l N i +W T i 1,1 l T i ⇤ . (109)
• Impulsive:

We remind of the free velocity (77) and of ( 6) in its row-by-row description (78). Then, we write also [START_REF] Arnold | Non-stiff integrators for differential algebraic systems of index 2[END_REF] on velocity level row-by-row:

L T i = 8 > > > > > < > > > > > : 0 if g N i 1,1 > 0 proj C T (L N i ) 2 6 6 6 4
L T i r(ġ + T i + e T ġ T i | {z } ġ+ T i ) 3 7 7 7 5 else . (110)
Elimination of ġ+ N i and ġ+ T i for active contacts on position level includes friction contributions:

L N i ,I i 1,1 1 = proj IR + 0 h L N i ,I i 1,1 1 r ġN i ,I i 1,1 1 i , (111)
L T i ,I i 1,1 1 = proj C T (L N i ,I i 1, 1 1
)

h L T i ,I i 1,1 1 r ġT i ,I i 1,1 1 i . ( 112
)
We solve for L N,I 1 , L T,I 1 as roots of a summarizing function

f : IR |I 1 | ⇥ IR 2|I 1 | ! IR |I 1 | ⇥ IR 2|I 1 | (113)
with a semi-smooth Gauss-Newton method.

Inria

State increment

Finally, we update the generalized velocity by

v + i = v i + M 1 i (W N i L N i +W T i L T i ) (114)
and the generalized position by (83).

Overview

Specify characteristics, start time t = 0, end time T , time step size Dt

Initialize i = 1 For t < T Evaluate left-hand side Evaluate q i 1,0 = q i 1 (50) Evaluate v i 1,0 = v + i 1 (53) Evaluate l + N i 1 , l + T i 1 on velocity level (I i 1,0 1 ) Predict right-hand side Compute v i 1,1 = v + i 1 + Dt i M 1 i 1 h h + i 1 +W N i 1 l + N i 1 +W T i 1 l + T i 1 i (54)
Compute

q i 1,1 = q i 1 + Dt i 2 {v i 1,0 + v i 1,1 } (51) Decide strategy if @k ⇤ : g N k ⇤ q i 1,0 > 0 ^gN k ⇤ q i 1,1  0 (14) Set L N i = 0, L T i = 0 Evaluate l N i , l T i on velocity level (I i 1, 1 1
)

Compute v i = v + i 1 + Dt i 2 M 1 i 1 h h + i 1 +W N i 1 l + N i 1 +W T i 1 l + T i 1 i + Dt i 2 M 1 i 1,1 ⇥ h i 1,1 +W N i 1,1 l N i +W T i 1,1 l T i ⇤ (109) if 9k ⇤ : g N k ⇤ q i 1,0 > 0 ^gN k ⇤ q i 1,1  0 (14) Set l + N i 1 = 0, l + T i 1 = 0 and l N i = 0, l T i = 0 Compute v i = v + i 1 + Dt i 2 n M 1 i 1 h + i 1 + M 1 i 1,1 h i 1,1 o (77) Compute L N i , L T i on velocity level (I i 1 ) State increment Compute v + i = v i + M 1 i (W N i L N i +W T i L T i ) (114) Compute q i = q i 1 + Dt i 2 {v i 1,0 + v i 1,1 } (83
) Write result of time step, e.g. q i and v + i Update t by Dt i and increment i by 1 Postprocessing Figure 7: Flowchart of the explicit timestepping scheme on velocity level.

The overall algorithm can be summarized as shown in Fig. 7.

Second analysis of simulation results

Simulation results of the slider-crank mechanism in Fig. 1 of Sect. 2.1 with e N = 0.4, Dt = 10 4 s and h ⌘ 0 are depicted in Fig. 8. The drift-off effect is quite strong, even stronger than for the half-explicit acceleration level approach. The choice h ⌘ h additionally induces bad oscillatory behavior into the nonlinear slider-crank mechanism, i.e., loss of contact which even cannot be resolved by keeping the activation rule from the left to the right interval end (I i 1,1

1 = I i 1,0 1
, not shown). An implicit calculation of the contact forces is proposed as a remedy.

Half-explicit timestepping schemes on velocity level

This section introduces the final version of the proposed schemes by correctly applying velocity level discretizations in an implicit way and hence reducing the degree of nonlinearity and the drift-off effect from half-explicit timestepping schemes on acceleration level (Sect. 3). Using an implicit evaluation of constraints [START_REF] Brasey | Hem 5 user's guide[END_REF][START_REF] Brasey | Methodes demi-explicites pour equations differentielles algebraiques semi-explicites d'index 2[END_REF][START_REF] Murua | Phem 56 user's guide[END_REF][START_REF] Murua | Partitioned half-explicit runge-kutta methods for differential-algebraicsystems of index 2[END_REF][START_REF] Arnold | Non-stiff integrators for differential algebraic systems of index 2[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations I: Nonstiff Problems[END_REF][START_REF] Hairer | Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic-Problems[END_REF], we do not have a projection anymore as in (98) and (99) for calculating both l + N i 1 , l + T i 1 and l N i , l T i . Instead, we solve the following equations for the second stage of the velocity

0 = l + N i 1 ,I i 1,0 1 proj IR + 0 ✓ l + N i 1 ,I i 1,0 1 rW T N,I i 1,0 1
(q i 1 + Dt i v + i 1 )v i 1,1 ◆ , ( 115
) 0 = l + T i 1 ,I i 1,0 1 proj C T (l + N i 1 ,I i 1,0 1 ) ✓ l + T i 1 ,I i 1,0 1 rW T T,I i 1,0 1
(q i 1 + Dt i v + i 1 )v i 1,1 ◆ (116)
and the following equations assuming that the output equation of the trapezoidal rule defines a virtual third stage:

0 = l N i ,I i 1,1 1 proj IR + 0 ✓ l N i ,I i 1,1 1 rW T N i ,I i 1,1 1 v i ◆ , ( 117
) 0 = l T i ,I i 1,1 1 proj C T (l N i ,I i 1,1 1 ) ✓ l T i ,I i 1,1 1 rW T T i ,I i 1,1 1 v i ◆ . ( 118
)
It is important to evaluate W on the right-hand side of the interval. As we do not know its value for the second stage of the velocity, we have to calculate a prediction, which, however, is not used for the Inria calculation of any index set. The virtual third stage of the velocity coincides with the output equation of the trapezoidal rule. Hence, the velocity level constraint is automatically satisfied also for the beginning of the next time step. The overall algorithm can be summarized as shown in Fig. 9.

Specify characteristics, start time t = 0, end time T , time step size Dt 53) Predict right-hand side Compute Lagrange multipliers on velocity level

Initialize i = 1 For t < T Evaluate left-hand side Evaluate q i 1,0 = q i 1 (50) Evaluate v i 1,0 = v + i 1 (
l + N i 1 ,I i 1,0 1 proj IR + 0 ✓ l + N i 1 ,I i 1,0 1 rW T N,I i 1,0 1
(q i 1 + Dt i v + i 1 )v i 1,1 ◆ = 0 l + T i 1 ,I i 1,0 1 proj C T (l + N i 1 ,I i 1,0 1 ) ✓ l + T i 1 ,I i 1,0 1 rW T T,I i 1,0 1
(q i 1 + Dt i v + i 1 )v i 1,1 ◆ = 0 with v i 1,1 = v + i 1 + Dt i M 1 i 1 h h + i 1 +W N i 1 l + N i 1 +W T i 1 l + T i 1 i Compute q i 1,1 = q i 1 + Dt i 2 {v i 1,0 + v i 1,1 } (51) Decide strategy if @k ⇤ : g N k ⇤ q i 1,0 > 0 ^gN k ⇤ q i 1,1  0 (14) Set L N i = 0, L T i = 0 Compute Lagrange multipliers on velocity level l N i ,I i 1,1 1 proj IR + 0 ✓ l N i ,I i 1,1 1 rW T N i ,I i 1,1 1 v i ◆ = 0 l T i ,I i 1,1 1 proj C T (l N i ,I i 1, 1 1
)

✓ l T i ,I i 1,1 1 rW T T i ,I i 1,1 1 v i ◆ = 0 with v i = v + i 1 + Dt i 2 M 1 i 1 h h + i 1 +W N i 1 l + N i 1 +W T i 1 l + T i 1 i + Dt i 2 M 1 i 1,1 ⇥ h i 1,1 +W N i 1,1 l N i +W T i 1,1 l T i ⇤ if 9k ⇤ : g N k ⇤ q i 1,0 > 0 ^gN k ⇤ q i 1,1  0 (14) Set l + N i 1 = 0, l + T i 1 = 0 and l N i = 0, l T i = 0 Compute v i = v + i 1 + Dt i 2 n M 1 i 1 h + i 1 + M 1 i 1,1 h i 1,1 o (77) Compute L N i , L T i on velocity level (I i 1,1 1 ) State increment Compute v + i = v i + M 1 i (W N i L N i +W T i L T i ) (82) Compute q i = q i 1 + Dt i 2 {v i 1,0 + v i 1,1 } (
83) Write result of time step, e.g. q i and v + i Update t by Dt i and increment i by 1 Postprocessing Figure 9: Flowchart of the half-explicit timestepping scheme on velocity level.

RR n°8623

Third analysis of simulation results

With the half-explicit evaluation of the constraints, we get the typical results of a velocity level discretization in comparison to the acceleration level discretization (Fig. 10) discussing the slider-crank mechanism in Fig. 1

Half-explicit timestepping schemes on velocity level in the sense of [9]

For the half-explicit trapezoidal rule, we just adapt the interpretation of how to evaluate the constraints. Hence, we keep the interpretation as a time-discontinuous Galerkin method. However, we can also interpret the scheme with the notation of [START_REF] Brasey | Hem 5 user's guide[END_REF]:

q i 1,0 = q i 1 , (119)
q i 1,1 = q i 1 + Dt i v + i 1 , (120)
v i 1,0 = v + i 1 , (121)
v i 1,1 = v + i 1 + Dt i M 1 i 1 h h + i 1 +W N i 1 l + N i 1 +W T i 1 l + T i 1 i , ( 122
) 0 = l + N i 1 ,I i 1,0 1 proj IR + 0 ✓ l + N i 1 ,I i 1,0 1 rW T N,I i 1,0 1 (q i 1,1 )v i 1,1 ◆ , ( 123
) 0 = l + T i 1 ,I i 1,0 1 proj C T (l + N i 1 ,I i 1,0 1 ) ✓ l + T i 1 ,I i 1,0 1 rW T T,I i 1,0 1 (q i 1,1 )v i 1,1 ◆ (124)

Inria

Half-explicit timestepping schemes on velocity level based on time-discontinuous Galerkin methods 27

for the first stage unknowns v i 1,1 , l + N i 1 and l + T i 1 . For this interpretation, we have to change (120), e.g. to a forward Euler step. Then, we proceed:

q i = q i 1 + Dt i 2 {v i 1,0 + v i 1,1 } , ( 125
)
v i = v + i 1 + Dt i 2 M 1 i 1 h h + i 1 +W N i 1 l + N i 1 +W T i 1 l + T i 1 i + Dt i 2 M 1 i ⇥ h(q i , v i 1,1 ) +W N i l N i +W T i l T i ⇤ , ( 126
) 0 = l N i ,I i 1 proj IR + 0 ⇣ l N i ,I i 1 rW T N i ,I i 1 v i ⌘ , ( 127
) 0 = l T i ,I i 1 proj C T (l N i ,I i 1 ) ⇣ l T i ,I i 1 rW T T i ,I i 1 v i ⌘ . ( 128
)
The unknowns for the second stage are v i , l N i and l T i . As an option, we can evaluate the impact equations without taking care, because they are also stated on velocity level and, as the case may be, automatically satisfied. Concerning efficient evaluations, however, we implement a robust procedure as shown in Fig. 9.

v + i = v i + M 1 i [W N i L N i +W T i L T i ] , ( 129
) 0 = L N i ,I i 1 proj IR + 0 h L N i r ġN i ,I i 1 i , ( 130
) 0 = L T i ,I i 1 proj C T (L N i ) h L T i r ġT i ,I i 1 i . (131)
The unknowns are v + i , L N i and L T i . The "forecasting trapezoidal rule" is the easiest member of D + or D timestepping schemes [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook[END_REF]. For the other representatives the internal stages for position and velocity as well as their corresponding constraint equations are coupled. For the efficient implementation on velocity level, adaptations are necessary which may result in the loss of the interpretation as time-discontinuous Galerkin method. The trapezoidal rule as a locally second order scheme, however, suffices for most practical considerations and industrial simulations.

Multi-contact examples

In this section, we discuss several examples to show the benefits of the proposed scheme. The slidercrank mechanism is a nonlinear example and the bouncing ball in a box is a linear example. We compare computational results, convergence and computing time of the half-explicit timestepping schemes on acceleration and velocity level, as well as of the classic explicit Moreau-Jean timestepping scheme. The computations have been done with Matlab. Thereby index sets are calculated without any additional tolerances. The constraints are calculated with a tolerance of 10 12 for satisfying the constraint equations.

Slider-crank mechanism

For the slider-crank mechanism in Sect. 2.1, we give some supplementary impressions without regarding friction. The calculations have been performed with the half-explicit timestepping scheme on velocity level. Figure 11 shows the slider's movement in the notch for different coefficients of restitution. In Fig. 12, the trajectories of the three generalized coordinates and generalized velocities are depicted, respectively. Figure 13 demonstrates nicely the splitting of the different Lagrange multipliers and therefore the possibility to achieve locally a higher order discretization. (d) e N = 0.9

Figure 11: Movement of the center of gravity of the slider for different coefficients of restitution and half-explicit timestepping on velocity level.

Convergence and computing time for the bilateral case

We compare the half-explicit timestepping scheme on velocity level with the classic explicit Moreau-Jean timestepping scheme regarding a bilateral slider-crank mechanism without friction, i.e., c = 0 m. Thereby, we analyze computing time and convergence of the schemes in comparison with a corresponding Simpack model, which is integrated with the SODASRT2 solver and a high tolerance 10 12 . Independent of different time step-sizes, the relative overhead of the half-explicit timestepping scheme on velocity level in comparison with the classic explicit Moreau-Jean timestepping scheme is about 1.2. Hence, the computation time per time step is about 1.2 times larger for the half-explicit timestepping scheme on velocity level. We calculate the error with respect to the reference Simpack solution for different time step-sizes. Thereby, we compute the differences of the generalized coordinates for the considered time instances and arrange them in a matrix. The 2-norm of this matrix results in Table 2. Concerning these evaluations, this is a trend to an order of convergence of 1 for the classic explicit Moreau-Jean timestepping scheme and to an order of convergence of 2 for the half-explicit timestepping scheme on velocity level. The overall relative cost to achieve a given tolerance, e.g. 10 2 , is given by the relation of the fractions of the relative overhead and the necessary time step-size, i.e., 1.2 • 10 1 for the half-explicit timestepping scheme on velocity level in comparison to the classic explicit Moreau-Jean timestepping scheme. Thus, it makes sense to apply higher-order approximations for appropriate

Convergence and computing time for the unilateral case

We compare the half-explicit timestepping scheme on velocity level with the classic explicit Moreau-Jean timestepping regarding a unilateral slider-crank mechanism with friction. Thereby, we analyze computing time and convergence of the schemes assuming a reference solution given by a simulation of the halfexplicit timestepping on velocity level with Dt = 10 7 s.

The relative overhead of the half-explicit timestepping scheme on velocity level is about 1.15. For the error, we get Table 3. The order of convergence concerning these evaluations drops down to 1 for both integration schemes. However, the results for the half-explicit timestepping scheme on velocity level are better perhaps because of the possibility to achieve a locally higher order. The overall relative cost for a given tolerance, e.g. 10 2 , is 1.15 • 10 1 for the half-explicit timestepping scheme on velocity level in comparison to the classic explicit Moreau-Jean timestepping scheme. Hence like in the bilateral case, the relative overhead could be compensated by an adaptive time-step size and a control of the accuracy.

Dt [s]

Bouncing ball in a box

The linear bouncing ball in a box example explains the difference between the half-explicit timestepping schemes on acceleration and velocity level, i.e., the drift-off effect, and the energy behavior of halfexplicit timestepping schemes based on the trapezoidal rule.

A bouncing ball of radius R in a planar box of width a and height b (Fig. 14) is described by the two translational coordinates, x and y, of its center of gravity. The rigid ball with mass m is subject to gravitation g in cos(p/6) sin(p/6) T -direction. The normal gap functions g N are shown in Fig. 14.

The vector containing the coefficients of restitution is e N . Friction is not considered. Fixed characteristics are given in Table 4.

The generalized mass matrix satisfies

M : IR 2 ! IR 2,2 , q 7 ! M (q) = ✓ m 0 0 m ◆ , (132)
and the linear generalized force can be obtained by The normal gap functions are

h : IR 2 ⇥ IR 2 ! IR 2 , (q, v) 7 ! h (q, v) = g ✓ cos(p/6) sin(p/6) ◆ . ( 133
) Inria y x R a b g N 3 g N 1 g N 4 g N 2
g N 1 (q) = y R , (134)
g N 2 (q) = a x R , (135)
g N 3 (q) = b y R , (136)
g N 4 (q) = x R , (137)
resulting in a Boolean matrix for the generalized force directions. Fig. 15 shows the trajectory of the ball starting from its initial position and ending in the lower left corner. The position and velocity curves are depicted in Fig. 16 and clearly show the nonsmooth, i.e., impulsive behavior. For the application of the half-explicit timestepping methods based on the trapezoidal rule to the slider-crank mechanism, we have mentioned a not always decreasing energy trend (Fig. 6). As we see in Fig. 17, the energy behavior for the linear bouncing ball in a box example actually shows the expected curve without any increasing periods. The previous non-decreasing trend is traced back to the nonlinear nature of the example. In Fig. 18, we can convince ourselves from the drift-off effect for the acceleration level discretization. We have to remark, that the half-explicit timestepping on acceleration level just evaluates the contact forces similar to the explicit timestepping on velocity level. This procedure is more a calculation of equivalent forces involving a projection than a computation of contact forces while adapting the unknowns of the system. That is why, we see the drift-off effect also for the linear bouncing ball in a box example. All calculations have been performed with Dt = 5 • 10 3 s.

Flexible examples

The elastic bar is a linear but flexible example and the rotor is a linear and flexible example from practice.

In this section, we discuss the behavior of waves and damping for half-explicit timestepping schemes on velocity level.

Elastic bar

In this section, we consider the classical example of a linear elastic bar that impacts a rigid obstacle at constant velocity v 0 . The example is depicted in Figure 19. The bar has a length L, a section area S, a density r and a Young modulus E. If we denote by c 0 = p E/r the longitudinal wave speed, the time of contact is twice the time of the traveling of elastic waves in the bar, that is T = 2L/c 0 . Within the contact time, the contact force is constant and equals t = ESv 0 /c 0 . The parameters for this example are summarized in Table 5 stiffness matrices are

M e = 1 6 rSl e  2 1 1 2 , K e = 1 l e ES  1 1 1 1 , (138)
where l e = L/N is the length of an element. The first simulation results with the half-explicit timestepping scheme on velocity level for a time-step size Dt 10 7 s and 50 finite elements are reported in Fig. 20. This simulation has been carried out with Siconos, an open-source software for the modelling and the simulation of nonsmooth dynamical systems [START_REF] Acary | An introduction to Siconos[END_REF]. In this simulation, the expected solution for the elastic bar is found exaclty with Lemke's algorithm. The time of contact of the bar is well-approximated and the velocity within the contact interval is maintained at zero. The violation of the constraint is equal to 5 • 10 9 m. The presence of the impulse is a consequence of the space discretization by finite elements. The finite mass of the node in contact is instantaneously stopped at the impact. This generates an impulse. This impulse associated with a plastic impact induces a loss of the total energy at impact which can be observed in the energy plot. Finally, we remark that the reaction force is well approximated in the first half of the contact period before the development of an instability which is inherent to standard finite element approximations of a travelling velocity jump trough the mesh. Nevertheless, it is noteworthy that the contact is never lost due to artificial numerical chattering.

In Figure 21, a similar simulation is reported but the number of elements is increased to 1000. We observe that the reaction impulse is reduced due to the decreased mass of the node in contact. The energy loss is also reduced. The instability is nevertheless still visible but its frequency has increased with the inverse of the element size. The travelling of the impulse excites the spurious high frequency mode associated with the discretization. One of the usual ways to circumvent this problem is to introduce an artificial damping in the simulation. Since the discretization family presented in this paper does not include numerical damping-even if it seems possible, it is left as further work-, we have added a small amount of Rayleigh damping (proportional to the stiffness C e = d K e ) to damp out this numerical artifact. In Figure 22, the results with the viscous damping are plotted. The stiffness proportional damping coefficient is equal to d = 10 3 /E s = 4.761 • 10 6 s. We observe that the instability is nearly damped out and the post velocity is smooth. Only a remaining impulse travels trough the bar. The decay of the total energy can be If we compare our results with those in [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF] for the standard Moreau-Jean timestepping scheme for the example of the elastic bar, we observe the same type of instability of the reaction forces when there is no numerical damping in the Moreau-Jean algorithm. Nevertheless, the proposed scheme is mainly interesting for two reasons. First, it enables a splitting between the reaction force and the artificial reaction impulse due to the space discretization. The value of the reaction force and the stress in the bar can be used for design as a good approximation of the forces in the space continuous bar. The magnitude of the reaction impulse gives a measure of the quality of the mesh with respect to the contact representation. Second, the scheme is half-explicit and thus requires less computational effort for the evaluation of the constitutive behavior in the nonlinear setting. However, we have to notice that a CFL type condition has to be satisfied since the internal forces are explicitly evaluated.

Rotor

A rotor test rig at the Institute of Applied Mechanics of the Technische Universität München is modelled in the multibody simulation framework MBSim [mbsim, 52] (Fig. 23). The idea of this example is the implementation of a variant of the half-explicit timestepping scheme on velocity level in MBSim and the application to an industrial example. Thereby, we test a simple and heuristic adaptation of the time-step size.

The rotor in Fig. 23 consists of drive system with a driving torque about the rotor's longitudinal axis T D and a failure given by a torque about the vertical axis T F . The drive system is mounted by a springdamper element with constants c D , d D . The rotor's axle is modelled with N A beam elements in a floating frame of reference framework. Each finite element node has five degrees of freedom, two small out-off plane translational deflections, two small out-off plane rotational deflections and the rotation about the rotor's longitudinal axis. The length of the axle is l A , its radius is r A , the density is r A , the Young's modulus is E A , and the shear modulus is G A . A mass proportional damping d 1 A and a torsional damping d 2 A are used. The fly-wheel is a rigid body with originally six degrees of freedom, which is mounted bilaterally at position x F from the axle's left end on the rotor's axle. It has the mass m F and the rotational inertia about the rotor's longitudinal axis J F . The journal bearing consists of a rigid frustum with six degrees of freedom mounted bilaterally in its center on the left end of the rotor's axle but with negligible inertia values. This frustum of radius r JF and half width b JF moves freely in a circular bearing with radius r JB , i.e., backlash and dry friction µ. The bearing is a rigid body with six degrees of freedom, mass m JB and negligible rotational inertia values. The bearing is mounted in its center by a spring-damper element with constants c JB , d JB . The actual contact may only occur at position x JB from the axle's left end. We consider the bearing, e.g., as an auxiliary bearing occurring in many rotor systems with active magnetic bearings [START_REF] Ginzinger | Model-based condition monitoring of an auxiliary bearing following contact events[END_REF] or as an example for a safety bearing. The whole system is subject to gravitation g in negative vertical direction. The characteristics are summarized in Table 6. Damping is added because of the experience in Sect. 7.1.

For an implementation of a variant of the half-explicit timestepping scheme on velocity level in MB-Sim, especially we have to discretize

q = Y (q)v (139)
instead of q = v in the sense of [START_REF] Brasey | Hem 5 user's guide[END_REF]. At the same time, we adapt the interpretation of the contact activity for the first stage unknowns. Instead of (123)-(124), we solve

0 = l + N i 1 ,I i 1,1 1 proj IR + 0 ✓ l + N i 1 ,I i 1,1 1 rW T N,I i 1,1 1
(q i 1,1 )v i 1,1 ◆ , 0 = l + T i 1 ,I i 1,1 1 proj C T (l + N i 1 ,I i 1,1 1 ) ✓ l + T i 1 ,I i 1,1 1 rW T T,I i 1,1 1
(q i 1,1 )v i 1,1 ◆ (141)
in a more natural way. An overview is depicted in Fig. 24. The time-step size for impulsive periods Dt is heuristically adapted concerning the experiences with the slider-crank example (Table 3): Dt = Dt 10 . As we do not know the error constant, we prefer this heuristic in comparison to Dt = Dt 2 , which uses the non-impulsive order of the schemes.

The typical forward whirl phenomenon because of the external loadings and the dry friction contact can be seen in Fig. 25. We are interested in the general behavior of the variant of the half-explicit timestepping scheme on velocity level. In MBSim, the tolerances define tubes around the corner laws [START_REF] Förg | Mehrkörpersysteme mit mengenwertigen Kraftgesetzen : Theorie und Numerik[END_REF]. Starting from the impulsive tolerances tol L = 10 10 Ns, tol ġ = 10 10 m/s, we use tol l = tol L Dt , tol g = tol ġ

Dt for the calculation of the constraint equations with a single-step fixed-point method [START_REF] Förg | Mehrkörpersysteme mit mengenwertigen Kraftgesetzen : Theorie und Numerik[END_REF]. The index sets are detected with a tolerance of 10 8 m. The variant of the half-explicit timestepping scheme on velocity level uses a fraction of impulsive integration steps of 1.3%. The overall computing time of the variant of the half-explicit timestepping scheme on velocity level in comparison with the half-explicit variant of the Moreau-Jean timestepping scheme in [mbsim] is 0.32. Hence, the new proposed scheme offers huge Inria Half-explicit timestepping schemes on velocity level based on time-discontinuous Galerkin methods 37

Summary and Conclusion

We have presented a time-discretization scheme for the simulation of nonsmooth mechanical systems with friction and impacts as well as rigid and flexible bodies with all possible interactions. Without switching between impulsive and non-impulsive periods, the time-discretization scheme is both consistent and efficient. It represents impulses, if necessary, and automatic local order elevation for state variables, if possible, at the same time and is unique in the literature concerning these criteria. Constraints are formulated on velocity level in an implicit way using an augmented Lagrangian technique with semi-smooth Newton schemes without penetration; all other evaluations are explicit, which yields a half-explicit method. The scheme is derived and explained in detail concerning different criteria and it is applied to various examples with multi-contact, flexibility and industrial relevance. As a compromise between classic timestepping schemes and event-driven integration in the context of time-discontinuous Galerkin methods finally, two half-explicit timestepping schemes on velocity level distinguish impulsive and non-impulsive periods for all examples in a geometrically consistent way and reduce the computing time significantly. For flexible examples, we recognize the drawback of explicit schemes, i.e., stability issues, which can be solved by applying schemes from the Newmark family as base integration schemes

Figure 2 :

 2 Figure 2: Definition of the gap functions for the slider-crank mechanism with unilateral constraints and friction [22].

1 = 7 . 4 • 10 5 kgm 2 J 2 = 5 . 9 • 10 4 kgm 2 J 3 = 2 . 7 •

 17422592327 An event-driven scheme discretizes the non-impulsive equations (1)-(4) by a classic numerical rule [3, Chap. 8]. At the same time, it observes the gap functions for open (g N > 0)-close (g N  0) transitions Geometrical characteristics l 1 = 0.1530 m l 2 = 0.3060 m a = 0.0500 m b = 0.0250 m c = 0.0010 m Inertia properties m 1 = 0.0380 kg m 2 = 0.0380 kg m 3 = 0.0760 kg J 10 6 kgm 2 Force elements g = 9.81 m/s 2 Contact parameters e N 1 = e N 2 = e N 3 = e N 4 = 0.4 e T 1 = e T 2 = e T 3 = e T 4 = 0.0 µ 1 = µ 2 = µ 3 = µ 4 = 0.01 Initial conditions

1 PostprocessingFigure 4 :

 14 Figure 4: Flowchart of the half-explicit timestepping scheme on acceleration level.

4 Figure 5 :

 45 Figure 5: Normal gap functions of the slider for half-explicit timestepping on acceleration level.

Figure 6 :

 6 Figure 6: Total energy of the slider-crank mechanism for half-explicit timestepping on acceleration level.

4 Figure 8 :

 48 Figure 8: Normal gap functions of the slider for explicit timestepping on velocity level.

4 Figure 10 :

 410 Figure 10: Normal gap functions of the slider for half-explicit timestepping on acceleration (dotted) and velocity level (solid).

3 Figure 12 : 4 Table 2 :

 31242 Figure 12: Angles and angular velocities of the slider-crank mechanism for half-explicit timestepping on velocity level.

3 Table 3 :

 33 classic explicit Moreau-Jean timestepping scheme 7.1 • 10 1 1.5 • 10 1 4.5 • 10 2 error of the half-explicit timestepping scheme on velocity level 4.7 • 10 1 1.8 • 10 2 3.3 • 10 Unilateral slider-crank example: comparison of the error for different time step-sizes.

4 Figure 13 :

 413 Figure 13: Lagrange multipliers of the slider-crank mechanism for half-explicit timestepping on velocity level.

Figure 14 :Table 4 :

 144 Figure 14: Bouncing ball in a box. Geometrical characteristics a = 4.0m b = 4.0m R = 1.0m Inertia properties m = 1.0 kg Force elements g = 9.81 m/s 2 Contact parameters e N 1 = e N 2 = e N 3 = e N 4 = 0.3 Initial conditions x 0 = 2.0m y 0 = 2.0m ẋ0 = 0.0m/s ẏ0 = 0.0m/sTable 4: Characteristics of the bouncing ball example.

Figure 15 :

 15 Figure 15: Trajectory of the bouncing ball in a box for half-explicit timestepping on velocity level.

Table 5 :Figure 16 :

 516 Figure 16: Position and velocity of the bouncing ball in a box for half-explicit timestepping on velocity level.

Figure 17 :

 17 Figure 17: Total energy of the bouncing ball in a box for half-explicit timestepping on velocity level.

Inria

4 Figure 18 :

 418 Figure 18: Normal gap functions of the bouncing ball in a box for half-explicit timestepping on acceleration (dotted) and velocity level (solid).

Figure 19 :

 19 Figure 19: Linear elastic impacting bar on a rigid obstacle.

  kinetic, potential and total energy

Figure 20 :

 20 Figure 20: Elastic bar simulation with the half-explicit timestepping scheme on velocity level Dt = 1 • 10 7 s, N = 50.

Figure 21 :

 21 Figure 21: Elastic bar simulation with the half-explicit timestepping scheme on velocity level Dt = 1 • 10 7 s, N = 1000.

Figure 22 :

 22 Figure 22: Damped elastic bar simulation with the half-explicit timestepping scheme on velocity level Dt = 1 • 10 7 s, N = 1000.

Figure 23 :

 23 Figure 23: Model of a rotor test rig.

Figure 25 :

 25 Figure 25: Curve of the journal frustum in the journal bearing.

  )

		7.5						Total energy
												energy
		7.48
		7.46
	Kinetic and potential energy [J]	7.38 7.4 7.42 7.44
		7.36
		7.34
		7.32	0	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.1
								Time [s]

Table 6 :

 6 Characteristics of the rotor example.

RR n°8623

http://www.netlib.org/odepack/

http://www.netlib.org/ode/ Inria

Acknowledgement

JK gratefully acknowledges the UK's Engineering and Physical Sciences Research Council (EPSRC) for funding through a studentship.

Specify characteristics, start time t = 0, end time T , time step size Dt Initialize i = 1 For t < T Evaluate left-hand side Evaluate q i 1,0 = q i 1 Evaluate v i 1,0 = v + i 1

Predict right-hand side

Compute Lagrange multipliers on velocity level

o Compute Lagrange multipliers on velocity level

Write result of time step, e.g. q i and v + i Update t by Dt i respectively Dt i and increment i by 1 Postprocessing Figure 24: Flowchart of the half-explicit timestepping scheme on velocity level, variant for MBSim.

[22] Paulo Flores, Remco Leine, and Christoph Glocker. Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach. Multibody