A nonsmooth generalized‐ alpha method for mechanical systems with frictional contact

In this article, the existing nonsmooth generalized‐ α method for the simulation of mechanical systems with frictionless contacts, modeled as unilateral constraints, is extended to systems with frictional contacts. On that account, we complement the unilateral constraints with set‐valued Coulomb‐type friction laws. Moreover, we devise a set of benchmark systems, which can be used to validate numerical schemes for mechanical systems with frictional contacts. Finally, this set of benchmarks is used to numerically assert the properties striven for during the derivation of the presented scheme. Specifically, we show that the presented scheme can reproduce the dynamics of the frictional contact adequately and no numerical penetration of the contacting bodies arises—a big issue for most popular time‐stepping schemes such as the one of Moreau. Moreover, we demonstrate that the presented scheme performs well for multibody systems containing flexible parts and that it allows general parametrizations such as the use of unit quaternions for the rotation of rigid bodies.


Introduction
In this paper, we derive a nonsmooth generalized-α method for the simulation of mechanical systems with frictional contact. Moreover, we introduce a set of benchmark systems to validate the presented scheme.
In many engineering applications, systems are modeled through rigid and flexible bodies, which are interconnected by joints and can come into contact with each other or their surroundings. Prominent examples are automotive and robotic systems. The dynamics of mechanical systems with frictional contact can be described within the theory of nonsmooth mechanics [21,33], where the velocities of the system are allowed to jump. This is particularly important for the description of contact between rigid bodies, where at the time instant when the contacting bodies touch, an impact might occur, which due to rigidity instantaneously changes the velocity of the body.
For the simulation of nonsmooth mechanical systems, two approaches can be distinguished, the event-driven and the event-capturing schemes. The event-driven schemes use standard ordinary differential equations (ODE) solvers, or differential algebraic equations (DAE) solvers, to compute the impact-free motion. Every time an impact event is detected, the algebraic impact equations are solved to find the post-impact velocities, which are then used to continue the integration with the ODE solver. The main strength of event-driven integration is that ODE solvers with high-order accuracy can be used. However, since every impact is resolved, these schemes are not suitable to find motions with accumulation points, i.e., motions with an infinite number of impacts occurring in a finite interval of time.
In contrast to event-driven schemes, the event-capturing schemes, also called time-stepping schemes, can overcome accumulation points, because they smear the effects of impacts over a time step. However, the most widespread schemes, such as Moreau's time-stepping scheme [33] and variants thereof, are of first-order accuracy, show a high level of numerical dissipation and allow penetration of the contacting bodies. These properties are problematic, especially for the simulation of mechanical systems containing elastic parts. Several publications present improved event-capturing schemes by addressing at least one of these drawbacks. In [10,31], integrators similar to Moreau's time-stepping scheme with improved long-term energy behavior have been derived. To arrive at event-capturing schemes with higher accuracy order for the impact-free motion, [1] combines high-order Runge-Kutta methods with Moreau's time-stepping scheme, [40] relies on extrapolation methods and [38] uses a discontinuous Galerkin method to discretize the dynamics of the mechanical system. The constraint drift, which is also responsible for contact penetration, is generally solved by a stabilization in the sense of Gear-Gupta-Leimkuhler (GGL) [18] and/or a projection approach, see for example [2,39].
In [8,9,11,14], the nonsmooth generalized-α methods were introduced, which alleviate many of the problems of the Moreau-type time-stepping schemes. In particular, the generalized-α schemes are second-order accurate for the impact-free motion, use the GGL stabilization to avoid penetration of the contacting bodies and it is known from structural mechanics applications, see [7,12,26], that generalized-α schemes perform well for flexible multibody systems. However, the nonsmooth generalized-α methods of [8,9,11,14] are only applicable to multibody systems with frictionless contacts, i.e., unilateral constraints without Coulomb friction. Furthermore, these nonsmooth generalized-α methods are restricted to mechanical systems for which the velocity of the system corresponds to the time derivative of the position coordinates. In [13,17], the schemes of [8,9,11,14] have been extended to cope with systems with frictional contact and a more general kinematic equation.
In [17], the augmented Lagrangian approach together with discrete set-valued Coulomb friction laws on position and velocity level have been used to arrive at a nonsmooth generalized-α scheme, which can describe systems with frictional contact. The scheme of [17] has been extended in [13] to allow for different rotation parametrizations for rigid bodies. This represents a particular case of a general kinematic equation since the angular velocity is not the time derivative of the rotation matrix describing the orientation of a rigid body. The nonsmooth generalized-α methods of [8,9,11,13,14,17] are derived from a splitting strategy, where the motion is artificially split into a nonsmooth and a smooth part. Furthermore, the contact laws describing the impenetrability of the contacting bodies (Signorini condition) as well as the impact between these bodies (Newton's impact law) are introduced either as inequality complementarity conditions or equivalently as normal cone inclusions. The discrete counterparts of the contact laws are given as an active set formulation, i.e., depending on the state of the system, conditions are established that decide which set of nonlinear equations describes the contact. This set of nonlinear equations is finally solved by a semi-smooth Newton method in every time step.
In nonsmooth mechanics, the set-valued Coulomb friction law is naturally stated as a force law on velocity level. Since in the case of sticking this friction law acts like a bilateral constraint, it can be brought to acceleration level through differentiation. Similar to [13,17], it is the aim of this paper to extend the nonsmooth generalized-α schemes [8,9,11,14] to account for friction as well as a general kinematic equation. However, the explained link between the friction law on velocity and acceleration level will be exploited instead of formulating the discrete friction laws on position and velocity level as done in [13,17].
We use the same generalized-α time discretization for the dynamics of the mechanical system as the schemes of [8,9,11,14], but include the general kinematic equation. In doing so, the presented scheme shares the beneficial properties of the existing generalized-α schemes, such as the second-order accuracy and the absence of numerical contact penetration. To account for friction, we invoke set-valued Coulomb-type friction laws introduced as normal cone inclusions. We exploit the fact that also the Signorini condition and Newton's impact law can be formulated as normal cone inclusions, such that the contact laws are given as a set of normal cone inclusions. We directly discretize these contact laws to arrive at their discrete counterparts, which are again a set of normal cone inclusions. This opens the possibility to numerically solve the time step using any strategy suitable for the solution of normal cone inclusion problems, see for example [3]. We present two such strategies, which both rely on a reformulation of the normal cone inclusion invoking the proximal point function to arrive at implicit nonlinear equations. For the first strategy, the proximal point function is used to achieve an active set formulation of the contact laws, which are then solved by a semi-smooth Newton method. The second strategy solves the nonlinear equations including the proximal point functions with fixed point iterations. Finally, we introduce a set of benchmark systems, which can be used to validate numerical schemes for mechanical systems with frictional contacts. Therefore, we devise the systems in such a way, that each benchmark can be used to test the performance of the scheme with respect to a specific feature. These are: overcoming accumulation points, slip-stick transitions, handling the presence of linearly dependent force directions, contact penetration, the Painlevé paradox, the suitability of the scheme for the simulation of flexible multibody systems, combined spatial friction laws and general kinematic equations. We use these benchmark systems, to validate the performance of the presented scheme.
The mathematical concepts used for the discretization of the contact laws as well as their numerical treatment is presented in Section 2. In Section 3, the equations of motion of a mechanical system subjected to ideal bilateral constraints as well as frictional contacts are established. To avoid that the resulting generalized-α scheme exhibits numerical constraint drift, and hence penetration of the contacting bodies, a constraint stabilization in the sense of Gear-Gupta-Leimkuhler is introduced in Section 4. Subsequently, the generalized-α discretization of the dynamics is presented in Section 5, the contact laws in normal direction of the contacting surfaces of the bodies are discretized in Section 6 and the discretization of the friction laws is conducted in Section 7. The extension to general relations between velocities and position is discussed in Section 8. The reformulation of the discrete contact laws as nonlinear equations and how these are solved with a semi-smooth Newton method or fixed point iterations is shown in Section 9. Finally, a set of mechanical benchmark systems is used to validate the derived nonsmooth generalized-α scheme in Section 10.

Convex sets and cones
A set C ⊆ R n is convex if and only if α a + (1 − α)b ∈ C for all a, b ∈ C and for all α ∈ [0, 1]. For λ ∈ R, we introduce the scalar multiple of a set C ⊆ R n as λC = {λx | x ∈ C}. If C is convex then it holds that for α 1 , α 2 ≥ 0. More generally, it holds for a closed convex set C that where α i ≥ 0. Let K ⊆ R n be a cone, i.e., the set K fulfills the cone property For a closed cone K it holds that when a ∈ K also α a ∈ K for all α ≥ 0. Moreover, if K is a closed convex cone, then it holds that and, likewise, where I ⊆ R and the Lebesgue measure dt may be generalized to any positive measure.

Properties of normal cones
Let C ⊆ R n be a closed convex non-empty set. The normal cone of C at x ∈ C is defined as whereas N C (x) is empty if x / ∈ C. One can show that the normal cone N C (x) is a closed convex cone, see [36]. It can easily be verified, that N C (x) = {0} if x is an interior point of the set C. Indeed, for an interior point x there exists a ball around x contained in C. This allows to choose x * from this ball such that any direction x * − x can be produced in (6), which leaves only y = 0 to fulfill the inequality.
The normal cone has a number of less known scaling properties needed in this paper. These may conveniently be proved using topics from Convex Analysis such as the subdifferential and the support function. However, here we will try to introduce as little machinery as possible and derive these properties by using only the definition (6). Proposition 1. Let α > 0 and C be a closed convex non-empty set. It then holds that Proof. Consider x ∈ αC. The definition (6) applied on the inclusion y ∈ N αC (x) implies that y T (x * − x) ≤ 0, ∀x * ∈ αC. We divide by α > 0 giving Furthermore, x / ∈ αC implies N αC (x) = ∅, which in turn yields N C ( 1 α x) = ∅. The other direction of the implication follows by taking the reciprocal value of α. Proposition 2. Let C be a closed convex non-empty set and C(α) = αC for all α ≥ 0. If where α i ≥ 0, then it holds that Proof. Retaining the case α i = 0 for the end of the proof, we first consider α i > 0 for all i. By using Proposition 1 together with the homogeneity Without loss of generality, consider i = 1, 2 and use the definition (6) for either, that is, We now sum the inequalities, but multiply by α 1 and α 2 respectively to arrive at where we used property (1) of a convex set C. Repeated summation gives which is equivalent to If α i = 0 for some i, then it holds that C(α i ) = C(0) = {0} and x i = 0 whereas N {0} (0) = R n . One easily verifies that the proposition still holds.
Proposition 2 may also be written in integral form as where α(t) ≥ 0 and where y is a fixed value, i.e., does not depend on t.

Reformulations and numerical solution of normal cone inclusion problems
The generalized-α discretization of the dynamics of mechanical systems with frictional contacts results in a normal cone inclusion problem, that is, a set of nonlinear equations and normal cone inclusions, describing a time step of the scheme. Hence, finding a numerical solution to a time step of the presented scheme in essence corresponds to finding a numerical solution to where R s : R n × R n → R n is an implicit relationship between y and z. An extensive collection of solutions strategies for (18) can be found in [3]. In what follows, we limit ourselves to two strategies based on the reformulation of the normal cone inclusion as an equation including the proximal point function.
For a closed convex non-empty set C ⊆ R n , we define the proximal point function which maps a point p to the closest point q ∈ C, where the distance between the points is measured by the Euclidean norm || . ||. It is immediately clear from (19) that prox C (p) = p if and only if p ∈ C. Using the just defined proximal point function, it can be shown that two points x and y fulfill the normal cone inclusion if and only if they fulfill the equation for any r > 0, see [29].
Using the equivalence of (20) and (21) allows to reformulate (18) as which reduces the problem of finding the solution of (18) to numerically finding the solution of a nonlinear equation.
The equation (22) can for example be used to define a residual R = Due to the presence of prox C , the residual R is continuous but nonsmooth, hence (23) can be solved using the semi-smooth (nonsmooth) Newton method [3,4], where in the Newton update any regular element ∇R(x ν ) of the generalized Jacobian ∂R(x ν ) can be used. The semi-smooth Newton method reduces to the Newton method whenever R is differentiable. In that case, ∇R(x ν ) is the Jacobian matrix of R at x ν . Solving (23) with the semi-smooth Newton method has two drawbacks if ∇R(x ν ) is computed by finite differences. Firstly, the accuracy of the Jacobian strongly depends on the parameter r. Moreover, using finite differences to compute ∇R(x ν ) at a point where R is not differentiable generally leads to an arbitrary element in ∂R(x ν ). Hence, we have no control over which element of ∂R(x ν ) is ultimately used for the Newton update.
As a remedy, we can use an active set formulation to treat (23), which relies on the specific knowledge of the set C. In our case, two sets are of interest. Either C = R − 0 := {z ∈ R | z ≤ 0 } or C = B(R) := {z ∈ R n | ||z|| ≤ R}, for which the respective proximal point functions are piecewise given as It follows immediately from using (25) in (23) that for C = R − 0 the residual in (23) is equivalent to where we abstain from the bold notation as all variables are scalar. Proceeding all the same for C = B(R), the residual in (23) is equivalent to where we have used that z + R ry−z ||ry−z|| = 0 is solved by z = −R y ||y|| . Furthermore, ry−z / ∈ B(R) implies with z ∈ B(R) that ||y|| > 0. In these active set formulations of (23) the parameter r is only present in the activation condition of the residuals and does not affect the accuracy of the computation of the Jacobian by finite differences. Moreover, since the active set formulation gives direct access to the non-differentiable points of the residual, we can choose to not switch between the two pieces of the residual during the numerical differentiation process, which leads to a well-behaved Jacobian and a more robust scheme.
Another popular strategy to solve (22) is the use of fixed point iterations, see [4,40]. For that, we use the implicit function theorem on R s (y, z) = 0 guaranteeing the existence of a function F such that y = F(z). Hence, the fixed point iterations follow as where generally F is not known analytically and y ν = F(z ν ) must be found numerically by solving R s (y ν , z ν ) = 0 for y ν while treating z ν as a constant.

Mechanical systems with frictional contacts
Consider a finite-dimensional mechanical system whose motion is described by the generalized coordinates q(t) ∈ R n which are considered to be functions of time t. We introduce the generalized velocities u corresponding to the time derivativeq of the generalized coordinates for almost all t, i.e.,q = u for almost all t. Equivalently, we can write We assume the generalized velocities to be special functions of locally bounded variation, i.e., functions of locally bounded variations with no singular part, see [6,15]. This implies that the left and right limits of u, respectively denoted as u − and u + , exist and are bounded at every time instant, and that the discontinuity points of u are countable. In order to have the notion of velocity for every time instant, we set u = u + , i.e., we consider the velocities to be right-continuous. It is well known, see [5,21], that the differential measure du can be decomposed into the sum of an absolutely continuous measure and a singular measure with respect to the Lebesgue measure dt.
To state the decomposition, we define the atomic measure η as a finite sum of Dirac point measures δ t k . Specifically, With that, the velocity measure can be decomposed as where we call the density a the generalized acceleration of the mechanical system and where dη is concentrated on the set of discontinuities t k of u.
To complete the description of the dynamics of the mechanical system, the link between the change in velocity du and the forces acting on the system is established by the equality of measures where M = M T denotes the mass matrix of the system and the forces are represented by the force measure dF. We assume that the mechanical system is subjected to n g + n γ ideal bilateral constraints formulated at position and velocity level, respectively, as g(t, q) = 0 ∈ R ng and γ(t, q, u) = 0 ∈ R nγ .
By ideality of the constraint forces [21], the corresponding force directions are given by such that the constraint forces are where dP g and dP γ denote the constraint percussion measures. Assuming that percussions P 2 are special functions of locally bounded variation, the force measures dP 2 are composed by a nonimpulsive force λ 2 and an impulsive force Λ 2 , that is Hereby the box 2 is used as a placeholder for any subscript, e.g., for g or γ.
To model the contacts occurring in the mechanical system, we assume that they can be described by n N ideal unilateral constraints at position level where the inequality holds component-wise and g N (t, q) ∈ R n N are the gap functions describing the distance between the tangent planes of the pairs of contact points on either contacting bodies, see [21,30]. The corresponding constraint forces are For the contact (pair) k we say that the contact is , the contacting bodies are separated • closed if g k N = 0, i.e., the contacting bodies are in contact • penetrating if g k N < 0, i.e., the contacting bodies penetrate each other. We define the set of active contacts as It is clear from this definition, that the set of inactive contacts, that is, the complementĀ = {1, . . . , n N } \ A of the set A, is the set of open contacts. The constraint force laws describing a contact in normal direction are formulated separately for the nonimpulsive forces λ N and the impulsive forces Λ N composing the normal contact percussions dP N by (36). For the nonimpulsive contact force λ k N of the k-th contact we choose Signorini's law which assures contact impenetrability and is also known as the Signorini condition. Hereby, we assume that the contact surfaces can only exert compressive normal contact forces λ k N on each other. The physical interpretation of Signorini's law is easiest understood by looking at (40) as an inequality complementarity condition Indeed, if λ k N > 0, the argument of the normal cone is an interior point of the set R − 0 and the normal cone is zero, implying g k N = 0. On the other hand, if λ k N = 0, it follows from the definition of the normal cone (6) that N R − 0 (0) is the set of positive numbers including zero and therefore g k N ≥ 0. It is easy to see from Signorini's law in the form (41) that for open contacts, the normal contact force is zero and if the contact is closed, only forces are allowed which push the contact surfaces apart.
For the impulsive contact forces Λ N , we use the gap velocity defined bẏ to formulate the Newton-type impact law as Herein, we have introduced the kinematic quantity for the k-th contact with restitution coefficient e k N . The impact law (43) implies that whenever an impact takes place, i.e., Λ k N > 0, the post-impact velocityġ k+ N =ġ k N (t, q, u + ) is related to the similarly defined pre-impact velocity by Newton's impact laẇ More details on the intricacies of the generalized Newton's impact law (43) can be found in [22,41]. Friction between the surfaces of the k-th contact is described by a set-valued force law. Therefore, n k F velocity parameters γ k F (t, q, u) ∈ R n k F describing the relative motion of the surfaces are typically introduced, where n k F = 1 for planar friction, n k F = 2 for spatial friction and n k F ≥ 3 for combined friction laws such as Coulomb-Contensou friction [28]. The friction forces have the form With C F denoting the set of admissible (negative) friction forces, the constitutive laws for the nonimpulsive and impulsive friction forces of the k-th contact are whenever the contact k is active, i.e., k ∈ A. Moreover, the set C F depends on the normal contact forces and we have introduced with restitution coefficient e k F . Otherwise, if the contact k is open, the friction forces are zero.
For the sake of simplicity, in this paper we restrict ourselves to sets of admissible (negative) friction forces having the form with friction coefficient µ k and mention how the general case can be treated wherever it seems appropriate. The most prominent example having the form (49) is isotropic spatial Coulomb friction. As described in [30], in that case we have n k F = 2 since the velocity parameters correspond to two orthogonal velocities which are tangent to the contact surfaces.
Under the premise that all forces which are not constraint or contact forces, such as spring forces, gyroscopic terms and dashpot forces, are nonimpulsive and can therefore be represented by a Lebesgue-density h(t, q, u), the totality of forces acting on the mechanical system is represented by Consequently, the equality of measure (32) takes the form With (31) and (36) in mind, we conclude from (51) that the acceleration of the system is characterized by the equations of motion holding for dt-almost everywhere in time and the velocity jumps are given by the impact equations

Acceleration level constraints and stabilization
In this section, we formulate the bilateral constraints as well as Signorini's law on acceleration level. Moreover, a stabilization on velocity and position level in the sense of Gear-Gupta-Leimkuhler is introduced to avoid constraint drift in the numerical scheme.
In order to formulate the bilateral constraints (33) on acceleration level, we introduce the constraint velocitẏ as well as the constraint accelerations With those, the bilateral constraints (33) are equivalently formulated by demandingg (t, q, u, a) = 0 andγ(t, q, u, a) = 0 (56) for dt-almost everywhere as well aṡ g(t, q, u + ) = 0 and γ(t, q, u + ) = 0 (57) whenever a velocity jump occurs. Clearly, these conditions are equivalent to the original constraints only if the initial conditions are chosen appropriately, i.e., the initial conditions must fulfill the original constraints (33). As described in [21, p.138], we can use the gap velocity (42) and the gap acceleration to formulate the Signorini condition (40) at velocity and acceleration level, respectively. Specifically, on velocity level Signorini's law reads as whereas on acceleration level we have where we have introduced the set characterizing the contacts that are active on position as well as on velocity level. The complement of B is again denoted asB = {1, . . . , n N } \ B and by definition includesĀ as a subset, i.e.,Ā ⊆B. It is important to point out that if the motion fulfills Signorini's law on one kinematic level, it does so also on all other kinematic levels provided that the initial conditions are compatible with the other kinematic levels.
It is well known that the constraints of a mechanical system can be formulated on acceleration level dt-almost everywhere without changing its motions. Moreover, the acceleration level constraints come with favorable mathematical properties. During impact free time intervals for example, the presence of position level constraints leads to a DAE of index three. A formulation of the system with constraints on acceleration level reduces to a DAE of index one, which are often easier to solve numerically than higher index DAEs. However, the described index reduction by differentiation is prone to numerical drift, meaning that although the constraint is satisfied on acceleration level, the corresponding position and velocity level constraints are violated due to numerical integration errors. As a remedy, the position and velocity level constraints can be stabilized by introducing additional Lagrange multipliers. This stabilization, initially proposed by [18] and hence known by the name Gear-Gupta-Leimkuhler (GGL) method, can analogously be extended to unilateral constraints, see [9].
To stabilize the constraints, we formally extend the kinematics of the system to dq = (u + u S ) dt du = (a + a S ) dt where we have added the velocity u S and the acceleration a S caused by the stabilization to (29) and (31), respectively. The multipliers ν g , ν γ and ν N then take care of the stabilization of (56) and (60) on velocity level by demanding Note, that the stabilization (63) demandsġ = 0 for almost all time instants, which combined to (57) yields the conditionġ = 0 for all time instants. The same reasoning holds for γ = 0. We use the multipliers µ g and µ N to stabilize the constraints (56) and (60) on position level by for all contacts k. It can be shown that, in absence of numerical errors, the solution of the equations of motion formulated with stabilized acceleration level constraints have vanishing Lagrange multipliers ν 2 and µ 2 dt-almost everywhere and therefore a S = u S = 0 for almost all time instants. Moreover, it can be shown that the remaining quantities solve the equations of motion with constraints formulated on position and velocity level described in Section 3. This fact establishes the mechanical equivalence of the original equations of motion of Section 3 and the equations of motion with stabilized acceleration level constraints presented in this section.

Nonsmooth generalized-α discretization
To compute the motion of the mechanical system numerically, in this section we derive a time-stepping scheme from the family of generalized-α methods. The scheme is derived by integrating the equations of motion with stabilized acceleration level constraints over a time interval I = (t i , t i+1 ] and introducing appropriate discrete variables. Considering the velocity u(t) as a right-continuous function, the velocity of the system at a time t can be written as where we have used (62). Similarly, the position of the system at time t i+1 is which with the help of (65) can be reformulated to where we have denoted the last integral in (65) by To derive the position and velocity updates of the scheme, numerical approximations of the integrals in (65) and (67) using quadratures have to be introduced. As approximants for the position, velocity and acceleration at some time instant t i we introduce q i , u i and a i , respectively. Moreover, we define the discrete variables The integrals of the acceleration a in (65) and (67) are discretized in the fashion of a generalized-α method [9] using the quadratures where the auxiliary acceleration variablesā are linked to the approximants of the acceleration by and the time step of the scheme is introduced as ∆t = t i+1 − t i . The coefficients α f , α m , β and γ can be chosen according to Newmark [34], Hilber-Hughes-Taylor [26] or Chung and Hulbert [12]. We choose the last option, which results in a second-order scheme with an adjustable level of numerical dissipation in the high-frequency range. We introduce the spectral radius at infinite frequencies ρ ∞ ∈ [0, 1], which controls the dissipation in the highfrequency range. The coefficients of the scheme are then given by For ρ ∞ = 1 the scheme shows the minimal and for ρ ∞ = 0 the maximal amount of dissipation in the high-frequency regime. Finally, the position and velocity updates are obtained by using (70) and (69) in (65) and (67), respectively. The equations of motion (52) are discretized as where the subscript i + 1 indicates that the quantity is evaluated at t i+1 , q i+1 and when applicable at u i+1 , e.g., . Moreover, we have introduced λ 2,i+1 to approximate the forces λ 2 .
To find the discrete equations for U i+1 , consider the approximation which is exact for a constant mass matrix. Furthermore, we approximate which are exact if W 2 is constant and where we have introduced the discrete variables In view of (53) and (63), the above discretizations (75) and (76) yield the discrete impact equations From the definition of the atomic measure (30), we see that the discrete variables Λ K,i+1 and Λ F,i+1 consist of the sum of impulsive forces k Λ K (t k ), respectively k Λ F (t k ), corresponding to collisions in the time interval I. In addition, Λ K,i+1 contains a contribution due to the stabilization of the constraints.
The discrete equation for Q i+1 is found from the approximation Finally, it is by combining the impact equation (53) with the stabilizing conditions (63) and (64) that we can motivate the discrete equation Hereby, we have introduced the discrete variables κ K,i+1 with K ∈ {g, N } as and have used the approximation similar to (76). Moreover, we have approximated the remaining double integrals by (83) Similar to standard DAE solvers [24], the bilateral constraints on all kinematic levels are discretized by just evaluating them at the end of the time step. For the constraints originating from a position level constraint that is whereas for constraints originating from a velocity level constraint, we have

Discrete normal contact laws
Since we aim at an event-capturing time-stepping scheme, we do not resolve the contact dynamics during a time step I = (t i , t i+1 ] in all detail, but rather derive discrete contact laws capturing the contact dynamics occurring during the time step. More precisely, we derive discrete normal contact laws such that at the end of the time step impenetrability is satisfied on all kinematic levels while capturing the effects of Newton's impact law. We start the discretization of the normal contact laws from the acceleration level Signorini conditions (60). Since the velocity of the system is continuous between velocity jumps, implyingġ k+ N =ġ k− N =ġ k N dt-a.e., we have that where we have used the definition (44) of ξ k N . Consequently, in view of e k N ≥ 0, we can replace the conditionġ k N ≤ 0 with ξ k N ≤ 0 in (61). The discrete Signorini conditions then result by evaluating all quantities in (60) at the end of the time step. Thus, we have where we have used the notation and defined the discrete version of (61) as To formulate the normal contact law on velocity level, we combine Signorini's law on velocity level, the impact law and the stabilization condition. For that, consider the case of active normal contact, i.e., k ∈ A. Then, Signorini's law and the stabilization conditioṅ hold dt-almost everywhere, respectively. Using the cone property of the normal cone, we may write (1 + e k N ) in front ofġ k N in (90), which in view of (86) takes the form The proposed reformulation of the combined Signorini and stabilization condition has the same form as the impact law (43), which reads as With this preparatory work, we can finally proceed toward a discrete law. Since we are interested in the end of the time step, we use the set A i+1 introduced in (89) to determine whether the contact k is active on position level at the end of the time step. If this contact is active, that is k ∈ A i+1 , we assume the contact to be active during the whole time step. Furthermore, we assume ξ k N to be constant within a time step I = (t i , t i+1 ] and to correspond to ξ k N,i+1 , which allows to combine (91) and (92) in integral form as where we have used Proposition 2. If k ∈Ā i+1 , the integral in (93) is zero. At this point it is appropriate to qualitatively discuss the assumptions leading to (93). For that, we look at two cases. Either there is no collision in the time interval I and only non-impulsive motion takes place on I (case i) or there is impulsive motion during the time interval I (case ii).
(i) Purely non-impulsive motion on I: Since the contact is closed at the end of the time step and there is no collision, the contact must either have been closed during the whole time step, or it must have closed at some t c ∈ I = (t i , t i+1 ] withġ k N (t c ) = 0. In the first case, the contact velocityġ k N vanishes on I and by (86) we have ξ k N = 0 on I justifying the approximation. In the second case, for which the contact closes during the time step we haveġ k N (t i ) < 0 andġ k N (t i+1 ) ≥ 0, which is at best approximated by ξ k N,i+1 = 0. This allows the integral over the contact forces in (93) to be non-zero, which can capture the exact dynamics.
(ii) Impulsive motion on I: There are one or more collisions or other events causing impulsive motion (a dynamic catastrophe). Then we have to admit that the impulsive part of the motion is dominating the dynamics on I and we may neglect the non-impulsive motion. The error which we then make in the non-impulsive description by falsely considering ξ k N to be constant is then small and of the order of the time step. If the collision takes place at t c ∈ I = (t i , t i+1 ], then we approximate ). This approximation becomes exact in the limit of ∆t ↓ 0.
Finally, equation (93) suggests a discrete normal contact law on velocity level given by where the discrete percussion is defined by Hereby, we have introduced the auxiliary force variablesλ N linked to the approximants of the contact forces by To motivate (94) as an approximation of (93), we first state that the second part of the discrete percussion (95) is a generalized-α discretization of the integral of λ N over I, see (70) and (71). This, in combination with (77), indeed shows that the discrete percussion approximates the integral appearing in (93), that is, We use the terminology "percussion", even though the discrete variable includes the stabilizing Lagrange multipliers.
In order to state the contact law on position level, we observe that (93) also holds if we only integrate over a time span (t i , t] ⊆ I. Moreover, the inclusion (93) implies that the integral is non-negative if the contact k is closed at the end of the time step. Since, in addition, the integral is zero if the contact k is open at the end of the time step, it follows that for all t ∈ I. Moreover, we approximate the stabilization on position level (64) by for all times t in I. Finally, we use the Proposition 2 to combine (98) and (99) to This inclusion suggests a discrete normal contact law on position level given by together withκ Hereby, the second part in (102) approximates the double integral of λ N using the generalized-α method, see (70) and (71), and it therefore follows from (81) that which confirms that (101) is indeed an approximation of (100).

