Benchmark cases for robust explicit time integrators in non-smooth transient dynamics

This article introduces benchmark cases for time integrators devoted to non-smooth impact dynamics. It focuses on numerical properties of explicit integrators. Each case tests one necessary numerical property in computational impact dynamics: energy behaviour at impact, angular momentum conservation, non-linear behaviour. The cases are easy to implement and analyse, providing a benchmark well-suited to first numerical studies. We rewrite explicit schemes for non-smooth impact dynamics with unified notations, and analyse them with the benchmark cases.


Introduction
Direct time-integration schemes are a main issue in structural transient dynamics. They have been developed and improved to meet some crucial numerical properties. The first ones are stability and accuracy in linear regime, where Newmark's schemes are the most used [1]. They were improved by adding a controlled high frequency numerical dissipation to get α-generalized schemes [2,3]. But these schemes do not keep their properties for non-linear problems. Symplectic (energy-momentum conserving) or variational time integrators have been then developed to meet non-linear issues: energy-decaying property, overshoot or numerical integration of internal forces. For variational integrators, a review can be found in [4]. For symplectic schemes, a pioneer work is provided by Simo and Tarnow in [5]. They present a symplectic scheme for non-linear dynamic, and show the Central Difference Method (CD method) is the only symplectic scheme of Newmark's family.
An other crucial issue in computational mechanic is to enforce contact constraints. The most common ways are penalty methods, augmented Lagrangian, Lagrange multipliers and Nitsche methods [6,7]. Penalty methods are widely used in finite element simulations. Relating the contact reaction to penetration, they brings stability problems especially for explicit schemes. On the contrary, Lagrange multiplier methods enforce exactly the con-tact constraints but induce an implicit contact problem. Contact problems are divided into smooth and non-smooth ones. A smooth contact problem, like persistent contact, keeps the regularity of velocities, forces and accelerations. For non-smooth contact problems, impacts imply velocity jumps which make forces and accelerations locally non-defined. In this paper, we focus only on impact dynamics. According to [8,9], two approaches coexist to deal with impact in a discrete time: 'event-driven' and 'time-stepping' schemes. Eventdriven schemes precisely locate impact times. They are suited for small system with few impact times, but have no convergence proof in case of infinite number of impacts in a finite time. Time-stepping schemes do not require the location of non-smooth events, and have convergence proof even for infinite number of impacts. The reader will find more information about general computational contact dynamics in [6,7], and scheme reviews in [10,11]. We just evoke here some standard schemes in impact dynamics.
The first one is Laursen, Chawla and Love's algorithm [12,13]. It follows the pioneer work of Simo, Tarnow and Wong on symplectic schemes [5,14] by extending its symplectic properties to contact part. More recently, Moreau-Jean's scheme [15,16] is also energy conservative for contact [17]. It moves the imposition of contact constraints from position to velocity level. This brings new numerical properties well-suited to non smooth impact dynamics. Based on a θ -method [16], Moreau-Jean's scheme is only first order. Acary presents an higher order Moreau-Jean scheme in [8] by adjusting the time-step with impact time. And finally Chen et al. manage to adapt extended Newmark's schemes (αgeneralized) in Moreau's framework in [18,19]. In [20,21] an other energy conservative approach is described. The authors redistribute the mass on the contact boundary in order to get a conservative contact. This approach stays valid whatever the scheme, provided it is implicit.
In an explicit framework the available schemes are less numerous. Carpenter et al. in [22] adapt the Central Difference method for contact with Lagrange multipliers. The contact constraints are imposed at position level. Paoli and Schatzman present an other adaptation of CD method in [23,24] close to Carpenter's one. They modify the writing of contact conditions in order to integrate an impact law. Then Cirak and West [9] develop an explicit integrator in a variational framework with Lagrange multiplier. And finally Fekak et al. [25] gather Moreau's contact formulation and Cirak and West approach in the CD-Lagrange scheme, in order to make a bridge between variational integrators and Moreau's formalism. Based on an asynchronous version of the central difference method, the CD-Lagrange scheme keeps the symplectic properties suitable to non-smooth impact dynamics. Recently, in [26], an other interesting explicit approach is provided using Nitsche's forms for contact. We have to mention also the work of Casadei on implementations of explicit time integrators in finite elements problems, see e.g. [27].
The goal of this article is to present benchmark cases in impact dynamics mainly designed for explicit schemes. It is divided into two parts: the first one presents the theoretical framework and some schemes in unified notations, the second one presents the test cases intended to provide a benchmark. We try to focus on one main property in each test case: impact behavior, angular momentum conservation, spurious high frequencies in finite elements problems, treatment of a non-linear damping term, or mass-scaling ability. The presented test cases are suitable for implicit schemes, even if the benchmark is not designed for.

