Élie Bretin
email: elie.bretin@insa-lyon.fr

Yves Renard
email: yves.renard@insa-lyon.fr

Stable IMEX schemes for a Nitsche-based approximation of elastodynamic contact problems. Selective mass scaling interpretation

Keywords: unilateral contact, elastodynamics, Nitsche's method, IMEX schemes, stability, finite element method, selective mass scaling

We introduce some IMEX schemes (implicit-explicit schemes with an implicit term being linear) for approximating elastodynamic contact problems when the contact condition is taken into account with a Nitsche method. We develop a theoretical and numerical study of the properties of the schemes, especially in terms of stability, provide some numerical comparisons with standard explicit and implicit scheme and propose some improvements to obtain a more reliable approximation of motion for large time steps. We also show how selective mass scaling techniques can be interpreted as IMEX schemes.

Introduction

This paper concerns the construction of implicit-explicit time integration schemes for the dynamics of deformable solids that can impact each other. The main addressed issue is to build schemes close to the computational cost of explicit ones but allowing not to be under the constraint of a CFL condition, i.e. allowing the use of large or very large time steps. One characteristic of the dynamic with impact of deformable bodies is the very low regularity of the solutions and the potential ill-posedness of the semi-discretized problem (see the analysis in [START_REF] Houari | Mass redistribution method for finite element contact problems in elastodynamics[END_REF] and the discussion in [START_REF] Chouly | Explicit Verlet time-integration for a Nitsche-based approximation of elastodynamic contact problems[END_REF] for example). A consequence is that time integration schemes must be chosen carefully, since most schemes, even the implicit ones, are subject to instabilities (see [START_REF] Houari | Mass redistribution method for finite element contact problems in elastodynamics[END_REF][START_REF] Hager | A stable energy-conserving approach for frictional contact problems based on quadrature formulas[END_REF][START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF]), except the most dissipative ones such as implicit Euler scheme (see [START_REF] Boumediène | Problèmes de contact unilatéral avec frottement de Coulomb en élastostatique et élastodynamique[END_REF]).