Discrete friction laws
In (47) we stated the friction law of the k-th contact as for the case where the contact is active, i.e., k ∈ A. Using the cone property of the normal cone, we may write (1 + e k F ) in front of γ k F in (104). Furthermore, for almost all t it holds that u + (t) = u − (t) = u(t). Hence, it follows from (48) that we may equivalently write the friction law as This reformulation of the friction law brings it in a similar form as the frictional impact law (47) In the following discretization process, we exploit that the set of admissible (negative) friction forces has the homogeneity property C F (α) = αC F (1). Moreover, similar as for the laws in normal direction, we consider ξ k F to be constant on the short time-interval I and assume it to take the value Specifically, this translates to the friction law as Because ξ k F,i+1 is constant on I, we can invoke (17), i.e., Proposition 2 in integral form, which directly gives For the impulsive part of the motion on I, we similarly have the approximation and using (17) cast the impact law in integral form Finally, we combine (108) and (110) by using Proposition 2 to which motivates as approximation the discrete friction law where Hereby, we have introduced the auxiliary discrete friction forces linked to the actual forces by Clearly, because of (77) and interpreting the second part of (113) as the generalized-α approximation of the integral of the nonimpulsive friction forces, we have Using the discrete percussion P k N,i+1 in (112) is finally justified by (97) since the stabilizing Lagrange multipliers are zero for the exact solution.
Since the discrete friction law (112) combines the effects of nonimpulsive and impulsive friction forces, an additional friction law is needed to distinguish these effects and compute values for both λ k F,i+1 and Λ k F,i+1 . For that, we first state that the friction law (112) basically consists of three cases. Either the k-th contact is open (k ∈Ā i+1 ) and P k F,i+1 = 0 or the contact is active (k ∈ A i+1 ) and one of the following two cases holds. In the first case, the negative discrete percussion lies in the interior of the set of admissible negative friction forces C F (P k N,i+1 ), which by the normal cone inclusion in (112) implies ξ k F,i+1 = 0. Hence, this can be seen as the case of sticking contact and we introduce the set of sticking contacts The second case is the case of slipping contact, where ξ k F,i+1 is nonzero and the negative discrete percussion must lie on the boundary of the set C F (P k N,i+1 ) by (112). We write k ∈ D sl i+1 in that case, where obviously the set of slipping contacts is just the complement of the set of sticking contacts, that is, We are now ready to state the remaining discrete friction law as where λ k,sl F,i+1 denotes the element on the boundary of andγ k F,i+1 denotes the evaluation at the end of the time step of the accelerationγ k F (t, q, u, a) = (W k F ) T (t, q) a + In our case, the set of admissible friction forces is given by (49) and one can easily verify that In essence, the discrete friction law (117) corresponds to the evaluation of (104) at the end of the time step, with the difference that in the sticking case, the friction law is formulated on acceleration level. We refer to [21,Sect. 10.4] for more details on this matter. We conclude this section with the remark, that the formulation of the sticking case on acceleration level is strictly necessary, because otherwise the two discrete friction laws would lead to an ambiguity. Specifically, in the sticking case, (112) would imply ξ k F,i+1 = 0 and the evaluation of (104) at the end of the time step would imply γ k F,i+1 = 0, which are the same condition in the case of e k F = 0. Hence, in that case the two conditions collapse and create an ambiguity which does not allow to compute λ k F,i+1 and Λ k F,i+1 independently.