Discrete model and notations
We present here a frictionless contact-impact formulation between a deformable body and a rigid obstacle. This framework is adequate to spotlight the main ideas of impact dynamics. Switching to frictional contact or deformable-deformable impact does not bring new algorithmic difficulties (related to the impact law), and will be detailed in specific remarks.
The problem is depicted in Fig. 1. We consider a deformable body impacting a rigid boundary R which represents the obstacle. The boundary of is separated into three distinct regions: u for Dirichlet's boundary conditions, F for Neumann's, and c which gathers all potential contact points.
x(t), u(t) andu(t) correspond respectively to position, displacement and velocity. is made with an elastic (potentially non-linear) material obeying a stress-strain law σ = f (u), with σ the Cauchy stress field. u d is a prescribed displacement on u , and F d a prescribed load on F . The external outer normal to is n .
As a starting point, the contact constraints are those of a unilateral contact law [28,29]: x M is the closest point of x on R , and n is the outer-pointing normal of R at x M . When contact occurs, n = −n . We separate n and n in order to define g N and λ N even for free-of-contact case, and λ N is defined along −n to use a complementary form (see below).
By impenetrability condition (1), solids can not penetrate each others. The intensility condition (2) implies that the contact stress is only compressive (i.e. no adhesion occurs). The complementary condition (3) reduces contact states to only two cases: active contact (g N = 0 and λ N > 0), or non-contact (g N > 0 and λ N = 0). (3) implies also that λ N does not produce any mechanical work. Although usually referred as Signorini's conditions, they were stated by Hertz, Signorini and Moreau as reported in [28]. We thus call them HSM conditions.  (1), (2) and (3) can be rewritten in a single equation by a complementary formulation:

Remark 1
The HSM conditions rule only the normal components on c . In case of frictional contact (4) stays valid, but a friction law is added for tangential components [6,7]. We consider e.g. a Coulomb law with a friction coefficient μ: The tangential velocity is defined by v T =u −u · n, in the orthogonal plane to n.
As shown by Moreau in [15] the stress-displacement formulation (4) is equivalent to an impulse-velocity one: With v N =u · n the normal part of velocity, and r N the impulse associated to λ N ( r N = dλ n ). The impulse-velocity formulation is a conditional complementary formulation. Such a formulation is well suited to an algorithmic framework.
The previous formulation (6) corresponds to a 'plastic' impact. For other contact behaviours, we introduce Newton's impact law: t i is an impact time, e c ∈ [0, 1] the restitution coefficient, andṽ N the formal velocity. This law (7) controls the given back energy at impact. With e c = 1 no energy is absorbed ('elastic' impact), with e c = 0 the impact is dissipative (and corresponds to 'plastic' case). The impulse-velocity formulation (6) is naturally gathered with Newton's impact law (7). And the final formulation for HSM conditions and Newton's impact law is, at each point of C : Finally one summarizes the strong form as follows: HSM conditions + Impact law on c (Normal contact law) This problem is discretized in space by the finite elements method based on a weak form of (9). Several approaches exist for formulating the HSM conditions under a weak form. The common methods are Lagrange multipliers, penalty method, augmented Lagrangian, or recently Nitsche methods. The reader will find more details in [6,7]. We focus only on Lagrange multipliers methods, because they enforce exactly the HSM conditions. After a space discretization, nodes are numbered by an index k ∈ {1, . . . , N nodes } with N nodes the number of nodes. We note the vector of nodal displacements U = u k , k ∈ {1, . . . , N nodes } , with u k the displacement of node k. More generally · k denotes the considered field at node k.U andÜ are the time derivatives of U. X is the vector of nodal positions. It depends on displacement: X = X 0 + U with X 0 the vector of initial nodal positions. M is the mass matrix, F int the vector corresponding to material response, and F ext the vector of external loads.
The nodes of c are specified by indexes {1, . . . , p}. We note respectively n k and g k N , the outer-pointing normal of R and normal gap associated to node k. We build a projection operator L on C : L = (n k ) t , ∀k ∈ {1, . . . , p} , and we note g N = g k N , ∀k ∈ {1, . . . , p} . As L and g N depending on the displacement, we note L(U) and g N (U). We note respectively λ and r, the vectors of normal nodal contact forces and impulses. In the following in order to simplify notations, we omit the subscript N being the contact frictionless.
In [9] the authors start from the discrete Lagrangian of an impact dynamics problem. They obtain the following semi-discrete equations by a variationally consistent approach. Over the time interval [t 0 , t f ] only one impact time t i is considered for simplicity, but several nodes can impact C at t i . And the generalisation at several impact times over [t 0 , t f ] is easy (see [8,9]).
Equation (10) corresponds to the smooth dynamics equilibrium, completed by nonsmooth dynamics Eq. (11). Equation (12) is linked to the kinetic energy balance during impact. Following Moreau, Eqs. (10) and (11) can be gathered in one Eq. [15,16,30]. This is developed and used in [8,18,19,25,31]. The following measures on velocity are defined: dU s deals with the smooth part of velocity, dU I with velocity jumps. From (10), and (11): dt is the standard Lebesgue measure for time, and: The smooth equilibrium (10) and impact part (11) can be joined, and augmented with HSM conditions (8) to get the non-smooth contact dynamics equation (NSCD equation): The NSCD Eq. (13) does not contain the Eq. (12) on kinetic energy conservation during impact. But this energy balance is satisfied thanks to HSM conditions with Newton's impact law (8) for e c = 1.

Time integration schemes for impact
In the following, we rewrite schemes used in the next benchmarks with unified notations. We focus on Carpenter's, Paoli-Schatzman's and CD-Lagrange for explicit schemes with Lagrange multipliers. And we choose Moreau-Jean's scheme, although implicit, as a reference to non-smooth dynamics. All are time-stepping schemes. As explained in Remark 2, only time-stepping schemes are robust in a finite elements framework with multiple impacts.

Carpenter's explicit scheme
In [22] Carpenter presents an explicit scheme with contact treated by Lagrange multiplier: the forward increment Lagrange multiplier method. The time integration scheme is a multi-step central difference method, the HSM conditions are considered in a displacement formulation, and contact loads added directly to equilibrium equation. Contact forces at t n are computed in order to assure HSM conditions at t n+1 . This offset is necessary to get an explicit scheme compatible with Lagrange multipliers. The discrete equations are: We note h the time-step. The mass matrix is lumped [29,32] to get an explicit scheme. L n+1 is defined below. The final scheme can be summarized as follows: (18) The displacement at t n+1 is computed by an predictor-corrector approach. U * n+1 is firstly updated without contact contribution. It is used to predict the gap at t n+1 and construct L n+1 . Then U c n+1 is computed to correct negative gap and added to U * n+1 to get total displacement U n+1 . This one satisfies exactly the HSM conditions: ∀k ∈ {1, . . . , p}, g k N (U n+1 ) = 0. Here it can be noticed that Eq. (16) stay implicit, and form a Linear Complementary Problem (LCP). The Delassus or Steklov-Poincaré operator defined by H = L n+1 M −1 lump L t n+1 is generally non-diagonal. In [22] an effective way is proposed to solve this LCP. If H is diagonal, the problem is then fully explicit.

