Keywords: AMS Subject Classification : 74M15, 74M20, 74M10, 74B20, 74H15, 74S30, 49M15, 90C53, 70F40, 70-08, 70E55, 35Q70 Hyperelasticity, Unilateral Contact, Normal Compliance, Moreau-Yosida Regularization, Coulomb Friction, Dynamics, Semi-Smooth Newton method, Primal-Dual Active Set, Numerical simulations

The purpose of this work is to present an improved energy conservation method for hyperelastodynamic contact problems based on specific normal compliance conditions. In order to determine this Improved Normal Compliance (INC) law, we use a Moreau-Yosida α-regularization to approximate the unilateral contact law. Then, based on the work of Hauret-LeTallec [1], we propose in the discrete framework a specific approach allowing to respect the energy conservation of the system in adequacy with the continuous case. This strategy (INC) is characterized by a conserving behavior for frictionless impacts and admissible dissipation for friction phenomena while limiting penetration. Then, we detail the numerical treatment within the framework of the semi-smooth Newton method and primal-dual active set strategy for the normal compliance conditions with friction. We finally provide some numerical experiments to bring into light the energy conservation and the efficiency of the INC method by comparing with different classical methods from the literature throught representative contact problems.

Introduction

Dynamic problems involving contact between deformable bodies or between a body and a foundation still remain today an important subject of study for mathematical and numerical analysis. In the literature, many references offer different approaches for the usual contact conditions with friction [START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF][START_REF] Chabrand | Various numerical methods for solving unilateral contact problems with friction[END_REF][START_REF] Chouly | An overview of recent results on Nitsche's method for contact problems[END_REF][START_REF] Khenous | Hybrid discretization of the Signorini problem with Coulomb friction. theoretical aspects and comparison of some numerical solvers[END_REF][START_REF] Lebon | Contact problems with friction: models and simulations[END_REF][START_REF] Pietrzak | Large deformation frictional contact mechanics: continuum formulation and augmented lagrangian treatment[END_REF][START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF][START_REF] Raous | Numerical methods for solving unilateral contact problem with friction[END_REF][START_REF] Sofonea | Mathematical models in contact mechanics[END_REF][START_REF] Wriggers | Computational contact mechanics[END_REF]. As a matter of fact, the equations resulting from frictional contact problems are difficult to solve both mathematically and numerically due to their inherent non-linearity and the non-smoothness. In order to handle these issues, several methods have already been successfully tested. They rely on handling the original conditions by normal compliance methods [START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF][START_REF] Oden | Interior penalty methods for finite element approximations of the Signorini problem in elastostatics[END_REF][START_REF] Martins | Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws[END_REF] or by other relevant methods such as the quasi-augmented lagrangian method [START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF], the bi-potential method [START_REF] De Saxcé | New inequality and functional for contact with friction: the implicit standard material approach[END_REF][START_REF] Dumont | On enhanced descent algorithms for solving frictional multicontact problems: application to the discrete element method[END_REF], the conjugate gradient method [START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF][START_REF] Wriggers | Computational contact mechanics[END_REF], Uzawa method [START_REF] Raous | Numerical methods for solving unilateral contact problem with friction[END_REF][START_REF] Joli | Uzawa and newton algorithms to solve frictional contact problems within the bi-potential framework[END_REF] and Nitsche finite element method [START_REF] Chouly | An adaptation of Nitsche's method to the tresca friction problem[END_REF][START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. space semi-discretization and time-marching schemes[END_REF][START_REF] Chouly | An overview of recent results on Nitsche's method for contact problems[END_REF] and references therein. Moreover, Newton's semi-smooth method with Primal-Dual Active Set (PDAS) procedure appears as one of the most relevant methods to solve contact problems with friction because of their efficiency and their simplicity of implementation [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF][START_REF] Hüeber | A primal-dual active set algorithm for threedimensional contact problems with coulomb friction[END_REF][START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF]. Some works have been dedicated to studying the efficiency of PDAS methods, as well as to solve linear elasticity problems with unilateral boundary conditions [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF][START_REF] Hintermüller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF][START_REF] Abide | Analysis of two active set type methods to solve unilateral contact problems[END_REF], or even to solve non-linear multi-body contact problems in elastodynamics [START_REF] Hüeber | A primal-dual active set algorithm for threedimensional contact problems with coulomb friction[END_REF][START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF][START_REF] Abide | Inexact primal-dual active set method for solving elastodynamic frictional contact problems[END_REF].