Extension to general velocities
We generalize the kinematic equation (29) to dq =q(t, q, u) dt , whereq(t, q, u) = B(t, q)u + β(t, q) with u(t) ∈ R m and B(t, q(t)) ∈ R n×m . Such a generalization is for example needed in rigid body dynamics, when the orientations of the bodies are described by unit quaternions and the components of the angular velocities with respect to a body fixed frame are chosen as velocity parameters of the system. Another prominent example is the use of minimal coordinates and minimal velocities for a nonholonomic system, where typically m < n. It follows immediately from (121) that Gathering the densities with respect to dt after using (31) and (121) allows to rewrite the differential measure ofq as dq =q (t, q, u, a)dt + B(t, q) where we have introduced the function q(t, q, u, a) = B(t, q)a + ∂q(t, q, u) ∂qq (t, q, u) + ∂q(t, q, u) ∂t .
It is straightforward to see, that after introducing the stabilization as in (62), we have dq =q(t, q, u + u S ) dt dq =q (t, q, u, a + a S )dt + B(t, q) By a similar reasoning as in Section 5, we can arrive at the corresponding position update formula given as which generalizes the update formula (73).

Numerical implementation
In this section we apply the ideas of Section 2.3 to the nonsmooth generalizedα method. As initial conditions of the stepping scheme, we assume that we know the values of the kinematic quantities q 0 , u 0 , a 0 , U 0 and Q 0 as well as all discrete forces at the initial time t 0 . The initial conditions must be compatible in the sense that they solve the discrete equations of motion as well as satisfy all constraints at the initial time t 0 . As initial value for the auxiliary variables we choose the first order approximation 1ā 0 = a 0 and λ 2,0 = λ 2,0 .
We assume the quantities at the beginning of the time step to be known and formulate the presented scheme as a system of nonlinear equations R(x) = 0, where The nonlinear equations are then solved by a semi-smooth Newton method.
The computed x can subsequently be used to find the remaining quantities at the end of the time step. Hence, the positions, velocities and percussions at t i+1 are regarded as dependent on x. In fact, we can solve (71) forā i+1 and use it to compute q i+1 and u i+1 from (73) or its generalized counterparts of Section 8. Similarly, the auxiliary contact forces can be computed from (96) and (114), respectively. In turn, these can be directly inserted into (102), (95) and (113) to retrieve the respective values ofκ N,i+1 , P N,i+1 and P F,i+1 .
Having in mind that the just mentioned quantities depend on x, the first part of the residual R is stated as where we assumed a splitting of the residual R = (R T s R T c ) T into a part R s containing all equations except the discrete contact laws, which are contemplated in R c . Hereby, we have chosen the subscript "s" for R s to indicate the smooth part of the residual. In order to state the remaining part of the residual, we have to formulate the normal cone inclusions arising in the discrete contact laws as equations. This is done using the results from Section 2.3 and gives a piecewise smooth residual which we subsequently set up part by part.
We start with the normal contact law on position level (101) and use (22) to restate the law aŝ Equivalently, we can by (26) restate (101) as the residual where we have implicitly defined the k-th component R k κ N of R κ N and where we have introduced the set Since for contacts in A i+1 we are demanding the gap to be closed at the end of the time step, see (129), we have that A i+1 contains the same contacts as A i+1 because due to the stabilization no contact is penetrated at the end of the time step. We can use this fact to state the residual for the contact law on velocity level.
We use again (26) to formulate the normal contact law (94) as Similar to before, the k-th component R k Λ N of R Λ N follows from (26) as where we have introduced the set as well as its complementB i+1 = {1, . . . , n N } \ B i+1 . Since, B i+1 is the subset of closed contacts A i+1 for which by (132) ξ k N,i+1 = 0, the set B i+1 and B i+1 contain the same contacts. The reasoning here is similar as on position level.
Finally, we can proceed in the same manner with the acceleration level constraint (87) and state it as With this intermediate step, where we have introduced the set as well as its complementC i+1 = {1, . . . , n N } \ C i+1 . Again, due to the stabilization and (135), the sets C i+1 and C i+1 contain the same contacts along the discrete motion. For the discrete friction laws the procedure is the same and we start with reformulating (112) as where we used (22). We can now use (27) to implicitly define the k-th contribution of the residual R T where we have introduced the sets and corresponds the set of sticking contacts D st i+1 and consequently D sl i+1 to the set of slipping contacts D sl i+1 . This, together with (22), allows to rewrite (117) as We gather the sticking contacts which are also sticking on acceleration level by and define E sl Consequently, (27) allows to reformulate (140) as where we implicitly defined the k-th contribution of the residual . Now that we have defined the residual R, we can see that it is piecewise smooth, i.e., it is smooth (differentiable) in x if the index sets A i+1 , B i+1 , . . . do not change in the vicinity of x. This enables us to use a semi-smooth Newton method to solve R(x) = 0. Hence, a time step of the scheme can be summarized as follows: Time step with semi-smooth Newton method: 1) As starting value (ν = 0) we use the vector x 0 constructed like (126) but by choosing the known values at the beginning of the time step.
2) While ||R(x ν )|| ∞ ≤ TOL n and ν ≤ MAXITER n do the Newton update and increase ν by one. Hereby, ∇R(x ν ) denotes the Jacobian matrix of R evaluated at x ν , where the index sets A i+1 , B i+1 , . . . are held constant (equal to the sets arising in the computation of R(x ν )) while taking the partial derivatives.
For the sake of completeness, we remark that step 3) is only needed for output purposes and is basically the first part of the computation of R in the subsequent time step. Moreover, it is clear from (143) that the Jacobian matrix ∇R must have full rank.
For mechanical systems without friction, the here presented scheme is very similar to the generalized-α method presented in [9]. The only difference lies in the residuals describing the normal contact law on position and velocity level. Specifically, in [9] the authors use κ k N,i+1 = 0 instead ofκ k N,i+1 = 0 in (129) as well as Λ k N,i+1 = 0 instead of P k N,i+1 = 0 in (132). Even though by our experience, this small difference produces if at all minimal differences in the motion and the forces of the simulated system, it is more than just a subtlety. In fact, replacing the contact force quantities in (129) and (132) destroys the equivalence of the contact law formulated using the sets A i+1 and B i+1 with the corresponding discrete contact laws formulated as normal cone inclusions.
In cases where multiple contacts with linearly dependent generalized force directions are present in the system, the Jacobian matrix ∇R is singular. Hence, for this kind of system the equation R(x) = 0 must be solved by a method that does not invert this Jacobian matrix. A popular choice is to reformulate the system R(x) = 0 such that it can be solved by fixed point iterations, see Section 2.3 and [4,40]. For the aforementioned reformulation, we split the vector of unknowns as x T = (y T z T ), where Since y does not contain any contact forces, ∂Rs ∂y can be regular despite the presence of linearly dependent contact force directions. Assuming that this is the case, the implicit function theorem applied to R s (y, z) = 0 states the existence of a function F such that y = F(z), i.e., we can regard y as depending on the contact forces at the end of the time step. Moreover, it can be seen from the derivation of the residual R c above, that it can be equivalently formulated as z = p(y, z), where p makes use of the proximal point function. Indeed, see for example that (129) is equivalent to (128) or that (138) is equivalent to (137). This leaves us with an equation z = p (F(z), z), which we can solve using fixed point iterations.
Numerically, we can find the value y = F(z) by solving R s (y, z) = 0 for a fixed value of z using Newton's method. The computed value y in turn is then used in the fixed point iteration. Therefore, a time step of the scheme can be implemented as follows: Time step with fixed point iterations: 1) As starting value (ν = µ = 0) we use the vectors y 0 and z 0 constructed like (144) but by choosing the known values at the beginning of the time step.
2) While ||z µ − p(y 0 , z µ )|| ∞ ≤ TOL fp and µ ≤ MAXITER fp do the fixed point update (i) While ||R s (y ν , z µ )|| ∞ ≤ TOL n and ν ≤ MAXITER n , perform a Newton step and increase ν by one. Hereby, ∇ y R s (y ν , z µ ) denotes the Jacobian matrix of R s with respect to its y dependence.
(ii) Use the converged solution y ν of step (i) to perform the fixed point update Subsequently increase µ by one and set ν = 0 as well as y 0 = y ν .
Last but not least, it has to be mentioned that the parameter r used in the prox-equations as well as in the index sets A i+1 , B i+1 , . . . could be chosen differently for every contact and even differently for the normal contact law and the friction law of the same contact.