Paoli-Schatzman's scheme
Paoli and Schatzman in [23,24] present a scheme close from Carpenter's but with modified HSM conditions in order to integrate Newton's impact law. The following version of Paoli-Schatzman's scheme is slightly modified for velocity calculation and proposed by Acary in [31]. We rewrite it in a form close from Carpenter's scheme.
The discrete HSM conditions is modified into: and finally we get the following scheme: Carpenter's and Paoli-Schatzman's schemes differ only in the gap calculation. And Carpenter's scheme is a particular case of Paoli-Schatzman's with e c = 0. The modified HSM conditions (19) impose Newton's impact law over three time-steps. If the node k stays in contact during three time-steps meaning that g k With n k the outer-pointing normal of R associated to node k over the three time-steps. A persistent contact over three time-steps is a strong assumption not always met in practice.

CD-Lagrange scheme
CD-Lagrange scheme [25] is the only explicit one built on the NSCD Eq. (13). This equation is integrated between t n+ 1 2 and t n+ 3 2 to obtain the discrete equilibrium equation on velocity. The estimation of time integrals is done by a midpoint rule, and treatment of impact done at the end of time-step. A CD method links velocities and displacements. The discrete equations are: With: By multiplying the discrete equilibrium (26) by L n+1 M −1 lump , we obtain the following equation for contact: , the Delassus operator. The final algorithm can be summarized as follows: A LCP solver is necessary to solve the contact problem (31) unless H is diagonal. This problem is similar to that of Carpenter's scheme, and thus the solver of [22] can be used.
The LCP theory and other resolution algorithms are presented in [33]. If H is diagonal, according to [25], the previous algorithm can be rewrite in a fully explicit form: see algorithm (1).
In [25], Fekak et al. write this scheme differently in order to link it with Newmark's CD method and the approach of Chen [18,19]. This other writing is also more suitable to discrete energy balance. For the smooth part, it is equivalent to a CD method keeping its properties. The scheme is therefore conservative for a discrete modified energy [25], and symplectic. CD-Lagrange is a fully explicit scheme for the smooth part if the internal stresses do not depend on velocity. In this case, Belytschko [29] proposes velocities at t n+ 1 2 to compute explicitly the internal stresses. However the contact problem stays generally implicit like in Carpenter's and Paoli-Schatzman's because of non-diagonal Delassus operator H.

Moreau-Jean scheme
From the NSCD framework we compare the explicit scheme CD-Lagrange to an implicit one, Moreau-Jean's scheme. It is introduced in [15,16].
This scheme is obtained by a time integration of NSCD Eq. (13) between t n and t n+1 and the θ-method [16]. The contact problem (34) is build with a explicit predicted position U n+ 1 2 . It ensures discrete energy balance for θ = 1/2 [17]. And in this case, it is equivalent to Newmark's or Crank-Nicholson constant average acceleration method (CAA method) for smooth part.
Remark 2 Time integration schemes for impact split up between event-driven and timestepping schemes. Event driven schemes detect impact times and compute the solution between them. They keep the computed solution smooth avoiding impact velocity jumps. Impact conditions become initial conditions when the computation restarts after an impact. Time-stepping schemes do not adapt their time-step dealing directly with impact velocity jumps. They consider impact contribution at the end of time-step allowing multiple impacts in one. According to [8], time-stepping schemes are more efficient than event-driven ones in non-smooth dynamics. In case of infinite number of impacts in a finite time, they have a convergence proof unlike event-driven ones. They do not require a detection of impact events, and stay efficient in case of multiple impacts in one time-step (no refinement of time-step). Moreover, in time-stepping schemes, a derivation of contact constraints is not required like in event-driven ones for initial conditions after impact. The discrete HSM conditions are written in term of the main unknowns: displacement for Carpenter's and Paoli-Schatzman's, velocity for Moreau-Jean's and CD-Lagrange. This brings stability to the system, and allows to prove convergence.

Remark 3
Only the velocity formulation of HSM conditions integrates a discrete Newton impact law. But, as highlighted above, this law is necessary for a energy consistent impact formulation. Therefore in presented schemes, only these using the velocity formulation of HSM conditions can ensure the discrete energy balance at impact. This form combined with a time-stepping approach requires to authorize a residual penetration. The displacement constraint has to be relaxed in order to impose the velocity form of HSM conditions only at the end of time-step. The residual penetration decreases with time-step and stays small for usual ones.

