Explicit Verlet time-integration for a Nitsche-based approximation of elastodynamic contact problems

The aim of the present paper is to study theoretically and numerically the Verlet scheme for the explicit time-integration of elastodynamic problems with a contact condition approximated by Nitsche’s method. This is a continuation of papers (Chouly et al. ESAIM Math Model Numer Anal 49(2), 481–502, 2015; Chouly et al. ESAIM Math Model Numer Anal 49(2), 503–528, 2015) where some implicit schemes (theta-scheme, Newmark and a new hybrid scheme) were proposed and proved to be well-posed and stable under appropriate conditions. A theoretical study of stability is carried out and then illustrated with both numerical experiments and numerical comparison to other existing discretizations of contact problems.


Introduction and problem setting
Explicit time-marching schemes for the dynamics of deformable solids with impact has already been the subject of an abundant literature (see, e.g., [1][2][3] for some recent contributions). They are appealing since they can be easy to implement, fast and adapted to parallel architectures. Nevertheless, there still remains important difficulties to design robust explicit methods and to obtain reliable numerical simulations in this context (see, e.g., [4]). Among these difficulties, numerical stability and energy conservation remains one of the most important ones. Another one is to preserve the quality and the accuracy of the numerical solution, which can present spurious oscillations in the displacement, the velocity or the contact stress. A last one is to enforce properly the contact condition, particularly the non-penetration condition.
A precursory method is the one developed by Taylor and Flanagan [5] in the framework of PRONTO3D software (see also the description in [6]). Nevertheless, the method is not fully explicit, except in a node-to-node contact approximation, in the sense that the contact pressure is computed in an iterative process on the whole contact surface. To mention some other of the most important contributions, we can say that a widely resumed theoretical work in dynamic impact problems is due to Moreau [7,8] for the impact of rigid body systems. The (implicit) schemes proposed by Moreau have been extended quite naturally to the elasticity case through finite element semi-discretization in space (for instance in [9]) which transforms the continuous impact problem into a discrete one very close to a rigid body system. These discrete impact problems, governed by a so-called measure differential inclusion are notoriously ill-posed and of very low regularity.
The ill-posedness can be fixed (for the most part) by the addition of an impact law with a restitution coefficient. As a matter of fact, standard schemes, such as the commonly used ones of Newmark's family [10], have an erratic behavior when they are applied to dynamic contact problems. This is mainly because they select a solution corresponding to an arbitrary (and potentially very large) restitution coefficient (see [11]). Alternatively, a valuable scheme in this context is that of Paoli and Schatzman [12,13] who implicitly takes into account this restitution coefficient. However, the addition of a restitution coefficient can be considered as artificial in the context of deformable solids. This does not diminish the interest for the Paoli-Schatzman scheme which will be a point of comparison with our proposed approach. The implicit inclusion of a restitution coefficient has also been considered in [14] to develop a wide range of schemes based on a time discontinuous Galerkin framework.
As noticed in [11], even in the case where the continuous problem is well-posed (see, e.g., [15,16] for well-posedness results), the ill-posed measure differential inclusion that results from finite element semi-discretization in space has an infinite number of solutions, depending on the choice of a restitution coefficient on each node of the contact boundary. Moreover, it is not possible to decide which solution is more suitable than other. Indeed, the two most remarkable solutions are, on the one hand, the one for a unitary restitution coefficient which ensures conserving energy but which causes very important spurious oscillations of the contact nodes and unexploitable contact stress, and, on the other hand, the solution for a vanishing restitution coefficient which ensures stability and a better approximation of the contact stress but is energy dissipative, while the continuous problem is not. This resulted in [11] to design the mass redistribution method (generalized in [17,18]) which allows a compromise in this context, i.e. well-posedness of the space semidiscretized problem, conservation of the energy and an improved quality of the contact stress. However, and this is also the case for the Paoli-Schatzman scheme, it introduces a global problem to be solved (at least on the contact nodes) when an explicit time-marching scheme is used. In the same spirit, a time-marching scheme has been designed in [19] for dynamic fracture problems, in which the cohesive forces are treated implicitly, while an explicit scheme is used for the dynamics of interior nodes.
For explicit time-integration, primal formulations of contact conditions are better suited. Indeed, since no additional unknown such as a Lagrange multiplier are introduced, they allow to enforce the contact conditions at the previous time-step, instead of the current one, so that the contact term appears at the right-hand side and does not require global (and non-linear) solving. A first possibility is to penalize / regularize the contact conditions (see, e.g., [20,21]): the resulting penalty method is simple to implement and only an inversion of the mass-matrix is needed at each time-step to solve the resulting fully discretized problem (and the scheme becomes fully explicit when the mass matrix is lumped). Nevertheless, the penalty method is not consistent and the choice of the penalty parameter remains a difficulty (see, e.g., [22]). The alternative we explore in this paper is a Nitsche treatment of contact conditions, which is still a primal method, with the same advantages as penalty, but that remains consistent with the original prob-lem, and more robust with respect to the Nitsche parameter. Nitsche's method, originally designed to enforce weakly Dirichlet boundary conditions [23,24], was adapted to unilateral contact in [25,26] (see also [27] for an overview of recent results on this topic).
We studied previously in [28,29] the behavior of Nitsche's method for contact in elastodynamics, when combined to various implicit time-marching schemes. Particularly, when applied to contact-impact in elastodynamics, Nitsche's method has the good property of leading to a well-posed semi-discrete problem in time (i.e., a system of Lipschitz differential equations) as it is shown in [28]. This feature is shared also by the penalty method and modified mass methods. Moreover the symmetric variant of Nitsche's space semi-discretization conserves an augmented energy [28], as does the penalty method [30]. We studied as well theoretically the well-posedness, the stability and energy conservation properties of fully discrete schemes based on space semidiscretization with Nitsche's method combined with the theta-scheme, the Newmark scheme and a new Hybrid scheme. This study was illustrated with some numerical experiments.
The aim of this paper is to study mathematically and numerically the approximation of contact problems in elastodynamics by Nitsche's method combined with the explicit Verlet time-marching scheme. The choice of the Verlet scheme is motivated both by its simplicity and its attractive theoretical properties (symplecticity) [31]. We will also make comparisons with some of the existing methods mentioned above and with the approximation by penalized contact. The numerical comparison will be mainly performed on the one-dimensional problem introduced in [15] whose advantage is to present a known periodic solution and to make clear the occurrence of parasitic oscillations, the convergence and energy conservation properties. Comparisons for 2D and 3D problems will also be presented.
Let us introduce some useful notations. In what follows, bold letters like u, v, indicate vector or tensor valued quantities, while the capital ones (e.g., V, K . . .) represent functional sets involving vector fields. As usual, we denote by (H s (.)) d , s ∈ R, d = 1, 2, 3 the Sobolev spaces in one, two or three space dimensions (see [32]). The usual scalar product of (H s (D)) d is denoted by (·, ·) s,D and the corresponding norm is denoted by · s,D -we keep the same notation when d = 1 or d > 1. The letter C stands for a generic constant, independent of the discretization parameters.
We consider an elastic body in R d with d = 1, 2, 3. Small strain assumptions are made (as well as plane strain when d = 2). The boundary ∂ of is polygonal (d = 2) or polyhedral (d = 3). The normal unit outward vector on ∂ is denoted n. We suppose that ∂ consists in three nonoverlapping parts D , N and the contact boundary C , with meas( D ) > 0 and meas( C ) > 0. In its initial stage, the body is in contact on C with a rigid foundation and we suppose that the unknown contact zone during deformation is included into C . The body is clamped on D for the sake of simplicity. It is subjected to volume forces f in and to surface loads g on N .
We deal with the unilateral contact problem in linear elastodynamics during a period of time where the notationẋ is used for the time-derivative of a vector field x on T , so thatu is the velocity of the elastic body andü its acceleration; u 0 andu 0 are initial displacement and velocity. The density of the elastic material is denoted by ρ, and is supposed to be a constant to simplify the presentation (this is not restrictive and the results can be extended straightforwardly for a variable density). The notation σ = (σ ij ), 1 ≤ i, j ≤ d, stands for the stress tensor field and div denotes the divergence operator of tensor valued functions. The notation ε(v) = (∇v + ∇v T )/2 represents the linearized strain tensor field and A is the fourth order symmetric elasticity tensor having the usual uniform ellipticity and boundedness property. For any displacement field v and for any density of surface forces σ(v)n defined on ∂ we adopt the following notation where v t (resp. σ t (v)) is the tangential component of v (resp. σ(v)n). The conditions describing unilateral contact without friction on CT are: Note additionally that the initial displacement u 0 should satisfy the compatibility Condition u 0n ≤ 0 on C . To our knowledge, the well-posedness of Problems (1), (2) is still an open issue. The few available existence results concern simplified model problems involving the (scalar) wave equation with Signorini's conditions (see, e.g., [16,[33][34][35][36]) or thin structures like membranes, beams (see [37]) or plates (see [38]). Even in these simplified cases, obtaining uniqueness and energy conservation still involves difficulties in 2D or 3D. For a review on some of these results, one can refer to the book [39].
We introduce the Hilbert space and the following forms: for any u and v in V, for all t ∈ [0, T ]. The (total) mechanical energy associated with the solution u of the dynamic contact problem (1, 2) is: Let us take t ∈ [0, T ]. Formally, we get from (1), after multiplication byu(t), integration by parts, with the boundary conditions on DT , NT and the absence of friction: Moreover, if we assume that the contact force does not dissipate any energy, i.e. satisfies the so called persistency Condition σ n (u(t))u n (t) = 0 (see, e.g., [30,40,41]) then we end up with: Notably, when L vanishes, we get energy conservation: Note that, even if it is expected that solutions to Problems (1), (2) satisfy the persistency Condition σ n (u(t))u n (t) = 0 in order to respect the non-dissipative character of the frictionless contact condition, it has only been rigorously proved in a one dimensional framework (elastic bar) for instance in [36,Lemma 2.5].
The rest of our paper is outlined as follows. The first section is dedicated to the description of the fully discrete formulation for dynamic contact with Nitsche and Verlet explicit time-integration. Then, a stability analysis is carried out, and finally, some numerical comparisons with other classical methods are investigated and analysed.