Examples
In this section we use the presented nonsmooth generalized-α scheme to obtain the time evolution of some benchmark systems, which are all chosen such that particular features of the scheme can be validated separately.

Rotating bouncing ball
Following [20,35], we look at a homogeneous rigid sphere of radius R = 0.1 and mass m = 1 which is constrained to move in the (e I x -e I y )-plane and which under the influence of gravity with gravitational acceleration g = 9.81 falls on a horizontal plane, see with the rotational inertia θ S = 2 5 mR 2 of the sphere. Finally, the contact with the plane is described by the gap function which is the relative horizontal velocity of the contact point with respect to the plane and with (49) describes planar Coulomb friction. We choose the friction coefficient µ = 0.2 and the restitution coefficient e F = 0. To validate the presented scheme, three simulations are instructive. For all of them we choose q(0) = (0 1 0) T as initial configuration and use u(0) = (0 0 ω) T , such that the ball has an initial rotational velocity ω.
The first case starts from rest, that is ω = 0, and we choose e N = 0.5. This results in the typical bouncing motion, which exhibits the Zeno phenomenon.
The simulation result 2 is shown in Figure 1 and asserts that the proposed scheme can overcome accumulation points. The subsequent two cases are used to test the behavior of the scheme with respect to friction forces. For both cases, we set e N = 0 implying that once the contact closes it remains closed, i.e., the post impact velocity is u + y = 0, which allows us to validate friction. At the closing time instant a frictional impact occurs. Moreover, because the ball is constantly accelerated by gravitation, the pre-impact velocity is u − y = − 2g(y(0) − R). It then follows from the impact equations (53), that the impulsive normal contact force is Λ N = m 2g(y(0) − R) ≈ 4.2. Furthermore, after the impact the non-impulsive normal contact force compensates the gravitational force and therefore takes the value λ N = mg = 9.81. Depending on the value of ω, two cases arise.
In the first case, the rotational velocity is high, we choose ω = 50, and the contact slides after the impact, implying that the friction forces attain the maximally allowed values Λ F = µΛ N ≈ 0.84 and λ F = µλ N ≈ 1.96. After a period of sliding contact, finally, the ball has slowed down enough such that a slip-stick transition takes place and the ball begins a pure rolling motion described by the kinematic condition γ F = 0. Since the rolling motion is described by constant velocities, no net forces occur, implying λ F = 0. Hence, at the slip-stick transition the non-impulsive friction force instantly 2 The generalized-α scheme with the following parameters was used: r = 0.3; ∆t = 2 · 10 −3 ; ρ ∞ = 0.5; TOL n = 10 −6 . The semi-smooth Newton method was used for all time steps and the Jacobian matrix ∇R was computed by finite differences with a step size of ε = 10 −6 . The maximal number of Newton steps (143) encountered in a time step was ν max = 1 and the average was ν avg = 0.01. jumps to zero. Figure 2 shows that the described behavior of the friction forces is perfectly reproduced by the presented generalized-α scheme 2 . In the second case, when the rotational velocity is small enough, the contact sticks at impact and the ball exhibits a rolling motion directly after. As no friction force is needed for rolling, we have λ F = 0 for all times. The pre-impact velocities are u − x = 0 as well as u − ϕ = ω. Since the contact sticks directly after the impact, i.e., γ + F = 0, from the impact equations (53) we can deduce that the impulsive friction force takes the value Λ F = − 2 7 mRω ≈ −0.29. Again this behavior is in perfect accordance with the simulation 2 , see Figure 2.