Remark 4
The above-mentioned explicit schemes are fully explicit only if H = LM −1 lump L t is diagonal. As the mass matrix is lumped, the diagonal feature depends only on L. This operator links global quantities to local contact ones. If the nodal contact condition depends only on its own quantities, H is diagonal. But if the nodal contact condition links several nodes, it is not. For example, in case of rigid-deformable contact, H is diagonal when the contact boundary is analytically described. The contact condition is easily expressed for each node with its displacement. In case of deformable-deformable contact, both sides of contact have unknown displacement. H is then diagonal only for node-tonode contact, i.e. for conforming meshes. Non-conforming meshes lead to a non-diagonal Delassus operator (see [22]).

Bouncing ball: one DOF-linear system with rigid impact
This test case focuses on impact. A rigid ball falls and impacts the ground, an analytic boundary. This simple problem is quite common in literature about non-smooth impact dynamics [18,31]. Impact times are analytically known, which facilitates result analysis.
When contact occurs, the energy conservation depends on restitution coefficient. If e c = 1 the ball bounces without any dissipation to its initial height. If e c < 1 the contact absorbs energy at each impact, and the time between two impacts decreases. The end of movement is thus an infinite accumulation of bounces in a finite time, sometimes called 'Zeno's paradox'.
If event-driven schemes seem compelling by keeping a high regularity, the simple case of bouncing ball eliminates them. They are obviously not able to represent Zeno's paradox: an infinite accumulation of impacts in a finite time. On the contrary time-stepping schemes deal efficiently with this paradox. The same problem appears in a finite elements framework, when multiple nodes impact nearly at same time. An adaptive time-step leads to time-consuming simulations.
The problem is depicted in Fig. 2, and corresponds to equation: We set problem parameters and initial values at: n indicates the index of impact times, t i 1 is the first impact time. By tending n towards infinity for e c < 1, the stopping time is: We test firstly the schemes based on HSM conditions in displacement formulation. Figure 3 shows the results obtained with the explicit and implicit Carpenter's schemes described in [22]. Being designed with HSM condition on position, the shock law is restricted to e c = 0 for these schemes. The explicit scheme follows the analytic solution, but the implicit one is unstable at impact. Chen in [18] finds similar results for others implicit Newmark's schemes with a displacement formulation for HSM conditions. Figure 4 focus on Paoli-Schatzman scheme results for e c = 1. Sometimes the ball reaches a height higher than initial one, meaning that impact injects energy. Even if this behaviour decreases with time-step, Paoli-Schatzman is therefore not satisfying from an energy point of view like Carpenter's scheme.
We focus now on schemes based on HSM conditions in velocity formulation: Moreau-Jean's and CD-Lagrange scheme. On Fig. 5 we plot position got for e c = 1 with these two schemes. Despite penetration at impact the ball always reaches the same height: the energy at impact is conserved. Acary shows in [17] the energy conservation at impact for Moreau-Jean's scheme with θ = 1/2. For CD-Lagrange, we demonstrate it by considering a node k which impacts in [t n+ 1 2 , t n+ 3 2 ]. The contact is active only in this time interval.
We have r k n+ 3 2 = 0 and r k n+ 1 2 = r k n+ 5 2 = 0. According to [25], the work of contact force at node k is:    The analytic solution eases convergence studies on this test case. We use the norm described in [8,25]. In [8] the author observed the same convergence order with this norm or Hausdorff norm. The norm is: We observe on Fig. 8 a numerical error of computer precision order for the phase before the first impact. For smooth part, both schemes are therefore second order. But they are only first order when an impact occurs, as visible on Fig. 9. Indeed these schemes keep the properties of the original one for smooth part (CD method for CD-Lagrange or CAA method for Moreau-Jean), but the contact part is only first order. In [8] these properties were demonstrated for Moreau-Jean's scheme.
This first problem evaluates the contact behaviour. Only schemes using a velocity formulation of HSM conditions have a energy consistent impact behaviour. This test case eliminates also directly event-driven schemes by its practical application of Zeno's paradox. In the following we only focus on CD-Lagrange and no more on Paoli-Schatzman's scheme. Indeed it does not conserve exactly the discrete energy at impact. The CD-Lagrange scheme will be compared to Moreau-Jean's scheme as a standard in non-smooth impact dynamics with HSM conditions in velocity formulation.