Discrete setting: Nitsche's method with Verlet scheme
We begin this section with preliminary notations and results. Then, we introduce our Nitsche-based finite element semi-discretization in space, and we recall its main properties of well-posedness and energy conservation. Finally we describe the fully discretized problem based on the Verlet explicit time-marching scheme.

Preliminary notations and results
We make use of the notation [·] R − , that stands for the projection onto R − ([x] R − = 1 2 (x − |x|) for x ∈ R). The notation H(·) will stand for the Heaviside function 1 2 if x = 0, and 0 if x < 0, which satisfies H(x) + H(−x) = 1, ∀x ∈ R. Moreover we will make use of the equality H(−x)[x] R − = [x] R − , ∀x ∈ R, and the following property of projection: Let V h ⊂ V be a family of finite dimensional vector spaces (see [42]) indexed by h coming from a family T h of triangulations of the domain (h = max K ∈T h h K where h K is the diameter of the triangle K ). The family of triangulations is supposed: • Regular, i.e., there exists σ > 0 such that ∀K ∈ T h , h K /ρ K ≤ σ where ρ K denotes the radius of the inscribed ball in K , • Conformal to the subdivision of the boundary into D , N and C , which means that a face of an element K ∈ T h is not allowed to have simultaneous non-empty intersection with more than one part of the subdivision, • Quasi-uniform, i.e., there exists c > 0, such that, ∀h > 0, ∀K ∈ T h , h K ≥ ch.
To fix ideas, we choose a standard Lagrange finite element method of degree k with k = 1 or k = 2, i.e.: However, our results would be similar for any C 0 -conforming finite element method.
We consider in what follows γ h , a positive piecewise constant function on the contact interface C which satisfies for every K that has a non-empty intersection of dimension where γ 0 is a positive given constant (the Nitsche parameter). Note that the value of γ h on element intersections has no influence. We next define convenient mesh-dependent norms, in fact weighted L 2 ( C )-norm