Ball in corner
A ball making frictional contact with a corner is a simple case where the Newton step (143) fails due to the singularity of the Jacobian matrix ∇R. Consider the ball described in Section 10.1 and assume it can get in contact with two inclined planes with inclination angles α = 45 • and β = 45 • , respectively, see Figure 3. The gap functions and the friction velocities are (149) For the simulation 3 shown in Figure 3, we assumed the restitution coefficients to be e 1 N = 0.5, e 2 N = 0 and e 1 F = e 2 F = 0. For the friction coefficients we set µ 1 = µ 2 = 0.3. Starting at rest with q(0) = (−0.5 1 0) T , the ball will eventually come to rest with both contacts closed. It is exactly in that situation that the Jacobian matrix ∇R becomes singular and the simulation can only be continued with fixed point iterations as solution strategy. In Figure 3, this is the case for t ≈ 1.36.

Ball in cylinder
The importance of the stabilization of the unilateral constraint at position level can impressively be shown by simulating the ball of Section 10.1 rolling inside a cylinder of radius R c = 1. We assume that the centerline of the cylinder is orthogonal to the e I x -e I y -plane, such that it can be identified with the point P , see Figure 4. The normal contact with the cylinder is described by the gap function with e N = 0. Using I n = (n x n y 0) T to denote the components of the inward normal n = r SP /||r SP || of the cylinder, we define the tangent unit vector I t = (n y − n x 0) T . With that, we can define the tangent velocity where the velocity of S is I v S = (u x u y 0) T . We choose the friction coefficient µ = 0.1 and the restitution coefficient e F = 0. Looking at the trajectory of  S simulated 4 with the generalized-α scheme shows that the contact does not penetrate. To show that contact penetration is a big issue for this system, we compared the trajectory of S resulting from the generalized-α scheme to the one gotten from a simulation 5 with the widely used time-stepping scheme of Moreau [23,33,40], which does not stabilize the unilateral constraint on position level. The comparison is shown in Figure 4. For the simulation we assume that the ball is initially at rest with q(0) = (−0.9 1 0) T .