Van der Pol oscillator: one DOF-Non-linear system with rigid impact
The Van der Pol oscillator does not represent a mechanical situation. It is an zero degree of freedom oscillator with a non-linear damping term. This test case evaluates the nonlinear behaviour of schemes, especially for long time simulations and explicit treatment of non-linear damping term.
The Van der Pol oscillator tends to a limit cycle in phase space. Only energy conservative schemes can represent it especially for long time simulations. This problem requires also to treat explicitly a term which depends on velocity in internal stresses. A solution is to use the last known velocity, this modifies slightly the dynamic equilibrium but converges to the same solution [29].
The equation of Van der Pol oscillator is: The two parameters ξ and ω 0 set the behaviour of the system. ξ rules the non-linear term preponderance: for ξ ω 0 the system tends to the linear oscillator. For following simulations we use as parameters: The discretizations of non-smooth dynamic equilibrium (36) are respectively for CD-Lagrange and Moreau-Jean:   Figure 11 shows the computed solutions in phase space at times t n over 300 s (only a few points are drawn). The arrow represents the velocity jump at impact. The limit cycle is visible, and both schemes stay on it for more than forty cycles. For the Van der Pol oscillator, the CD-Lagrange scheme has results as good as Moreau-Jean with the same time-step with a smaller number of iterations (no Newton algorithm).
We compute a refined solution with a small time step h = 10 −7 s for a convergence study. We use the norm (35) described previously. And we add a modified CD-Lagrange scheme for comparison purpose, with a damping term evaluate withẋ n+1 : This scheme is implicit, we consider it just for comparison. In the following we call it 'CD-Lagrange Implicit'. On Fig. 12 the convergence order is studied without impact. Moreau-Jean scheme (θ = 1/2) keeps the convergence order of Newmark's CAA method. But the CD-Lagrange is only first order, whereas its implicit version finds back the second order. Taking velocity This second test case evaluates the scheme in a non-linear framework with impact. It verifies the energy conservation for long time simulation by means of the limit cycle in phase space. And it investigates also the error induced by explicit treatment of damping term.

Rotating spring: two DOF-angular momentum conservation with impact
The rotating spring is a mass-spring system in rotation only subjected to the internal spring force. This system evaluates the scheme behaviour with large rotations. It tests particularly the conservation of angular momentum. This property is ensured by symplectic algorithms [5,14] as CD method.
In order to add a frictional contact, we consider now the following discrete contact problem (written for CD-Lagrange scheme): With L T a projection operator on the orthogonal plane to n, r T and v T the tangential components for impulse and velocity, and H T = L t T M −1 lump L T the Delassus operator for tangential contact problem.
Equations (39) and (41) are the preceding ones for normal part, and (40) and (42) are those for tangential contact problem ruled by the Coulomb law. The algorithm below (2) describes the CD-Lagrange algorithm for a frictional contact, with a diagonal H. This algorithm is very close from the preceding one (1) with frictionless contact, which corresponds to the case μ = 0. It is fully explicit.