However, when considering impact problems whether in small or large deformations, even the standard implicit schemes (θ-method, Newmark schemes, midpoints or Hilber-Hughes-Taylor type methods [START_REF] Hilber | Improved numerical disspation for time integration algorithms in structural dynamics[END_REF]) lose their unconditional stability, which leads either to an energy explosion or to a numerical dissipation which is neither realistic nor acceptable mechanically, as explained in [START_REF] Hilber | Improved numerical disspation for time integration algorithms in structural dynamics[END_REF][START_REF] Simo | The discrete energy-momentum method. part i: Conserving algorithms for nonlinear elastodynamics[END_REF][START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in non linear elasticity[END_REF][START_REF] Armero | On the formulation of high-frequency dissipative time-stepping algorithms for nonlinear dynamics. part ii: second-order methods[END_REF]. Therefore, it is necessary to use appropriate implicit schemes providing energy conservation type properties. To this purpose, many references [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF][START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF][START_REF] Simo | The discrete energy-momentum method. part i: Conserving algorithms for nonlinear elastodynamics[END_REF][START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in non linear elasticity[END_REF][START_REF] Armero | On the formulation of high-frequency dissipative time-stepping algorithms for nonlinear dynamics. part ii: second-order methods[END_REF][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] Laursen | Improved implicit integrators for transient impact problems: dynamic frictional dissipation within an admissible conserving framework[END_REF][START_REF] Ayyad | Formulation and analysis of two energy-consistent methods for nonlinear elastodynamic frictional contact problems[END_REF][START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and coulomb's friction[END_REF][START_REF] Acary | Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact[END_REF][START_REF] Barboteu | A hyperelastic dynamic frictional contact model with energy-consistent properties[END_REF] propose relevant approaches to solve these impact problems with balance energy properties. This work proposes two traits of novelty. The first one arises from the improved frictional contact model we consider, which provides intrinsic energy-controlling properties. Contact is modeled with a general Moreau-Yosida regularization [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF][START_REF] Jourani | Moreau-yosida regularization of state-dependent sweeping processes with nonregular sets[END_REF][START_REF] Nacry | Regularization of sweeping process: old and new[END_REF] of the unilateral condition. The Moreau-Yosida regularization seems to be an appropriate tool to find a regular contact model (Normal Compliance) which respects the kinetic energy of the system while preserving the non-interpenetration of contact. The discrete Improved Normal Compliance (INC) will be well suited to respect energy conservation in adequation with the continuous framework. The second trait of novelty consists in the analysis and the implementation of an energy-controlling method, based on a semi-smooth Newton method combined with an active set method via the complementarity functions for normal compliance with friction models. In summary, the major novelty comes from the combination of the Active Set approach with the Improved Normal conformance method and the fact that such a combination is also useful and relevant for the conservation of the energy in hyper-elastodynamics. Based on representative examples from the literature [START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF][START_REF] Ayyad | Formulation and analysis of two energy-consistent methods for nonlinear elastodynamic frictional contact problems[END_REF][START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF], we study and analyze numerically this Improved Normal Compliance method for dynamic hyperelastic problems with the main objective of respecting the conservation of energy and the non-interpenetration condition during impacts.

The remainder of the article is organized as follows. In Section 2, we present and explain the physical framework and the mathematical model studied, then we recall the formulation of hyperelastic problems with frictional contact. We present in detail various contact models as unilateral contact, persistent contact, and normal compliance conditions. After formulating the strong and variational problems, we detail the energy conservation and dissipation properties in the continuous case using the specific properties of the normal compliance conditions. Section 3 is devoted to the discretization of the hyperelastodynamic problem with contact and the approach (INC) adapted to respect the energy conservation in adequation to the continuous framework. In Sections 4 and 5, we propose an innovative, fast and efficient Primal Dual Active Set (PDAS) method to solve a hyperelastodynamic problem with Improved Normal Compliance and Coulomb friction. Contact and friction conditions are realized by applying an Active Set strategy via a non-linear complementarity function based on the semi-smooth Newton iterative scheme. In the last section, we provide the numerical experiments of the hyperelastodynamic problems in contact with and without friction carried out with a dynamic elastic ball, then a hyperelastic ring both launched in the direction of a rigid foundation. We present comparative studies between different numerical energy conservation methods during the impacts of the systems [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF][START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF][START_REF] Armero | Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems[END_REF][START_REF] Laursen | Improved implicit integrators for transient impact problems: dynamic frictional dissipation within an admissible conserving framework[END_REF].

2. Hyperlastic problems for low velocity impact with friction

Hyperlastic framework

A hyperelastic body occupies a bounded domain Ω ⊂ R d (d = 1, 2, 3) with a continuous Lipschitz boundary Γ, divided into three disjoint measurable parts Γ 1 , Γ 2 , and Γ 3 . We denote by x = (x i ) the point in Ω ∪ Γ used and we designate by n = (n i ) the unit outward normal over Γ. The indices i, j, k, l vary between 1 and d (d is the space dimension), and unless otherwise specified, the summation convention on repeated indices is employed. We denote by M d the space of second-order tensors on R d or, equivalently, the space of square matrices of order d. The scalar product and the norm on R d and M d are defined by

u • v = u i v i , u = (u • u) 1 2 ∀ u, v ∈ R d , σ : γ = σ ij γ ij , σ = (σ : σ) 1 2 ∀ σ, γ ∈ M d .
Let u and Π denote the displacement field and the first Piola-Kirchhoff stress tensor, respectively, and the normal and tangential components of u on Γ, which are given by

u n = u • n, u τ = u -u n n,
where n is the unit normal outside Γ. We consider that an index following a comma represents the partial derivative with respect to the corresponding spatial variable of x, u i,j = ∂u i ∂x j . Dots above a function represent partial derivatives with respect to time, i.e. u = ∂u ∂t and ü = ∂ 2 u ∂t 2 . Moreover, we recall that the divergence operator acting on a tensor field τ is Div τ = (τ ij,j ). In the problems studied later, the behavior of the material is described by a hyperelastic constitutive law. We recall that hyperelastic constitutive laws are characterized by the first Piola-Kirchhoff tensor Π, which derives from a deformation energy density W :

Ω × M d + → R, Π(x, F) = ∂ ∂F W (x, F) = ∂ F W (x, F), for all x ∈ Ω and F ∈ M d + , where M d + = {F ∈ M d : det F > 0}.
Here F is the deformation gradient defined by F = I + ∇u and ∂ F represents the differential with respect to the variable F, for more details on hyperelasticity see [START_REF] Ciarlet | Mathematical Elasticity: Volume I: three-dimensional elasticity[END_REF][START_REF] Tallec | Numerical methods for nonlinear three-dimensional elasticity[END_REF][START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF]. In what follows, we consider a dynamic problem with contact and friction in which the hyperelastic body comes into contact with a perfectly rigid obstacle Ω obs ⊂ R d (d = 1, 2, 3) which is a rigid bounded domain with a continuous Lipschitz boundary Γ obs , the so-called foundation (see Figure 1). The hyperelastic body is subjected to the action of volumetric forces with density f 0 and surface tractions with density f 2 which act on Γ 2 .

In the rest of the paper, we consider the time interval of interest [0, T ] with T > 0. We denote by t ∈ [0, T ] the time variable.

Frictional contact conditions

Suppose the body is fixed on Γ 1 and can come into contact at Γ 3 with the foundation. In the following, the frictional contact conditions are based on the combination of normal compliance conditions with a Coulomb law of dry friction on Γ 3 . We denote by ϕ : Ω × [0, T ] → R d the deformation field, with ϕ(x, t) = x + u(x, t) the position at time t ∈ [0, T ] of point x ∈ Ω. For any point x ∈ Γ 3 , we define the point y(x, t) of the foundation closest to x: y(x, t) = arg min y∈Γ obs ϕ(x, t) -y 2 .

In this way, one can define the minimal allowed contact distance (gap) between a point of Γ 3 and its orthogonal projection on the rigid foundation as follows:

g ν = (ϕ(x, t) -y(x, t)) • ν, ∀x ∈ Γ 3 ,
where ν is the inner unit normal vector to the rigid foundation. The normal force of contact Π ν , assumed to be negative, can be written in the direction ν:

Π ν = ν • Πn.
In the same way, the tangential force of contact can also be expressed according to the first tensor of Piola-Kirchhoff:

Π τ = Πn -Π ν ν,
With these definitions in place, the tangential contact velocity ġτ of a point x ∈ Γ 3 , relative to the opposite surface of the foundation, is given by

ġτ = [I d -ν ⊗ ν] u(x, t),
The conditions of contact and friction are posed on the boundary Γ 3 or Γ 3 × [0, T ] in a timedependent problem and pertain to the normal and tangential components Π ν and Π τ of the surface contact force, respectively, and the displacements u ν , u τ . With these considerations, the unilateral conditions on the boundary of contact Γ 3 are given in the following section.

Law of unilateral contact with friction

The law of unilateral contact of a solid on a rigid obstacle was proposed in 1933 by Signorini [START_REF] Signorini | Sopra alcune questioni di elastostatica[END_REF] and is written in the form of three conditions: a condition of non-penetration, a condition of compression and a condition of complementarity. Thereafter and for convenience, we will change notation for the contact variables: δ ν = g ν and λ ν = -Π ν . Therefore, the unilateral contact conditions read Non-penetration condition:

δ ν ≤ 0, Compression condition: λ ν ≥ 0,
Complementarity condition:

λ ν δ ν = 0. (2.1)
From a mechanical viewpoint, this amounts to considering a perfectly rigid foundation; no matter how much compressive force is applied, no penetration occurs. Coulomb's law of friction involves the tangential friction stress Π τ , the normal contact pressure Π ν , and the tangential contact velocity uτ that we now denote δτ = uτ , as well as the tangential stress λ τ = -Π τ , as follows:

λ τ ≤ µ |λ ν | stick status, λ τ = µ |λ ν | δτ δτ if δτ = 0 slip status. (2.2)
where µ ≥ 0 is the coefficient of friction.

If the norm of the tangential stress λ τ is less than the friction threshold µ |λ ν |, then there is sticking between the body and the foundation. If, on the other hand, this threshold is reached, then the body slides on the foundation while the tangential stress is constant and depends on the unit tangent τ = δτ / δτ .

Note that it is possible to write conditions (2.1) as the following subdifferential inclusion:

λ ν ∈ ∂Ψ R -(δ ν ) on Γ 3 × (0, T ), (2.3)
where ∂ represents the sub-differential operator in the sense of convex analysis and Ψ A denotes the indicator function of the set A ⊂ R. A similar consideration for frictional stress leads to

λ τ ∈ µλ ν ∂ δτ on Γ 3 × (0, T ). (2.4)

Law of persistent contact with friction

Such a law is slightly on the margins of the previous one insofar as its main interest resides in its natural properties of energy conservation. The persistency condition is expressed as a complementarity condition between the normal stress λ ν and the tangential velocity δν , namely:

λ ν δν = 0. (2.5)
This condition alone is sufficient to guarantee that the work of the contact normal force Γ3 λ ν δν = 0 vanishes. By combining with the unilateral contact law, it takes the following form:

if δ ν < 0, λ ν = 0, if δ ν = 0, λ ν ∈ -∂Ψ R + ( δν ). (2.6)
We refer the reader to [START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF][START_REF] Ayyad | Formulation and analysis of two energy-consistent methods for nonlinear elastodynamic frictional contact problems[END_REF] for more details.

Normal compliance law via a α-Moreau-Yosida regularization

We briefly present the α-Moreau-Yosida regularization [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF][START_REF] Jourani | Moreau-yosida regularization of state-dependent sweeping processes with nonregular sets[END_REF] of the unilateral condition of Signorini and begin with a reminder of the concepts of variational analysis. Let S be a subset of a Hilbert space H endowed with the norm • . Definition 2.1. Let f be a lower bounded semicontinuous function defined by f : H → R ∪ {+∞}. For all r > 0, the α-Moreau-Yosida envelope [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF][START_REF] Jourani | Moreau-yosida regularization of state-dependent sweeping processes with nonregular sets[END_REF][START_REF] Vilches | Regularization of perturbed state-dependent sweeping processes with nonregular sets[END_REF][START_REF] Nacry | Regularization of sweeping process: old and new[END_REF] of f , with α ≥ 2, is defined by:

f α r (z) = inf y∈H f (y) + 1 r y -z α , ∀z ∈ H. (2.7)
Theorem 2.2. Let f : H → R ∪ {-∞}, the regularization of ∂f for all r > 0 is the gradient ∇f α r associated with the envelope f α r .

For the details of proof, see [START_REF] Moreau | Proximité et dualité dans un espace hilbertien[END_REF][START_REF] Jourani | Moreau-yosida regularization of state-dependent sweeping processes with nonregular sets[END_REF][START_REF] Vilches | Regularization of perturbed state-dependent sweeping processes with nonregular sets[END_REF][START_REF] Nacry | Regularization of sweeping process: old and new[END_REF]. If we take f = Ψ R -, for r > 0 and α ≥ 2, we have

(Ψ R -) α r (z) = inf y∈H Ψ R -(y) + 1 r y -z α = inf y∈R - 1 r y -z α =: 1 r dist α R -(z). Proposition 2.3. Let α ≥ 2 on; ∇(Ψ R -) α r (z) = ∇ 1 r dist α R -(z) = 1 r proj R + (z) proj R + (z) α-2 . (2.8)
The contact with a deformable foundation is modeled by the normal compliance condition. It attributes a reactive normal pressure depending on the interpenetration of the foundation; this means that the normal stress λ ν is a function of the normal displacement δ ν . A general expression for the normal compliance condition is then given by

λ ν = p(δ ν ) on Γ 3 × (0, T ) (2.9)
where p(•) vanishes for negative arguments. The normal compliance condition was first introduced in [START_REF] Martins | Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws[END_REF][START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF]. This standard condition can be considered as a Moreau-Yosida regularization of the unilateral contact condition of Signorini with α = 2 and z = δ ν : Standard Normal Compliance (SNC) with friction

λ ν = c ν α 2 ([δ ν ] + ) α-1 on Γ 3 × (0, T ), (2.10)
where [x] + is the positive part of x ∈ R, and c ν = 1 r can be assimilated to the stiffness coefficient of the foundation. This law of normal compliance (2.10) has two particularities: first, it allows to reduce the interpenetration while the second aspect lies in its natural quasi-energy-conservation properties (see energy balance in section 2.5), which is critical from a physical point of view. In adequacy with this regularization process, it is assumed that the law of friction can expressed as

     if ||λ τ || < µ|λ ν | λ τ = c τ δτ , stick if ||λ τ || = µ|λ ν | λ τ = µλ ν δτ δτ , slip on Γ 3 × (0, T ), (2.11)
where c τ > 0 is the tangential compliance parameter.

Strong formulation of the mechanical problem

With the preceding notation, the strong formulation of the problem is the following one. Problem P. Find displacement field u : Ω × [0, T ] → R d and the stress field

Π : Ω × [0, T ] → M d such that Π = ∂ F W (F)
in Ω × (0, T ), (2.12) ) and (2.17) respectively represent the conditions of contact with normal compliance and friction described in the preceding section. Finally, (2.18) represents the initial conditions in which u 0 and u 1 are respectively the initial displacement and velocity.

Div Π + f 0 = ρü in Ω × (0, T ), (2.13)
u = 0 on Γ 1 × (0, T ), (2.14)
Πν = f 2 on Γ 2 × (0, T ), (2.15)
λ ν = c ν α 2 [δ ν ] α-1 + on Γ 3 × (0, T ), (2.16)
     ||λ τ || < µ|λ ν | λ τ = c τ δτ ||λ τ || = µ|λ ν | λ τ = µλ ν δτ δτ on Γ 3 × (0, T ), (2.17)
u(0) = u 0 , u(0) = u 1 in Ω. ( 2

Variational formulation of the problem

In order to derive the variational formulation of Problem P, additional notation and some preliminary elements are necessary. The classical notation for the Sobolev and Lebesgue spaces associated with Ω and Γ is used. We consider a closed subspace of H 1 as follows: Following the approach of Duvaut and Lions [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF], we introduce the following Hilbert spaces:

V = {v ∈ H 1 (Ω; R d ) : v = 0 on Γ 1 } and H = L 2 (Ω; R d )
they are Hilbert spaces endowed with the scalar products (u, v) V and (Π, τ ) H and their associated norms • V and • H respectively. Note that the Lagrange multipliers λ ν and λ τ are taken equal to -Π ν and -Π τ , respectively. In order to introduce a variational formulation of the mechanical problem, we consider the following spaces [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF][START_REF] Khenous | Hybrid discretization of the Signorini problem with coulomb friction. theoretical aspects and comparison of some numerical solvers[END_REF]:

X = { v| Γ3 : v ∈ V } ⊂ H 1/2 (Γ 3 ; R d ) X ν = { v ν | Γ3 : v ∈ V } , X τ = { v τ | Γ3 : v ∈ V }
and their topological dual spaces X , X ν and X τ . We assume that Γ 3 is sufficiently regular so that

X ν ⊂ H 1/2 (Γ 3 ; R), X τ ⊂ H 1/2 (Γ 3 ; R d-1 ), X ν ⊂ H -1/2 (Γ 3 ; R) and X τ ⊂ H -1/2 (Γ 3 ; R d-1
). Moreover, we denote by •, • X ν ,Xν and •, • X τ ,Xτ the corresponding duality products. For more details on the operators traces, we can refer to [START_REF] Khenous | On the discretization of contact problems in elastodynamics[END_REF][START_REF] Khenous | Hybrid discretization of the Signorini problem with coulomb friction. theoretical aspects and comparison of some numerical solvers[END_REF].

We move on to the variational formulation (or weak formulation) for Problem P. Multiplying (2.13) by any virtual velocity v and applying Green's formula, we get

- Ω Π(t) : ∇v dx + ∂Ω Π(t)n • v da + Ω f 0 (t) • v dx = Ω ρü(t) • v dx ∀v ∈ V. (2.19)
To establish the variational formulation of Problem P (2.12)-(2.18), we need some additional notation. Thus, we consider the function f : (0, T ) → V * where the exterior volume and surface force densities are assumed to be such that

f 0 ∈ L 2 (0, T ; L 2 (Ω)), f 2 ∈ L 2 (0, T ; L 2 (Γ 2 )), (2.20) so that (f (t), v) V = (f 0 (t), v) H + (f 2 (t), v) L 2 (Γ2) , ∀v ∈ V. (2.21)
Using the duality between V * and V , Green's formula and the contact conditions with friction on the boundary Γ 3 , the equation (2.19) becomes:

ρü(t), v V * ×V + Π(t), ∇v V * ×V = (f (t), v) V + Γ3 Π ν (t)v ν da + Γ3 Π τ (t) • v τ da.
Finally, we obtain the variational formulation of the contact problem with friction P using the Lagrange multiplier λ ν , related to the normal contact stress Π ν , and the Lagrange multiplier λ τ , related to the tangential contact stress Π τ , in terms of two unknown fields (here and in the following, we drop explicit mention of time dependence for ease of presentation).

Problem P V : Find the displacement field u ∈ L ∞ (0, T ; V ), with u ∈ L 2 (0, T ; V ) and ü ∈ L 2 (0, T ; V * ), the normal stress field λ ν : (0, T ) → X ν and the tangential stress field λ τ : (0, T ) → X τ such that, ∀ v ∈ V , ρü, v V * ×V + Π, ∇v V * ×V + λ ν , v ν X ν ,Xν + λ τ , v τ X τ ,Xτ = (f , v) V , in (0, T ), (2.22
)

λ ν = c ν α 2 [δ ν ] α-1 + on Γ 3 × (0, T ), (2.23)      ||λ τ || < µ|λ ν | λ τ = c τ δτ ||λ τ || = µ|λ ν | λ τ = µλ ν δτ δτ on Γ 3 × (0, T ), (2.24)
and, moreover,

u(0) = u 0 , u(0) = u 1 . (2.25)
According to Hadamard's definition of the meaning of well-posed problems, it is known that the variational problem P V is ill-posed. Indeed, this problem is very complex due to the multiplicity of its solutions for several reasons: hyperelastic medium, dynamic process, and friction conditions. For further details on the existence of solution in particular hyperelastic cases and with more regularity, one can refer to the following references: [START_REF] Ciarlet | Mathematical Elasticity: Volume I: three-dimensional elasticity[END_REF][START_REF] Tallec | Numerical methods for nonlinear three-dimensional elasticity[END_REF][START_REF] Barboteu | Analysis of a dynamic frictional contact problem for hyperviscoelastic material with non-convex energy density[END_REF].

Energy conservation properties in the continuous case

From a physical point of view, the solution of a hyperelastic problem should satisfy some conservation properties like conservation of energy, conservation of kinematic momentum, and conservation of linear momentum. We are now particularly interested in the energy conservation properties. In the absence of contact and friction, the conservation of energy can be written as follows:

t 0 Ė(s)ds = E(t) -E(0) = t 0 Ω f 0 • u dx ds + t 0 Γ2 f 2 • u da ds, (2.26)
where (s, t) ∈ (0, T ) × (0, T ) and E(t) denotes the internal energy of the system at time t, defined as

E(t) = 1 2 Ω ρ| u| 2 dx + Ω W (C)dx (2.27)
where C = F T F is the Cauchy-Green tensor and W (C) := W (F) . Moreover, the conservation of energy for hyperelastic phenomena with frictional contact is written in the following way:

E(t) -E(0) = t 0 Ω f 0 • u dxds + t 0 Γ2 f 2 • u dads - t 0 Γ3 λ ν δν dads - t 0 Γ3 λ τ • δτ dads (2.28)
Using the variational formulation of the problem by taking v = u(t, x), we obtain the frictional contact reaction work given by

W c+f = Γ3 (λ ν δν + λ τ • δτ ) da. (2.29)
where δν and δτ represent the time derivatives of δ ν and δ τ respectively. When the volume and surface forces are not taken into account, and with the presence of the normal contact compliance law without friction, the expressions of the energy allow us to obtain

E(t) -E(0) = - t 0 Γ3
λ ν δν da ds.

(2.30)

Using the general normal compliance condition (2.23), we get:

E(t) -E(0) = - t 0 Γ3 c ν α 2 (δ ν α-1 + δν da ds = - c ν 2 Γ3 δ ν α + (t) -δ ν α + (0) da (2.31)
The difference δ ν α

+ (t) -δ ν α + (0) in (2.31
) is very small since the penetrations δ ν + (t) at any time t are also small as long as the compliance parameter c ν is sufficiently large. Therefore, the energy of the system is "almost" conserved: E(t) ≈ E(0). Considering now the friction, we obtain

λ τ • δτ ≥ 0 ⇒ -W c+f ≤ 0 ⇒ E(0) ≥ E(t) (2.32)
Concerning the expression (2.32), we note a dissipation of energy between the instants 0 and t because of the friction, which is physically acceptable due to the dissipative nature of this phenomenon. When we take these energy expressions with the "persistent" conditions (2.5), we get the following results:

Case without friction :

λ ν δν = 0, λ τ • δτ = 0 ⇒ W c+f = 0 ⇒ E(0) = E(t), (2.33)
Case with friction :

λ ν δν = 0, λ τ • δτ ≥ 0 ⇒ W c+f ≥ 0 ⇒ E(0) ≥ E(t). (2.34)
Here, (2.33) expresses the conservation of total energy when the persistence condition is applied, whereas by (2.34), because of friction, we observe a dissipation of energy between the instants 0 and t.

Discrete formulation of the frictional contact problems

Variational approximation

In this section, we introduce a discrete approximation in time and space of the problem P V , based on arguments similar to those used in [START_REF] Ayyad | Formulation and analysis of two energy-consistent methods for nonlinear elastodynamic frictional contact problems[END_REF][START_REF] Ayyad | A frictionless viscoelastodynamic contact problem with energy consistent properties: Numerical analysis and computational aspects[END_REF][START_REF] Barboteu | A hyperelastic dynamic frictional contact model with energy-consistent properties[END_REF][START_REF] Barboteu | A dynamic viscoelastic contact problem with normal compliance, finite penetration and nonmonotone slip rate dependent friction[END_REF][START_REF] Barboteu | Analysis of a dynamic frictional contact problem for hyperviscoelastic material with non-convex energy density[END_REF][START_REF] Khenous | On the discretization of contact problems in elastodynamics[END_REF][START_REF] Khenous | Hybrid discretization of the Signorini problem with Coulomb friction. theoretical aspects and comparison of some numerical solvers[END_REF]. First of all, we recall some preliminary elements concerning the time discretization step. Let N be an integer, and ∆t = T N a time step. For a continuous function f with respect to time, we will use the notation

f j = f (t j ) for 0 ≤ j ≤ N . In what follows, we consider a collection of discrete times {t n } N n=0 which defines a uniform partition of the time interval [0, T ] = N n=1 [t n-1 , t n ]
, with t 0 = 0, t n = t n-1 + ∆t, and t N = T . Using truncated Taylor expansions, we find

u n = u n-1 + ∆t 2 ( un + un-1 ) + θ(∆t 2 ), (3.1)
Finally, for a finite sequence {u n } N n=1 , we denote the divided differences of the midpoints by: un-

1 2 = u n -u n-1 ∆t = 1 2 ( un + un-1 ), (3.2)
Then we use the notation n-1 2 = 1 2 ( n + n-1 ), where n represents the approximation of (t n ). For regularity in time of the solution, see [START_REF] Hilber | Improved numerical disspation for time integration algorithms in structural dynamics[END_REF][START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF][START_REF] Wriggers | Computational contact mechanics[END_REF]. Let us note that the time integration scheme employed is based on an implicit scheme of order 2 which one finds in (3.2). We now present some elements concerning the spatial discretization step. Let Ω be a polyhedral domain. Consider a regular partition {T h } of triangular finite elements of Ω which are compatible with the boundary decomposition Γ = Γ 1 ∪ Γ 2 ∪ Γ 3 , i.e., if one side of an element T ∈ T h has more than one point on Γ, then the side lies entirely on Γ 1 , Γ 2 or Γ 3 . The space V is approximated by the finite dimensional space V h ⊂ V of continuous and piecewise affine functions, that is:

V h = { v h ∈ [C(Ω)] d : v h | T ∈ [P 1 (T)] d ∀ T ∈ T h , v h = 0 at nodes on Γ 1 },
where P 1 (T) represents the space of polynomials of degree less than or equal to 1 in T and h > 0 is the spatial discretization parameter. Let us consider the following spaces [START_REF] Khenous | On the discretization of contact problems in elastodynamics[END_REF][START_REF] Khenous | Hybrid discretization of the Signorini problem with coulomb friction. theoretical aspects and comparison of some numerical solvers[END_REF]:

X h ν = v h ν | Γ3 : v h ∈ V h , X h τ = v h τ | Γ3 : v h ∈ V h
and their topological dual spaces X h ν and X h τ . We assume that

X h ν ⊂ X ν ∩ L 2 (Γ 3 ; R) and X h τ ⊂ X ν ∩ L 2 (Γ 3 ; R d-1
). For the discretization of the Lagrange multiplier spaces X h ν and X h τ , we use piecewise constant functions as done in [START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF][START_REF] Pietrzak | Large deformation frictional contact mechanics: continuum formulation and augmented lagrangian treatment[END_REF][START_REF] Khenous | On the discretization of contact problems in elastodynamics[END_REF][START_REF] Barboteu | A hyperelastic dynamic frictional contact model with energy-consistent properties[END_REF][START_REF] Bonaldi | Energy-stable discretization of two-phase flows in deformable porous media with frictional contact at matrix-fracture interfaces[END_REF][START_REF] Bonaldi | Numerical analysis of a mixed-dimensional poromechanical model with frictionless contact at matrix-fracture interfaces[END_REF]. The discrete Lagrange multiplier spaces, denoted by X h ν and X h τ , are related to the discretization of the normal stress λ ν and the discretization of the friction stress λ τ , respectively.

With these previous notation and the midpoint scheme (3.2), we have the following discrete approximation for the problem P V at time t n- 1 2 :

Problem P h V . Find a discrete displacement field u h = {u h n } N n=0 ⊂ V h , a discrete normal stress field λ h nu = {λ ν h n } N n=0 ⊂ X h ν , and a discrete tangential stress field λ h τ = {λ τ h n } N n=0 ⊂ X h ν such that, for all n = 1, . . . , N , ρü h n-1 2 + B(u h n-1 2 ) + λ ν h n-1 2 ν n-1 2 + λ τ h n-1 2 -f = 0 (3.3) λ ν h n-1 2 = c ν α 2 [δ h ν n-1 2 ] α-1 + , (3.4)
         λ τ h n-1 2 < µ|λ ν h n-1 2 | λ τ h n-1 2 = c τ δτ h n-1 2 λ τ h n-1 2 = µ|λ ν h n-1 2 | λ τ h n-1 2 = µλ ν h n-1 2 δτ h n-1 δτ h n-1 2 (3.5) u h 0 = ūh 0 , uh 0 = ūh 1 .
where B(•) is the internal stress operator representing the first Piola-Kirchhoff tensor Π. The initial values ūh 0 ∈ V h and ūh 1 ∈ V h are discrete values resulting from the finite element approximation of u 0 and u 1 , respectively.

Usual discrete framework of energy conservation

In the rest of the section, to simplify notation and readability, we do not indicate the dependence of the different variables on the discretization parameter h, i.e., we write u instead of u h . We begin by recalling some preliminaries regarding discrete energy conservation in the hyperelastodynamic contactless framework. In order to solve a hyperelastic dynamic problem, we have to use adapted time integration schemes. When considering nonlinear dynamical problems, standard implicit schemes (θ-method, Newmark schemes, midpoint methods or HHT methods) lose their unconditional stability, as explained in [START_REF] Hilber | Improved numerical disspation for time integration algorithms in structural dynamics[END_REF][START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in non linear elasticity[END_REF][START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF]. Therefore, it is necessary to use implicit energy-conserving schemes like those used in [START_REF] Armero | On the formulation of high-frequency dissipative time-stepping algorithms for nonlinear dynamics. part ii: second-order methods[END_REF][START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in non linear elasticity[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF][START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF][START_REF] Simo | The discrete energy-momentum method. part i: Conserving algorithms for nonlinear elastodynamics[END_REF][START_REF] Ayyad | Formulation and analysis of two energy-consistent methods for nonlinear elastodynamic frictional contact problems[END_REF][START_REF] Barboteu | A hyperelastic dynamic frictional contact model with energy-consistent properties[END_REF][START_REF] Acary | Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact[END_REF] which are appropriate due to their long-term time integration accuracy and stability. In all these methods, the corresponding discrete mechanical conservation properties are satisfied. By taking into account the implicit second-order temporal integration scheme of midpoint (3.2), we obtain the weak form of a nonlinear hyperelastodynamic problem integrated between times t n-1 and t n :

   Find u n ∈ U such that 1 ∆t Ω ρ( un -un-1 ) • v dx + Ω Π algo n-1 2 : ∇v dx - Γ2 f 2n-1 2 • v dx - Ω f 0n-1 2 • v da = 0. (3.6)
There, the discrete tensor Π algo is introduced in order to satisfy the exact properties of the discrete energy. This tensor, defined by Gonzalez in [START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in non linear elasticity[END_REF], takes the form:

     Π algo n-1 2 = F n-1 2 Σ algo Σ algo = 2 ∂ W ∂C (C n-1 2 ) + 2[ W (C n ) -W (C n-1 ) - ∂ W ∂C (C n-1 2 ) : ∆C n-1 ] ∆C n-1 ∆C n-1 : ∆C n-1 (3.7) with ∆C n-1 = C n -C n-1 and C n-1 = F T n-1 F n-1 .
Using the arguments of [START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in non linear elasticity[END_REF] and as shown by the axiom of material indifference which implies that W (F ) = W (C), it follows that (3.7) verifies the exact conservation of energy characterized by the following condition:

Π algo n-1 2 : (F n -F n-1 ) = Π algo n-1 2 : (∇u n -∇u n-1 ) = W (C n ) -W (C n-1 ). (3.8)
For more details on the standard energy conservation framework, we refer the reader to [START_REF] Armero | On the formulation of high-frequency dissipative time-stepping algorithms for nonlinear dynamics. part ii: second-order methods[END_REF][START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in non linear elasticity[END_REF][START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF][START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF][START_REF] Simo | The discrete energy-momentum method. part i: Conserving algorithms for nonlinear elastodynamics[END_REF]. Much work has been devoted to extending the conservative properties of the frictionless contact; more precisely, Lauren and Chawla [START_REF] Laursen | Design of energy-conserving algorithms for frictionless dynamic contact problems[END_REF] and Armero and Petocz [START_REF] Armero | Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems[END_REF] showed the advantage of the persistence condition to conserve energy in the discrete setting. Nevertheless, in all these works, the numerical method shows that the interpenetration only disappears when the time step tends towards zero. In order to overcome this drawback, Laursen and Love [START_REF] Laursen | Improved implicit integrators for transient impact problems: dynamic frictional dissipation within an admissible conserving framework[END_REF] have developed an efficient method, by introducing a discrete jump in velocity; however, this method requires the solution of an auxiliary system in order to compute the velocity update results. Furthermore, Hauret and Le Tallec [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF] have considered a specific penalized enforcement of the contact conditions which allows to provide energy conservation properties. Then, Khenous, Laborde and Renard [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF] have introduced the Equivalent Mass Matrix method (EMM), based on a procedure of redistribution of the mass matrix. Interpretations and extensions of this method can be found in [START_REF] Hauret | Mixed interpretation and extensions of the equivalent mass matrix approach for elastodynamics with contact[END_REF]. The resulting problem exhibits Lipschitz regularity in time and achieves good energy evolution properties, due to the fact that the persistency condition is automatically satisfied. The EMM approach was studied and used in many works; for instance, theoretical and computational aspects related to this model can be found in [START_REF] Hauret | Energy consistent co-rotational schemes for frictional contact problems[END_REF][START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF].

Improved approach to "almost" conserve the energy

In what follows, based on the papers [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF][START_REF] Barboteu | A hyperelastic dynamic frictional contact model with energy-consistent properties[END_REF], we present a improved energy conservation method for hyperelastodynamic contact problems with its extension to dissipation phenomena with friction. This method allows to enforce the general normal compliance law during each time step with "minimal" contact penetrations and with conservation properties which respect "almost" the energy. In order to take into account contact at time t n- 1 2 , we choose to implicitly approach the term of contact with friction; thus the weak formulation integrated between the times t n-1 and t n is

           Find u n ∈ V such that 1 ∆t Ω ρ( un -un-1 ) • v dx + Ω Π algo n-1 2 : ∇v dx - Γ2 f 2n-1 2 • v dx - Ω f 0n-1 2 • v da + Γ3 [λ ν n-1 2 δν n-1 2 + λ τ n-1 2 • δτ n-1 2
] da = 0 (3.9)

In order to obtain energy conservation properties, we propose to change the normal constraint λ ν n-1 2 by an improved discrete value (3.11) which will respect the energy balance of the continuous case (2.31).

Discrete form of Improved Normal Compliance (INC) conditions

In order to "almost" conserve the energy and to respect the energy balance of the continuous case, we replace the normal contact distance (δ n ν ) α-1 by

δ n ν := δ n ν + α -δ n-1 ν ) + α α δ n ν -δ n-1 ν (3.10)
We then obtain the normal stress value λ ν n-1 2 of the improved normal compliance condition in the discrete case:

λ ν n-1 2 = c ν α 2 δ n ν . (3.11)

Discrete Energy Evolution Analysis

This part is devoted to establishing the energy conservation properties induced by the improved normal compliance described in the previous paragraph. We use above the notation E n and E n-1 for the energy E of the hyperelastic system of contact with friction evaluated at times t n and t n-1 respectively. For example, the discrete energy at time t n can be written as follows:

E n = 1 2 Ω ρ| u2 n | dx + Ω W (C n ) dx. (3.12)
The general evaluation of the discrete energy of the contact problem with friction between times t n and t n-1 is based on the following proposition.

Proposition 3.1. The following discrete energy balance holds between times t n and t n-1 :

E n -E n-1 = ∆t f n-1 2 , u n-1 2 V * ×V -∆t Γ3 λ ν n-1 2 δν n-1 2 + λ τ n-1 2 • δτ n-1 2 da (3.13)
Proof. Using the variational formulation (3.9) with

v = un-1 2 = u n -u n-1 ∆t = un + un-1 2
we get the following equality

1 2∆t Ω ρ( un -un-1 ) • ( un + un-1 ) dx + 1 ∆t Ω Π algo : ∇(u n -u n-1 ) dx = f n-1 2 , u n-1 2 V * ×V + Γ3 λ ν n-1 2 δν n-1 2 + λ τn • δτ n-1 2 da (3.14)
Also, using the identity ( unun-1 )

• ( un + un-1 ) = [ un ] 2 -[ un-1 ] 2
and the conservation property of the Gonzalez scheme given in equality (3.8), we get that

1 2∆t Ω ρ [ un ] 2 -[ un-1 ] 2 dx + 1 ∆t Ω W (C n ) -W (C n-1 ) dx = f n-1 2 , u n-1 2 V * ×V - Γ3 λ ν n-1 2 δν n-1 2 + λ τ n-1 2 • δτ n-1 2 da (3.15)
By using the definition (3.12) of the discrete energy at times t n-1 and t n , we obtain the assertion.

Using the previous proposition, we can give an estimate of the discrete energy balance for the contact law with Improved Normal Compliance (3.11). When the external forces are assumed to be zero, by using

δ n ν -δ n-1 ν ∆t = δν n-1 2
and considering the formula (3.10), the energy balance is

E n -E n-1 = - Γ3 c ν 2 [(δ n ν ) + ] α -[(δ n-1 ν ) + ] α da -∆t Γ3 λ τ n-1 2 • δτ n-1 2 da. (3.16)
We notice that the condition of Improved Normal Compliance allows in the case without friction an evaluation of the discrete energy which is in agreement with the continuous case (2.31).

Case without friction

The difference [(

δ n ν ) + ] α -[(δ n-1 ν
) + ] α is very small since the penetrations (δ n ν ) + and (δ n-1 ν ) + are also small. So the energy of the system is "almost" conserved, i.e. E n ≈ E n-1 .

Case with friction

In this case the product

λ τ n-1 2 • δτ n-1 2
is always positive due to the friction law, so we observe a dissipation of the energy: E n ≤ E n-1 . In other words this strategy limits the dissipation of energy between times t n and t n-1 . In summary, the INC strategy respects the dissipation in the case of friction and "almost" conserves the energy in the case of contact without friction, and this is achieved by limiting the penetration. In the following sections, we will propose both semi-smooth Newton's method as well as the Primal Dual Active Set (PDAS) algorithm in the case of Normal Compliance conditions.

Semi-Smooth Newton approach for solving normal compliance conditions

The semi-smooth Newton/PDAS methods appear to be one of the most relevant methods for solving frictional contact problems (cf [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF][START_REF] Hüeber | A primal-dual active set algorithm for threedimensional contact problems with coulomb friction[END_REF][START_REF] Hüeber | A primal-dual active set strategy for non-linear multibody contact problems[END_REF]). These methods are based on the following principle: the conditions of contact and friction are reformulated in terms of non-linear complementarity equations whose solution is provided by the semi-smooth iterative method of Newton [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF][START_REF] Hintermüller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF]. To this end, we need the generalized derivative of complementary functions for contact and friction. In practice, the conditions of contact with Coulomb's friction can be formulated in terms of a fixed point problem related to a quasi-optimization one. From a purely algorithmic point of view, the main goal of these methods is to separate the nodes potentially in contact into two subsets (active and inactive) and to find the correct subset of all the nodes actually in active contact (subset A), as opposed to those that are inactive (subset I). In practice, the semi-smooth Newton/PDAS methods do not require the use of Lagrange multipliers. In fact, the boundary conditions on the subsets A and I are directly enforced thanks to the fixed points found by the semi-smooth Newton method, and consequently, their implementation can be achieved without much effort.

Standard Normal Compliance (SNC) conditions

Complementary function

The standard normal compliance contact conditions (2.10) with α = 2 can be formulated from the following non-linear complementary function:

C λ ν (δ ν , λ ν ) = λ ν -c ν [δ ν ] + , (4.1)
where we have dropped, for the sake of readability, the time index n + 1.

Generalized derivative of complementary functions

We provide the generalized derivative of the complementary functions in the gap and contact cases.

• Gap case: δ ν ≤ 0.

According to the complementary function C λ ν (δ ν , λ ν ) = λ ν we have

d δν C λ ν = 0, (4.2)
d λν C λ ν = dλ ν . (4.3)
• Contact case: δ ν > 0.

Given the complementary function

C λ ν (δ ν , λ ν ) = λ ν -c ν δ ν , we have d δν C λ ν = -c ν dδ ν , (4.4)
d λν C λ ν = dλ ν . (4.5)
By combining (4.2)-(4.5), with the generalized derivative D C λ ν of C λ ν , we obtain

D C λ ν (δ ν , λ ν )(dδ ν , dλ ν ) = -c ν (1 Contact )dδ ν + dλ ν , (4.6)
where

1 Contact = 0 if δ ν ≤ 0, 1 Contact = 1 if δ ν > 0.
D C λ ν (δ k ν , λ k ν )(∆δ k+1 ν , ∆λ k+1 ν ) = -C λ ν (δ k ν , λ k ν ), (4.7) (δ k+1 ν , λ k+1 ν ) = (δ k ν , λ k ν ) + (∆δ k+1 ν , ∆λ k+1 ν
).

• Gap case: 1 Contact = 0.

From the equations (4.7) we have

λ k+1 ν -λ k ν = -λ k ν . (4.8)
Next, the gap conditions of the semi-smooth Newton formalism are as follows λ k+1 ν = 0. (4.9)

• Contact case:

1 Contact = 1.
From the equations (4.7) we have

-c ν (δ k+1 ν -δ k ν ) + (λ k+1 ν -λ k ν ) = -λ k ν + c ν δ k ν . (4.10)
Next,

λ k+1 ν = c ν δ k+1 ν . (4.11)

Improved Normal Compliance (INC) conditions

Complementary function

The improved normal compliance contact condition (3.11) can be formulated from the following non-linear complementary function:

C λ ν (δ n+1 ν , λ n+1 ν ) = λ n+1 ν -[c ν α 2 δ n+1 ν ] + . (4.12)
For the sake of readability, from now on, we use this equation:

C λ ν (δ ν , λ ν ) = λ ν -[c ν α 2 δ n+1 ν ] + .

Generalized derivative of complementary functions

Now, we provide the generalized derivative of the complementary functions in the gap and contact cases.

• Gap case:

δ n+1 ν ≤ 0.
According to the complementary function C λ ν (δ ν , λ ν ) = λ ν we have the following derivative

d δν C λ ν = 0, (4.13)
d λν C λ ν = dλ ν . (4.14) • Contact case: δ n+1 ν > 0.
Given the complementary functions C λ ν (δ ν , λ ν ) = λ ν -c ν δ ν , we have

d δν C λ ν = -c ν α(([δ n+1 ν ] + ) α-1 -δ n+1 ν ) 2(δ n+1 ν -δ n ν ) dδ ν , (4.15)
d λν C λ ν = dλ ν . ( 4
D C λ ν (δ ν , λ ν )(dδ ν , dλ ν ) = -c ν α(([δ n+1 ν ] + ) α-1 -δ n+1 ν ) 2(δ n+1 ν -δ n ν ) (1 Contact )dδ ν + dλ ν , (4.17)
where

1 Contact = 0 if δ n+1 ν ≤ 0, 1 Contact = 1 if δ n+1 ν > 0.
D C λ ν (δ k ν , λ k ν )(∆δ k+1 ν , ∆λ k+1 ν ) = -C λ ν (δ k ν , λ k ν ), (4.18)
(δ k+1 ν , λ k+1 ν ) = (δ k ν , λ k ν ) + (∆δ k+1 ν , ∆λ k+1 ν
).

• Gap case:

1 Contact = 0
From the equations (4.18) we have

λ k+1 ν -λ k ν = -λ k ν . (4.19)
Next, the gap conditions of the semi-smooth Newton formalism are as follows

λ k+1 ν = 0. (4.20)
• Contact case:

1 Contact = 1
From the equations (4.18) we have

-c ν α(([δ k,n+1 ν ] + ) α-1 -δ k,n+1 ν ) 2(δ k,n+1 ν -δ n ν ) (δ k+1,n+1 ν -δ k,n+1 ν ) + (λ k+1 ν -λ k ν ) = -λ k ν + c ν α 2 δ k,n+1 ν . (4.21)
Next,

λ k+1 ν = c ν α 2 δ k,n+1 ν + c ν α(([δ k,n+1 ν ] + ) α-1 -δ k,n+1 ν ) 2(δ k,n+1 ν -δ n ν ) (δ k+1,n+1 ν -δ k,n+1 ν ). ( 4
(δ n+1 ν , δn+1 τ , λ n+1 ν , λ n+1 τ ) C λ τ (δ n+1 ν , δn+1 τ , λ n+1 ν , λ n+1 τ ) = max(µλ n+1 ν , c τ δn+1 τ )λ n+1 τ -µλ n+1 ν (c τ δn+1 τ ). (4.23)
For the sake of readability, in the following we use this equation:

C λ τ (δ ν , δτ , λ ν , λ τ ) = max(µλ ν , c τ δτ )λ τ -µλ ν (c τ δτ ).

Generalized derivative of complementary functions

Now, we provide the generalized derivative of the complementary function in the gap and friction cases.

• Gap case:

λ τ = 0, C λ τ (δ ν , δτ , λ ν , λ τ ) = 0. • Stick case: λ τ < µλ ν : C λ τ (δ ν , δτ , λ ν , λ τ ) = µλ ν λ τ -µc τ λ ν δτ .
Then

d δν C λ ν = 0, (4.24) d δτ C λ τ = -µc τ λ ν d δτ , (4.25
)

d λν C λ τ = (µλ τ -µc τ δτ )dλ ν , (4.26
)

d λτ C λ τ = µλ ν dλ τ . (4.27) • Slip case: λ τ ≥ µλ ν C λ τ (δ ν , δτ , λ ν , λ τ ) = c τ δτ λ τ -µλ ν (c τ δτ ).
Then

d δν C λ τ = 0, (4.28)
d δτ C λ τ = c τ λ τ ( δτ ) T δτ -µc τ λ ν I 2 d δτ , (4.29)
d λν C λ τ = -µc τ δτ dλ ν , (4.30)
d λτ C λ τ = c τ δτ dλ τ . ( 4
G C λ τ (δ ν , δτ , λ ν , λ τ )(∆δ ν , ∆ δτ , ∆λ ν , ∆λ τ ) = 1 Stick µλ τ ∆λ ν + 1 Slip c τ λ τ ( δτ ) T δτ ∆ δτ -µc τ δτ ∆λ ν -µc τ λ ν ∆ δτ + 1 Stick µλ ν ∆λ τ + 1 Slip c τ δτ ∆λ τ
where

1 Stick = 1, 1 Slip = 0 if λ τ < µλ ν , 1 Stick = 0, 1 Slip = 1 if λ τ ≥ µλ ν .
Using now the semi-smooth Newton formalism at the current iterate (δ

(k) ν , δ(k) τ , λ (k) ν , λ (k) τ ), one can derive the new iterate (δ (k+1) ν , δ(k+1) τ , λ (k+1) ν , λ (k+1) τ
) as follows:

G C λ τ (δ (k) ν , δ(k) τ , λ (k) ν , λ (k) τ )(∆δ (k+1) ν , ∆ δ(k+1) τ , ∆λ (k+1) ν , ∆λ (k+1) τ ) = -C λ τ (δ (k) ν , δ(k) τ , λ (k) ν , λ (k) τ ), (δ (k+1) ν , δ(k+1) τ , λ (k+1) ν , λ (k+1) τ ) = (δ (k) ν , δ(k) τ , λ (k) ν , λ (k) τ ) + (∆δ (k+1) ν , ∆ δ(k+1) τ , ∆λ (k+1) ν , ∆λ (k+1) τ
).

• Stick case:

1 Stick = 1, 1 Slip = 0.
We have

-µc τ λ (k) ν ( δ(k+1) τ - δ(k) τ ) + µλ (k) ν (λ (k+1) τ -λ (k) τ ) + (µλ (k) τ -µc τ δ(k) τ )(λ (k+1) ν -λ (k) ν ) = -µλ (k) ν (λ (k) τ + µc τ λ (k) ν δ(k) τ .
Next, with λ (k) τ = c τ δ(k) τ , we obtain

λ (k+1) τ = c τ δ(k+1) τ .
(4.32)

• Slip case:

1 Stick = 0, 1 Slip = 1.
We obtain

c τ λ (k) τ ( δ(k) τ ) T δ(k) τ -µc τ λ (k) ν I 2 ( δ(k+1) τ - δ(k) τ ) (4.33) -µ(c τ δ(k) τ )(λ (k+1) ν -λ (k) ν ) + c τ δ(k) τ (λ (k+1) τ -λ (k) τ ) = -c τ δ(k) τ λ (k) τ + µλ (k) ν (c τ δ(k) τ ).
Therefore, after an elementary computation, with

λ (k) τ = µλ (k) ν δ(k) τ δ(k) τ we have λ (k+1) τ = µλ (k+1) ν δ(k) τ δ(k) τ -λ (k) τ ( δ(k) τ ) T δ(k) τ 2 -µλ (k) ν I 2 δ(k) τ ( δ(k+1) τ - δ(k) τ ).
For a two dimensional case, one obtains a simplified version of the previous condition:

λ (k+1) τ = µλ (k+1) ν δ(k) τ δ(k) τ .
For more details regarding the obtention of such a simplified version, we refer to the proof provided in [START_REF] Abide | Inexact primal-dual active set method for solving elastodynamic frictional contact problems[END_REF].

Primal-Dual Active Set methods

This section is devoted to the numerical treatment of the contact conditions using a Primal-Dual Active Set method within the framework of dynamic contact problems. After defining the active and inactive subsets of all nodes that are potentially in contact, we compute the contact conditions on each subset only in terms of contact reaction, using the local general equations of motion.

Primal-Dual Active Set method for standard normal compliance

Let us denote by S the set of potential contact and γ a potential contact node belonging to S. The standard normal contact condition (2.10) with α = 2 is enforced by applying an active set strategy which derives directly from the computation of the fixed point on the non-linear complementary functions C λ ν and C λ τ based on the Newton semi-smooth scheme. The active and inactive sets are defined as follows

A k+1 ν = {γ ∈ S : δ γ,k ν ≥ 0}, I k+1 ν = {γ ∈ S : δ γ,k ν < 0}, A k+1 τ = {p ∈ S : λ γ,k τ < µλ γ,k ν }, I k+1 τ = {p ∈ S : λ γ,k τ ≥ µλ γ,k ν }.
The status of a given potential γ at the non-linear iteration k depends on the set it belongs to. It can be either in the non-contact or frictional contact status (either stick or slip status). It yields the following Algorithm 1

Algorithm 1 PDAS for standard normal compliance (i) Choose (δ (0) , λ (0) ), c ν > 0, c τ > 0 and set k = 0.

(ii) Compute:

τ γ ν = δ γ,k ν and τ γ τ = -λ γ,k τ + µλ γ,k ν
for each γ ∈ S. (iii) Set the active and inactive sets:

A k+1 ν = {γ ∈ S : τ γ ν ≥ 0}, I k+1 ν = S \ A k+1 ν , A k+1 τ = {p ∈ S : τ γ τ > 0}, I k+1 τ = S \ A k+1 τ . (iv) Find (δ γ,k+1 , λ γ,k+1 ) such that λ γ,k+1 ν = 0, λ k+1 τ,p = 0 ∀γ ∈ I k+1 ν , (5.1)
λ γ,k+1 ν = c ν δ γ,k+1 ν ∀γ ∈ A k+1 ν , (5.2)
λ γ,k+1 τ = c τ δγ,k+1 τ ∀γ ∈ A k+1 τ ∩ A k+1 ν , (5.3) λ γ,k+1 τ = µλ γ,k+1 ν δγ,k τ δγ,k τ -λ γ,k τ ( δγ,k τ ) T δγ,k τ 2 -µλ γ,k ν I2 δγ,k τ ( δγ,k+1 τ - δγ,k τ ) ∀γ ∈ I k+1 τ ∩ A k+1 ν . (5.4) (v) If (δ γ,k+1 , λ γ,k+1 ) -(δ γ,k , λ γ,k ) ≤ , A k+1 ν = A k ν and A k+1 τ = A k τ stop, else go to (ii).

Primal-Dual Active Set method for improved normal compliance

Likewise, using similar notation, the improved contact condition (3.11) are enforced by applying another active set strategy yielding a different Algorithm 2, described below.

Note that the main difference between Algorithms 1 and 2 lies in step (iv), regarding the normal compliance handling (respectively for 5.2 and 5.6), as one would expect. For the convergence of this type of Semi-Smooth Newton algorithm, we can refer to the following papers: [START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF], [START_REF] Hintermüller | The primal-dual active set strategy as a semismooth Newton method[END_REF], [START_REF] Hintermüller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF].

Numerical experiments

In what follows, we carry out a comparative study of methods of the energy conservation type, focusing more particularly on the behavior of the discrete energy of the system during and after the impact. To do this, we consider two numerical examples, the first concerns the impact of an elastic ball and the second represents the impact of a hyperelastic ring. The goal is to demonstrate that the Active Set-Improved Normal Compliance method respects the conservation of energy after impact and practically guarantees non-penetration.

Algorithm 2 PDAS for improved normal compliance (i) Choose (δ (0) , λ (0) ), c ν > 0, c τ > 0 and set k = 0. (ii) Compute: τ γ ν = δ γ,k,n+1 ν and τ γ τ = -λ γ,k τ + µλ γ,k ν
for each γ ∈ S. (iii) Set the active and inactive sets:

A k+1 ν = {γ ∈ S : τ γ ν ≥ 0}, I k+1 ν = S \ A k+1 ν , A k+1 τ = {p ∈ S : τ γ τ > 0}, I k+1 τ = S \ A k+1 τ . (iv) Find (δ γ,k+1 , λ γ,k+1 ) such that λ γ,k+1 ν = 0, λ γ,k+1 τ,p = 0 ∀γ ∈ I k+1 n , (5.5)
λ γ,k+1 ν = c ν α 2 δ γ,k,n+1 ν + c ν α(([δ γ,k,n+1 ν ]+) α-1 -δ γ,k,n+1 ν ) 2(δ γ,k,n+1 ν -δ γ,n ν ) (δ γ,k+1,n+1 ν -δ γ,k,n+1 ν ) ∀γ ∈ A k+1 ν , (5.6) λ γ,k+1 τ = c δγ,k+1 τ ∀γ ∈ A k+1 τ ∩ A k+1 ν , (5.7)
λ γ,k+1 τ = µλ γ,k+1 ν δγ,k τ δγ,k τ -λ γ,k τ ( δγ,k τ ) T δγ,k τ 2 -µλ γ,k ν I2 δγ,k τ ( δγ,k+1 τ - δγ,k τ ) ∀γ ∈ I k+1 τ ∩ A k+1 ν . (5.8) (v) If (δ γ,k+1 , λ γ,k+1 ) -(δ γ,k , λ γ,k ) ≤ , A k+1 ν = A k ν and A k+1 τ = A k τ stop, else goto (ii).

Impact of a linear elastic ball against a foundation

This representative benchmark problem describes the frictionless impact of a linear elastic ball against a foundation. The elastic ball is launched with an initial velocity (u 1 = (0, -10) m/s) towards the foundation (x 1 , x 2 ) ∈ R 2 : x 2 ≤ 0 . The domain Ω represents the cross-section of the ball, under the assumption of plane stress. The behavior of the material is described by a linear elastic constitutive law defined by the energy density

W e (ε) = λ * 2 (tr ε) 2 + µ tr(ε 2 ), ∀ε ∈ M d . (6.1)
with

λ * = 2λµ λ + 2µ , µ = Eκ 2(1 + κ) and λ = Eκ 2(1 + κ)(1 -2κ)
Here, E and κ are respectively the Young modulus and the Poisson ratio of the material and tr(•) is the trace operator. Note that ε = 1 2 (∇u T + ∇u) represents the linearized strain tensor within the framework of the theory of small strains ( u 1 and ∇u 1 in Ω). The physical setting is shown in Figure 2. Here:

Ω = (x 1 , x 2 ) ∈ R 2 : (x 1 -100) 2 + (x 2 -100) 2 ≤ 100 , Γ 1 = ∅, Γ 2 = ∅, Γ 3 = (x 1 , x 2 ) ∈ R 2 : (x 1 -100) 2 + (x 2 -100) 2 = 100 .
We assume that the volume forces do not act on the body during the process. For the discretization of the problem of contact represented in Figure 2, we use 7820 elastic nodes. For numerical experiments, we use the following data: ρ = 1000 kg/m 3 , T = 2 s, ∆t = 0.001 s, u 0 = (0, 0) m, u 1 = (0, -10) m/s, E = 100 GP a, κ = 0.35, f 0 = (0, 0) P a, g = 50 m, µ = 0. In Figure 3, the successive positions of the deformed ball as well as the contact forces are represented before, during, and after the impact. The interest of this representative example lies in the comparison of the numerical results of the Improved Normal Compliance-PDAS methods (α = 2) with other classical numerical methods. For this, we consider five existing methods:

-The classical quasi-Lagrangian method with the Signorini contact condition.

-The penalty method with a standard normal compliance condition of the form λ ν = c ν (δ ν ) + .

-The Equivalent Mass Matrix (EMM) method proposed by Khenous [START_REF] Laursen | Improved implicit integrators for transient impact problems: dynamic frictional dissipation within an admissible conserving framework[END_REF], which represents a specific distribution of the mass matrix without any inertia of the contact nodes. This method is characterized by relevant stability properties of the contact stress.

-Newton's adapted continuity method, developed by Ayyad and Barboteu [START_REF] Alart | A mixed formulation for frictional contact problems prone to newton like solution methods[END_REF], which is characterized by enforcing, after two steps, the unilateral contact law and the persistence condition during each time increment.

-The Active Set method with the persistent contact law.

In what follows, we analyze the methods in terms of discrete energy evolution. For this we introduce the total discrete energy at time t n which is given by the following formula:

E n = 1 2 Ω ρ| u2 n |dx + Ω σ e n : ε(u n )dx,
where σ e = ∂W e ∂ε denotes the stress tensor for infinitesimal strains. Figure 4 represents the evolution of the total discrete energy of the dynamical system. We note that after the impact (i.e. for t ≥ 1.52 s) and for the considered time step ∆t = 1.0 e -3 s, the classical quasi-Lagrangian method with Signorini's law (curve ) as well as the method with standard normal compliance condition (curve ) are characterized by non-conservation of energy, which is not realistic from a physical point of view. We also notice that the EMM method (curve ) strongly reduces the energy dissipation, without obtaining the exact conservation. It turns out the schemes developed by Ayyad and Barboteu [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF] (curve •) and the Improved Normal Compliance-PDAS method (curve ) do conserve energy after impact. However, for the penalized method, we find energy fluctuations which disappear after the impact. Besides. for the penalized method and the method used in [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF], the unilateral contact is not exactly satisfied, see Table 1. Indeed, the penalized method (standard normal compliance) generates a maximum error on the displacement of normal contact of 1.4e -4 m and 5.1e -3 m for the method of Ayyad and Barboteu [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF]. The Improved Normal Compliance-PDAS method allows to obtain a better energy conservation and in addition to limit the penetration: 5.7e -4 m. The Active Set method for persistent contact (shown by the curve ) enforces exactly the energy conservation without any fluctuation. Due to the "leapfrog" time step predictor, this Active Set method generates a maximum error on the normal contact displacement of 1.54e -2 m, larger than all the other methods. In order to overcome these difficulties, we consider the Active Set Methods Maximum error on the δ ν (m) Quasi-Lagrangian with Signorini law 0. Scheme with Signorini and persistent (Ayyad- Barboteu) 5.1e -3 Active set with persistent conditions 1.54e -2 Active set with classical normal compliance 1.4e -4 Active set with improved normal compliance (α = 2) 5.7e -4 Active set with improved normal compliance (α = 3 ) 8.05e -3

Table 1: Maximum error on normal contact displacement (∆t = 0.001, cν = 1e 3 ).

method with improved normal compliance, by analyzing the behavior of discrete energy with respect to several parameters (c ν , α, ∆t). Figure 5 allows to evaluate the influence of the normal compliance parameter c ν . We see that, for c ν > 1e 4 , this method respects the conservation of energy after impact and practically guarantees the non-penetration for α = 2. Nevertheless, this method generates slight fluctuations of the discrete energy during the impact. For c ν < 1e 4 , the Active Set method with improved normal compliance produces more fluctuations which can be explained by energy dissipation during the impact, but the energy of the system is recovered and conserved after the impact. In Figure 6, we analyze the discrete energy behavior of the Active Set scheme for the improved normal compliance condition for different time steps. For ∆t varying from 1e -2 s to 1e -4 s, we can notice that the obtained numerical results are similar and this method is characterized by a dissipation of energy during the impact, but after the impact the energy of the system is conserved. For ∆t = 1e -2 s, the method generates even more fluctuations during impact. So taking ∆t below 1e -3 s does not guarantee the minimization of fluctuations.

In Figure 7, we observe the behavior of the discrete energy given by the Active Set scheme with classic normal and improved compliance with respects to the parameter α and mesh refinement related to the contact boundary, that is to say the number of contact points nbc. We notice that the numerical results obtained by using the classical normal compliance with a different number of contact points show a strong energy dissipation (between 24% and 28%). Such a phenomenon is expected in this configuration and can be explained by the discrete energy balance; we refer to [START_REF] Ayyad | Formulation and analysis of two energy-consistent methods for nonlinear elastodynamic frictional contact problems[END_REF][START_REF] Ayyad | A frictionless viscoelastodynamic contact problem with energy consistent properties: Numerical analysis and computational aspects[END_REF] for more details. However, using the Improved Normal Compliance method (α = 3) with the same contact boundary discretization, we see a remarkable improvement regarding the conservation of energy in each case after the impact and with less fluctuation during the impact. It is interesting, albeit expected, to observe to what extent the mesh refinement related to the contact boundary impacts the energy conservation properties in both cases. In Figure 8 with a time step of 1e -2 s, we notice that the Active Set scheme produces a significant energy dissipation (28%) with the standard normal compliance method. For the improved normal compliance, we notice some fluctuations, but the system recovers exactly the energy after the impact. By comparing with Figure 9 (∆t = 1e -3 s) and that of Figure 8, we notice large fluctuations for different values of α greater than 4 for the improved normal compliance method. Note that, while it seems at first sight that the classical scheme with normal compliance conserves the energy, that is not actually the case as there is a slight energy dissipation (1%). Such a phenomena is expected and can be explained by the smallness of ∆t considered.

Regarding the fluctuations, the results observed seem consistent with the theory. As a matter of fact, let us consider the following configurations:

• Let α 1 and α 2 be two values of the α coefficient used in INC, such that α 2 > α 1 ;

• Let the associated normal stress values λ 1 ν-1 2 , λ 2 ν-1 2 be given by (3.11).

As the penetrations are very small for this physical setting, for the same penetration we have λ

1 ν-1 2 > λ 2 ν-1 2
. It means that the penetrations observed with α = α 2 are to be larger than with α = α 1 (see Table 1). Therefore, the energy dissipated is expected to be larger with α = α 2 and it turns out to be the case with Figure 9, as the fluctuation amplitude is strictly increasing with respect to the value of α. Results of the Active Set method and the quasi-Lagrangian method in comparison with the number of degrees of freedom (dof), the number of contact nodes (nbc) and the total CPU time (CPU) in seconds (∆t = 0.01 s, α = 3, cν = 1e 3 ).

In Table 2, we study the convergence in CPU time of the Active Set method compared to different methods (persistent contact, classic normal compliance and the quasi-Lagrangian approach for unilateral contact), depending on the number of degrees of freedom (dof) and the number of contact points (nbc) on the boundary Γ 3 . When nbc < 64, the quasi-Lagrangian method is slightly more efficient than the Active Set method. When we increase the number of degrees of freedom (dof) we notice that the Active Set method is much faster in terms of CPU time than the quasi-Lagrangian method; this can be explained by the fact that the Active Set method does not require the use of Lagrange multipliers, and this method requires less Newton iterations; for more details we refer for instance to [START_REF] Abide | Analysis of two active set type methods to solve unilateral contact problems[END_REF][START_REF] Abide | Inexact primal-dual active set method for solving elastodynamic frictional contact problems[END_REF].

Impact of a hyperelastic ring on a foundation

The interest of this example is to give a validation of the Active Set method with the improved normal compliance conditions on a hyperelastodynamic problem. This non trivial example, introduced by Laursen [START_REF] Laursen | Computational contact and impact mechanics: fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF] concerns an academic problem of frictional impact of a hyperelastic ring against a foundation. The details on the physical setting of the problem are given below:

Ω , = (x 1 , x 2 ) ∈ R 2 : 81 ≤ (x 1 -100) 2 + (x 2 -100) 2 ≤ 100 , ∂ 0 Ω = ∅, ∂ g Ω = (x 1 , x 2 ) ∈ R 2 : (x 1 -100) 2 + (x 2 -100) 2 = 81 , ∂ c Ω = (x 1 , x 2 ) ∈ R 2 : (x 1 -100) 2 + (x 2 -100) 2 = 100 .
The domain Ω represents the cross-section of a three-dimensional deformable body under the assumption of plane stress. The ring is launched with an initial velocity towards a foundation, as shown in Figure 10. The foundation is given by (x 1 , x 2 ) ∈ R 2 : x 2 ≤ 0 . For the discretization, we use 1664 elastic nodes. For numerical experiments, the data are: The response of the compressible material, considered for Ogden's constitutive law [START_REF] Ciarlet | Sur les lois de comportement en élasticité non-linéaire compressible[END_REF] is characterized by the following energy density: In Figure 11, we observe the evolution of the total discrete energy of the frictionless dynamic system depending on selected time-integration schemes, i.e, leapfrog scheme with persistent conditions, classical normal compliance, quasi-Lagrangian method, two-step scheme [START_REF] Ayyad | Formulation and analysis of two energy-consistent methods for nonlinear elastodynamic frictional contact problems[END_REF] and improved normal compliance. We find that the leapfrog scheme and the two-step scheme (Signorini and persistent condition) conserve the energy of the system unlike the quasi-Lagrangian method with unilateral contact (9% of energy dissipation) and the classical normal compliance (9% of energy dissipation) which generates maximum error on the normal contact displacement of 1e -5 m (c ν = 1e 4 ). Due to the increment step of the leapfrog scheme, this method generates a maximum error on the normal contact displacement of 8.9e -2 m (see Table 3), therefore higher than the other methods. Moreover, the improved normal compliance method is characterized by an exact conservation of energy after the impact. This method also allows to minimize interpenetration with α = 2 (3.4e -4 m, see Table 3).

W (F) = c 1 (I 1 -3) + c 2 (I 2 -3) + d(I 3 -1) -(c 1 + 2c 2 + d) ln I 3 ,
In Figure 13, with a time step of 1e -2 s, we present the evolution of the total discrete energy of the frictionless dynamic system with the Active Set scheme with improved normal compliance conditions according to the parameter α. We observe an exact conservation of energy after the impact. However, the method generates some fluctuations during impact. Based on these observations, the choice α = 3 seems to perform at best.

In Table 4, we study the convergence in CPU time of the Active Set-INC method compared with different methods (persistent contact, standard normal compliance and quasi-Lagrangian method for unilateral contact) depending on the number of degrees of freedom (dof) and the number of contact points (nbc) on the boundary Γ 3 . For different values of nbc, the Active Set-INC method is better than the quasi-Lagrangian method in terms of CPU time; as in the elastic ball numerical example, this can be explained by the fact that the Active Set-INC method does not involve Lagrange multipliers and that it generates fewer nonlinear Newton iterations. In Figure 14, we are interested in the Active Set-INC method (α = 3) depending on the time step. For ∆t = 1e -3 s, we notice that the method is characterized by a better conservation of energy during the impact. On the other hand, for a time step ∆t ranging from 1e -1 s to 1e -2 s, we notice that the method generates some energy fluctuations during the impact, but it conserves energy after the impact.

In Figure 15, we observe the behavior of the discrete energy using an Active Set scheme with normal compliance condition with friction (µ = 0.2). This method is characterized by a slightly dissipative behavior of friction; this is consistent with the physical dissipative nature of the friction phenomenon. The difference of the dissipative behavior comes from the model of frictional contact according to the parameter α.

Conclusion and perspectives

In this work, we investigated a new energy conservation method for hyperelastodynamic contact problems, from a theoretical and numerical point of view. We first introduced the mathematical framework for general frictional contact problems in large deformations, and demonstrated how the Improved Normal Compliance (INC) condition can be applied in order to obtain energy conservation properties during impact. Next, we assessed these conservation properties in the continuous case and in its space-time discrete counterpart. Then, we derived two semi-smooth Newton Primal Dual algorithms to handle respectively the Standard Normal Compliance (SNC) and INC conditions. The following section presents numerical results on two classical representative academic test cases, namely the impact of an elastic ball and the impact of a hyperelastic ring on a foundation. The aim was to compare the performances of the INC with respect to the SNC and the classical methods used within the literature. It turned out that the methods implemented are at least just as much relevant as several other methods, regarding energy conservation properties after impact, as they display physically realistic behaviors. In the second test case, even more challenging from a numerical point of view, the energy is even almost perfectly conserved during and after the impact. Also, one of the interesting and well-known advantages when using the Active Set approach is the numerical efficiency. Indeed, in both test cases, this method is much more efficient (CPU time-wise, between twice and five times as fast for large problems) than the classical quasi-Lagrangian approach.

At this stage, we outline some future perspectives. As the results in 2D cases seem quite promising, an extension to 3D ought to be considered in order to assess the behavior of this approach on more complex cases. We can reasonably expect the performances with respect to the quasi-Lagrangian approach in 3D to be of an even greater order of magnitude than in 2D. Regarding an underlying practical use case, it could be interesting to investigate a potential application to granular media. Finally, this method was only applied for hyperelastic and elastic material so far, and provided quite successful results, which suggests to consider alternative constitutive laws. Amongst them, in particular, an extension to viscosity and plasticity would make sense as such materials involves energy dissipation phenomena.

Figure 1 :

 1 Figure 1: Reference (a) and deformed (b) configurations of a body.

4. 1 . 3 .

 13 Fixed point conditions from Newton's Semi-Smooth approach Using now the semi-smooth Newton formalism (indexed by the superscript k) at the current fixed point iterate (δ k ν , λ k ν ) of the complementary functions C λ ν , one can derive the new iterate (δ k+1 ν

4. 2 . 3 .

 23 Fixed point conditions from Newton's Semi-Smooth approach Using now the semi-smooth Newton formalism (indexed by the superscript k) at the current fixed point iterate (δ k ν , λ k ν ) of the complementary functions C λ ν , one can derive the new iterate (δ k+1 ν

.31) 4 . 3 . 3 .

 433 Fixed point conditions from Newton's Semi-Smooth approach By combining (4.24)-(4.31), with G C λ τ the generalized derivative of C λ τ , respectively, we obtain

Figure 2 :

 2 Figure 2: Discretization of the elastic ball in contact with a foundation.

Figure 3 :

 3 Figure 3: Sequence of deformed ball and contact forces before, during, and after impact.

Figure 4 :

 4 Figure 4: Discrete energy behavior of selected time integration schemes during impact (∆t = 0.001, cν = 1 1000 ).

Figure 5 :

 5 Figure 5: Discrete energy behavior of the Active Set scheme with Improved Normal Compliance (α = 3) as a function of cν (∆t = 0.01).

Figure 6 :

 6 Figure 6: Discrete energy behavior of the Active Set scheme with Improved Normal compliance (α = 3) as a function of ∆t (cν = 1e 3 ).

Figure 7 :

 7 Figure 7: Discrete energy behavior of Active Set schemes with classic Normal Compliance and with Improved Normal compliance (α = 3) as a function of the number of contact points (∆t = 0.01, cν = 1e 3 ).

Figure 8 :

 8 Figure 8: Discrete energy behavior of the Active Set scheme with Improved Normal compliance (∆t = 0.01 s) as a function of α (cν = 1e 3 ).

Figure 9 :

 9 Figure 9: Discrete energy behavior of the Active Set scheme with Improved Normal compliance (∆t = 0.001) as a function of α (cν = 1e 3 ).

  nbc

Figure 10 :

 10 Figure 10: Sequence of the deformed hyperelastic ring before, during, and after impact.

ρ

  = 1000 kg/m 3 , T = 10 s, ∆t = 0.01 s,u 0 = (0, 0) m, u 1 = (10, -10) m/s, f 0 = (0, 0) N/m 2 , f 1 = (0, 0) N/m,c 1 = 0.5 MP a, c 2 = 5.0e -3 MP a, c 3 = 5.0e -5 MP a, D = 100 MP a, c ν = 1000, a = 1.5, b = 0.5, α = 100.

with c 1

 1 = 0.5 MP a, c 2 = 0.5 • 10 -2 MP a and d = 0.35 MP a and, using the Green-Lagrange tensor defined by C = F T F, the invariants I 1 , I 2 and I 3 are defined by I 1 (C) = tr(C), I 2 (C) = (tr(C)) 2 -tr(C 2 ) 2 , I 3 (C) = det(C).

Figure 11 :

 11 Figure 11: Discrete energy behavior of selected time integration schemes during frictionless impact (∆t = 0.01, cν = 1e 3 ).

Figure 12 :

 12 Figure 12: CPU time of the Active Set method and the quasi-Lagrangian method in comparison with the number of degrees of freedom (∆t = 0.01, cν = 1e 3 ).

Figure 13 :

 13 Figure 13: Discrete energy behavior of the Active Set scheme with Improved Normal compliance (∆t = 0.01) as a function of α without friction (cν = 1e 3 ).

Figure 14 :

 14 Figure 14: Discrete energy behavior of the Active Set scheme with Improved Normal compliance (α = 3) as a function of ∆t (cν = 1e 3 ).

Figure 15 :

 15 Figure 15: Discrete energy behavior of the Active Set scheme with Improved Normal compliance (∆t = 1e -2 s) as a function of α with friction (µ = 0.2, cν = 1e 3 , cτ = 1e 3 ).

Table 2 :

 2

Table 3 :

 3 Maximum error on normal contact displacement (∆t = 1e -2 , cν = 1e 3 ).

	nbc

Table 4 :

 4 Results of the Active Set method and the quasi-Lagrangian method in comparison with the number of

degrees of freedom (dof), the number of contact nodes (nbc) and the total CPU time (CPU) in seconds (∆t = 0.01, cν = 1e 3 ).