Painlevé rod
It is well known that during the sliding motion of rigid bodies over a rough surface, a "frictional dynamic catastrophe" can occur, i.e., impulsive motion which is not the result of a collision. This phenomenon, called Painlevé paradox, is for example studied in [19,20,27], where the following benchmark system is used. Consider a rigid homogeneous slender rod of mass m = 1, length 2l and rotational inertia θ S = 1 3 ml 2 , where l = 1. As shown in Figure 5, the rod moves in the e I x -e I y -plane and is under the influence of gravity with gravitational acceleration g = 10. Describing the orientation of the rod by 4 The generalized-α scheme with the following parameters was used: r = 0.3; ∆t = 10 −2 ; ρ ∞ = 0.5; TOL n = 10 −6 . The semi-smooth Newton method was used for all time steps and the Jacobian matrix ∇R was computed by finite differences with a step size of ε = 10 −6 . The maximal number of Newton steps (143) encountered in a time step was ν max = 74 and the average was ν avg = 0.37.
5 Parameters for Moreau's time-stepping scheme: r = 0.3; ∆t = 10 −2 ; TOL fp = 10 −6 . the angle ϕ, we choose the minimal coordinates q = (x y ϕ) T , where the center of mass S of the rod is addressed by I r OS = (x y 0) T . Using the natural velocity coordinates u corresponding toq almost everywhere, the mass matrix and the force vector of the system have the form (147). The contact of the rod's tip with the ground is described by The friction coefficient is set to µ = 5/3 and e N = e F = 0 is used to model inelastic impact behavior. With the initial conditions q(0) = (0 l sin(ϕ 0 ) ϕ 0 ) T and u(0) = (v 0 0) T , the rod's tip is initially in contact and the rod slides with an initial inclination of ϕ 0 = 31 • and an initial horizontal velocity v = 30. As analyzed in [20], this sliding motion eventually results in detaching of the contact. This detaching comes with a blowup of the accelerations and nonimpulsive contact forces, which is seen as an impact. The blowup of the accelerations leads to vertical asymptotes in the velocities, which constitute problems for integration schemes. However, since the asymptote is just an isolated singularity, event-capturing schemes with constant time step might overcome these blowups, as they do not try to fully resolve it. The simulation 6 result in Figure 5 shows that the presented generalize-α scheme can indeed deal with the Painlevé paradox. 6 The generalized-α scheme with the following parameters was used: r = 0.1; ∆t = 8 · 10 −4 ; ρ ∞ = 0.9; TOL n = 10 −8 . The semi-smooth Newton method was used for all time steps and the Jacobian matrix ∇R was computed by finite differences with a step size of ε = 10 −6 . The maximal number of Newton steps (143) encountered in a time step was ν max = 2 and the average was ν avg = 0.07.