Algorithm 2 CD-Lagrange-frictional contact
The system and notations are depicted in Fig. 14, and the equation is: We use the following parameters: We add a concentric contact boundary. Its radius of 1.4 m is larger than relaxed spring length. The time-step is taken to h = 10 −1 s to be in order of the critical one. We present two cases: frictionless and then frictional contact. We consider a Newton impact law for normal part, and Coulomb friction for tangential part. For the frictionless case, we set e c = 1 and μ = 0 for an energy conservative contact. In the second one, we set e c = 0 and μ = 0.2 for a dissipative contact. As no analytic solution exists, a reference solution is computed by CD-Lagrange scheme with refined time step at h = 10 −5 s. This test case is geometrically non-linear needing a Newton algorithm for implicit schemes. The case of frictionless contact is depicted in Figs. 15 and 16. On Fig. 15 we plot the mass position over the first 5 s. We denote a large penetration at each impact for both schemes due to the large time-step. In this frictionless case, the impact impulse should not change  the angular momentum being parallel to position vector. But for Moreau-Jean's scheme the angular momentum oscillates and decreases (more than 10% of its initial value over 100 s) as shown by Fig. 16. The conservation of angular momentum is achieved only at convergence. For CD-Lagrange the angular momentum is exactly conserved to its initial value illustrating the symplectic aspect.  Fig. 18 the angular momentum over 100 s. At each impact the angular momentum is decreased by the tangential contact force. The system reaches a state where it tangents the contact boundary but without contact. The final angular momentums for both Moreau-Jean and CD-lagrange are smaller than the refined one due to numerical errors especially at contact. These discrete solutions include more contact phase because of numerical errors at contact detection. It can be noticed when Moreau-Jean scheme reaches its tangential state, the angular momentum is still oscillating.
This test case with frictionless contact focuses on angular momentum conservation, a crucial property of symplectic schemes. Moreau-Jean's scheme e.g. does not conserve discrete angular momentum. The induced error can become large for long time simulations. On the contrary the discrete angular momentum is exactly conserved at its analytical value with symplectic CD-Lagrange scheme. For the frictional case, this test case shows also the ability of schemes to deal with frictional contact with no additional algorithmic complexity.

Impacting bar: multiple DOF-FE discretization and impact
The impacting bar problem is a frequent test case in literature [9,19,22,25,31] (see [10] for a review of contact schemes on it). Indeed it is interesting from mathematical and numerical point of view. The analytical solution is known and unique, which is rare for elastodynamics with unilateral constraint. And from numerical point of view, the discretization introduces spurious high frequency oscillations at contact node when impact occurs.

Fig. 19 Impacting bar
We consider a multiple degrees-of-freedom, linear and elastic bar, discretized in space by linear finite elements. One of its ends impacts a wall with a contact area reduced to one node. Figure 19 depicts this problem.
For the bar, we consider a density ρ, cross section S, Young Modulus E, length L and initial velocity v 0 . It is discretized by N linear finite elements of equal length l = L/N . The mass matrix is lumped by summing the terms of rows into the diagonal. We use the lumped mass matrix only for explicit scheme (CD-Lagrange). The elementary mass and stiffness matrices are: The material parameters are those of [22,25]: The equation of discrete problem is then: with M the global mass matrix (lumped or not), and K the global stiffness matrix. For discretization: The time step is equal at 0.7 × h crit with h crit = 9.82 × 10 −7 s the larger stable time-step for explicit schemes (based on CFL condition). We set e c = 0 to be closer to analytical solution: for e c = 1, oscillations on displacement appear at contact node after impact. The contact is therefore dissipative, and energy conservation is only achieved at space and time convergence. The lost energy corresponds to the kinetic energy of contact node, which is tiny with an accurate discretization. On Fig. 20 the position of contact node is represented for both schemes. The discrete solution is close from the analytical one with a persistent contact. But velocity oscillations appear after the contact release at contact node. They are visible in Fig. 21 for velocity,   and also for impulsion in Fig. 22. Both schemes are quite equivalent with respect to these oscillations. As mentioned in [25] the amplitude of oscillations for CD-Lagrange decreases with a refinement in space and time, but their frequency increases. Dissipative schemes could damp these oscillations like EG-α scheme [34] or Noh and Bathe scheme [35]. But in non-smooth dynamic framework, no explicit ones are reported in literature. Chen and Acary present implicit dissipative schemes for non-smooth impact dynamic in [18,19] based on implicit α-generalized schemes. The Time Discontinuous Galerkin schemes are also particularly effective for damping spurious mode. In [36], a test case close from impacting bar is presented without spurious oscillations. These schemes are implicit too.
The impacting bar tests the behaviour of schemes in a finite element framework with impact. It brings the difficulties linked to discretization, in particular spurious high frequencies.