Definition 1 For any
Additionally, it will be convenient to endow V h with the following mesh-and parameterdependent scalar product: and note · γ h := (·, ·) 1 2 γ h the corresponding norm. Remark that the two norms · γ h and · 1, are equivalent on V h , in the following sense (for a quasi-uniform mesh T h ): for any v h ∈ V h . The positive constant C comes from the trace inequality and the constant of quasi-uniformity of the mesh T h . For a mesh T h that is not quasi-uniform, the same relationship holds, replacing h by (min K ∈T h h K ).
We end this section with the following statement: a discrete trace inequality (see, e.g., [43]), that is a key ingredient for the whole mathematical analysis of Nitsche's based methods.

Lemma 3
There exists C > 0, independent of the parameter γ 0 and of the mesh size h,

Semi-discrete problem in space
Our Nitsche-based discretization of the contact condition comes from the following result (see [44] and as well [25] for a detailed formal proof).

Proposition 4
Let γ be a positive function defined on C . The contact Condition (2) can be reformulated as follows: As in [28,29] we will consider a family of methods indexed by a parameter ∈ R (with, in general, = −1, 0, 1, see, e.g., [26]). Let us introduce the discrete linear operator P n ,γ h : Define as well the bilinear form: The space semi-discretized Nitsche-based method for unilateral contact problems in elastodynamics then reads (see, e.g, [27,28]): where u h 0 (resp.u h 0 ) is an approximation in V h of the initial displacement u 0 (resp. the initial velocityu 0 ), for instance the Lagrange interpolant or the L 2 ( ) projection of u 0 (resp.u 0 ).

Remark 5
Note that, as in [27], we adopted in this presentation a different convention for notations compared to previous works [28,29]. This is in order to get closer to the formulations provided in most of the papers on Nitsche's method and on the augmented Lagrangian method.
We can reformulate (8) as a system of (non-linear) second-order differential equations. To this purpose, using Riesz's representation theorem in (V h , (·, ·) γ h ) we first introduce the mass operator M h : Still using Riesz's representation theorem, we define the (non-linear) operator B h : V h → V h , by means of the formula Moreover, we recall the results of well-posedness and the energy estimate for the semidiscrete problem in space, that were established in [28]. First, the following theorem together with the boundedness of (M h ) −1 γ h [see [28] show that Problem (8) or equivalently Problem (9)] is well-posed.

Theorem 6
The operator B h is Lipschitz-continuous in the following sense: there exists a constant C > 0, independent of h, and γ 0 such that, for As a consequence, for every value of ∈ R and γ 0 > 0, Problem Remark 7 Note that, conversely to the static case (see [25,26,45]) and the fully-discrete case there is no condition on γ 0 for the space (semi-)discretization, which remains wellposed even if γ 0 is arbitrarily small. However, this does not imply that the solution remains consistent when γ 0 becomes small (see Remark 19 and Fig. 4 in the sequel).
We recall that the standard (mixed) finite element semi-discretization for elastodynamics with unilateral contact leads to ill-posed problems (see, e.g., [11,22]), which is not the case of Nitsche's formulation that leads to a well-posed (Lipschitz) system of differential equations. This feature is shared with the standard penalty method, the difference being that Nitsche's method remains consistent in a strong sense (see [28]). Note that the standard (mixed) finite element semi-discretization is consistent as well as the singular dynamic method introduced in [18]. The mass redistribution method introduced in [11] is asymptotically consistent when h vanishes. Now, we consider the energy estimates which are counterparts of the Eq. (3), in the semi-discretized case. Let us define the discrete energy as follows: which is associated to the solution u h (t) to Problem (8). This is the direct transposition of the mechanical energy E(t) from the continuous system. As in [28], we define also a modified energy more suited to Nitsche's method in which a consistent term is added. This term denoted R h (t) represents, roughly speaking, the non-fulfillment of the contact Condition (7) by u h . The relationship between E h (t) and E h (t) is provided in the Lemma below: Proposition 8 For ≥ 0 and γ 0 large enough, there exists C > 0 independent of h, of γ 0 and of the solution to Problem (8), such that, for all t ∈ [0, T ]: Proof This result is obtained using the coercivity of a(·, ·) and applying Lemma 3.

Remark 9 Proposition 8 states that the energy E h (t) remains always positive (if E h (0) is)
for ≥ 0 and γ 0 large enough. For γ 0 small, the existence of zero energy spurious modes cannot be excluded.
Remark 10 For < 0, such a result with a constant independent of the mesh parameter h cannot be obtained. As a consequence, for < 0, the quantity E h (t) cannot be used for optimal energy evolution estimates and might become even negative for h small.
Still in [28], the following evolution of E h is obtained: Theorem 11 Suppose that the system associated to (1,2) is conservative, i.e., that L(t) ≡ 0 for all t ∈ [0, T ]. The solution u h to (8) then satisfies the following identity: Notably, when = 1, we get for any t This result links the non-satisfaction of the energy conservation to the non-satisfaction of the so-called persistency condition. However, it appears in the present study that it would be preferable to use E h 1 (t) even for the variants = 1 (see Remark 10), for which the following result can be established: The solution u h to (8) then satisfies the following identity: (8): which we reformulate as: We split the second term, use P n The result is obtained by re-ordering the terms, using the property d