Guided hopper
The suitability of the presented nonsmooth generalized-α scheme for the simulation of flexible multibody systems is demonstrated by simulating a guided hopper. It consists of a vertically guided main body of mass M = 3, which is addressed by the coordinate y. At the hip H with I r OH = (x 0 y 0) T a rigid homogeneous rod of mass m = 1.56, length L = 0.2 and rotational inertia θ S = 8.5 · 10 −3 around the rod's center of mass is attached to the main body. x 0 is the arbitrary horizontal position of the guidance. The orientation of the rod is prescribed by the angle α(t) which is a given function of time. A straight planar Euler-Bernoulli beam [16] with undeformed length L is connected to the knee K of the rod by an actuated rotational joint with prescribed actuation angle β(t) as a given function of time. For the linear elastic beam we choose the axial stiffness EA = 1.89·10 7 , the bending stiffness EI = 14.2 and the mass line density = 0.71. Moreover, we follow [25] and discretize the centerline of the beam with B-Spline shape functions. Using the parameter ξ ∈ [0, 1] to parametrize the beam and denoting the generalized coordinates of the beam by q b , a point C on the centerline of the beam is addressed by r OC (ξ, q b ). We introduce the generalized coordinates q = (y q T b ) T of the hopper as well as the auxiliary quantities d = r OK −r OC (0, q b ) and ϕ = (r KH , r OC (0, q b )), where (·) denotes the derivative with respect to ξ. With these quantities, the bilateral constraints composing the knee joint are Finally, the contact of the endpoint P of the beam with the horizontal plane is described by where v C is the velocity of the centerline and corresponds to the time derivative of r OC almost everywhere. The contact parameters are set to e N = e F = 0 and µ = 0.2. For the simulation, 7 the system is initially assumed to be at rest with y(0) = 0.31. Moreover, the beam is undeformed initially. For the actuation 7 The generalized-α scheme with the following parameters was used: r = 0.15; ∆t = 5 · 10 −5 ; ρ ∞ = 0; TOL n = 10 −6 . The beam was discretized by two elements of polynomial degree 2 and 5 Gauss quadrature points were used. The semi-smooth Newton method was used for all time steps and the Jacobian matrix ∇R was computed analytically. The maximal number of Newton steps (143) encountered in a time step was ν max = 5 and the average was ν avg = 1.23. angles have been chosen and gravity is contemplated by the gravitational acceleration g = 9.81. The percussions, plotted in Figure 6, show that the generalized-α scheme can cope with the complex contact dynamics arising in multibody systems containing flexible parts and time dependent bilateral constraints. This makes the presented scheme well suited for engineering applications.