Discrete arch: multiple DOF-mass scaling feature with impact
This test case shows the ability of an explicit scheme to be used in a 'mass-scaling' solving strategy. We are not interested in the dynamical behaviour of the structure, but only with the final equilibrium configuration. The dynamic simulation is therefore used to allow a cheap determination of the static configuration. This test corresponds to a limit point instability of a discrete shallow (circular) arch which 'snaps' into its buckled configuration. It is quite similar e.g. to those of [37,38], but with an impact. In this problem neither the load, the mass, nor the viscous dissipation have an influence on the final solution They are therefore only numerical artefacts allowing an explicit code to provide the static solution. The physical model is restricted to the computation of the internal nodal forces, due to large displacement elasticity problem in a given configuration (corresponding to the Cauchy internal force description). The test configuration is depicted in Fig. 23. To avoid using Lagrangian beam theory, the test is built on a discrete chain of elements that encounter traction-compression solicitations similar as springs, with concentrated rotational springs at internal nodes. The two end nodes are simply supported, without rotational springs. A temporary vertical nodal load F d (t) is applied on this structure to allow it to pass the buckling limit point. After that instability point it snaps to a new equilibrium configuration with contact at distance d c below horizontal line. The internal stresses are composed of spring-like forces − k(l − l 0 ) and torsion springs torques − κθ. k is an equivalent stiffness, l the current length of element, and l 0 the initial one. κ is an equivalent rotational stiffness, and θ is the difference between the initial angle between two contiguous bars and the current angle. Each node is submitted to actions of closest neighbours. We note the vector of internal stresses F int which only depends on current nodal displacements. The arch has a We consider a lumped mass matrix. We set the restitution coefficient to e c = 0 for a dissipative contact. We add a viscous term to dissipate kinetic energy and reach the static final state. It is a low frequency dissipation term, proportional to lumped mass. For CD-Lagrange scheme its expression is −μM lumpU n+ 1 2 (the velocity is taken at t n+ 1 2 to make this contribution explicit).
The discrete problem (with CD-Lagrange scheme) is summarized as follows: On Fig. 24, we show the initial, final and an intermediate configurations.

Remark 5
Multiple impacts can occur in the same time-step, and for several nodes. But the Delassus operator H is here diagonal as explained in Remark 4. In fact the nodal contact condition is simply expressed with the sign of nodal position.

Termination criteria
To end the simulation, i.e. when a static solution is obtained, a termination criteria is required. We do not rely on a sufficiently small nodal velocity since it could depend on the level of viscous damping. We therefore prefer to compare the vertical inertial impulse k −μhM k lumpU k n+ 1 2 · y to the contact impulse k r k n+ 3 2 . When the ratio is sufficiently small, the simulation can be stopped and the static configuration is obtained. Figure 25 shows the evolution of these quantities over time.
The discrete arch problem illustrates the mass scaling feature of explicit schemes in a non-linear problem with impact. And we propose a termination criteria based on impact reaction.

Concluding remarks
These five test cases form a benchmark for explicit schemes dedicated to impact dynamics. Each case focuses on one critical property: impact behaviour and Zeno's paradox for bouncing ball, non-linear behaviour with a damping term for Van der Pol oscillator, conservation of angular momentum for rotating spring, spurious high frequencies behaviour with finite elements discretization for impacting bar, and mass scaling ability for discrete arch. These test cases are light and easy to implement, forming a benchmark well suited to development phase. A missing issue is the solving of contact problem for non-diagonal Delassus operator. A problem with non-conforming meshes should be added in order to complete this benchmark.
With this benchmark, we find that the discrete velocity formulation of HSM conditions is more suited to non-smooth impact dynamics than the displacement formulation. The CD-Lagrange scheme uses this formulation, being explicit and symplectic. It is therefore well-suited to non-linear impact dynamics.
Contact in FE problems forces to consider a dissipative form of HSM conditions, and causes spurious oscillations on velocity for both Moreau-Jean and CD-Lagrange schemes. For spurious oscillations, a dissipative scheme could damp them. Efficient explicit dissipative schemes exist [34,35], but no one for non-smooth dynamic. Such a scheme would be an improvement in this framework. An other improvement would be an discrete formulation for impact FE problems, that would be energy conservative as the implicit ones [20,21] and suitable to explicit schemes.