Remark 13
The above result still states that E h 1 (t) is conserved for the symmetric variant = 1, and, for = 1 the variations of E h 1 (t) come from the non-fulfillment of the contact Condition (7) by u h .
The discretization of Problem (9) with the velocity-Verlet scheme reads: This scheme corresponds to the variant of the Newmark scheme with γ = 1 2 and β = 0 (see, e.g., [22,28]). As a result, this is a second order consistent scheme in τ . Note that, for practical implementation, the acceleration can be eliminated using the relationship M hüh,n = L h,n − B h u h,n . This result into the following reformulation, where the only unknowns are the displacement and the velocity: Since this scheme only requires the inversion of the mass matrix M h at each time-step, it becomes then fully explicit when the mass matrix M h is lumped.
We can transform the Scheme (13) into a two-steps scheme. Indeed, the first equation of (13) applied for time-steps n and n + 1 reads: So, using the above relationships, the second equation of (13), at time-step n, can be written as which can be simplified as: We add the above equation to the second relationship in (15) and get Using finally the third equation in (13), combined with the above relationship, we obtain: This is a two-steps scheme, so called Störmer-Verlet scheme or central difference scheme, that involves only the displacement as an unknown (the first step u h,1 is classically obtained via the first equation of (14) at n = 0). Note finally that Leapfrog scheme is also equivalent to Verlet one (see e.g. [31,46]): it suffices to define an intermediate velocityu with initial Conditions u h,0 = u h 0 ,u h, [1/2]

Stability properties of Verlet scheme
First, we present different energies associated to the solution to Problem (13), and make explicit their relationships. Then, we derive energy estimates associated to the fully discrete Problem (13), and a (non-optimal) stability result is deduced. We conclude with some comments and arguments that a better result may be expected.

Discrete energies
We next define the following energy: which is associated with the solution u h,n to Problem (13). Set also Note that the energies E h,n and E h,n are the fully discrete counterparts of the semi-discrete energies E h (t) and E h (t). Note additionally that a result similar to Proposition 8 can be established, that allows to bound the physical energy by the augmented energy under appropriate conditions on the numerical parameters (and the statements of Remark 9 and Remark 10 still applies):

Proposition 14
For ≥ 0 and γ 0 large enough, there exists C > 0 independent of h, of γ 0 and of the solution to Problem (13), such that, for all n = 0, . . . , N : To simplify slightly the notations in the energy estimates below, we will make use of the convention: P n := P n 1,γ h (u h,n ) for any n ∈ N. We will denote as well by [·] R + the projection onto R + . Additionally, we define a modified energy adapted to the Verlet scheme: Of course, this variant of the energy definition makes also sense and is usable for stability analysis only if it can be used to bound the physical energy E h,n . This is the aim of the following result:

Proposition 15
Suppose that L n ≡ 0 for all n = 0, . . . , N and that the mesh T h is quasiuniform. Then, for γ 0 large enough and τ h small enough, there exists C > 0, independent of h, of γ 0 and of the solution to Problem (13), such that, for all n = 0, . . . , N : Proof We suppose L n ≡ 0 and take v h =ü h,n in Problem (13): We assume that γ 0 is large enough, then use the Cauchy-Schwarz inequality, the Lemma 3 and the inverse inequality [47, Corollary 1.141, Remark 1.143] to bound the first right term: with C > 0 independent of γ 0 and h. We use the Cauchy-Schwarz inequality for the second term: Using once more Lemma 3 we bound: We insert the above bounds into (18), take into account that γ 0 is large, apply once again the inverse inequality and get: As a result, we obtain which allows to conclude that E h,n 1 ≤ CE h,τ ,n for τ h small enough using the coercivity of a(·, ·). The estimate E h,n ≤ CE h,τ ,n is then deduced from Proposition 14.

Proposition 16
Suppose that L n ≡ 0 for all n ≥ 0. The following energy identity holds for all n ≥ 0: This result can be easily adapted as follows when the energy E h,n 1 is considered instead, even for = 1:

Proposition 17
Suppose that L n ≡ 0 for all n ≥ 0. The following energy identity holds for all n ≥ 0: Proof This identity is deduced from combined with the following rearrangement (where we use the identity As an interesting consequence, we obtain the following result for the discrete energy E h,τ ,n by simplifying the previous one: Proposition 18 Suppose that L n ≡ 0 for all n ≥ 0. The following energy identity holds for all n ≥ 0: Proof We use Eq. (13) to rewrite: Then we just use the above identity in (21).

Remark 19
For γ 0 small, the property of Proposition 15 can be lost and the energy E h,n 1 may become negative. In that case, some deformation corresponding to a negative energy may exist, which is of course a non-physical situation. This highlights that, even thought the semi-discrete problem (9) has a unique solution for γ 0 small, the reliability of the discretization is guaranteed only for γ 0 large enough.
Remark 20 Still referring to [29,Proposition 3.4], and instead of Verlet scheme, if we consider the explicit Newmark scheme γ = 1 and β = 0 and = 1 as Nitsche's variant, the pending energy evolution corresponding to Proposition 18 in that case involves the sole . This term being non-positive, the stability of this explicit scheme can be established thanks to Proposition 15 for τ h small enough.

Stability analysis in the case = 1
The main result of this section is the following (non-optimal) stability result for the Scheme (13) in the case = 1:

Proposition 21
Suppose that L n ≡ 0 for all n ≥ 0 and that the mesh T h is quasi-uniform. Then, for = 1, for γ 0 large enough and for with C > 0 independent of γ 0 , h and τ , the energies E h,τ ,n and E h,n remain bounded.
Proof We already know from Proposition 18 that, for = 1: We note that P n+1 Since the mesh is quasi-uniform, we then use Lemma 3, the trace inequality of [48, Following exactly the same path as above, but using · 0, ≤ · 1, after the trace inequality, we bound also: We combine the above result with the bound (19). This yields: This results into: , from which we can deduce, using Proposition 15: This means that, still with N = T τ , which remains bounded under the assumption (23). The boundedness of E h,N is then deduced from Proposition 15.

Comments on the stability analysis
The stability result given by Proposition 21 is submitted to a CFL Condition τ = O(h 2 ). Of course, this is not the expected one which would be τ = O(h) in accordance with the result of Proposition 15 and the stability analysis of Verlet scheme for a linear evolution equation (see, e.g., [46]). The reason of the non-optimality of Proposition 21 is that we did not succeed to optimally bound the involved contact term P n+1 However, an important remark is that this term vanishes unless the terms P n and P n+1 are of opposite signs, which occurs only when the contact status changes. Moreover it is positive only when the status changes from contact to non-contact. If we assume that the number of such transitions is finite during the simulation, a stability result with a Condition τ = O(h) may be recovered. However, an infinite number of changes of the contact status cannot be excluded. Another argument in favor of such a stability result is obtained via the definition of the following linear (but depending on u h,n ) operator B h,n : V h → V h : x for x ∈ R, there holds: The following bound on the norm of B h,n can be established, following an argument similar to [28, Theorem 2.8]:

Lemma 22
Let us suppose that γ 0 is large enough. There exists a constant C > 0, independent of , γ 0 and h, such that where · γ h is the operator norm induced by the norm · γ h on V h .
Proof Let us take v h and w h in V h . First, using Lemma 3 there holds and the same bound holds for P n 1,γ h (v h ) 0, C , replacing | | by 1. Then, using Lemma 3 and the above result, we bound: which proves the assertion (24).
Using B h,n we can rewrite the operator associated to the Scheme (14) as: where I is the identity operator. It links the successive values u h,n+1 , u h,n and u h,n−1 thanks to u h,n+1 u h,n = C h,n u h,n u h,n−1 , when neglecting the source term L h,n . For a linear problem this would correspond to the amplification matrix. The following result can be established for C h,n :

Proposition 23
Let us suppose that the mesh T h is quasi-uniform, that γ 0 is large and where C > 0 is a constant independent of γ 0 , h and τ . Then, C h,n is diagonalizable and its spectral radius ρ(C h,n ) is equal to 1. Furthermore, the same conclusion holds for C h,n C h,n+1 , which is diagonalizable with a spectral radius ρ(C h,n C h,n+1 ) also equal to 1.
Proof Let us consider λ ∈ C an eigenvalue of the operator C h,n . Denoting With help of the second equation, the first one can be written λA h,n w h − w h = λ 2 w h , or equivalently We then use [28, Lemma A.1] and Lemma 22 to bound We now consider γ 1 2 0 τ /h small enough so that the eigenvalues of A h,n = (2I − τ 2 (M h ) −1 B h,n ) are positive. We call α j these eigenvalues and w h j the corresponding eigenvectors. Taking w h = w h j in Eq. (25), we deduce that, for each index j, we can compute two eigenvalues λ ± j which are solution to λ 2 − α j λ + 1 = 0.
Since the eigenvalues of (M h ) −1 B h,n are all positive, we infer α j < 2, and j = α 2 j − 4 < 0. Therefore the above algebraic equation has two imaginary solutions Remark that these eigenvalues are such that This allows to conclude that C h,n is diagonalizable and ρ(C h,n ) = 1. We can now make a similar computation for two successive iterations. Since If we denote by β j the eigenvalues of A h,n+1 A h,n which are close to 4 for γ 1 2 0 τ /h small enough, and v h j the corresponding eigenvectors, we conclude that the eigenvalues of C h,n+1 C h,n are Since |λ ± j | = 1, the operator C h,n+1 C h,n is diagonalizable with a unit spectral radius.

Remark 24
For a linear problem, we would conclude that the scheme is stable, under the Condition τ h small enough. However, in a nonlinear framework, the conclusion cannot be drawn since the matrix C h,n changes from an iteration to another. Moreover, it seems difficult to pursue the reasoning made on two iterations to an arbitrary number of iterations.

Numerical experiments
We first carry out numerical experiments in 1D, where we can compare our results with an exact solution. Then, numerical experiments in 2D/3D will be described. These numerical tests are performed with the help of our freely available finite element library GetFEM++ (see http://getfem.org).

1D numerical experiments: multiple impacts of an elastic bar
We first present the setting, and then the results obtained by combination of Nitsche's contact discretization and Verlet scheme. These results are also compared with computations using other methods: the Paoli-Schatzman scheme, the Taylor-Flanagan scheme, the mass redistribution method and the penalty method. This section is ended with numerical convergence tests.

Setting
We first deal with the one-dimensional case d = 1 with a single contact point, namely an elastic bar = (0, L) with C = {0}, D = {L} and N = ∅. The elastodynamic equation is then reduced to find u : T = (0, T ] × (0, L) → R such that: where E is the Young modulus and the Cauchy stress tensor is given by σ (u) = E(∂u/∂x). Note that σ n (u) = (σ (u)n) · n = σ (u) on C . In this case, Problem (1, 2) admits one unique solution (see e.g. [36]) for which the following energy conservation equation holds, for t a.e. in (0, T ]: We consider a finite element space using linear (k = 1) or quadratic (k = 2) finite elements and a uniform subdivision of [0, L] with M segments (so L = Mh). We denote the vector which contains all the nodal values of u h,n (resp.u h,n andü h,n ) by U n := [U n 0 , . . . , U n N ] T (resp.U n ,Ü n ). The component of index 0 corresponds to the node at the contact point C . We also note M, resp. K, the mass, resp. the stiffness, matrix that results from the finite element discretization. We introduce the Courant number which is defined as: where c 0 is the wave speed associated to the bar. For each simulation, we compute and plot the following time-dependent quantities: 1. The displacement u at the contact point C , given at time t n by u h,n (0)(= U n 0 ). 2. The contact pressure σ C , which, in the discrete case, is different from σ (u). If a standard (mixed) method is used for the treatment of contact, it is directly given by the Lagrange multiplier, i.e., σ n C := λ h,n at time t n . In the case of the Nitsche-based formulation, it can be computed as follows at time t n : which comes from the contact Condition (7). 3. The discrete energy E h which is at time t n E h,n = 1 2 (U n ) T MU n + (U n ) T KU n , and the discrete augmented energy E h 1 : We propose a benchmark associated to multiple impacts. This allows to check both the presence of spurious oscillations and the long term energetic behavior of the method. In the absence of external volume forces, the bar is initially compressed. Then, it is released without initial velocity. It impacts first the rigid ground, located at x = 0, and then gets compressed again. We take the following values for the parameters: This problem admits a closed-form solution u which derivation and expression is detailed in [15]. Notably, it has a periodic motion of period 3. At each period, the bar stays in contact with the rigid ground during one time unit (see Fig. 1). The chosen simulation time is T = 12, so that we can observe 4 successive impacts.

Numerical results for Nitsche's method
We discretize the bar with M = 20 linear finite elements (k = 1, h = 0.05) and take τ = 0.01. The resulting Courant number is ν C = 0.2. Note that almost all the parameters have been taken identical as in [29] for comparison purposes. The number of element is smaller (M = 20 instead of 100 in [29]) and the time-step τ is 0.01 for stability reasons. We first investigate the variant = 0 with a parameter γ 0 = 1. The mass matrix is computed in a standard fashion. The choice γ 0 equal to 1 is guided by the concern to obtain a Fig. 1 Multiple impacts of an elastic bar. The bar is clamped at x = L and the contact node is located at the bottom. The closed-form solution is periodic of period 3, with one impact during each period (here between t = 1 and t = 2) The results are depicted in Fig. 2 where the approximated solution corresponds to the solid blue line and the red dotted line is the exact solution. Note also that the dashed energy is E h 1 , the modified energy given by expression (11). For an element size which remains relatively rough, one notes the good approximation of the displacement. Significant oscillations can be deplored on the velocity of the contact point, but unfortunately they are very difficult to avoid. Indeed, since a velocity shock is propagating without attenuation or dissipation in the proposed test case, the Gibbs phenomena are inevitable in the finite element approximation. It would be the case even without a contact condition. The approximation of the contact stress is, although polluted by oscillations too, of good quality given the discontinuous character of the exact solution and the relatively coarse mesh size (one can compare with the results obtained in [29] for elements five times smaller). The evolution of the energy reveals some variations which are far from being negligible, but remain moderate, and without appearance of instabilities. Moreover these variations tend to decrease for a finer element size.
The calculation for the variant = −1 and for Nitsche's parameter γ 0 still equal to 1 is presented on Fig. 3. It can be seen that the non-penetration condition is slightly better  Fig. 4, one gets a non convergent approximation of the displacement. This is due to the loss of coercivity arising when the assumption of Proposition 14 is not satisfied and probably to the existence of zero-energy modes. Indeed, the variant = 1 is the most restrictive from this point of view and needs a larger parameter γ 0 . Figure 5 represents the simulation for γ 0 = 2 which allows to recover the coercivity. We observe the very good conservation of the discrete energy together with a good approximation of the nonpenetration condition. The level of oscillations on the contact point velocity and on the contact stress is similar with the case = 0.  Figure 6 displays a simulation still for ( = 1, γ 0 = 2) but using a lumped mass matrix. Following a standard strategy for P 1 Lagrange finite elements, on each row, the extra-diagonal components of the mass matrix are set to zero and added to the diagonal component. The comparison with Fig. 5 allows to notice more pronounced oscillations on the displacement, the velocity and the contact stress, but still with a very good energy conservation.
Finally, Figs. 7 and 8 show the evolution of the solution for decreasing discretization parameters, and for the variant = 1, γ 0 = 2, and the standard mass matrix. We note in Fig. 7 a rapid decrease of the oscillations on the displacement with the refinement of the discretization. Conversely, the convergence of the contact stress as depicted in Fig. 8 is rather slow, as it could be expected from the very low regularity of the solution. Indeed we observe a very gradual decrease in the amplitude of the oscillations.

Comparison with Paoli-Schatzman scheme
The Paoli-Schatzman scheme directly addresses the measure differential inclusion that results from the finite element semi-discretization of the dynamic contact problem. Following [12,13,[49][50][51], it can be summarized as follows. Let U h adm be the approximated set of admissible displacements such that U ∈ U h adm means that the vector of degrees of freedom U satisfies a chosen approximated non-penetration condition. In our onedimensional test case, this can be simply written U 0 ≥ 0, where U 0 is the displacement of the contact point. Then, still denoting M, resp. K, the mass, resp. stiffness, matrices that result from the finite element discretization, the (generally ill-posed) measure differential inclusion resulting from the semi-discretization of the dynamic contact Problem (1), (2) can be written: Of course, one easily recognize the Verlet scheme except for the contact constraint which is taken in an implicit manner and prescribed to the intermediate displacement U n+1,e = U n+1 + eU n−1 1 + e .
When implementation is considered, one usually introduces a multiplier to prescribe the constraint. Denoting by B the matrix having N c lines such that the non-penetration condition reads (BU) i ≤ 0, i = 1 . . . N c , the scheme may be re-written or, as a one-step scheme still with the addition of the complementarity conditions (29). The proof that the restitution coefficient e is asymptotically reached for a vanishing time-step is detailed in [12,13]. Note that, even though Verlet scheme is an explicit scheme, the implicitation of the contact force in Paoli-Schatzman scheme results in a global problem to be solved on the contact surface at each time-step. Of course, this corresponds to a scalar algebraic equation which can be explicitly solved in the one-dimensional test.
The numerical tests for h = 0.05 and τ = 0.01 are presented in Figs. 9, 10 and 11 for a restitution coefficient equal to 0, 1/2 and 1, respectively. The results for e = 0 and e = 1/2 are very similar to each other despite the difference between the restitution coefficients,

Comparison with Taylor-Flanagan scheme
Taylor and Flanagan [5] developed an explicit scheme in the framework of PRONTO3D software which rapidly became a reference for explicit integration of contact and impact problems. To summarize the principle of the method, it is based on the Leapfrog form of Verlet scheme (17). When contact occurs, the persistency condition is prescribed at the half time-step by enforcing the relative velocity to vanish (see Eq. (21) in [6], for instance). To this aim, a Lagrange multiplier is introduced which is taken into account in an implicit way. The equation associated to the Lagrange multiplier can be solved locally only in a node-to-node contact approximation. For a more general contact condition, the Lagrange multiplier is obtained by solving a global problem on the contact surface. However, Taylor and Flanagan propose an iterative method to compute the Lagrange multiplier which needs only a few iterations.
Since the Taylor-Flanagan scheme prescribes the contact condition with an implicited Lagrange multiplier and enforces the persistency condition, it is very close to the Paoli-Schatzman scheme with a restitution coefficient e = 0 even if the contact condition is prescribed in a slightly different way. The consequence is that the results of the simulations shown on Fig. 12 for the Taylor-Flanagan scheme are almost identical to the results shown on Fig. 9 for the Paoli-Schatzman scheme with e = 0. Particularly, a loss of energy occurs at each impact.

Comparison with the mass redistribution method
The mass redistribution method, introduced in [11], considers a semi-discretization that comes from the finite element approximation of the dynamic contact problem combined with a Lagrange multiplier method to enforce the contact condition: for a.e. t ∈ (0, T ], still with K the stiffness matrix, with B the matrix representing the discrete normal trace operator on the contact boundary, and with M r a modified mass matrix with a vanishing contribution on the contact boundary. The matrix M r is simply built from the standard mass matrix M by setting to zero the lines and columns corresponding to the degrees of freedom on the contact boundary and redistributing the removed mass on the internal degrees of freedom (see a possible redistribution algorithm in [11] and other strategies in [17,52,53]). In the one-dimensional test-case, this just means setting to zero the first column and first row and adding the removed mass on the first degree of freedom, which as been proved to be the most optimal strategy in [52]. It is proved in [11] that the mass redistribution method allows to recover the well-posedness of the discretization and that the solution to the approximated problem is energy conserving. Some convergence results can be found in [15,17,52,54,55]. Note also that such singular mass matrices can be obtained with the method introduced in [18] by considering different approximations for the displacement and the velocity instead of using a post-modification of the mass matrix.
Since the mass matrix admits a kernel containing the vectors being only nonzero on the contact boundary, the system (32), (33) consists in an algebraic variational inequality when reduced on this kernel. Due to the Lipschitz continuity, with respect to the data, of the solution to this variational inequality, Problems (32), (33) reduces to a system of ordinary differential equations on the orthogonal of the kernel. This property, detailed in [11] allows to use quite arbitrary time-marching schemes to approximate (32), (33), among others the Verlet scheme. Of course, the method is not strictly an explicit one since a global solving has to be done on the kernel of the modified mass matrix. However, in the one-dimensional test case, this kernel is one-dimensional which allows an explicit solving.
The corresponding simulations can be seen on Fig. 13. One characteristic of the mass redistribution method is to produce low oscillating velocity and contact stress compared to other discretizations. One can see that the energy is conserved, but slightly modified compared to the standard energy.

Comparison with the penalty method
The penalty method is one of the simplest and most popular way to approximate the contact condition (see for instance [21]). With the notations used to write system (9) for Nitsche's method and with the non-linear operator B h p : the dynamic contact problem approximated by a penalty method reads: This problem corresponds to a nonlinear system of ordinary differential equations on which a Verlet scheme can be applied. The parameter γ h is now the penalty parameter and, following for instance the analysis in the static case of [56] and similarly to Nitsche's method, we will still consider that γ h | K ∩ C = γ 0 /h K for each finite element K .
Figures 14, 15 and 16 depict three simulations for the one-dimensional test case for γ 0 = 5, γ 0 = 1 and γ 0 = 0.25, respectively. The augmented energy associated to penalty is the following one [30]: For γ 0 = 0.25 the contact interface is clearly too soft and a large penetration occurs, which makes the approximated solution being far from the exact one. Conversely, for γ 0 = 5, the non-penetration condition is better approximated, but some important oscillations on the velocity of the contact point and on the contact stress occur. Similarly to Nitsche's method, an acceptable compromise seems to set γ 0 = 1, which corresponds to comparable stiffnesses on the contact point due to both the penalty term, on one side, and to the interior elasticity terms, on the other.
It is worth comparing Fig. 15 to Figs. 2, 3 and 4 for Nitsche's method and the same value of the parameter γ 0 . The non-penetration condition is better satisfied with Nitsche's method, which highlights its consistency. However, energy conservation is better preserved by the penalty method except when the variant = 1 of Nitsche's method is used.

Numerical convergence
Simulations in the previous sections allow a qualitative comparison of the different studied methods on the one-dimensional test-case. The aim of this section is to complete this  comparison with a convergence study, still for the same test-case. This is done for both linear finite elements (P 1 Lagrange) on Fig. 17 and quadratic finite elements (P 2 Lagrange) on Fig. 18. For M the number of elements, h = 1/M is the element size and the time-step is chosen to be τ = h/10 for P 1 elements and τ = h/20 for P 2 elements.
A comparison of Figs. 17 and 18 leads to the conclusion that despite the very low regularity of the exact solution (velocity and stress are discontinuous), there is a substantial gain in using quadratic elements. It even improves the convergence rate for the L 2 (0, T, L 2 ( ))norm of the error of the displacement. Globally, the mass redistribution with quadratic elements provides the best compromise. However, as the Paoli-Schatzman scheme, it necessitates to solve a global problem on the contact surface.
So, if we limit the comparison to primal discretizations, which do not require to solve such a global problem (except the inversion of the mass matrix if the mass matrix is not lumped), we can compare only Nitsche and penalty methods. What can be observed, espe- cially on the L 2 (0, T, L 2 ( ))-norm of the error in displacement, is that Nitsche's method is less sensitive to the parameter γ 0 due to its consistency. The difficulty to achieve a good compromise for penalty method is illustrated on Fig. 17: a low γ 0 allows a good approximation of the contact stress, but the worst approximation of the L 2 (0, T, H 1 ( ))-norm of the displacement for coarse meshes. Conversely, a large γ 0 leads to a better approximation for the L 2 (0, T, H 1 ( ))-norm but a too much oscillatory solution which prevents the L 2 (0, T ) convergence of the contact stress. There is, however, also a constraint for the choice of γ 0 with Nitsche's method because it has to be chosen sufficiently large to preserve the coercivity of the formulation (see Remark 19 and Fig. 4).

2D/3D numerical experiments: multiple impacts of a disc / a sphere
Numerical experiments are then carried out in 2D and 3D, to assess the behavior of Nitsche's method in a more realistic situation. We study the impact of a disc and a sphere on a rigid support. The physical parameters are the following: the diameter of the disc is D = 40, the Lamé coefficients are λ = 30 and μ = 30, and the material density is ρ = 1. The total simulation time is T = 120.
The volume load in the vertical direction is set to f = 0.1 (gravity, oriented towards the support). On the upper part of the boundary is applied a homogeneous Neumann condition g = 0 and the lower part of the boundary is the contact zone C . There an initial vertical displacement (u 0 = (0, 4)) and no initial velocity (u 0 = 0). In such a  For space semi-discretization, Lagrange isoparametric finite elements of order k = 2 have been used. The mesh size is h = 4 for the ball and h = 8 for the sphere (see Fig. 19). Integrals of the non-linear term on C are computed with standard quadrature formulas  The comparison of the simulations for the different methods is depicted Fig. 22 for the two-dimensional case and Fig. 23 for the three-dimensional case. For the sake of shortness, only the penalty, the singular mass and Nitsche methods are compared. First of all, a conclusion that can be drawn from these numerical experiments is that the tested methods are all capable of reliably approximating two and three-dimensional dynamic contact problems. An important difference between simulations in dimension 2 and 3 is a much smaller oscillation of the contact stress in dimension 3, except for the mass redistribution method which is not subjected to spurious oscillations. The energy is conserved more strictly with the penalty method, the mass redistribution method and the variant = 1 of Nitsche's method, the other two variants presenting significant disturbances in the energy evolution. The mass redistribution method appears to give the best compromise between energy conservation and the low level of oscillation on the contact boundary. Note however that it produces a weakening of the rebound, mainly in dimension 3, which we do not explain. The lake of consistency of the penalty method is illustrated on the normal displacement graph where we can note a larger interpenetration compared to the other methods. Finally, among the variants of Nitsche's method, the symmetric variant = 1 is the one that achieves the best compromise between energy

Concluding remarks
In this paper, we studied the application of an explicit Verlet scheme for the approximation of elastodynamic contact problems with Nitsche's method. The explicit method being commonly used in elastodynamic contact problems, it seemed important to complete the study that had been performed in [28,29] for implicit schemes. We tried to characterize the stability properties of the different variants of Nitsche's method for the Verlet scheme and we introduced a number of necessary tools for this analysis. Of course, we are aware that the stability result we establish is very partial (only for = 1) and certainly suboptimal: a Fig. 23 Comparison of the penalty, singular mass and Nitsche methods in the 3D case stability condition such as τ = O(h) would be more satisfactory and would correspond to what we noted in numerical tests. This result remains to be refined. Moreover, it would certainly be possible to prove a convergence result in dimension one, as in [15], because in this context the existence and uniqueness of the solution is theoretically proven. We numerically compared the Nitsche method to the main existing methods that can support an explicit scheme: the Paoli-Schatzman scheme, the Taylor-Flanagan scheme, the mass redistribution method and the penalty method.
We first performed this comparison on a one-dimensional test case whose exact solution consists of a shock wave indefinitely travelling between the two ends of a bar. We can see globally, by comparing the Figs. 2,3,4,5,6,7,8,9,10,11,12,13,14,15 and 16 that Nitsche's method, especially the variant corresponding to = 1, yields an approximation of a comparable quality as the one obtained with the schemes using an implicitation of the contact force (Paoli-Schatzman scheme, Taylor-Flanagan scheme, mass redistribution). Only the mass redistribution method results in lower oscillation levels. Compared to the penalty method, the oscillation levels are of the same magnitude, but contact penetration is more limited. This is reflected in Figs. 17 and 18 by quite similar convergence rates for all the methods, except for penalty. For this latter, a compromise remains difficult to find between a large penalization coefficient, which corresponds to a good approximation of the displacement but a poor approximation of the contact stress, and a small penalization coefficient, for which the interpenetration becomes too large. The decisive advantage of Nitsche's method over the other methods, with the exception of the penalty method, is to be a primal method for which there is no need for an implicit resolution of the contact force. This allows a really explicit resolution in case of lumped mass matrix. The 2D and 3D test cases we performed also confirm the good behavior of Nitsche's method. We can see in Figs. 22 and 23 the advantage in comparison to the penalty method in term of interpenetration, which is less. Still some better approximation results are obtained for the variant = 1.
We can thus conclude that, among the variants of Nitsche's method, the symmetric variant = 1 seems to be the most suitable for solving dynamic contact problems mainly because of its energy conservation properties. For the other variants a gain of energy can be observed, especially for low values of Nitsche's parameter γ 0 . Some perspectives of this work could be to gain further insight into the properties of energy conservation, for instance using other definitions of the discrete energies, such as in [57, Theorem 4.1, Remark 4.1] in which an energy that remains positive irrespectively of the value of numerical parameters is introduced. Also some new explicit time-marching schemes endowed with appealing properties of energy conservation could be considered (see, e.g., [58]). Moreover, further study of the effect of the mass matrix lumping, particularly on the stability of the method, and of the proper choice of the Nitsche's parameter γ 0 are other perspectives of this work.