Tippetop
The tippetop consists of a spherical main body at which a stick is attached. Starting from the standing position with stick pointing up and with high spinning velocity around its symmetry axis, the top inverts and spins on the stick. Since the spinning velocity decreases due to dissipation, the top tumbles back to the standing position after a while. Let R 1 = 1.5 · 10 −2 be the radius of the sphere with center C 1 characterizing the main body of the top. The rounded end of the stick is described by a sphere of radius R 2 = 5 · 10 −3 and midpoint C 2 . The two midpoints as well as the top's center of mass S lie on the axis of symmetry of the top and the distances between S and the points C 1 and C 2 are a 1 = 3 · 10 −3 and a 2 = 1.6 · 10 −2 , respectively, see Figure 7. We describe the position of the top by the components I r OS of the position vector of S with respect to the resting basis I. To characterize the orientation of the top, we introduce the body fixed K-frame such that e K z lies on the symmetry axis of the top and points towards the stick. The transformation matrix A IK = ( I e K x I e K y I e K z ) is then parametrized using a unit quaternion p. Hence, the configuration of the top is described by q = ( I r T OS p T ) T . As generalized velocities we choose u = ( I v T S K Ω T ) T composed by the representations of the velocity v S of S and the angular velocity of the top Ω with respect to the bases I and K, respectively. As mentioned in Section 8, this choice leads to a model with the generalized kinematic equation (121). For the relevant quantities B, β, M and h describing a so parameterized rigid body under the influence of gravity, we refer to [37].
The mass of the top is m = 6 · 10 −3 and the representation of the inertia tensor is the diagonal matrix K Θ S = diag(I 1 , I 1 , I 3 ) with I 1 = 8 · 10 −7 and I 3 = 7 · 10 −7 . To describe the contact between the ground and either the main body (k = 1) or the stick (k = 2), we use the gap The friction between the top and the ground is modeled as Coulomb-Contensou friction [28], which uses the n k F = 3 velocity parameters where v P k = v C k − R k Ω × e I z denotes the velocity of the respective contact point and R = 5 · 10 −4 denotes the assumed contact radius. This friction model assumes that the normal contact force λ k N stems from a parabolic force distribution over a circular contact area of radius R. A closer look at the velocity parameters (157) reveals that the first two velocities correspond to tangential contact velocities and capture translational Coulomb friction. Moreover, the third component is a representative radial contact velocity, which accounts for drilling friction. It is shown in [32] that using the friction velocity (157) together with the set of admissible friction forces (49) results in a good approximation for the Coulomb-Contensou friction law used for the tippetop simulation in [28].
For the simulation, 8 the top is initially at the position I r OS (0) = (0 0 z 0 ) T with z 0 = 1.2015 · 10 −2 and is inclined by the angle θ(0) = 0.1, where θ denotes the angle between e I z and e K z , i.e., e I z · e K z = cos θ. The center of mass S is assumed at rest and the top spins with initial angular velocity K Ω = (0 0 180) T . To assure that the quaternion p remains a unit-quaternion, p is normalized after every step.
It is apparent from Figure 7 that the simulation results using the presented scheme are in line with the results of [32] 9 . This shows that the presented scheme is well suited for mechanical systems with spatial friction as well as models with a general kinematic equation (121).

Conclusion
We presented a nonsmooth generalized-α method for the simulation of mechanical systems with frictional contact. The dynamics of such systems, which additionally to the contacts can be subjected to ideal bilateral constraints, is described within the theory of nonsmooth mechanics. We have modeled the frictional contact as unilateral constraints described by Signorini's law, the generalized Newton's impact law and a Coulomb-type friction law. All these laws are formulated as normal cone inclusions, allowing for a coherent discretization of the contact laws. All constraints are formulated on acceleration level and a numerical constraint drift, and hence contact penetration, is avoided by stabilizing the constraints on velocity and position level using the Gear-Gupta-Leimkuhler approach. The stabilized dynamics is discretized using generalized-α quadratures, which leads to a second-order accurate scheme during impact-free time intervals. Exploit- 8 The generalized-α scheme with the following parameters was used: r = 0.001; ∆t = 1 · 10 −3 ; ρ ∞ = 0.5; TOL n = 10 −6 . The semi-smooth Newton method was used for all time steps and the Jacobian matrix ∇R was computed by finite differences with a step size of ε = 10 −6 . The maximal number of Newton steps (143) encountered in a time step was ν max = 3 and the average was ν avg = 0.36. 9 To reproduce the result of [32] with Moreau's time-stepping scheme, we used: r = 0.001; ∆t = 1 · 10 −4 ; TOL fp = 10 −6 .
ing particular properties of the normal cone enabled to directly discretize the contact laws leading to their discrete counterparts, which depend on the discrete kinematic and kinetic quantities arising from the generalized-α discretization of the dynamics. The resulting time-stepping scheme still contains normal cone inclusions, which we have proposed to numerically treat by either a semi-smooth Newton method or by fixed point iterations. Finally, the derived nonsmooth generalized-α method has been validated using a set of mechanical benchmark systems.
The main contribution of this paper is the extension of the nonsmooth generalized-α schemes [8,9,11,14] to systems with frictional contact and general kinematic equations. This has been achieved by formulating the setvalued Coulomb-type friction on velocity and acceleration level. Moreover, we have introduced the discrete contact laws as normal cone inclusions, which are shown to be a convenient interface to the numerical implementation of the scheme. Last but not least, we have devised a collection of benchmark systems suitable for validating the performance of any numerical scheme for mechanical systems with frictional contact, where every benchmark lies the focus on a different characteristics of such systems or of the numerical scheme. These characteristics are: accumulation points, slip-stick transitions, the presence of linearly dependent force directions, contact penetration, the Painlevé paradox, the suitability of the scheme for the simulation of flexible multibody systems, combined spatial friction laws and general kinematic equations. We showed that the presented nonsmooth generalized-α scheme performs very well for all benchmark systems and hence qualifies for a broad range of engineering applications.
In fact, the presented scheme can reproduce the dynamics of a frictional contact adequately. Moreover, it performs well for multibody systems containing flexible parts and allows general parametrizations such as the use of unit quaternions for the rotation of rigid bodies.