This ill-posed character can then be adressed by adding an impact law, which is a classical approach in the rigid body case, but lacks a clear physical interpretation in the context of deformable bodies. Most of the stable schemes that have been developed so far correspond to a vanishing restitution coefficient, which implies that they generally dissipate energy at each impact, regardless of the size of the time step (this dissipation decreases, however, when the mesh size decreases). Among the first stable schemes proposed are those of LM Taylor and DP Flanagan [START_REF] Taylor | PRONTO3D: A three-dimensionnal transient solid dynamics program[END_REF] (see also [START_REF] Heinstein | Contact-impact modeling in explicit transient dynamics[END_REF]), which are explicit schemes based on Verlet's scheme (also called Leapfrog or central difference scheme) for which the contact force is treated implicitly. The scheme is not fully explicit, in the sense that it remains a non-linear problem to solve at each time step, restricted to the contact boundary, which, in [START_REF] Heinstein | Contact-impact modeling in explicit transient dynamics[END_REF], is solved iteratively. J.J. Moreau's work on the sweeping process and its numerical approximation [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF][START_REF] Moreau | Numerical aspects of the sweeping process[END_REF] has led to many developments, mainly in the context of rigid bodies but also in the context of deformable ones (for instance in [START_REF] Vola | Friction and instability of steady sliding: squeal of a rubber/glass contact[END_REF]). The developed schemes are implicit ones, often based on an expression of the contact condition in term of sliding velocity and an implicit consideration of the coefficient of restitution (see also [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: definition and outlook[END_REF] for a generalization using time discontinuous Galerkin schemes). In parallel, L. Paoli and M. Schatzman also developed and mathematically analyzed in [START_REF] Paoli | A numerical scheme for impact problems. I. The one-dimensional case[END_REF][START_REF] Paoli | A numerical scheme for impact problems. II. The multidimensional case[END_REF] central difference schemes with an implicitation of contact condition which also implicitly takes into account the restitution coefficient.

Notice that the ill-posed character of the finite element semi-discretization is not present in the case of the approximation of the contact condition by a penalty method [START_REF] Belhytschko | Contact-impact by the pinball algorithm with penaly and lagrangian methods[END_REF][START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF][START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF]. The penalty method is however not consistent in the strong sense and induces an additional approximation. The interest of Nitsche's methods in this context (see [START_REF] Chouly | An overview of recent results on Nitsche's method for contact problems[END_REF][START_REF] Annavarapu | A Nitsche stabilized finite element method for frictional sliding on embedded interfaces. Part I: Single interface[END_REF]) is that it combines the fact of being strongly consistent and the well-posed character of semi-discretization (see [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and time-marching schemes[END_REF][START_REF] Chouly | A Nitsche finite element method for dynamic contact: 2. Stability of the schemes and numerical experiments[END_REF]).In [START_REF] Chouly | Explicit Verlet time-integration for a Nitsche-based approximation of elastodynamic contact problems[END_REF] fully explicit schemes based on Verlet's scheme have been introduced, analysed and compared to other schemes previously introduced for impact dynamics. Of course, the disadvantage of schemes based on an explicit time integration is their conditional stability, which makes it necessary to consider a time step that can be extremely small.

In this context, using an implicit-explicit (IMEX) scheme can be advantageous. Indeed, following D.J. Eyre approach [START_REF] David | Unconditionally gradient stable time marching the cahn-hilliard equation[END_REF], by breaking down the operators in the difference of two monotonous parts, it is possible to build unconditionally stable schemes, i.e. without constraint on the time step size, with the cost of a single linear system resolution per time step. This cost is not so far from an explicit time integration scheme when the mass matrix is not lumped. These schemes therefore present an interesting compromise by avoiding the resolution of a non-linear problem at each time step and having stability properties close to the implicit schemes. We refer to [START_REF] Rodolfo R Rosales | Unconditional stability for multistep imex schemes: Theory[END_REF] for a comprehensive theory of unconditional stability of such IMEX schemes.

It should also be remembered that certain techniques has been developed in the context of explicit schemes to allow the use of a larger time step than the critical one imposed by the CFL condition. The simplest technique, the mass scaling, consists in adding some mass to the structure in order to obtain a greater critical time step. However, for rapid transient dynamics, this additional inertia may fundamentally change the solution. An alternative, called selective mass scaling has been proposed in [START_REF] Olovsson | Selective mass scaling for thin walled structures modeled with tri-linear solid elements[END_REF][START_REF] Olovsson | Selective mass scaling for explicit finite element analyses[END_REF] and further developed in [START_REF] Cocchetti | Selective mass scaling and critical time-step estimate for explicit dynamics analyses with solid-shell elements[END_REF][START_REF] Tkachuk | Local and global strategies for optimal selective mass scaling[END_REF] for instance. It consists in perturbing the mass matrix using the stiffness matrix which has the advantage not to modify too much the lowest eigenmodes of the structure. We will show that this latter strategy is equivalent to the use of an IMEX scheme, and that it may shed light on the choice of the terms that can be used to perturb the mass matrix.

In the continuity of the work presented in [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and time-marching schemes[END_REF][START_REF] Chouly | A Nitsche finite element method for dynamic contact: 2. Stability of the schemes and numerical experiments[END_REF][START_REF] Chouly | Explicit Verlet time-integration for a Nitsche-based approximation of elastodynamic contact problems[END_REF] respectively for implicit and explicit schemes, we present in this paper some IMEX schemes adapted to the contact problem approximated by Nitsche's method. The rest of the paper is described as follows. Section 2 is dedicated to the description of the semi-discrete formulation. Then, our IMEX schemes are proposed and their energy conservation properties are analyzed in Section 3. Some numerical tests are presented in Section 4 that confirm the presented properties but reveal a slowing down of the motion for large time steps. The section 5 provides an analysis of this phenomenon and interprets IMEX schemes in terms of selective mass scaling. We then propose different techniques to improve the approximation for large time step and end the document with a conclusion.

2 Problem setting and Nitsche's formulation

Problem setting

Let Ω in R d with d = 1, 2, 3 be the reference configuration of a linearly elastic body (with plain strain assumption for d = 2). Let us describe its dynamic evolution submitted to a contact condition with a rigid obstacle. We suppose that ∂Ω consists in three non-overlapping parts Γ D , Γ N and the contact boundary Γ C , with meas(Γ D ) > 0 and meas(Γ C ) > 0. 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 . The body is in potential contact on Γ C with a rigid foundation. We assume for simplicity a vanishing gap in the reference configuration. Considering T > 0 the final time, the evolution of the displacement field u : [0, T ] × Ω → R d satisfies the equations and conditions (1)-( 2):

ρü -div σ(u) = f , σ(u) = A ε(u) in (0, T ] × Ω, u = 0 on (0, T ] × Γ D , σ(u)n = g on (0, T ] × Γ N , u(0, •) = u 0 u(0, •) = u0 in Ω, (1)
where the following notations have been used: the time derivative of a quantity x is denoted ẋ, ρ is the density which is assumed to be constant for simplicity, u 0 and u0 are initial displacement and velocity, σ = (σ ij ), 1 ≤ i, j ≤ d, is the Cauchy stress tensor field, div denotes the divergence operator of tensor valued functions, ε(v) = (∇v + ∇v T )/2 represents the linearized strain tensor field, A is the fourth order symmetric elasticity tensor and n is the outward normal unit vector on ∂Ω. We consider the following decomposition into normal and tangential components

v = v n n + v t and σ(v)n = σ n (v)n + σ t (v),
for any displacement field v and density of surface forces σ(v)n defined on ∂Ω. This allows to express the frictionless unilateral contact condition on Γ C as follows:

u n ≤ 0 σ n (u) ≤ 0 σ n (u) u n = 0 σ t (u) = 0. (2)
with u 0 satisfying the compatibility condition u 0n ≤ 0 on Γ C . We refer to [START_REF] Eck | Unilateral contact problems[END_REF] for the mathematical analysis of elastodynamic contact problems. Apart for the one-dimensional case, the well-posedness of Problem (1)-( 2) is still an open issue. A few existence results has been proposed for simplified models (scalar wave equations, thin structures, one-dimensional case) in [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: the vibrating string with a concave obstacle[END_REF][START_REF] Schatzman | Un problème hyperbolique du 2ème ordre avec contrainte unilatérale: la corde vibrante avec obstacle ponctuel[END_REF][START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF][START_REF] Uhn | A boundary thin obstacle problem for a wave equation[END_REF][START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF][START_REF] Ahn | Existence of solutions for a class of impact problems without viscosity[END_REF][START_REF] Pozzolini | Vibro-impact of a plate on rigid obstacles: existence theorem, convergence of a scheme and numerical simulations[END_REF].

We note the Hilbert space

V := v ∈ H 1 (Ω) d : v = 0 on Γ D ,
where H s (D), s ∈ R stands for the classical Sobolev space (see [START_REF] Adams | Sobolev spaces[END_REF]) on the domain D. The usual scalar product of V will be denoted (•, •) s,D and the corresponding norm • s,D . We consider the following forms, for any u and v in V, for all t ∈ [0, T ]:

a(u, v) := Ω σ(u) : ε(v) dΩ, L(t)(v) := Ω f (t) • v dΩ + Γ N g(t) • v dΓ. (3)
Note that, introducing the mechanical energy

E(t) := 1 2 ρ u(t) 2 0,Ω + 1 2 a(u(t), u(t)), ∀t ∈ [0, T ],
and assuming that the contact force does not dissipate any energy (i.e. that the so-called persistency condition σ n (u(t)) un (t) = 0 is satisfied, which is expected from a mechanical viewpoint but difficult to prove mathematically, see, e.g., [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF][START_REF] Armero | Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF]) then the solution to the dynamic contact problem (1)-( 2) is such that

d dt E(t) = L(t)( u(t)). (4)

A finite element Nitsche approach

Let us now present the Nitsche-based finite element semi-discretization of the dynamic contact problem (1)-( 2) which was introduced in [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and time-marching schemes[END_REF][START_REF] Chouly | A Nitsche finite element method for dynamic contact: 2. Stability of the schemes and numerical experiments[END_REF] together with some basic properties of wellposedness and energy conservation. Let V h ⊂ V be a family of finite dimensional vector spaces (see [START_REF] Ciarlet | The finite element method for elliptic problems, volume II of Handbook of Numerical Analysis[END_REF]) indexed by h coming from a finite element method on a family T h of triangulations, supposed regular in Ciarlet's sense, of the domain Ω (h = max K∈T h h K where h K is the diameter of the element K). For instance, for a standard Lagrange finite element method of degree k > 0, we have

V h := v h ∈ (C 0 (Ω)) d : v h |K ∈ (P k (K)) d , ∀K ∈ T h , v h = 0 on Γ D .
However, any C 0 -conforming finite element method would be convenient.

With the use of a piecewise constant parameter γ h > 0 defined on the contact boundary Γ C satisfying for every K ∈ T h having a face on Γ C

γ h | K∩Γ C = γ 0 h K , (5)
where γ 0 is a positive given constant (the so-called Nitsche's parameter), we use the following equivalent reformulation of the contact condition (2) (see [START_REF] Alart | A generalized Newton method for contact problems with friction[END_REF][START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF]):

σ n (u) = -(σ n (u) -γ h u n ) -= -(γ h u n -σ n (u)) + (6)
where (•) + , the positive part is defined by (x) + := (x + |x|)/2 and (•) -, the negative part by (x) -:= (-x) + .

As in [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and time-marching schemes[END_REF][START_REF] Chouly | A Nitsche finite element method for dynamic contact: 2. Stability of the schemes and numerical experiments[END_REF] we consider a family of methods indexed by an additional parameter Θ ∈ R (generally, Θ = -1, 0, 1, see, e.g., [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF]) which leads to the following expression of the space semi-discretized elastodynamic contact problem (see, e.g, [START_REF] Chouly | An overview of recent results on Nitsche's method for contact problems[END_REF][START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and time-marching schemes[END_REF]):

                   Find u h : [0, T ] → V h such that for t ∈ [0, T ] : (ρü h (t), v h ) 0,Ω + a(u h (t), v h ) - Γ C Θ γ h σ n (u h )σ n (v h ) dΓ - Γ C 1 γ h σ n (u h ) -γ h u h n (t) - Θσ n (v h ) -γ h v h n dΓ = L(t)(v h ), ∀ v h ∈ V h , u h (0, •) = u h 0 , uh (0, •) = uh 0 , (7)
where u h 0 (resp. uh 0 ) is an approximation in V h of the initial displacement u 0 (resp. the initial velocity u0 ).

Following [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and time-marching schemes[END_REF], we consider the following mesh-dependent norms for any v ∈ L 2 (Γ C ):

v -1 2 ,h,Γ C := (h K ) 1 2 v 0,Γ C , v 1 2 ,h,Γ C := (h K ) -1 2 v 0,Γ C , the scalar product for all v h , w h ∈ V h : (v h , w h ) γ h := (v h , w h ) 1,Ω + (γ h 1 2 v h n , γ h 1 2 w h n ) 0,Γ C ,
and

• γ h := (•, •) 1 2
γ h the corresponding norm. Then we reformulate [START_REF] Boiveau | A penalty-free Nitsche method for the weak imposition of boundary conditions in compressible and incompressible elasticity[END_REF] as a system of non-linear second-order differential equations, using Riesz's representation theorem in (V h , (•, •) γ h ). Let M h : V h → V h , be the mass operator defined for all v h , w h ∈ V h by (M h v h , w h ) γ h := (ρv h , w h ) 0,Ω , and B h : V h → V h , for all v h , w h ∈ V h the non-linear operator defined by

(B h v h , w h ) γ h := a(v h , w h ) - Γ C Θ γ h σ n (v h )σ n (w h ) dΓ - Γ C 1 γ h σ n (v h ) -γ h v h n - Θσ n (w h ) -γ h w h n dΓ. (8)
Finally, we denote by L h (t) the vector in V h such that, for all t ∈ [0, T ] and for every w h in V h : (L h (t), w h ) γ h := L(t)(w h ). Problem [START_REF] Boiveau | A penalty-free Nitsche method for the weak imposition of boundary conditions in compressible and incompressible elasticity[END_REF] then reads:

     Find u h : [0, T ] → V h such that for t ∈ [0, T ] : M h üh (t) + B h u h (t) = L h (t), u h (0, •) = u h 0 , uh (0, •) = uh 0 . (9)
The following theorem together with the boundedness of (M h ) -1 γ h (see [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and time-marching schemes[END_REF]) show that Problem (7) (or equivalently Problem ( 9)) is well-posed.

Theorem 1. 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 all

v h 1 , v h 2 ∈ V h : B h v h 1 -B h v h 2 γ h ≤ C(1 + γ -1 0 )(1 + |Θ|) v h 1 -v h 2 γ h . (10)
As a consequence, for every value of Θ ∈ R and γ 0 > 0, Problem (7) admits one unique solution

u h ∈ C 2 ([0, T ], V h ).
Concerning the energy evolution, and considering the discrete energy as follows:

E h (t) := 1 2 ρ uh (t) 2 0,Ω + 1 2 a(u h (t), u h (t)), ∀t ∈ [0, T ].
associated to the solution u h (t) to Problem [START_REF] Boiveau | A penalty-free Nitsche method for the weak imposition of boundary conditions in compressible and incompressible elasticity[END_REF]. we define also, as in [START_REF] Chouly | Explicit Verlet time-integration for a Nitsche-based approximation of elastodynamic contact problems[END_REF], the modified energy more suited to Nitsche's method

E h 1 (t) := E h (t) - 1 2γ 0   σ n (u h (t)) 2 -1 2 ,h,Γ C -σ n (u h (t)) -γ h u h n (t) - 2 -1 2 ,h,Γ C   . (11)
The two following results are stated in [START_REF] Chouly | Explicit Verlet time-integration for a Nitsche-based approximation of elastodynamic contact problems[END_REF]:

Proposition 2. For γ 0 large enough, there exists C > 0 independent of h, of γ 0 and of the solution to Problem [START_REF] Boiveau | A penalty-free Nitsche method for the weak imposition of boundary conditions in compressible and incompressible elasticity[END_REF], such that, for all t ∈ [0, T ]:

E h (t) ≤ CE h 1 (t).
Theorem 3. 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 (7) then satisfies the following identity:

d dt E h 1 (t) = -(1 -Θ) Γ C 1 γ h σ n (u h (t)) -γ h u h n (t) - + σ n (u h (t)) σ n ( uh (t)) dΓ.
Remark 4. As a result, 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 ( 6) by u h .

IMEX schemes

As mentioned in the introduction, most of the difficulties come from the treatment of the nonlinear contact term. Indeed, It raises some instability issues when using an explicit treatment while an implicit one leads to the resolution of a non-linear system at each time step. The main idea of the IMEX scheme is then to split the non-linear term B h as

B h = K h -A h ,

using

• an implicit integration of K h as a linear symmetric positive operator

• an explicit treatment of A h which contains all the non-linear contributions of B h .

In our context, a first idea is to treat implicitly the stiffness operator and explicitly the Nitsche contact term. It leads to define the following splitting:

(K h -1 v h , w h ) γ h := a(v h , w h ), (A h -1 v h , w h ) γ h := Γ C 1 γ h σ n (v h ) -γ h v h n - Θσ n (w h ) -γ h w h n dΓ + Γ C Θ γ h σ n (v h )σ n (w h ) dΓ.
However, the stability of the IMEX scheme requires that A h is a monotonous operator which is clearly not the case using this splitting. More precisely, this means that for all

v h , w h ∈ V h , v h = 0, (K h v h , v h ) γ h > 0, (A h v h -A h w h , v h -w h ) γ h ≥ 0.
Then, a more adequate choice of splitting decomposition for the dynamic contact problem (9) verifying these properties for γ 0 large enough and Θ ∈ [-1, 1] is the one defined by

(K h 0 v h , w h ) γ h := a(v h , w h ) - Γ C Θ γ h σ n (v h )σ n (w h ) dΓ + Γ C 1 γ h σ n (v h ) -γ h v h n σ n (w h ) -γ h w h n dΓ +(1 -Θ) Γ C 1 γ h σ n (v h )σ n (w h ) + γ h v h n w h n dΓ, (12)
(A h 0 v h , w h ) γ h := Γ C 1 γ h σ n (v h ) -γ h v h n + σ n (w h ) -γ h w h n dΓ -(1 -Θ) Γ C 1 γ h σ n (v h ) -γ h v h n - σ n (w h )dΓ +(1 -Θ) Γ C 1 γ h σ n (v h )σ n (w h ) + γ h v h n w h n dΓ, (13)
where the last term of the two operators has been added for Θ = 1 to ensure the monotonicity of A h 0 .

Remark 5. Such a splitting decomposition is not unique and the basic idea is of course to consider the simplest non-linear operator A h satisfying the monotonous assumption.

Proposition 6. The operator

A h 0 is a monotonous one for Θ ∈ [-1, 1].
Proof. Due to the monotonicity of the positive part, the first term of ( 13) is monotonous and then

(A h 0 v h -A h 0 w h , v h -w h ) γ h ≥ -(1 -Θ) Γ C 1 γ h -σ n (v h ) + γ h v h n + --σ n (w h ) + γ h w h n + σ n (v h -w h )dΓ +(1 -Θ) Γ C 1 γ h (σ n (v h -w h )) 2 + γ h (v h n -w h n ) 2 dΓ.
However, using again the monotonicity of the positive part and its Lipschitz-continuity, one obtains

--σ n (v h ) + γ h v h n + --σ n (w h ) + γ h w h n + σ n (v h -w h ) ≥ - 1 2 -σ n (v h ) + γ h v h n + --σ n (w h ) + γ h w h n + σ n (v h -w h ) - 1 2 -σ n (v h ) + γ h v h n + --σ n (w h ) + γ h w h n + γ h (v h n -w h n ) ≥ - 1 2 |σ n (v h -w h ) -γ h (v h n -w h n )| |σ n (v h -w h ) + γ h (v h n -w h n )| ≥ - 1 2 |(σ n (v h -w h )) 2 -γ h 2 (v h n -w h n ) 2 | ≥ - 1 2 (σ n (v h -w h )) 2 - γ h 2 2 (v h n -w h n ) 2 ,
which allows to conclude.

Remark 7. Note that operator K h 0 derives of course from the potential

ψ K h 0 (v h ) := 1 2 (K h 0 v h , v h ) γ h
and, in the case Θ = 1, the operator A h 0 also derives from the convex potential

ψ A h 0 (v h ) := Γ C 1 2γ h σ n (v h ) -γ h v h n 2 + dΓ.

A semi-implicit β-Newmark scheme

Let us now consider a uniform discretization of the time interval [0, T ]: (t 0 , . . . , t N ), with t n = nτ , n = 0, . . . , N where τ = T /N is the time step. In the following, we will use the notation

x h,n+α := (1 -α)x h,n + αx h,n+1 , x h,n-α := (1 -α)x h,n + αx h,n-1 for α > 0 and arbitrary quantities x h,n-1 , x h,n , x h,n+1 ∈ V h
. Moreover, we denote by u h,n (resp. uh,n and üh,n ) the discretized displacement (resp. velocity and acceleration) at time step t n . We present a family of schemes indexed by two parameters β ∈ [ 1 4 , 1 2 ] and α ∈ [0, 1 2 ]. We propose first the following two-step scheme:

       Find u h,n+1 ∈ V h such that: M h u h,n+1 -2u h,n + u h,n-1 τ 2 + K h βu h,n+1 + (1 -2β)u h,n + βu h,n-1 -A h u h,n-α = L h,n , (14)
with initial conditions u h,0 = u h 0 and u h,1 = u h 1 . Note that only a linear system corresponding to M h + βτ 2 K h is to be solved at each time step. This scheme is based on a β-Newmark scheme for the linear part K h and when α = 0, a central difference scheme (Störmer-Verlet scheme) for the non-linear part A h (a similar choice is presented in [START_REF] Stern | Implicit-explicit variational integration of highly oscillatory problems[END_REF], for instance). This scheme is second order accurate for α = 0 and only first order accurate otherwise.

From a numerical viewpoint, it is also useful to introduce the following two equivalent schemes:

• The one-step Newmark-like scheme:

                   Find u h,n+1 , uh,n+1 ∈ V h such that: M h u h,n+1 = M h u h,n + τ M h uh,n - τ 2 2 K h 2βu h,n+1 + (1 -2β)u h,n + τ 2 2 L h,n + A h u h,n , M h uh,n+1 = M h uh,n + τ 2 L h,n+1 + L h,n -K h (u h,n+1 + u h,n ) + τ 2 2A h u h,n+1-α -A h u h,n+1 + A h u h,n , (15)
associated with the initial conditions

u h,0 = u h 0 , v h,0 = v h 0 , • The Leapfrog-like scheme          Find u h,n+1 , uh,n+ 1 2 ∈ V h such that: M h uh,n+ 1 2 = M h uh,n-1 2 + βτ 2 K h uh,n-1 2 -uh,n+ 1 2 + τ L h,n -K h u h,n + A h u h,n-α , u h,n+1 = u h,n + τ uh,n+ 1 2 , ( 16
) associated with the initial conditions u h,0 = u h 0 , v h, 1 2 = v h 1 2
.

Note that the latter scheme is a one-step scheme only for α = 0, however, even for α = 0 starting from u h,0 and v h, 1 2 , the value of u h,1 is given by the second relation of ( 16) then the first relation gives v h, 3 2 from u h,1 , u h,0 and v h, 1 2 .

Proposition 8. Apart from the initial conditions, the three family of schemes ( 14), ( 15) and (16) are equivalent.

Proof. Subtracting the first relation of ( 15) with itself replacing n by n -1 and using the second relation of ( 15) to eliminate the velocity, one obtains directly [START_REF] Ciarlet | The finite element method for elliptic problems, volume II of Handbook of Numerical Analysis[END_REF]. In a similar manner, subtracting the second relation of ( 16) with itself replacing n by n -1, multiplying the result by M h + βτ 2 K h and using it in the first relation of ( 16) leads also to [START_REF] Ciarlet | The finite element method for elliptic problems, volume II of Handbook of Numerical Analysis[END_REF].

Discrete Energy evolution

We make the additional assumption that the non-linear operator A h derives from a positive convex potential ψ A h (u h ) in the sense

(A h u h , v h ) γ h = Dψ A h (u h )[v h ],
where

Dψ A h (u h )[v h ] := lim ε→0 ψ A h (u h + ε v h ) -ψ A h (u h ) ε is the directional derivative of ψ A h at u h in the direction v h . This means in particular that ψ B h (u h ) = 1 2 (K h u h , u h ) γ h -ψ A h (u h
) is the potential of B h and we assume also ψ B h (u h ) positive and convex.

For the sake of simplicity, we consider the case of a constant source term L h,n = L h and we integrate the source term potential to the energy. We consider the energy at mid-point

E h,n+ 1 2 := 1 2 (M h u h,n+1 -u h,n τ , u h,n+1 -u h,n τ ) γ h + ψ B h (u h,n+ 1 2 ) -(L h , u h,n+ 1 2 ) γ h . (17)
A more adapted energy for the β-Newmark scheme can then be deduced as

E h,n+ 1 2 imex1 := E h,n+ 1 2 + 1 2 β - 1 4 K h (u h,n+1 -u h,n ), u h,n+1 -u h,n γ h , (18)
which is ensured to be positive for β ≥ 1 4 and can be rewritten equivalently

E h,n+ 1 2 imex1 = 1 2 (M h u h,n+1 -u h,n τ , u h,n+1 -u h,n τ ) γ h + (1 -2β)(K h u h,n+ 1 2 , u h,n+ 1 2 ) γ h + β - 1 4 K h u h,n+1 , u h,n+1 γ h + β - 1 4 K h u h,n , u h,n γ h ( 19
)
-ψ A h (u h,n+ 1 2 ) -(L h , u h,n+ 1 2 ) γ h .
Proposition 9. When ψ A h and ψ B h are both convex potentials, for a constant source term

L h,n = L h and in the case α = 1 2 and 0 ≤ β ≤ 1 2 the energy E h,n+ 1 2
imex1 decreases with respect to n.

Proof. The evolution of this discrete energy can be computed for schemes ( 14), ( 15) and ( 16) taking the scalar product of the relation ( 14) with

u h,n+ 1 2 -u h,n-1 2 = u h,n+1 -u h,n-1 2 = u h,n+1 -u h,n 2 + u h,n -u h,n-1 2 ,
and using the decomposition

βu h,n+1 +(1-2β)u h,n +βu h,n-1 = 1 2 -β (u h,n+1 +2u h,n +u h,n-1 )+ 2β - 1 2 (u h,n+1 -u h,n-1 ).
This leads to

1 2 M h u h,n+1 -u h,n τ , u h,n+1 -u h,n τ γ h - 1 2 M h u h,n -u h,n-1 τ , u h,n -u h,n-1 τ γ h + (1 -2β) K h u h,n+ 1 2 , u h,n+ 1 2 γ h -(1 -2β) K h u h,n-1 2 , u h,n-1 2 γ h + β - 1 4 K h u h,n+1 , u h,n+1 γ h -β - 1 4 K h u h,n-1 , u h,n-1 γ h -(A h u h,n-α , u h,n+ 1 2 -u h,n-1 2 ) γ h = (L h , u h,n+ 1 2 -u h,n-1 2 ) γ h ,
which can be rewritten thanks to (19)

E h,n+ 1 2 imex1 = E h,n-1 2 imex1 -ψ A h (u h,n+ 1 2 ) + ψ A h (u h,n-1 2 ) + (A h u h,n-α , u h,n+ 1 2 -u h,n-1 2 ) γ h . (20)
The results follows directly from the convexity of ψ A h , since it implies

ψ A h (x) ≥ ψ A h (u h ) + A h v h , u h -v h γ h for all u h , v h ∈ V h .
This ensure the unconditional stability of the scheme α = 1 2 for 1 4 ≤ β ≤ 1 2 . However, α = 1 2 corresponds to first order dissipative schemes where the schemes with α = 0 are second order and reversible in time schemes. The reversibility in time of the schemes for α = 0 make it impossible to have a strict dissipativity property. Proving a strict conservation of energy, when this is possible, means to find the particular conserved discrete energy. Even in the very simple case where ψ A h is a quadratic potential (i.e. A h a linear symmetric monotonous operator), the energy E h,n+ 1 2 is not strictly conserved. In order to obtain a strict conservation in that case, it is necessary to consider the slightly modified energy

E h,n+ 1 2 imex2 = E h,n+ 1 2 imex1 + 1 8 A h u h,n+1 -A h u h,n , u h,n+1 -u h,n γ h . ( 21
)
Note that the additional term is a non-negative one whenever A h is monotonous. The following result can be deduced from basic simplifications in equation [START_REF] Glasner | Improving the accuracy of convexity splitting methods for gradient flow equations[END_REF] with the additional term considered in [START_REF] Hager | A stable energy-conserving approach for frictional contact problems based on quadrature formulas[END_REF].

Proposition 10. When ψ A h and ψ B h are both quadratic convex potentials, for a constant source term L h,n = L h and in the case α = 0 and 0

≤ β ≤ 1 2 the energy E h,n+ 1 2
imex2 is constant with respect to n.

This leads also to the unconditional stability of the corresponding scheme α = 0 for

1 4 ≤ β ≤ 1 2
in the case of quadratic potentials.

First numerical comparison

We now perform a numerical comparison of the IMEX schemes described in section 3 with some standard implicit and explicit schemes. The comparison is done first on a 1D case where an exact solution is available and then in a 2D Hertz-like case with the help of our freely available finite element library GetFEM++ (see [START_REF] Renard | Automated fe modeling of multiphysics problems based on a generic weak form language[END_REF] and http://getfem.org). The explicit scheme is Verlet's one which corresponds to scheme [START_REF] Ciarlet | The finite element method for elliptic problems, volume II of Handbook of Numerical Analysis[END_REF] for A h = -B h and K h = 0 (see also [START_REF] Chouly | Explicit Verlet time-integration for a Nitsche-based approximation of elastodynamic contact problems[END_REF]) and the implicit scheme is Crank-Nicolson one (also called trapezoidal rule) which corresponds conversely to A h = 0 and K h = B h when β = 1 4 . The Nitsche variant used is the energy conserving one Θ = 1 which gave the most satisfactory results for the approximation of elastodynamic contact problems in [START_REF] Chouly | Explicit Verlet time-integration for a Nitsche-based approximation of elastodynamic contact problems[END_REF].

1D numerical experiments: multiple impacts of an elastic bar

We first consider the one-dimensional case d = 1 described in [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF] corresponding to a single contact point. It consists in an elastic bar Ω = (0, L) with Γ C = {0}, Γ D = {L} and Γ N = ∅. The elastodynamic equation is then reduced to find u :

(0, T ] × (0, L) → R such that ρü -E ∂ 2 u ∂x 2 = f, in (0, T ] × (0, L), ( 22
)
where E is the Young modulus and the Cauchy stress tensor is given by σ(u) = E ∂u ∂x . We consider a finite element space using linear finite elements and a uniform subdivision of [0, L].

Let us denote U n := [U n 0 , . . . , U n N ] T (resp. Un , Ün ) the vector of all the nodal values of u h,n (resp. uh,n and üh,n ). The component of index 0 corresponds to the node at the contact point Γ C .

We take the following values for the parameters:

f = 0, E = 1, ρ = 1, L = 1, u 0 (x) = 1 2 -x x = 0 L t = 0 t 2 = 2 t 3 = 3 t 1 = 1
Figure 1: Multiple impacts of an elastic bar. The bar is clamped at x = L and the contact node is located at the bottom. The displacement is periodic of period 3, with one impact during each period (here between t = 1 and t = 2). impacts first the rigid ground, located at x = 0, and then gets compressed again (see Figure 1). This problem admits a closed-form solution u which derivation and expression is detailed in [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF]. 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. The chosen simulation time is T = 12, so that we can observe four successive impacts.

We discretize the bar with 50 linear finite elements (h = 0.02) and take τ = 0.01 (ν C = 0.5) and τ = 0.05 (ν C = 2.5), for

ν C := c 0 τ h = E ρ τ
h the Courant number and c 0 the wave speed. The numerical tests are presented for a time step smaller than the critical one on Figure 2 and a time step larger than the critical one on Figure 3. For these two figures, the left plots correspond to the displacement on the contact point (u h,n (0) = U n 0 ). The dotted red curve is the exact solution. The right plots correspond to the evolution of the discrete energy. The plotted energy is E h,n+ for IMEX α = 0. In addition to the plots for the explicit Verlet scheme and for the implicit Crank-Nicolson one, two versions of the IMEX scheme are considered: the version α = 1 2 which has been proved to be unconditionally stable and the version α = 0 which is second order an reversible in time. For both cases, we consider β = 1 4 and Nitsche's parameter is set to γ 0 = 5 for all experiments.

We see in Figure 2 that for a time step slightly lower and close to the critical time step, Verlet's, Crank-Nicolson and IMEX scheme for α = 0 give some accurate and similar results. This is not the case for the IMEX scheme for α = 1 2 . In addition to be first order while the other schemes are second order, it is exaggeratedly dissipative. Despite it guaranteed unconditional stability property, it seems to be unfeasible in that context since it would necessitate a too much small time step to recover an acceptable dissipation. This is confirmed on Figure 3 where for a time step ten times larger (and larger than the critical time step), the dissipation is so important that the motion do not look like the exact solution at all. Since one of the goal of the IMEX scheme is to be able to treat large time steps, this completely disqualifies this version of the IMEX scheme. On this figure, the fact that the time step is larger than the critical one make that Verlet's explicit scheme is no longer stable, which was the expected behavior. The implicit Crank-Nicolson scheme still gives a reasonable approximation of the motion and still has a good energy conserving property. The IMEX scheme for α = 0 gives a reasonable approximation on the two first periods with a good energy conservation. However, there is then a degradation of the solution, accompanied by a significant lengthening of the period of the solution.

2D numerical experiments: multiple impacts of a disc

Numerical experiments are then carried out in 2D on the impact of a disc on a rigid support at y = 0. The physical parameters are the following: the diameter of the disc is D = 40, the Lamé coefficients are λ = 3 × 10 4 and µ = 3 × 10 4 , and the material density is ρ = 1. The total simulation time is T = 30. The volume load in the vertical direction is set to f = 0.1 (gravity, oriented towards the support). The upper part of the boundary is a traction free boundary and the lower part is the contact zone Γ C . We consider an initial vertical displacement (u 0 = (0, 2)) and no initial velocity ( u0 = 0). We use a P 2 isoparametric finite element method, whose mesh is represented on Figure 4.

A first comparison is presented on Figure 5 for a time step τ = 0.0025 chosen smaller but close to the critical time step for the explicit scheme. For each scheme, the three plots correspond to the evolution of the normal displacement at the lowest point of the disc (the first point which enter into contact with the rigid support), the contact stress at this lowest point and the evolution of the discrete energy. The plotted energy is still E h,n+ 1 2 for Verlet and Crank-Nicolson schemes and E h,n+ 1 2 imex2 for IMEX α = 0. We can see a good accordance between the three schemes : Verlet's, Crank-Nicolson and IMEX α = 0 schemes. We no longer compare with the IMEX scheme with α = 1 2 which gives a too poor approximation. There is some overshoots of the discrete energy at impact when using the IMEX scheme, probably due to the non-regularity of the contact terms.

A comparison for a time step ten times larger is presented on Figure 6. The result for Verlet's scheme is not presented since the time step is larger than the critical value for stability. The implicit Crank-Nicolson still gives a good approximation of the solution. The IMEX scheme for α = 0 is stable and do not present some spurious oscillations, however, as in 1D, there is a significant increase of the time between two impacts.

At this point, we can conclude that the built IMEX schemes have very good stability properties for the contact problem. However, there is a slowing down of motion for large time steps that is not present with the Crank-Nicolson scheme. The next section is devoted to give an interpretation of this phenomenon.

Selective mass scaling effects, difficulty and improvements

In this section, we will try to understand the reason for the alteration of the observed motion when the IMEX scheme is used for large time steps and the link with the existing works on the selective mass scaling. We also propose some variants to try to reduce theses alterations. We remark first that the IMEX scheme in the central difference form ( 14) for α = 0 can easily be rewritten

       Find u h,n+1 ∈ V h such that: M h + βτ 2 K h u h,n+1 -2u h,n + u h,n-1 τ 2 + B h u h,n = L h,n , (23)
and the discrete energy ( 18) can be rewritten in a form more related to ( 23):

E h,n+ 1 2 imex1 = 1 2 M h + β - 1 4 τ 2 K h u h,n+1 -u h,n τ , u h,n+1 -u h,n τ γ h + ψ B h (u h,n+ 1 2 ).
It is particularly noteworthy that [START_REF] Heinstein | Contact-impact modeling in explicit transient dynamics[END_REF] corresponds exactly to a central difference scheme (Verlet's scheme) applied to the original problem but with a mass operator M h replaced by

M h + βτ 2 K h .
A perturbation of the mass matrix in order to increase the critical time step in an explicit time integration scheme is generally called a mass scaling. The perturbation of the mass matrix with a term proportional to the stiffness matrix has been proposed in [START_REF] Olovsson | Selective mass scaling for thin walled structures modeled with tri-linear solid elements[END_REF] and further studied in [START_REF] Olovsson | Selective mass scaling for explicit finite element analyses[END_REF][START_REF] Cocchetti | Selective mass scaling and critical time-step estimate for explicit dynamics analyses with solid-shell elements[END_REF][START_REF] Tkachuk | Local and global strategies for optimal selective mass scaling[END_REF] under the name "selective mass scaling". As far as we know, the interpretation in term of an implicit-explicit scheme has not been given yet, although this seems to be the most straightforward one.

The numerical study in [START_REF] Olovsson | Selective mass scaling for explicit finite element analyses[END_REF] (Figure 1 in this reference) reveals a thresholding of the highest eigenfrequencies and a moderate modification of the lowest eigenfrequencies, the rigid modes not being impacted. Our stability result imply that the thresholding of the highest eigenfrequencies is obtained not only for the stiffness matrix of the problem, but for any symmetric matrix K h such that both K h and A h = K h -B h are monotonous operators.

However, in our case, even if we recover the increase of the critical time step (even finding unconditionally stable methods), there is a consequent alteration of the rigid modes. These rigid lowest point modes are however not impacted by a standard stiffness matrix, of course, since the rigid modes lie in its kernel. Indeed, in our situation, the operator

(K h 0 v h , w h ) γ h = a(v h , w h ) - Γ C 1 γ h σ n (v h )σ n (w h ) dΓ + Γ C 1 γ h σ n (v h ) -γ h v h n σ n (w h ) -γ h w h n dΓ (24)
contains, additionally to the stiffness terms, the Nistche's contact terms in which we can extract the penalty term

Γ C γ h v h n w h n dΓ. ( 25
)
This term is in fact the only one being non-zero for a rigid motion. This is a term on the contact boundary which adds a significant mass on the boundary in the modified mass matrix M h + βτ 2 K h . This term is necessary to stabilize the explicit treatment of the non-linear Nitsche contact term but causes rigid body motion alterations even in the absence of contact.

Remark 11. About the implicit Crank-Nicolson scheme, it can also be rewritten in the form close to [START_REF] Heinstein | Contact-impact modeling in explicit transient dynamics[END_REF] for

K h = B h and β = 1 4 replacing K h (u h,n+1 -2u h,n + u h,n-1 ) by B h (u h,n+1 ) - 2B h (u h,n ) + B h (u h,n-1
). The reason why the rigid body motions are not perturbed using the Crank-Nicolson scheme is simply that the additional penalty term in B h is only present when contact occurs which has a limited influence since in that case the normal displacement is close to zero.

Unfortunately, keeping a constant K h , it is not possible to add the penalty term in K h only in the occurrence of contact and keeping a monotonous operator A h and so it is not possible to obtain the desired increase of the critical time step. We will however see in the next sections some possible ways to overcome this difficulty and improve the approximation for large time steps.

A non-constant splitting decomposition of B h

Since the perturbation of the additional penalty term in K h occurs mainly in absence of contact, a first idea is to add this term only when there is contact with a non-constant operator K h . This can be done, for instance, by considering the following splitting decomposition

(K h 1 (v h,0 )v h , w h ) γ h := a(v h , w h ) - Γ C 1 γ h σ n (v h )σ n (w h ) dΓ + Γ C H(γ h v h,0 n -σ n (v h,0 )) 1 γ h σ n (v h ) -γ h v h n σ n (w h ) -γ h w h n dΓ, (26)
and

(A h 1 v h , w h ) γ h := Γ C 1 γ h σ n (v h ) -γ h v h n + σ n (w h ) -γ h w h n dΓ - Γ C H(γ h v h,0 n -σ n (v h,0 )) 1 γ h σ n (v h ) -γ h v h n σ n (w h ) -γ h w h n dΓ. (27)
Here H(•) is the Heaviside function (H(x) = 1 for x ≥ 0, H(x) = 0 for x < 0) and v h,0 is the displacement at the previous time step. This is of course a clear drawback that the matrix of the implicit term is changing from an iteration to another. Moreover, the stability results of Section 3.2 no longer apply. Finally, this scheme can be viewed as a special implicit-explicit scheme where the contact status is taken at the previous time step. Some numerical tests in the two-dimensional case are shown on Figure 7. Compared to Figure 6 for our first IMEX scheme, it can be seen that a good approximation of the motion is recovered. The stability is preserved for the time step τ = 0.025 but this is not the case for a larger time step τ = 0.1. In conclusion, this first approach allows to overcome the problem of perturbation of the motion due to the additional contact term and allows the use of larger time steps than with explicit Verlet's scheme. However, the unconditional stability property is not preserved.

Small number of fixed points iterations

Another way to improve the accuracy of IMEX schemes, proposed for instance in [START_REF] Glasner | Improving the accuracy of convexity splitting methods for gradient flow equations[END_REF], is to add a small number of fixed points iterations taking the solution to the IMEX scheme as a predictor step. At time step n + 1, the quantities u h,n-1 and u h,n being known, we consider the following linear system which consists to find u h,n+1

solution to

M h + βτ 2 K h   u h,n+1 (i)
-2u h,n + u h,n-1

τ 2   + K h u h,n -G h (i) = L h,n . (28)
Notice that u h,n+1

, defined as the solution to [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF] for G h (0) = A h u h,n , is the solution to our IMEX scheme [START_REF] Ciarlet | The finite element method for elliptic problems, volume II of Handbook of Numerical Analysis[END_REF] (or ( 23)) for α = 0. Then, one obtains u h,n+1 (i+1) from u h,n+1

by the recurrence relation

u h,n+1 (
i+1) solution to [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF] for

G h (i+1) βA h u h,n-1 + (1 -2β)A h u h,n + βA h u h,n+1 (i)
.

In particular, it leads to the following fixed point iteration

u h,n+1 (i+1) = H n - M h βτ 2 + K h -1 A h u h,n+1 (i)
, which is a contraction for τ small enough. Here

H n = 2u h,n -u h,n-1 +τ 2 M h + βτ 2 K h -1 L h,n -K h u h,n -βA h u h,n-1 -(1 -2β)A h u h,n ,
and the iteration converges toward the solution of the implicit β-Newmark scheme.

The numerical results on the two-dimensional case of Section 4.2 is presented on Figure 8. Comparing with the results on Figure 6, one notices that even with only one fixed point iteration, there is an important correction of the period of the motion. However, the initial period is reached for ten fixed point iterations, which corresponds to a non-negligible computational cost. Compared for instance with a Newton method on the implicit scheme, note that the matrix of the linear system to be solved at each iteration do not vary, allowing for instance a unique factorization.

Compensate the penalty term in K h

Since the main difficulty comes from the penalty term

Γ C γ h v h n w h n dΓ
in the operator K h , another option is to try to substitute it by a stiffness term which do not perturbate the rigid body motion. When the bilinear form a(u, v) defined by ( 3) is coercive (i.e. in the present case when the Dirichlet boundary Γ D is of non-zero measure in ∂Ω), then, using additionally a trace inequality, there exists a constant α h > 0 such that

a(u, u) ≥ α h Γ C γ h u 2 n dx, ∀u ∈ V.
The fact that - Proposition 6, we can conclude, that in the coercive case, the following splitting corresponds to two monotonous operators for γ 0 large enough:

Γ C γ h u 2 n dx + 1 α h a(u,
(K h 2 v h , w h ) γ h := 1 + 2 α h a(v h , w h ) - Γ C 1 γ h σ n (v h )σ n (w h ) dΓ, (29)
(A h 2 v h , w h ) γ h := Γ C 1 γ h σ n (v h ) -γ h v h n - σ n (w h ) -γ h w h n dΓ + 2 α h a(v h , w h ). (30)
The numerical results are presented on Figure 9 for the one-dimensional case of Section 4.1 and on Figure 10 for the two-dimensional case of Section 4.2. The important difference between the two test-cases is that in the one-dimensional case the stiffness term is coercive due to the Dirichlet condition on the top of the bar while this is not the case in the two-dimensional case. Nevertheless, the numerical results are similar in both the two cases. For a time step approximately five time larger than the critical time step, a reasonable approximation is obtained without the augmentation of the period noticed in the initial IMEX scheme. However, some fluctuation of the energy are noted. For a larger time step, the motion is perturbed by this augmentation of the stiffness term (the approximation with Crank-Nicolson implicit scheme has been added for comparison). and on Figure 12 for a time step ten times larger make appear some instabilities for the IMEX scheme. We did not find any splitting that avoid these instabilities for large time steps. Note that the contact stress at lowest point in Figure 11 and 12 has smaller negative values during impact than in Figure 5 and 6 for instance. This means that the contact area is larger with this variant of Nitsche's method than with the classical one. This is due to the fact that a lower Nitsche parameter γ 0 as been considered. Unlike the classical Nitsche's method, the stability of this variant is obtained for small values of γ 0 .

Second order correction of the β-Newmark scheme

A last proposed improvement is to consider a perturbed version of β-Newmark scheme in order to minimize the influence of the penalty term in the implicit treatment of K h . Recall that the first splitting introduced here reads as

B h = K h -1 -A h -1
where only the stiffness part is treated implicitly. This then leads to the following β-Newmark scheme

u h,n+1 = 2u h,n -u h,n-1 + M h + βτ 2 K h -1 -1 (L h,n -B h u h,n ), (35)
which is unstable in practice. Indeed, the last part of the instabilities comes from the integration of the non-monotonous contact term A h -1 which corresponds to

M h + βτ 2 K h -1 -1 (A h -1 u h,n ).
We then proposed to use the splitting B h = K h 0 -A h 0 where A h 0 can be viewed as a monotonous relaxation of A h -1 . As expected, the new scheme is numerically stable but the presence of the penalty term

u h,n+1 = 2u h,n -u h,n-1 + M h + βτ 2 K h 0 -1 (L h,n -B h u h,n ), (36)
Γ C γ h v h n w h n dΓ in K h 0
gives some alterations of the rigid modes even when there is no contact. Notice that theses perturbations can be localized on both source and stiffness terms:

M h + βτ 2 K h 0 -1 (L h,n -K h -1 u h,n ).
Our idea here is then to combine the benefits of the two previous schemes by introducing a correction of the source term L h,n and the stiffness term K h -1 u h,n before applying the selective mass scaling effect M h + βτ 2 K h 0 -1 . We can then consider the following scheme

u h,n+1 = 2u h,n -u h,n-1 + M h + βτ 2 K h 0 -1 (C 1 τ L h,n -C 2 τ K h -1 u h,n -(B h -K h -1 )u h,n ).
where the consistence of the scheme requires that the two correction operators C 1 τ and C 2 τ satisfy the limits lim

τ →0 C 1 τ = lim τ →0 C 2 τ = I d .
In that case, it is not difficult to see that the scheme is still a β-Newmark one M h u h,n+1 -2u h,n + u h,n-1 τ 2 + K h βu h,n+1 + (1 -2β)u h,n + βu h,n-1 -Ãh u h,n-α = Lh,n (37) where the source and the explicit term are now given by

Ãh = A h 0 -(C 2 τ -I d )K h -1
and Lh,n = C 1 τ L h,n .

A second order correction can then be made by constructing two operators C 1 τ and C 2 τ satisfying

M h + βτ 2 K h 0 -1 [C 1 τ L h,n ] = (M h ) -1 L h,n + O(τ 2 ), (38)
and

M h + βτ 2 K h 0 -1 [C 2 τ K h -1 u h,n ] = M h + βτ 2 K h -1 -1 K h -1 u h,n + O(τ 2 ). ( 39
)
Here is an example of these operators

C 1 τ = M h + βτ 2 K h 0 (M h ) -1 and C 2 τ = M h + βτ 2 K h 0 M h + βτ 2 K h -1 -1 . (40)
However, recall that the stability of the associated β-Newmark scheme requires the monotonous property of the operator Ãh = A h 0 -(C 2 τ -I d )K h -1 which is satisfied only asymptotically when τ goes to zero.

Figure 13 shows the result of the simulation in the two-dimensional case for the correction [START_REF] Schatzman | Un problème hyperbolique du 2ème ordre avec contrainte unilatérale: la corde vibrante avec obstacle ponctuel[END_REF]. The numerical results shows a good correction of the slowing down noted in Figure 6 for IMEX scheme [START_REF] Pozzolini | Vibro-impact of a plate on rigid obstacles: existence theorem, convergence of a scheme and numerical simulations[END_REF]. The plotted energy is the discrete one [START_REF] Hager | A stable energy-conserving approach for frictional contact problems based on quadrature formulas[END_REF] adapted to the modified operators given in [START_REF] Renard | Automated fe modeling of multiphysics problems based on a generic weak form language[END_REF]. It is well conserved but its level is different to the uncorrected schemes.

Conclusion

The interest of our IMEX schemes is to allow an acceleration of the simulation compared to explicit schemes by using larger time steps and with a cost corresponding only to the resolution of a (constant) linear system at each time step. We exhibited some theoretical stability results in Section 3.2 and the numerical tests of Section 4 confirm the possibility of obtaining some unconditionally stable IMEX schemes.

The first order IMEX scheme (α = 1 2 ) for which we have a proof of unconditional stability is however too much dissipative in practice. The second order schemes do not raise any stability issues in our tests, despite we do not have such a proof of stability.

A less positive aspect is that our numerical tests have revealed a drawback of IMEX schemes in the treatment of the penalty contact term: a slowing down of the motion is observed for large time steps.

An analysis of this phenomenon in Section 5 shows that it is due to the addition of the contact penalty term to the mass matrix in the IMEX scheme (this also occurs for the implicit schemes but only when contact occurs).

In this section we also provide an interpretation of the IMEX scheme in terms of selective mass scaling, which has not been proposed before, as far as as we know. We then proposed five different techniques to try to overcome the difficulty of too large time steps.

In subsection 5.1, we first tested a non-constant implicit term where the idea is to apply the penalty term only when contact occurs. This effectively removes the slowing down of the motion, but a non-constant linear system must be solved at every time step, and the unconditional stability of the scheme is now loose. In subsection 5.2, we also tested to add a few fixed points iterations. It is then possible to recover a good approximation with this method for large time steps, but at the price of several (but constant) linear resolutions. In subsection 5.3, we tested to replace the penalty term with a stiffness one, exploiting the coercivity of the problem. This also partially overcome the difficulty but only for moderately large time steps. In subsection 5.4, we tested the penalty-free version of Nitsche's method proposed in [START_REF] Boiveau | A penalty-free Nitsche method for the weak imposition of boundary conditions in compressible and incompressible elasticity[END_REF]. However, we did not find any stable splitting using this approach. Finally, in subsection 5.5, we proposed a framework for a second-order correction of the β-Newmark scheme and we give an example leading to a second-order correction with a numerical example showing that adapted corrections can improve the approximation for large time steps.

Figure 2 :

 2 Figure 2: Simulation on the one-dimensional case for the time step τ = 0.01 (below CFL) and the other numerical parameters Θ = 1, γ 0 = 5, h = 0.02, β = 0.25 for the different schemes.

Figure 3 :

 3 Figure 3: Simulation on the one-dimensional case for the time step τ = 0.05 (above CFL) and the other numerical parameters Θ = 1, γ 0 = 15, h = 0.02, β = 0.25 for the different schemes.

Figure 4 :Figure 5 :

 45 Figure 4: P 2 mesh used for the disc.

Figure 6 :

 6 Figure 6: Simulation on the two-dimensional case for the time step τ = 0.025 (above CFL) and the other numerical parameters Θ = 1, γ 0 = 5 × 10 5 , β = 0.25 for the different schemes.

Figure 7 :

 7 Figure 7: IMEX scheme with non-constant implicit term. Simulation on the two-dimensional case for different time steps and Θ = 1, γ 0 = 5 × 10 5 , β = 0.25.

Figure 8 :

 8 Figure 8: IMEX scheme with fixed point iterations. Simulation on the two-dimensional case for the time step τ = 0.025 (above CFL) and the other numerical parameters Θ = 1, γ 0 = 5×10 5 , β = 0.25 for the different schemes.

Figure 11 :

 11 Figure 11: Penalty free approximation. Simulation on the two-dimensional case and a time step τ = 0.001 (below CFL) and the other numerical parameters Θ = 1, γ 0 = 5 × 10 4 , β = 0.25 for the different schemes.

Figure 12 :

 12 Figure 12: Penalty free approximation. Simulation on the two-dimensional case and a time step τ = 0.01 (above CFL) and the other numerical parameters Θ = 1, γ 0 = 5 × 10 4 , β = 0.25 for the different schemes.

Figure 13 :

 13 Figure 13: Second order correction. Simulation on the two-dimensional case and a time step τ = 0.025 (above CFL) and the other numerical parameters Θ = 1, γ 0 = 5 × 10 4 , β = 0.25.

and u0 (x) = 0. The bar is initially compressed. Then, it is released without initial velocity. It

Acknowledgements

We thank Lionel Morançay from ALTAIR Engineering for the constructive interactions and discussions on explicit dynamics for contact problems.

A penalty free variant of Nitsche's method

Another possibility to avoid the occurrence of the penalty term is to turn to the penalty-free variant of Nitsche's method. An alternative to the presented Nitsche-based approximation [START_REF] Boiveau | A penalty-free Nitsche method for the weak imposition of boundary conditions in compressible and incompressible elasticity[END_REF] which has been proposed in [START_REF] Boiveau | A penalty-free Nitsche method for the weak imposition of boundary conditions in compressible and incompressible elasticity[END_REF], is based on the following reformulation of the contact condition (2):

and can be written as follows:

There is indeed no penalty term of the form ( 25) in [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF]. Additionally, as it is noted in [START_REF] Seitz | Computational Methods for Thermo-Elasto-Plastic Contact[END_REF][START_REF] Seitz | Nitsche's method for finite deformation thermomechanical contact problems[END_REF], and since this formulation is equivalent to the following one still based on (6):

This formulation is not symmetric and consequently it does not derive from a potential. An example of splitting that both K h and A h are monotonous operators is the following (the monotonicity of A h 3 can be proved thanks to an adaptation of the proof of Proposition 6):

However, one can note in this splitting the appearance of a penalty term. We did not find a splitting verifying the monotonicity of the two operators that do not contain a penalty term. Moreover, the numerical results presented on Figure 11 for a time step lower than the critical one