Franz Chouly

Mathieu Fabre

Patrick Hild

Rabii Mlika

Jérôme Pousin

Yves Renard

An overview of recent results on Nitsche's method for contact problems

Keywords: contact, friction, finite elements, Nitsche's method AMS Subject Classification: 65N12, 65N15, 65N30, 74M10, 74M15, 74M20

We summarize recent achievements in applying Nitsche's method to some contact and friction problems. We recall the setting of Nitsche's method in the case of unilateral contact with Tresca friction in linear elasticity. Main results of the numerical analysis are detailed: consistency, well-posedness, fully optimal convergence in H 1 (Ω )-norm, residual-based a posteriori error estimation. Some numerics and some recent extensions to multibody contact, contact in large transformations and contact in elastodynamics are presented as well.

Introduction

For a wide range of systems in structural mechanics, it is crucial to take into account contact and friction between rigid or elastic bodies. Among numerous applications, let us mention foundations in civil engineering, metal forming processes, crash-tests of cars, design of car tires (see, e.g., [START_REF] Wriggers | Computational Contact Mechanics[END_REF]). Contact and friction conditions are usually formulated with a set of inequalities and non-linear equations on the boundary of each body, with corresponding unknowns that are displacements, velocities and surface stresses. Basically, contact conditions allow to enforce non-penetration on the whole candidate contact surface, and the actual contact surface is not known in advance. A friction law may be taken into account additionally, and various models exist that correspond to different surface properties, the most popular one being Coulomb's friction (see, e.g., [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF] and references therein).

Franz Chouly Laboratoire de Mathématiques de Besanc ¸on -UMR CNRS 6623, Université Bourgogne Franche-Comté, 16 route de Gray, 25030 Besanc ¸on Cedex, France. e-mail: franz.chouly@univ-fcomte.fr Mathieu Fabre EPFL SB MATHICSE (Bt. MA), Station 8, CH 1015 Lausanne, Switzerland. Istituto di Matematica Applicata e Tecnologie, Informatiche "E. Magenes" del CNR, via Ferrata 1, 27100 Pavia, Italy. e-mail: mathieu.fabre@epfl.ch Patrick Hild Institut de Mathématiques de Toulouse -UMR CNRS 5219, Université Paul Sabatier, 118 route de Narbonne, 31062 Toulouse Cedex 9, France. e-mail: patrick.hild@math.univ-toulouse.fr Rabii Mlika, Jérôme Pousin, Yves Renard Université de Lyon, CNRS, INSA-Lyon, ICJ UMR5208, F-69621, Villeurbanne, France. e-mail: rabii.mlika@insa-lyon.fr, jerome.pousin@insa-lyon.fr,yves.renard@insa-lyon.fr 1 Frictional contact problems can be formulated weakly within the framework of variational inequalities (see, e.g., [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF][START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF][START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]). Those are the very basis of most existing Finite Element Methods (FEM), see e.g. [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF][START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF][START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF][START_REF] Han | Quasistatic contact problems in viscoelasticity and viscoplasticity[END_REF][START_REF] Laursen | Computational contact and impact mechanics[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF][START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF]. For numerical computations with the FEM, various techniques have been devised to enforce contact and friction conditions at the discrete level, and the foremost are:

1. Penalty methods (see, e.g., [START_REF] Kikuchi | Penalty-finite-element approximation of a class of unilateral problems in linear elasticity[END_REF][START_REF] Oden | Interior penalty methods for finite element approximations of the Signorini problem in elastostatics[END_REF][START_REF] Oden | Finite element methods for constrained problems in elasticity[END_REF][START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF][START_REF] Chernov | A priori error estimates for hp penalty BEM for contact problems in elasticity[END_REF][START_REF] Chouly | On convergence of the penalty method for unilateral contact problems[END_REF]), where the set of inequations associated to contact is replaced with a non-linear inequation that approximates them. These methods remain primal, and are easy to implement. Nevertheless, consistency is lost, as a small amount of penetration, controlled by the penalty parameter, is allowed. Therefore the penalty parameter needs to be chosen with some care. Indeed, when the penalty parameter gets smaller to improve the approximation of contact conditions, the discrete problem gets stiffer and ill-conditioned, so iterative solvers such as generalized Newton may fail to converge. 2. Mixed methods (see, e.g., [START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF][START_REF] Hild | Numerical implementation of two nonconforming finite element methods for unilateral contact[END_REF][START_REF] Ben Belgacem | Hybrid finite element methods for the Signorini problem[END_REF][START_REF] Hüeber | An optimal a priori error estimate for nonlinear multibody contact problems[END_REF][START_REF] Laborde | Fixed point strategies for elastostatic frictional contact problems[END_REF][START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF]), where a Lagrange multiplier is introduced, that stands for the normal stress on the contact boundary (and for the tangential stress as well in case of frictional contact). The resulting weak form remains consistent, and characterizes the saddle-point of the corresponding Lagrangian.

Inf-sup compatibility between the primal space of displacements and the dual space of Lagrange multipliers must be satisfied at the discrete level to ensure well-posedness: see, e.g., [START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF] and references therein for different solutions that have been proposed to overcome this issue. As usual for mixed formulations, stabilized finite element methods, such as Barbosa and Hughes stabilization [START_REF] Barbosa | Circumventing the Babu ška-Brezzi condition in mixed finite element approximations of elliptic variational inequalities[END_REF], can be designed to circumvent the discrete inf-sup condition (see, e.g., [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF]).

Nitsche's method has hardly been considered to discretize contact and friction conditions, despite the fact that it has gained popularity for other boundary conditions. The Nitsche method orginally proposed in [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF][START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF] aims at treating the boundary or interface conditions in a weak sense, with appropriate consistent terms that involve only the primal variables. It differs in this aspect from standard penalization techniques which are generally non-consistent [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]. Moreover, no additional unknown (Lagrange multiplier) is needed and, therefore, no discrete inf-sup condition must be fulfilled, contrarily to mixed methods (see, e.g., [START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF][START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF]). Most of the applications of Nitsche's method during the last two decades involved linear conditions on the boundary of a domain or at the interface between subdomains: see, e.g,. [START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF] for the Dirichlet problem, [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF] for domain decomposition with non-matching meshes and [START_REF] Hansbo | Nitsche's method for interface problems in computational mechanics[END_REF] for a global review. In some recent works [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF][START_REF] Heintz | Stabilized Lagrange multiplier methods for bilateral elastic contact with friction[END_REF] it has been adapted for bilateral (persistent) contact, which still corresponds to linear boundary conditions on the contact zone. We remark furthermore that an algorithm for unilateral contact which makes use of Nitsche's method in its original form is presented and implemented in [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF], and an extension to large strain bilateral contact has been performed in [START_REF] Wriggers | A formulation for frictionless contact problems using a weak form introduced by Nitsche[END_REF].

In [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF][START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF] a new Nitsche-based FEM was proposed and analyzed for Signorini's problem, where a linear elastic body is in frictionless contact with a rigid foundation. Conversely to bilateral (persistent) contact, Signorini's problem involves non-linear boundary conditions associated to unilateral contact, with an unknown actual contact region. For this Nitsche-based FEM, optimal convergence in the H 1 (Ω )-norm of order O(h 1 2 +ν ) has been proved, provided the solution has a regularity H 3 2 +ν (Ω ), 0 < ν ≤ k -1/2 (k = 1, 2 is the polynomial degree of the Lagrange finite elements). To this purpose there is no need of additional assumption on the contact/friction zone, such as an increased regularity of the contact stress or a finite number of transition points between contact and non-contact. The proof applies in two-dimensional and three-dimensional cases, and for continuous affine and quadratic finite elements. Besides, the standard FEM for contact consists in a direct approximation of the variational inequality, with the elastic displacement as the only unknown. For this standard FEM and for many variants such as mixed/hybrid methods (e.g., [START_REF] Hild | Numerical implementation of two nonconforming finite element methods for unilateral contact[END_REF][START_REF] Ben Belgacem | Hybrid finite element methods for the Signorini problem[END_REF][START_REF] Laborde | Fixed point strategies for elastostatic frictional contact problems[END_REF]), stabilized mixed methods (e.g., [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF]), penalty methods (e.g., [START_REF] Chouly | On convergence of the penalty method for unilateral contact problems[END_REF]), it has been quite challenging to establish optimal convergence in the case the solution u belongs to H 3 2 +ν (Ω ) (0 < ν ≤ 1/2). As a matter of fact, the first fully optimal result, without extra assumptions, for the standard FEM has been achieved only recently, in 2015, see [START_REF] Drouet | Optimal convergence for discrete variational inequalities modelling Signorini contact in 2D and 3D without additional assumptions on the unknown contact set[END_REF]. The first analyses in the 1970s were indeed sub-optimal with a convergence in O(h [START_REF] Scarpini | Error estimates for the approximation of some unilateral problems[END_REF][START_REF] Haslinger | Finite element analysis for unilateral problems with obstacles on the boundary[END_REF][START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF] . In the 2000s, with additional assumptions on the finiteness of transition points between contact and non-contact optimality has been recovered (see [START_REF] Ben Belgacem | Numerical simulation of some variational inequalities arisen from unilateral contact problems by the finite element methods[END_REF] when 0 < ν < 1/2 and [START_REF] Hüeber | An optimal a priori error estimate for nonlinear multibody contact problems[END_REF] when ν = 1/2). We refer to, e.g., [START_REF] Hüeber | An optimal a priori error estimate for nonlinear multibody contact problems[END_REF][START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF][START_REF] Hild | An improved a priori error analysis for finite element approximations of Signorini's problem[END_REF][START_REF] Drouet | Optimal convergence for discrete variational inequalities modelling Signorini contact in 2D and 3D without additional assumptions on the unknown contact set[END_REF] for more detailed reviews on a priori error estimates for contact problems in elasticity. Moreover the Nitsche-based FEM encompasses symmetric and nonsymmetric variants depending upon a parameter called θ . The symmetric case of [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF] is recovered when θ = 1. When θ = 1 positivity of the contact term in the Nitsche variational formulation is generally lost. Nevertheless some other advantages are recovered, mostly from the numerical viewpoint. Namely, one of the variants (θ = 0) involves a reduced quantity of terms, which makes it easier to implement and to extend to contact problems involving non-linear elasticity. In addition, this nonsymmetric variant θ = 0 performs better in the sense it requires less Newton iterations to converge, for a wider range of the Nitsche parameter, than the variant θ = 1, see [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF]. Concerning the skew-symmetric variant θ = -1, the well-posedness of the discrete formulation and the optimal convergence are preserved irrespectively of the value of the Nitsche parameter. Note that for other boundary conditions, such as non-homogeneous Dirichlet, the symmetric variant (θ = 1) as originally proposed by Nitsche [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF] is the most widespread, since it preserves symmetry, and allows efficient solvers for linear systems with a symmetric matrix. However some nonsymmetric variants have been reconsidered recently, due to some remarkable robustness properties (see, e.g., [START_REF] Burman | A penalty-free nonsymmetric Nitsche-type method for the weak imposition of boundary conditions[END_REF][START_REF] Boiveau | A penalty-free Nitsche method for the weak imposition of boundary conditions in compressible and incompressible elasticity[END_REF]). In the context of discontinuous Galerkin methods, such nonsymmetric variants are well-known as well (see, e.g., [START_REF] Di Pietro | Mathematical aspects of discontinuous Galerkin methods[END_REF]Section 5.3.1,p.199]). Lately, various extensions of the method proposed in [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF][START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF] have been carried out:

1 2 + ν 2 )
• An extension to Coulomb's friction has been formulated in [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF] and tested numerically using a generalized Newton algorithm. • An extension to Tresca's friction has been studied in [START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF]. Optimal convergence in H 1 (Ω )-norm has been established as well, without any assumption other than usual Sobolev regularity. For the standard FEM, and other methods, technical assumptions on the contact/friction set are needed to recover optimal convergence in 2D, and in 3D, it remains an open issue [START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF]. • The case of contact in elastodynamics is dealt with in [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and timemarching schemes[END_REF][START_REF] Chouly | A Nitsche finite element method for dynamic contact: 2. Stability of the schemes and numerical experiments[END_REF]. As opposed to mixed methods, and identically as penalty and modified mass methods [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF], Nitsche's discretization yields a well-posed semi-discrete problem in space (system of Lipschitz differential equations). It can be combined with various time-marching schemes, such as the theta-scheme, Newmark or a new "hybrid" scheme. The papers [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and timemarching schemes[END_REF][START_REF] Chouly | A Nitsche finite element method for dynamic contact: 2. Stability of the schemes and numerical experiments[END_REF] present a theoretical study of well-posedness and stability of the discretized schemes, as well as numerical experiments. • The case of contact between two elastic bodies is addressed in [START_REF] Fabre | A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method[END_REF][START_REF] Chouly | An unbiased Nitsche's approximation of the frictional contact between two elastic structures[END_REF][START_REF] Mlika | An unbiased Nitsche's formulation of large deformation frictional contact and self-contact[END_REF]. In [START_REF] Fabre | A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method[END_REF] Nitsche's method is combined with a cut-FEM / fictitious domain discretization, in the small deformations framework. In [START_REF] Chouly | An unbiased Nitsche's approximation of the frictional contact between two elastic structures[END_REF][START_REF] Mlika | An unbiased Nitsche's formulation of large deformation frictional contact and self-contact[END_REF] an unbiased variant implements the contact between two elastic bodies without making any difference between master and slave contact surfaces. The contact condition is the same on each surface. This is an advantage for treatment of self-contact or multi-body contact. • Residual-based a posteriori error estimates are presented in [START_REF] Chouly | Residual-based a posteriori error estimation for contact problems approximated by Nitsche's method[END_REF]. Upper and lower bounds are proved under a saturation assumption, and the performance of the error estimates is investigated numerically. • The topic of small-sliding frictional contact on 3D interfaces is the object of [START_REF] Annavarapu | A weighted Nitsche stabilized method for small-sliding contact on frictional surfaces[END_REF], where a weighted-Nitsche method is designed, and tested numerically. • In [START_REF] Hansbo | Least-squares stabilized augmented Lagrangian multiplier method for elastic contact[END_REF] a least-square stabilized augmented lagrangian method, inspired by Nitsche's method, is described for unilateral contact. It shares some common features with Nitsche's method and allows increased flexibility on the discretization of the contact pressure. This has been followed recently by some papers [START_REF] Burman | Augmented Lagrangian finite element methods for contact problems[END_REF][START_REF] Burman | Galerkin least squares finite element method for the obstacle problem[END_REF][START_REF] Burman | Deriving robust unfitted finite element methods from augmented lagrangian formulations[END_REF] that explore further the link between Nitsche and the augmented lagrangian, for the contact problem, the obstacle problem and interface problems with adhesive contact. • A penalty-free Nitsche's method has been designed and studied in [START_REF] Burman | The Penalty Free Nitsche Method and Nonconforming Finite Elements for the Signorini Problem[END_REF] for scalar Signorini's problem, that is an extension of the method studied in [START_REF] Burman | A penalty-free nonsymmetric Nitsche-type method for the weak imposition of boundary conditions[END_REF] for the Dirichlet problem. It is combined with a non-conforming discretization based on Crouzeix-Raviart finite elements. Stability and optimal convergence rate in H 1 (Ω )-norm are established.

This paper is outlined as follows. In Section 2 we introduce a frictional contact problem and its Nitsche-based finite element approximation. The model problem that is focused on consists in unilateral contact between an elastic body and a rigid support (Signorini's problem), with Tresca's friction. We discuss as well the relationship with other methods, such as Barbosa and Hughes stabilized method, and the augmented lagrangian of Alart and Curnier. In Section 3 we summarize the main results that have been obtained regarding the consistency, well-posedness and convergence of the proposed FEM. Results of optimal convergence in H 1 (Ω )-norm are stated. Residual-based a posteriori error estimates are provided, as well as their main theoretical properties. Section 4 presents some numerical results of convergence for the H 1 (Ω )-norm of the displacement and for the contact condition, in the frictionless case, as well as some numerical results for the a posteriori error estimates. Some recent extensions are described in Section 5: contact between two elastic bodies, in small and large transformation frameworks, and contact in elastodynamics. Concluding remarks are given in Section 6.

Let us introduce some useful notations. In what follows, bold letters like u, v, indicate vector or tensor valued quantities, while the capital ones (e.g., V, K . . .) represent functional sets or operators involving vector fields. As usual, we denote by (H s (•)) d , s ∈ R, d ∈ N * the Sobolev spaces in d space dimensions (see [START_REF] Adams | Sobolev spaces[END_REF]). The standard scalar product (resp. norm) of (H s (D)) d is denoted by (•, •) s,D (resp. • s,D ) and we keep the same notation for all the values of d. We use the same notations as in [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF] for the Gâteaux derivative (or for the directional derivative) DF(v), w of a functional F at point v and in the direction w. The letter C stands for a generic constant, independent of the discretization parameters.

div σ (u) + f = 0 in Ω , σ (u) = C ε(u) in Ω , u = 0 on Γ D , σ (u)n = f N on Γ N , (1)
where σ = (σ i j ), 1 ≤ i, j ≤ d, stands for the stress tensor field and div denotes the divergence operator of tensor valued functions. The notation ε(v) = (∇v + ∇v T )/2 represents the linearized strain tensor field and C is the fourth order symmetric elasticity tensor having the usual uniform ellipticity and boundedness property. For any displacement field v and for any density of surface forces σ (v)n defined on ∂ Ω we adopt the following decomposition into normal and tangential components:

v = v n n + v t and σ (v)n = σ n (v)n + σ t (v).
The unilateral contact conditions (Kuhn-Tucker conditions) on Γ C are formulated as follows:

u n ≤ 0, σ n (u) ≤ 0, σ n (u) u n = 0. (2)
Let s ∈ L 2 (Γ C ), s ≥ 0 be a given threshold. The Tresca friction condition on Γ C reads:

   |σ t (u)| ≤ s if u t = 0, (i) σ t (u) = -s u t |u t | otherwise, (ii) (3)
where | • | stands for the euclidean norm in R d-1 . Note that conditions (3)-(i) and ( 3)-(ii) imply that |σ t (u)| ≤ s in all cases, and that if |σ t (u)| < s, we must have u t = 0.

Remark 1. The case of bilateral contact with Tresca friction can be considered too, simply substituting to equations (2) the following one on Γ C :

u n = 0. ( 4
)
The case of frictionless contact is recovered setting s = 0 in (3).

Remark 2. The conditions of (static) Coulomb friction can be written similarly as:

   |σ t (u)| ≤ -F σ n (u) if u t = 0, (i) σ t (u) = F σ n (u) u t |u t | otherwise, (ii) (5)
where F ≥ 0 is the friction coefficient. Conditions [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF] mean that, at a given contact point on Γ C , sliding can not occur while the magnitude of the tangential stress |σ t (u)| is strictly below the threshold -F σ n (u). When the threshold -F σ n (u) is reached, sliding may happen, in a direction opposite to σ t (u) (see, e.g., [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]Chapter 10]).

Remark that this static form of Coulomb friction is an adaptation of the quasi-static (or dynamic) Coulomb's law, in which the tangential velocity ut plays the same role as the displacement u t (see formulation [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF], and, e.g., [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF][START_REF] Eck | Unilateral contact problems[END_REF][START_REF] Ballard | Steady sliding frictional contact problems in linear elasticity[END_REF]). A formulation such as ( 5) is obtained for instance when the quasi-static Coulomb's law is discretized with a time-marching scheme (see, e.g., [START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF]).

Remark 3. In the Tresca friction model (3), it is assumed that the amplitude of the normal friction threshold is known (i.e., F |σ n (u)| = s, see, e.g., [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]Section 10.3]). The introduction of the Tresca friction model is motivated more by its mathematical simplicity than by physical reasons, though it can be of use in special situations. For instance:

• when a thin belt/tape is pressed against an elastic body, with a known pressure (see, e.g., [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]Chapter 10]),

• for persistent contact between solids with high intensity of contact pressures, such as it may occur in earth sciences (see, e.g., [START_REF] Pipping | Variational formulation of rate-and state-dependent friction problems[END_REF]).

Moreover the Tresca friction model can be useful for instance when Coulomb friction is approximated iteratively (see, e.g., [START_REF] Haslinger | Approximation of the Signorini problem with friction by a mixed finite element method[END_REF] and [67, Section 9, Theorem 7]).

We introduce the Hilbert space V and the convex cone K of admissible displacements which satisfy the nonpenetration condition on the contact zone Γ C :

V := v ∈ H 1 (Ω ) d : v = 0 on Γ D , K := {v ∈ V : v n = v • n ≤ 0 on Γ C } . Define a(u, v) := Ω σ (u) : ε(v) dΩ , L(v) := Ω f • v dΩ + Γ N f N • v dΓ , j(v) := Γ C s|v t | dΓ ,
for any u and v in V.

The weak formulation of Problem (1)-(3) as a variational inequality of the second kind is:

Find u ∈ K such that: a(u, v -u) + j(v) -j(u) ≥ L(v -u), ∀ v ∈ K. (6)
which admits a unique solution (see, e.g., [42, Theorem 5.1, Remark 5.2, p.69]). Moreover this solution is the unique minimizer on K of the functional

J : V v → 1 2 a(v, v) -L(v) + j(v) ∈ R. (7)
Remark 4. In the case of bilateral contact (condition (4) instead of (2)), the same weak formulation (6) holds, replacing the convex cone K by the space:

V b := {v ∈ V : v n = 0 on Γ C } .

The Nitsche-based finite element method

Let V h ⊂ V be a family of finite dimensional vector spaces (see [START_REF] Ciarlet | The finite element method for elliptic problems, Handbook of Numerical Analysis[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF][START_REF] Brenner | The mathematical theory of finite element methods[END_REF]) indexed by h coming from a family T h of triangulations of the domain Ω (h = max T ∈T h h T where h T is the diameter of T ). We suppose that the family of triangulations is regular, i.e., there exists σ > 0 such that ∀T ∈ T h , h T /ρ T ≤ σ where ρ T denotes the radius of the inscribed ball in T . Furthermore we suppose that this family is conformal to the subdivision of the boundary into Γ D , Γ N and Γ C (i.e., a face of an element T ∈ T h is not allowed to have simultaneous non-empty intersection with more than one part of the subdivision). We choose a standard Lagrange finite element method of degree k with k = 1 or k = 2, i.e.:

V h := v h ∈ (C 0 (Ω )) d : v h | T ∈ (P k (T )) d , ∀T ∈ T h , v h = 0 on Γ D . (8)
However, the analysis would be similar for any C 0 -conforming finite element method. We make use of the notation [•] R -, that stands for the projection onto R

-([x] R -= 1 2 (x -|x|) for x ∈ R).
Moreover, for any α ∈ R + , we introduce the notation [•] α for the orthogonal projection onto B(0, α) ⊂ R d-1 , where B(0, α) is the closed ball centered at the origin 0 and of radius α. This operation can be defined analytically, for x ∈ R d-1 by:

[x] α = x if |x| ≤ α, α x |x| otherwise.
The notation H(•) will stand for a "Heaviside" function: for any x ∈ R,

H(x) :=    1 if x > 0, 1 2 if x = 0, 0 if x < 0.
We adopt the convention H(0) = 1/2 to allow the property H(x) + H(-x) = 1, ∀x ∈ R. Moreover we will make use of the property

H(-x)[x] R -= [x] R -, ∀x ∈ R. (9)
The next properties of projections are useful in the mathematical analysis of the method:

(y -x)([y] R --[x] R -) ≥ ([y] R --[x] R -) 2 ∀x, y ∈ R, (10)
(y -x) • ([y] α -[x] α ) ≥ |[y] α -[x] α | 2 ∀x, y ∈ R d-1 , ( 11
)
where • is the euclidean scalar product in R d-1 .

The next result has been pointed out earlier in [START_REF] Alart | A generalized Newton method for contact problems with friction[END_REF] (see as well [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF][START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF] for detailed formal proofs).

Proposition 1. Let γ be a positive function defined on Γ C . The contact with Tresca friction conditions (2)-(3) can be reformulated as follows:

σ n (u) = [σ n (u) -γ u n ] R -, (12)
σ t (u) = [σ t (u) -γ u t ] s . (13)
Remark 5. Equation ( 12) is an example of nonlinear complementarity (NCP) function that allows to reformulate complementarity conditions such as expressed in (2) using a single nonlinear relationship (see, e.g., [START_REF] Galántai | Properties and construction of NCP functions[END_REF] and references therein). This is not the unique possible formulation, but is among the simplest ones.

Let now θ ∈ R be a fixed parameter and γ a positive function defined on Γ C . Let u be the solution of the unilateral contact problem in its strong form (1)-( 2)-(3). We assume that u is sufficiently regular so that all the following calculations make sense. From the Green formula and equation (1), we get for every v ∈ V:

a(u, v) - Γ C σ n (u) v n dΓ - Γ C σ t (u) • v t dΓ = L(v).
Note that

v n = - 1 γ (θ σ n (v) -γv n ) + θ γ σ n (v), v t = - 1 γ (θ σ t (v) -γv t ) + θ γ σ t (v).
So:

a(u, v) - Γ C θ γ σ n (u) σ n (v) dΓ + Γ C 1 γ σ n (u) (θ σ n (v) -γv n ) dΓ - Γ C θ γ σ t (u) • σ t (v) dΓ + Γ C 1 γ σ t (u) • (θ σ t (v) -γv t ) dΓ = L(v).
Finally, using conditions ( 12) and ( 13), we obtain:

a(u, v) - Γ C θ γ σ n (u) σ n (v) dΓ + Γ C 1 γ [σ n (u) -γ u n ] R -(θ σ n (v) -γv n ) dΓ - Γ C θ γ σ t (u) • σ t (v) dΓ + Γ C 1 γ [σ t (u) -γ u t ] s • (θ σ t (v) -γv t ) dΓ = L(v). (14)
Formula ( 14) is the starting point of the Nitsche-based formulation. We remark that it may have no sense at the continuous level if u lacks of regularity (the only assumption u ∈ V is not sufficient to justify the above calculations). Nevertheless we consider in what follows that γ is a positive piecewise constant function on the contact and friction interface Γ C which satisfies

γ| T ∩Γ C = γ 0 h T , (15)
for every T that has a non-empty intersection of dimension d -1 with Γ C , and where γ 0 is a positive given constant (the Nitsche parameter). Note that the value of γ on element intersections has no influence. Let us introduce the discrete linear operators P n θ ,γ :

V h → L 2 (Γ C ) v h → θ σ n (v h ) -γv h n and P t θ ,γ : V h → (L 2 (Γ C )) d-1 v h → θ σ t (v h ) -γv h t .
Define as well the bilinear form:

A θ γ (u h , v h ) := a(u h , v h ) - Γ C θ γ σ (u h )n • σ (v h )n dΓ .
The Nitsche-based method for unilateral contact with Tresca friction then reads:

   Find u h ∈ V h such that: A θ γ (u h , v h ) + Γ C 1 γ [P n 1,γ (u h )] R -P n θ ,γ (v h ) dΓ + Γ C 1 γ P t 1,γ (u h ) s • P t θ ,γ (v h ) dΓ = L(v h ), ∀ v h ∈ V h . ( 16
)
Note that we adopted in this presentation a different convention for notations compared to previous works [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF][START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF][START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF]. This is in order to get closer to the formulations provided in most of the papers on Nitsche's method and on the augmented lagrangian method (see Section 2.4). Furthermore, and as already stated in [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF] the parameter θ can be set to some particular values, namely:

1. for θ = 1 we recover a symmetric method for which the contact term

Γ C 1 γ P n 1,γ (u h ) R - P n 1,γ (v h ) dΓ
is positive when we set v h = u h . This method can be derived from an energy functional (see Section 2.3). The tangent matrix resulting from a linearization with generalized Newton is symmetric. 2. for θ = 0 we recover a simple method close to penalty and to the augmented lagrangian. It involves only a few terms and is of easiest implementation. 3. for θ = -1 the skew-symmetric method admits one unique solution and converges optimally irrespectively of the value of the Nitsche parameter γ 0 > 0.

Remark 6. The operators P n θ ,γ (•) and P t θ ,γ (•) have no special physical meaning but can be related to proximal augmented lagrangian formulations (see 2.4.3 for more details).

Remark 7. For frictionless unilateral contact (s = 0 in (3)) we can follow the same path as previously for the derivation of the method. With the solution u to (1)-( 2)-( 3) and the condition σ t (u) = 0 there holds first for every v ∈ V:

a(u, v) - Γ C σ n (u) v n dΓ = L(v),
and then

a(u, v) - Γ C θ γ σ n (u) σ n (v) dΓ + Γ C 1 γ σ n (u) (θ σ n (v) -γv n ) dΓ = L(v).
Using condition [START_REF] Brenner | The mathematical theory of finite element methods[END_REF] we obtain:

a(u, v) - Γ C θ γ σ n (u) σ n (v) dΓ + Γ C 1 γ [σ n (u) -γ u n ] R -(θ σ n (v) -γv n ) dΓ = L(v).
As a result, the counterpart of ( 16) reads:

   Find u h ∈ V h such that: A n θ γ (u h , v h ) + Γ C 1 γ [P n 1,γ (u h )] R -P n θ ,γ (v h ) dΓ = L(v h ), ∀ v h ∈ V h , (17)
where

A n θ γ (u h , v h ) := a(u h , v h ) - Γ C θ γ σ n (u h ) σ n (v h ) dΓ .
Remark that setting s = 0 in formulation [START_REF] Burman | Deriving robust unfitted finite element methods from augmented lagrangian formulations[END_REF] does not lead to formulation [START_REF] Burman | Augmented Lagrangian finite element methods for contact problems[END_REF] because of the difference between

A θ γ (•, •) and A n θ γ (•, •). Remark 8.
For bilateral contact with friction (equations ( 1)-( 3)-( 4)), the Nitsche-based formulation reads:

   Find u h ∈ V h b such that: A t θ γ (u h , v h ) + Γ C 1 γ P t 1,γ (u h ) s • P t θ ,γ (v h ) dΓ = L(v h ), ∀ v h ∈ V h b , (18)
where

A t θ γ (u h , v h ) := a(u h , v h ) - Γ C θ γ σ t (u h ) • σ t (v h ) dΓ and V h b := V h ∩ V b .
Remark 9. Following the same path as in Proposition 1 the Coulomb friction conditions (5) can be reformulated as:

σ t (u) = [σ t (u) -γu t ] (-F σ n (u)) = [σ t (u) -γu t ] (-F [σ n (u)-γu n ] R -)
. This motivates the introduction of the following Nitsche-based formulation for unilateral contact with Coulomb friction (equations ( 1)-( 2)-( 5)):

               Find u h ∈ V h such that: A θ γ (u h , v h ) + Γ C 1 γ [P n 1,γ (u h )] R -P n θ ,γ (v h ) dΓ + Γ C 1 γ P t 1,γ (u h ) -F P n 1,γ (u h ) R - • P t θ ,γ (v h ) dΓ = L(v h ), ∀ v h ∈ V h . (19)
We next define convenient mesh-dependent norms, in fact weighted

L 2 (Γ C )-norm (since (γ 0 /γ)| T = h T ). Definition 1. For any v ∈ L 2 (Γ C ), we set v -1 2 ,h,Γ C := (γ 0 /γ) 1 2 v 0,Γ C , v 1 2 ,h,Γ C := (γ/γ 0 ) 1 2 v 0,Γ C .
The same definitions extend straightforwardly to functions in

(L 2 (Γ C )) d-1 .
Additionally, it will be sometimes convenient to endow V h with the following mesh-and parameter-dependent scalar product:

Definition 2. For all v h , w h ∈ V h we set (v h , w h ) γ := (v h , w h ) 1,Ω + (γ 1 2 v h , γ 1 2 w h ) 0,Γ C ,
and note

• γ := (•, •) 1 2
γ the corresponding norm. Remark that the two norms • γ and • 1,Ω are equivalent on V h , in the following sense (for a quasi-uniform mesh T h ):

v h 1,Ω ≤ v h γ ≤ 1 +C γ 0 h v h 1,Ω ,
for any v h ∈ V h . The positive constant C comes from the trace inequality and the constant of quasi-uniformity of the mesh T h . For a mesh T h that is not quasi-uniform, the same relationship holds, replacing h by (min T ∈T h h T ).

We end this section with the following statement: a discrete trace inequality (see, e.g., [START_REF] Thomée | Galerkin finite element methods for parabolic problems[END_REF][START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF]), that is a key ingredient for the whole mathematical analysis of Nitsche's based methods.

Lemma 1. There exists C > 0, independent of the parameter γ 0 and of the mesh size h, such that, for all v h ∈ V h :

σ n (v h ) -1 2 ,h,Γ C + σ t (v h ) -1 2 ,h,Γ C ≤ C v h 1,Ω . (20)

Energy minimization for the symmetric variant

We show in this section that an energy functional can be associated to Problem [START_REF] Burman | Deriving robust unfitted finite element methods from augmented lagrangian formulations[END_REF] in the symmetric case (θ = 1), that is a discrete counterpart of J (•). Using Riesz's representation theorem, we identify (V h , (•, •) γ ) to its dual.

Let us first introduce a functional for the total potential energy, i.e. the strain energy and the potential energy of the external forces:

J E (v h ) := 1 2 a(v h , v h ) -L(v h ),
for any v h ∈ V h . For the contact condition [START_REF] Brenner | The mathematical theory of finite element methods[END_REF] we add the term:

J n γ (v h ) := 1 2 Γ C 1 γ [P n 1,γ (v h )] 2 R -dΓ - 1 2 Γ C 1 γ σ n (v h ) 2 dΓ .
And for the Tresca friction condition (13) we take:

J t γ (v h ) := - 1 2 Γ C 1 γ P t 1,γ (v h ) -P t 1,γ (v h ) s 2 dΓ + 1 2 Γ C 1 γ P t 1,γ (v h ) 2 dΓ - 1 2 Γ C 1 γ σ t (v h ) 2 dΓ .
The energy functional associated to Problem ( 16) is then:

J γ := J E + J n γ + J t γ .
Now, when θ = 1, we are able to characterize Problem ( 16) as the first-order optimality condition associated to the minimization of J γ (•) on V h : Proposition 2. Suppose that θ = 1 and that γ 0 is large enough. Then:

1. (J γ + L)(•) is non-negative. 2. J γ (•) is a Gâteaux-differentiable and convex functional on V h .
3. Any u h that minimizes J γ (•) on V h is solution to Problem [START_REF] Burman | Deriving robust unfitted finite element methods from augmented lagrangian formulations[END_REF], that can be written equivalently

DJ γ (u h ), v h = 0, ∀v h ∈ V h .
Proof. First, provided a large enough γ 0 , (J γ (v h ) + L(v h )) is a non-negative quantity due to the ellipticity of the elasticity tensor A, the discrete trace inequality [START_REF] Chernov | A priori error estimates for hp penalty BEM for contact problems in elasticity[END_REF] and the relationship

P t 1,γ (v h ) -P t 1,γ (v h ) s ≤ P t 1,γ (v h
) , this latter being a property of the projection onto a closed ball. Let us rewrite in a slightly different form the potential J γ (•):

J γ := JE + J n γ + J t γ , with, for v h ∈ V h : JE (v h ) := 1 2 A γ (v h , v h ) -L(v h ), J n γ (v h ) := J n γ (v h ) + 1 2 Γ C 1 γ σ n (v h ) 2 dΓ , J t γ (v h ) := J t γ (v h ) + 1 2 Γ C 1 γ σ t (v h ) 2 dΓ .
The potential JE (•) is Gâteaux-differentiable on V h and its derivative is:

D JE (v h ), w h = A γ (v h , w h ) -L(w h ), (21)
for all v h , w h ∈ V h .

Similarly we check that the potential J n γ (•) is Gâteaux-differentiable too on V h . Its derivative is obtained as

D J n γ (v h ), w h = Γ C 1 γ [P n 1,γ (v h )] R -D([P n 1,γ (v h )] R -), w h dΓ = Γ C 1 γ [P n 1,γ (v h )] R -H(-P n 1,γ (v h ))P n 1,γ (w h ) dΓ = Γ C 1 γ [P n 1,γ (v h )] R -P n 1,γ (w h ) dΓ (22)
where we used the property ( 9). For the last potential J t γ (•) let us consider the functional

J : R d-1 x → 1 2 |x -[x] s | 2 ∈ R.
After simple calculations we check that J(•) is Gâteaux-differentiable and that:

DJ(x), y = (x -[x] s ) • y, for all x, y ∈ R d-1 .
Using the above formula we obtain that J t γ (•) is Gâteaux-differentiable, with derivative

D J t γ (v h ), w h = - Γ C 1 γ P t 1,γ (v h ) -P t 1,γ (v h ) s • P t 1,γ (w h ) dΓ + Γ C 1 γ P t 1,γ (v h ) • P t 1,γ (w h ) dΓ .
This simplifies further into:

D J t γ (v h ), w h = Γ C 1 γ P t 1,γ (v h ) s • P t 1,γ (w h ) dΓ . ( 23
)
The convexity of JE (•) (resp. J n γ (•) and J t γ (•)) results from the ellipticity of A and the inequality (20) (resp. inequalities [START_REF] Ben Belgacem | Hybrid finite element methods for the Signorini problem[END_REF] and ( 11)) combined with the characterization of Gâteaux-differentiable convex functions that can be found, in, e.g., [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]Theorem 3.3,Chapter 3]. This ends the proof of the second point in the theorem. To prove the last Point 3, we apply, e.g., [64, Theorem 3.7, (v), Chapter 3], with the expression of DJ γ (•) that is the sum of the expressions ( 21)-( 22)-( 23).

Remark 10. Note that a same result can be obtained for the two-body contact problem discretized with an unbiased Nitsche method, see [START_REF] Chouly | An unbiased Nitsche's approximation of the frictional contact between two elastic structures[END_REF] for details.

Relationship with other methods

We discuss the relationship with other methods here: Nitsche for Dirichlet boundary conditions, Barbosa and Hughes stabilized method and the augmented lagrangian of Alart and Curnier.

Nitsche for contact and Nitsche for Dirichlet

We consider the case of frictionless contact (s = 0) and Nitsche's formulation [START_REF] Burman | Augmented Lagrangian finite element methods for contact problems[END_REF] in this situation. Let us split the contact boundary Γ C into two portions:

• Γ - C := {x ∈ Γ C : P n 1,γ (u h )(x) < 0}, • Γ + C := {x ∈ Γ C : P n 1,γ (u h )(x) ≥ 0}(= Γ C \Γ - C ).
In this case, we can rewrite formally the Nitsche-based method [START_REF] Burman | Augmented Lagrangian finite element methods for contact problems[END_REF]:

     Find u h ∈ V h such that: A n θ γ (u h , v h ) + Γ - C 1 γ P n 1,γ (u h )P n θ ,γ (v h ) dΓ = L(v h ), ∀ v h ∈ V h .
Note that this only a formal writing since in fact the splitting of Γ C into Γ + C and Γ - C is an unknown, that depends on u h . Then using the detailed expression of A n θ γ (•, •) and P n 1,γ (•), and after re-ordering of the terms we get:

a(u h , v h ) - Γ + C θ γ σ n (u h ) σ n (v h ) - Γ - C σ n (u h ) v h n dΓ -θ Γ - C u h n σ n (v h ) dΓ + Γ - C γ u h n v h n dΓ = L(v h ), ∀ v h ∈ V h .
We recognize on Γ - C Nitsche's method for imposition of Dirichlet boundary conditions on the normal component of the displacement u n (see, e.g., [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF][START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF]). It results that Γ - C can be viewed as a discrete approximation of the actual contact surface. On Γ + C a free Neumann boundary condition is imposed weakly, in the same fashion as in [START_REF] Juntunen | Nitsche's method for general boundary conditions[END_REF][START_REF] Lüthen | An improved a priori error analysis of Nitsche's method for Robin boundary conditions[END_REF]. Therefore Γ + C may represent a discrete approximation of the unsticked contact surface.

Link with Barbosa & Hughes stabilization

We still consider the case of frictionless contact (s = 0) and Nitsche's formulation [START_REF] Burman | Augmented Lagrangian finite element methods for contact problems[END_REF]. Let us introduce

L 2 -(Γ C ) := {µ ∈ L 2 (Γ C ) | µ ≤ 0 a.e. on Γ C }
and the discrete multiplier

λ h := [P n 1,γ (u h )] R -.
Following the same steps as in the symmetric case θ = 1 (see [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF]Section 2.3] for details as well as [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF] when θ ∈ R), we can rewrite the formulation (17) into an equivalent mixed form:

           Find (u h , λ h ) ∈ V h × L 2 -(Γ C ) such that: a(u h , v h ) - Γ C λ h v h n dΓ + Γ C θ γ -1 (λ h -σ n (u h )) σ n (v h ) dΓ = L(v h ), ∀ v h ∈ V h , Γ C (µ -λ h )u h n dΓ + Γ C γ -1 (µ -λ h )(λ h -σ n (u h )) dΓ ≥ 0, ∀ µ ∈ L 2 -(Γ C ).
The first line is simply expression [START_REF] Burman | Augmented Lagrangian finite element methods for contact problems[END_REF] recasted after introduction of the multiplier λ h , and the second line means that λ h is the projection of P n 1,γ (u h ) onto L 2 -(Γ C ). Note as well that the inverse of Nitsche parameter γ -1 0 can be interpreted as a stabilization parameter. We recover indeed a mixed form close to the stabilized method [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF], the only difference being that in [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF], the dual set L 2 -(Γ C ) is approximated by using finite elements on the contact boundary. The stabilized method of [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF] is an adaptation to unilateral contact of Barbosa & Hughes stabilization [START_REF] Barbosa | Circumventing the Babu ška-Brezzi condition in mixed finite element approximations of elliptic variational inequalities[END_REF]. The method of stabilized Lagrange multiplier at the boundary proposed by Barbosa & Hughes [START_REF] Barbosa | Circumventing the Babu ška-Brezzi condition in mixed finite element approximations of elliptic variational inequalities[END_REF] originates from a stream of works dedicated to the use of a penalization technique for recovering coercivity for the Lagrange multiplier in order to avoid handling the Babuska-Brezzi condition in the finite element context. At the beginning the proposed formulation was inconsistent [START_REF] Oden | Penalty-finite-element methods for constrained problems in elasticity[END_REF], then supplementary terms were added for ensuring consistency [START_REF] Barbosa | Circumventing the Babu ška-Brezzi condition in mixed finite element approximations of elliptic variational inequalities[END_REF]. This method of stabilized Lagrange multiplier has been adapted for the unilateral contact problem in the frictionless case [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF]. Optimal error estimates for the Lagrange multiplier have been obtained provided an extra regularity result for the Lagrange multiplier is satisfied, which in certain circumstances is not relevant, see remark 3.7 in [START_REF] Hild | A stabilized Lagrange multiplier method for the finite element approximation of contact problems in elastostatics[END_REF] (note that using the new results published in [START_REF] Drouet | Optimal convergence for discrete variational inequalities modelling Signorini contact in 2D and 3D without additional assumptions on the unknown contact set[END_REF] this analysis can be improved now). In the seminal paper [START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF] a simplified formulation of Barbosa & Hughes (where just the essential added terms are considered) has been proved equivalent to a Nitsche formulation for a Laplace problem with Dirichlet boundary conditions (for domain decomposition, see also the recent paper of Juntunen [START_REF] Juntunen | On the connection between the stabilized Lagrange multiplier and Nitsche's methods[END_REF]). In this context the Lagrange multiplier belongs to L 2 (Γ C ) and is approximated with discontinuous finite elements. Therefore this Lagrange multiplier can be eliminated and a Nitsche formulation is recovered.

Proximal augmented lagrangian and Nitsche

A popular formulation for solving contact problems is the augmented lagrangian (see, e.g., [START_REF] Alart | A generalized Newton method for contact problems with friction[END_REF][START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF][START_REF] Renard | Approximation des conditions de contact en élastostatique : méthode de Nitsche et lien avec le lagrangien augmenté proximal[END_REF]). In the case of unilateral contact problem without friction, its expression is:

L r (u h , λ H ) := 1 2 a(u h , u h ) -L(u h ) + Γ C 1 2r λ H -ru h n 2 R - -(λ H ) 2 dΓ , (24)
for u h ∈ V h and where the discrete multiplier λ H belongs to a finite element space W H of functions defined on the contact boundary Γ C . The multiplier λ H is a new unknown that approximates the normal stress σ n (u). We introduced as well r > 0, which is the augmentation parameter. We note that d dx

1 2 [x] 2 R - = [x] R -.
Therefore, using the notation

∂ ∂ u (•), • (resp. ∂ ∂ λ (•),
• ) for the directional derivative according to the first variable u (resp. to the second variable λ ):

∂ ∂ u Γ C 1 2r λ H -ru h n 2 R - dΓ , v h = Γ C 1 r λ H -ru h n R - (-rv h n ) dΓ = - Γ C λ H -ru h n R - v h n dΓ . Similarly ∂ ∂ λ Γ C 1 2r λ H -ru h n 2 R - dΓ , µ H = Γ C 1 r λ H -ru h n R - µ H dΓ .
Using the above expressions, let us write explicitly the optimality system associated to the augmented lagrangian [START_REF] Chouly | On convergence of the penalty method for unilateral contact problems[END_REF]:

0 = ∂ L r ∂ u (u h , λ H ), v H = a(u h , v h ) -L(v h ) - Γ C λ H -ru h n R - v h n dΓ , ∀v h ∈ V h , 0 = ∂ L r ∂ λ (u h , λ H ), µ H = Γ C 1 r λ H -ru h n R - -λ H µ H dΓ , ∀µ H ∈ W H .
This is an unconstrained formulation, that is more appropriate for numerical solving. Now, remark that the second equation is a way to enforce weakly, at the discrete level, the condition [START_REF] Brenner | The mathematical theory of finite element methods[END_REF]. Another way, straightforward, to enforce this condition, is to substitute σ n (u h ) to λ H in the first equation of the optimality system (and we forget about the second equation):

a(u h , v h ) - Γ C σ n (u h ) -ru h n R - v h n dΓ = L(v h ), ∀v h ∈ V h .
We now recognize Nitsche's method [START_REF] Burman | Augmented Lagrangian finite element methods for contact problems[END_REF] for θ = 0. Moreover the parameter γ can be identified with the augmentation parameter r. Note that some recent works are filling the gap between Nitsche and augmented lagrangian formulations in the case of contact and obstacle problems [START_REF] Burman | Augmented Lagrangian finite element methods for contact problems[END_REF][START_REF] Burman | Galerkin least squares finite element method for the obstacle problem[END_REF][START_REF] Hansbo | Least-squares stabilized augmented Lagrangian multiplier method for elastic contact[END_REF].

Analysis of the Nitsche-based method

This section sums up the main results about the numerical analysis of the Nitsche-based formulation [START_REF] Burman | Deriving robust unfitted finite element methods from augmented lagrangian formulations[END_REF]. For detailed proofs, the reader is refered to [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF][START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF][START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF][START_REF] Chouly | Residual-based a posteriori error estimation for contact problems approximated by Nitsche's method[END_REF]. First of all we recall the consistency of the method, that is a direct consequence of reformulations ( 12)-( 13) followed by integration-by-parts:

Lemma 2. The Nitsche-based method (16) is consistent: suppose that the solution u to (1)-(3) belongs to (H 3 2 +ν (Ω )) d , with ν > 0, then u is also solution to

A θ γ (u, v h ) + Γ C 1 γ P n 1,γ (u) R - P n θ ,γ (v h ) dΓ + Γ C 1 γ P t 1,γ (u) s • P t θ ,γ (v h ) dΓ = L(v h ), (25)
for any v h ∈ V h .

Note that for the same reasons the formulation (17) (resp. ( 18)) for frictionless unilateral contact (resp. for frictional bilateral contact) is consistent too.

Well-posedness

To show that Problem ( 16) is well-posed we use an argument by Brezis for M-type and pseudo-monotone operators [START_REF] Brezis | Équations et inéquations non linéaires dans les espaces vectoriels en dualité[END_REF] (see also [START_REF] Lions | Quelques méthodes de résolution des problèmes aux limites non linéaires[END_REF] and [START_REF] Kikuchi | Penalty-finite-element approximation of a class of unilateral problems in linear elasticity[END_REF]). We define a (non-linear) operator B h s : V h → V h , by using the Riesz representation theorem and by means of the formula:

(B h s v h , w h ) γ := A θ γ (v h , w h ) + Γ C 1 γ [P n 1,γ (v h )] R -P n θ ,γ (w h ) dΓ + Γ C 1 γ P t 1,γ (v h ) s • P t θ ,γ (w h ) dΓ . ( 26
)
Note that Problem ( 16) is well-posed if and only if B h s is one-to-one. The following result characterizes wellposedness:

Theorem 1. The operator B h s is hemicontinuous. Moreover there exist C,C > 0 such that, for all v h , w h ∈ V h :

(B h s v h -B h s w h , v h -w h ) γ ≥C 1 - C(1 + θ ) 2 2γ 0 v h -w h 2 γ . ( 27
)
As a result, when the condition below holds

γ 0 ≥ C(1 + θ ) 2 , (28)
Problem ( 16) admits one unique solution u h in V h .

Remark 11. In the symmetric case θ = 1, we remark that:

DJ γ (v h ), w h = (B h s v h , w h ) γ -L(w h ), for v h , w h ∈ V h , in other terms B h s is the gradient of (J γ + L)(•).
In this case the equation ( 27) means that J γ (•) is strongly convex under the condition [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF]. As a result, when θ = 1, well-posedness can alternatively be established using a minimization argument, such as [64, Theorem 3.4, Chapter 3], and the unique solution to ( 16) is also the unique minimizer of J γ (•) on V h .

A priori error estimates in H 1 (Ω )-norm

First we recall the abstract error estimate.

u -u h 1,Ω + γ -1 2 0   σ n (u)-P n 1,γ (u h ) R - -1 2 ,h,Γ C + σ t (u)-P t 1,γ (u h ) s -1 2 ,h,Γ C   ≤ C inf v h ∈V h u -v h 1,Ω + γ 1 2 0 u -v h 1 2 ,h,Γ C + γ -1 2 0 σ (u -v h )n -1 2 ,h,Γ C , ( 29
)
where C is a positive constant, independent of h, u and γ 0 .

2. Set θ = -1. Then for all values of γ 0 > 0, the solution u h to Problem (16) satisfies the abstract error estimate [START_REF] Chouly | An unbiased Nitsche's approximation of the frictional contact between two elastic structures[END_REF] where C is a positive constant, dependent of γ 0 but independent of h and u.

The optimal convergence of the method is stated below.

Theorem 3. Suppose that the solution u to Problem (6) belongs to (H

3 2 +ν (Ω )) d with 0 < ν ≤ k -1 2 (k = 1, 2
is the degree of the finite element method, given in (8)) and d = 2, 3. When θ = -1, suppose in addition that the parameter γ 0 is sufficiently large. The solution u h to Problem (16) satisfies the following error estimate:

u -u h 1,Ω + σ n (u)-P n 1,γ (u h ) R - -1 2 ,h,Γ C + σ t (u)-P t 1,γ (u h ) s -1 2 ,h,Γ C ≤ Ch 1 2 +ν u 3 2 +ν,Ω , (30)
where C is a positive constant, independent of h and u.

We can easily obtain the following error estimate on the Cauchy constraint σ (uu h )n

-1 2 ,h,Γ C
in the weighted

L 2 (Γ C )-norm (note that σ n (u h ) = P n 1,γ (u h ) R - and σ t (u h ) = P t 1,γ (u h )
s on Γ C conversely to the continuous case).

Corollary 1. Suppose that the solution u to Problem (6) belongs to (H 3 2 +ν (Ω )) d with 0 < ν ≤ k -1 2 and d = 2, 3. When θ = -1, suppose in addition that the parameter γ 0 is sufficiently large. The solution u h to Problem (16) satisfies the following error estimate:

σ (u -u h )n -1 2 ,h,Γ C ≤ Ch 1 2 +ν u 3 2 +ν,Ω , (31)
where C is a positive constant, independent of h and u.

Residual-based a posteriori error estimate

An explicit residual-based a posteriori error estimate can be derived for Problem [START_REF] Burman | Deriving robust unfitted finite element methods from augmented lagrangian formulations[END_REF], that is an extension of the one presented in [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF] (see also, e.g., [START_REF] Verfürth | A review of a posteriori error estimation techniques for elasticity problems[END_REF] for linear elasticity). We introduce standard notations for this purpose:

• We define E h the set of edges/faces of the triangulation and define E int h := {E ∈ E h : E ⊂ Ω } as the set of interior edges/faces of T h . We denote by E N h := {E ∈ E h : E ⊂ Γ N } the set of Neumann edges/faces and similarly

E C h := {E ∈ E h : E ⊂ Γ C }
is the set of contact edges/faces. • For an element T , we denote by E T the set of edges/faces of T and according to the above notation, we set

E int T := E T ∩ E int h , E N T := E T ∩ E N h , E C T := E T ∩ E C h .
• For an edge/face E of an element T , introduce ν T,E the unit outward normal vector to T along E. Furthermore, for each edge/face E, we fix one of the two normal vectors and denote it by ν E . The jump of some vector valued function v across an edge/face E ∈ E int h at a point y ∈ E is defined as

v E (y) := lim α→0 + v(y + αν E ) -v(y -αν E ).
• Let ω T be the union of all elements having a nonempty intersection with T . Similarly for a node x and an edge/face E, let ω x := ∪ T :x∈T T and ω E := ∪ x∈E ω x . • f T (resp. f N,E ) is a computable quantity that approximates f on the element T ∈ T h (resp. f N on the edge E ∈ E N h ). The a posteriori error estimator is defined below. Definition 3. The local error estimators η T and the the global estimator η are defined by

η T := 4 ∑ i=1 η 2 iT 1/2 , η 1T := h T div σ (u h ) + f T 0,T , η 2T := h 1/2 T   ∑ E∈E int T ∪E N T J E,n (u h ) 2 0,E   1/2 , η 3T := h 1/2 T   ∑ E∈E C T σ t (u h )-P t 1,γ (u h ) s 2 0,E   1/2 , η 4T := h 1/2 T   ∑ E∈E C T σ n (u h )-P n 1,γ (u h ) R - 2 0,E   1/2 , η := ∑ T ∈T h η 2 T 1/2 ,
where J E,n (u h ) means the constraint jump of u h in the normal direction, i.e.,

J E,n (u h ) := σ (u h )ν E E , ∀E ∈ E int h , σ (u h )ν E -f N,E , ∀E ∈ E N h . (32)
The local and global approximation terms are given by

ζ T :=   h 2 T ∑ T ⊂ω T f -f T 2 0,T + h E ∑ E⊂E N T f N -f N,E 2 0,E   1/2 , ζ := ∑ T ∈T h ζ 2 T 1/2 .
We need for the analysis a "saturation" assumption as in [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF] for Nitsche-based domain decomposition, and as in [START_REF] Wohlmuth | A residual based error estimator for mortar finite element discretizations[END_REF] for mortar methods.

Assumption 4

The solution u to (1)-( 2)-( 3) and the discrete solution u h to [START_REF] Burman | Deriving robust unfitted finite element methods from augmented lagrangian formulations[END_REF] are such that:

σ n (u -u h ) -1/2,h,Γ C + σ t (u -u h ) -1/2,h,Γ C ≤ C u -u h 1,Ω , (33)
where C is a positive constant independent of h.

The following statement guarantees the reliability of the a posteriori error estimator:

Theorem 5. Let u be the solution to (1)-( 2)-(3), with u ∈ (H 3 2 +ν (Ω )) d (ν > 0 and d = 2, 3), and let u h be the solution to the corresponding discrete problem [START_REF] Burman | Deriving robust unfitted finite element methods from augmented lagrangian formulations[END_REF]. Assume that, for θ = -1, γ 0 is sufficiently large. Assume that the saturation assumption (33) holds as well. Then we have

u -u h 1,Ω + σ n (u)-P n 1,γ (u h ) R - -1 2 ,h,Γ C + σ t (u)-P t 1,γ (u h ) s -1 2 ,h,Γ C + σ n (u) -σ n (u h ) -1/2,h,Γ C + σ t (u) -σ t (u h ) -1/2,h,Γ C ≤ C(1 + γ -1 0 )(η + ζ ),
where the positive constant C is independent of h and γ 0 .

The last result concerns the local lower error bounds of the discretization error terms: Theorem 6. For all elements T ∈ T h , the following local lower error bounds hold:

η 1T ≤ C u -u h 1,T + ζ T , (34)
η 2T ≤ C u -u h 1,ω T + ζ T . ( 35
)
For all elements T such that T ∩ E C h = / 0, the following local lower error bounds hold:

η 3T ≤ C ∑ E∈E C T h 1/2 T σ t (u)-P t 1,γ (u h ) s 0,E + σ t (u -u h ) 0,E , (36)
η 4T ≤ C ∑ E∈E C T h 1/2 T σ n (u)-P n 1,γ (u h ) R -0,E + σ n (u -u h ) 0,E , (37)
where the positive constant C is independent of h and γ 0 .

Remark 12. From Theorem 6, optimal convergence rates of order O(h min(k, 1 2 +ν) ) are expected for the estimator of Definition 3.

Numerical experiments

This section is devoted to numerical results that illustrate the theoretical analysis and the practical interest of the method, in the frictionless case. First, in 4.1 we provide practical details concerning the implementation. Then in 4.2 we assess numerically the a priori error estimates of Subsection 3.2. Numerical assessment of the a posteriori error estimator described in 3.3 and adaptive computations are object of Subsection 4.3.

Some implementation issues

The discrete contact problem is solved by using a generalized Newton method, which means that Problem ( 16) is derived with respect to u h to obtain the tangent system. The term "generalized Newton's method" comes from the fact that the operators such as [•] R -and [•] s are not Gâteaux-differentiable at some specific points. However, no special treatment is considered. If a point of non-differentiability is encountered, the tangent system corresponding to one of the two alternatives (x < 0 or x > 0 for [•] R -) is chosen arbitrarily. Integrals of the non-linear term on Γ C are computed with standard quadrature formulas. Note that, for frictionless contact, the situation where the solution is non-differentiable at an integration point is very rare and corresponds to what is called a "grazing contact" (both u n = 0 and σ n = 0). In [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF] one can find further details and references on generalized Newton's method, and especially a numerical study of its convergence when applied to contact problems discretized by Nitsche's method (for variants θ = 1 and θ = 0) as well as other methods. The finite element library Getfem++1 has been used for all the computations presented in this paper.

Numerics for Hertz's contact

The numerical results obtained in [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF] for frictionless Hertz's contact problems of a disk/sphere with a plane rigid foundation are summarized here. This slightly exceeds the scope defined in Section 2 since a non-zero initial gap between the elastic solid and the rigid foundation is considered in the computations. Moreover, the tests are performed with P 1 and isoparametric P 2 Lagrange finite elements on meshes which are approximations of the real domain. The numerical situation in two-dimensions is represented in Fig. 2. A disc of radius 20 cm is considered with a contact boundary Γ C which is restricted to the lower part (y < 20 cm) of the boundary. A homogeneous Neumann condition is applied on the remaining part of the boundary. Since no Dirichlet condition is considered, the problem is not fully coercive. To overcome the non-definiteness coming from the free rigid motions, the horizontal displacement is prescribed to be zero on the two points of coordinates (0 cm, 10 cm) and (0 cm, 30 cm) which blocks the horizontal translation and the rigid rotation. Homogeneous isotropic linear elasticity in plane strain approximation is considered with a Young modulus fixed at E = 25 MPa and a Poisson ratio P = 0.25. A vertical density of volume forces of 20 MN/m 3 is applied. The solution for mesh sizes h = 0.5 cm, 1 cm, 3 cm, 4.5 cm and h = 10 cm are compared with a reference solution on a very fine mesh (h = 0.15 cm) using quadratic isoparametric finite elements. Moreover, the reference solution is computed with a different discretization of the contact problem (Lagrange multipliers and Alart-Curnier augmented lagrangian, see [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF]). In complement to the relative error in H 1 (Ω )-norm we compute also the following relative error in L 2 (Γ C )-norm:

γ -1 2 P n 1,γ (u h ) R - -P n 1,γ (u h re f ) R - 0,Γ C P n 1,γ (u h re f ) R -0,Γ C
, where u h is the discrete solution and u h re f the reference solution. Note that

P n 1,γ (u h ) R -
is an approximation of the contact stress with a convergence of order 1 (see Theorem 3). The convergence rates computed thanks to this test are reported in Table 1 and Table 2, for P 1 and P 2 finite elements, respectively, and for different values of θ and γ 0 . 1 Computed convergence rates for Hertz's problem in 2D and P 1 finite elements. The column "H 1 (Ω )" stands for the relative error in H 1 (Ω )-norm on the displacement, the column "L 2 (Γ C )" stands for the relative error in L 2 (Γ C ) on the contact condition. 45 1 0.14 0.44 0.01 1.00 0.58 2 Computed convergence rates for Hertz's problem in 2D and P 2 finite elements. The column "H 1 (Ω )" stands for the relative error in H 1 (Ω )-norm on the displacement, the column "L 2 (Γ C )" stands for the relative error in L 2 (Γ C ) on the contact condition.

γ 0 /E H 1 (Ω ) L 2 (Γ C ) 100
γ 0 /E H 1 (Ω ) L 2 (Γ C ) 100 1.62 1.
γ 0 /E H 1 (Ω ) L 2 (Γ C ) 100 1.
For the 3D Hertz's problem and P 1 finite elements, convergence rates are reported in Table 3.

γ 0 /E H 1 (Ω ) L 2 (Γ C ) 100 1.62 1.12 1
-0.21 0.81 0.01 -0.47 0.51

γ 0 /E H 1 (Ω ) L 2 (Γ C ) 100 1.59 1.10 1
1.00 1.94 0.01 0.40 1.25 3 Computed convergence rates for Hertz's problem in 3D and P 1 finite elements. The column "H 1 (Ω )" stands for the relative error in H 1 (Ω )-norm on the displacement, the column "L 2 (Γ C )" stands for the relative error in L 2 (Γ C ) on the contact condition.

γ 0 /E H 1 (Ω ) L 2 (Γ C ) 100 1.56 1.07 1 1.44 1.85 0.01 1.42 1.41 θ = 1 θ = 0 θ = -1 Table
When θ = 1, optimal convergence is obtained for both H 1 (Ω ) and weighted L 2 (Γ C )-norms of the error, but only for the largest value of the parameter γ 0 (γ 0 = 100 E). This corroborates the theoretical result of Theorem 3 for which the optimal rate of convergence is obtained for a sufficiently large γ 0 . When θ = 0, for the smallest value of γ 0 the convergence remains sub-optimal. However, for the intermediate value of γ 0 (γ 0 = E) the optimal convergence is reached. Concerning the version with θ = -1, which corresponds to an unconditionally coercive problem, one can see that optimal convergence is reached for all values of γ 0 . Remark that the same conclusions hold both for 2D and 3D cases, and linear and quadratic finite elements, though the difference between the variants θ = 1 and the others can be greater for quadratic elements or in the 3D case. A strategy to guarantee an optimal convergence is of course to consider a sufficiently large γ 0 . However, the price to pay is an ill-conditioned discrete problem. The study presented in [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF] for the versions θ = 1 and θ = 0 shows that Newton's method has important difficulties to converge when γ 0 is large. When symmetry is not required, a better strategy seems to consider the version with θ = -1 or an intermediate value of θ = 0 which ensure both a optimal convergence rate and few iterations of Newton's method to converge.

A posteriori error estimation

We report the test case taken from [START_REF] Hild | Residual a posteriori error estimators for contact problems in elasticity[END_REF] (see also [START_REF] Hild | Residual error estimators for Coulomb friction[END_REF][START_REF] Lleras | Modélisation, analyse et simulation de problèmes de contact en mécanique des solides et des fluides[END_REF] in the frictional case). We consider the domain Ω = (0, 1) × (0, 1) with material characteristics E = 10 6 and P = 0.3. A homogeneous Dirichlet condition on Γ D = {0} × (0, 1) is prescribed to clamp the body. The body is potentially in contact on Γ C = {1} × (0, 1) with a rigid obstacle and Γ N = (0, 1) × ({0} ∪ {1}) is the location of a homogeneous Neumann condition. The body Ω is acted on by a vertical volume density of force f = (0, f 2 ) with f 2 = -76518 such that there is coexistence of a slip zone and a separation zone with a transition point between both zones. For error computations, since we do not have a closed-form solution, a reference solution is computed with Lagrange P 2 elements, h = 1/160, γ 0 = E and θ = -1. First of all we illustrate in Fig. 3 the difference between uniform and adaptive refinement. For the latter we refine only the mesh elements T in which the local estimator η T is below a given threshold s = 2.5 × 10 -3 . The minimal (respectively maximal) size of the adaptive mesh is equal to 1/160 (respectively h = 1/40). As expected the rate of convergence with respect to the number of degrees of freedom is far better in the case of adaptive refinement than with uniform refinement. The solution obtained with adaptive refinement and θ = -1 is depicted in Fig. 4. We observe that the error is concentrated at both left corners (transition between Dirichlet and Neumann conditions) and near the transition point between contact and separation.

A detailed numerical convergence study of the error estimator η for this test-case, as well as for 2D and 3D Hertz contact, is provided in [START_REF] Chouly | Residual-based a posteriori error estimation for contact problems approximated by Nitsche's method[END_REF].

Recent extensions

This section is about recent extensions. We describe numerical approximations of the contact between two elastic bodies, first with the master/slave (biased) formulation in 5.1, where the meshes are conformal to the geometry, and in 5.2 with cut-finite elements and fictitious domains. An unbiased formulation is provided in 5.3. Finally the extension to elastodynamics is given in 5.4.

Contact between two elastic bodies

Now, we consider two elastic bodies Ω 1 and Ω 2 expected to come into contact. To simplify notations, a general index i is used to represent indifferently each body (i = 1, 2). We denote by Γ i C a portion of the boundary of the body Ω i which is a candidate contact surface with an outward unit normal vector n i . For the contact surfaces, let us assume a sufficiently smooth one to one application (projection for instance) mapping each point of the first contact surface to a point of the second one (see also Fig. 5):

Π 1 : Γ 1 C → Γ 2 C .
Let J 1 be the Jacobian of the transformation Π 1 and J 2 = 1 J 1 the Jacobian of Π 2 = (Π 1 ) -1 . We suppose in the following that J 1 > 0. We define on each contact surface another "normal" vector ñi such that:

. Ω 2 n 1

x 1 = Π 2 (x 2 )

x 2 = Π 1 (x 1 ) Γ 2 C n 2 Γ 1 C Ω 1 .
Fig. 5 Two bodies Ω 1 and Ω 2 in contact: mapping

Π 1 from Γ 1 C to Γ 2 C and its inverse mapping Π 2 . ñi (x) =      Π i (x) -x Π i (x) -x if x = Π i (x)
,

n i if x = Π i (x).
Note that ñ1 = -ñ2 • Π 1 and ñ2 = -ñ1 • Π 2 . For any displacement field v i and for any density of surface forces σ i (v i )n i defined on ∂ Ω i , we adopt the following notation:

v i = v i n ñi + v i t and σ i (v i )n i = σ i n (v i ) ñi + σ i t (v i ).
The two-body contact problem, in the linear elastic framework, consists in finding the displacement field u = (u 1 , u 2 ) verifying the equations ( 38) and the contact conditions described hereafter:

div σ i (u i ) + f i = 0 in Ω i , (38a)
σ i (u i ) = C i ε(u i ) in Ω i , (38b)
u i = 0 on Γ i D , (38c)
σ i (u i )n i = f i N on Γ i N , (38d)
We consider a zero initial normal gap to simplify the notations (see [START_REF] Fabre | A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method[END_REF] for a nonzero one) and we define the relative normal displacements

u 1 n = (u 1 -u 2 • Π 1 ) • ñ1 and u 2 n = (u 2 -u 1 • Π 2 ) • ñ2 .
The classical master/slave (biased) formulation is obtained by selecting for instance Γ 1

C to be the slave surface and Γ 2 C to be the master one. Then the unilateral contact condition is written on the slave side Γ 1 C :

u 1 n ≤ 0, σ 1 n (u 1 ) ≤ 0, σ 1 n (u 1 ) u 1 n = 0. (39)
Let

s 1 ∈ L 2 (Γ 1 C ), s 1 ≥ 0, u 1 t = u 1 t -u 2 t • Π 1 . The Tresca friction condition on Γ 1 C reads:        σ 1 t (u 1 ) ≤ s 1 if u 1 t = 0, σ 1 t (u 1 ) = -s 1 u 1 t u 1 t otherwise. ( 40
)
As in Section 2, we reformulate the contact and friction conditions (39)-( 40) as follows:

σ 1 n (u 1 ) = σ 1 n (u 1 ) -γ 1 u 1 n R -, (41)
σ 1 t (u 1 ) = σ 1 t (u 1 ) -γ 1 u 1 t s 1 , (42)
where γ 1 is the counterpart of γ on the slave surface Γ 1 C . We note T h,i a triangulation of the domain Ω i and introduce the finite element spaces as in Section 2:

V h := V h,1 × V h,2 , with V h,i := v h,i ∈ C 0 (Ω i ) : v h,i |T ∈ (P k (T )) d , ∀T ∈ T h,i , v h,i = 0 on Γ i D .
We use the Green formula and equations (38a)-(38d) as well as the Nitsche's writing of the contact and friction conditions ( 41)-( 42) to get the following biased finite element approximation:

           Find u h ∈ V h such that, A 1 θ γ (u h , v h ) + Γ 1 C 1 γ 1 P 1,n 1,γ 1 (u h ) R - P 1,n θ ,γ 1 (v h ) dΓ + Γ 1 C 1 γ 1 P 1,t 1,γ 1 (u h ) s 1 • P 1,t θ ,γ 1 (v h ) dΓ = L(v h ), ∀v h ∈ V h , (43)
where:

P i,n θ ,γ i (v) := θ σ i n (v i ) -γ i v i n , , P i,t θ ,γ i (v) := θ σ i t (v i ) -γ i v i t ,
and

A 1 θ γ (u, v) := 2 ∑ i=1 Ω i σ (u i ) : ε(v i ) dΩ - Γ 1 C θ γ 1 σ 1 (u 1 )n • σ 1 (v 1 )n dΓ , L(v) := 2 ∑ i=1 Ω i f i • v i dΩ + Γ i N f i N • v i dΓ .
All the mathematical properties presented in Section 3 for the deformable/rigid case can be transposed to this biased deformable/deformable Nitsche formulation (see [START_REF] Fabre | A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method[END_REF]).

Contact between two elastic bodies with a fictitious domain

A complete analysis in a frictionless case is presented in [START_REF] Fabre | A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method[END_REF] in the framework of a fictitious domain discretization using cut-elements. A fictitious domain Ω that contains both Ω 1 and Ω 2 is considered, generally with a simple geometry such that a structured mesh T h of Ω can be used (see the example in Fig. 6). Then a single finite element space

W h := v h ∈ C 0 (Ω ) : v h,i |T ∈ (P k (T )) d , ∀T ∈ T h ,
is used to approximate the displacement of the two bodies, in the sense that we consider the following approximations space

V h := W h | Ω 1 × W h | Ω 2 . . Ω Γ 1 C Γ 2 C Γ 1 D Γ 1 N Γ 2 D Ω 2 Γ 1 N Ω 1 .
Fig. 6 Two bodies Ω 1 and Ω 2 in contact: a single mesh of the fictitious domain Ω .

in the discrete problem [START_REF] Glowinski | Augmented Lagrangian and operator-splitting methods in nonlinear mechanics[END_REF], where W h | Ω i is the space of restrictions to Ω i of functions of W h . Additionnaly, since the boundaries of the bodies are independent of the mesh element edges, the Dirichlet conditions are also prescribed with Nitsche's method. This yields the following approximation (in the frictionless case):

             Find u h ∈ V h such that, Ā1 θ γ (u h , v h ) + Γ 1 C 1 γ 1 P1,n 1,γ 1 (u h ) R - P1,n θ ,γ 1 (v h ) dΓ + ∑ i=1,2 Γ i D 1 γ 1 u h,i • v h,i -R ρ (u h,i ) • v h,i -θ u h,i • R ρ (v h,i ) dΓ = L(v h ), ∀v h ∈ V h , (44)
where now Pi,n θ ,γ i (v

) := θ R ρ (v i ) -γ i v i n , and
Ā1 θ γ (u, v) := 2 ∑ i=1 Ω i σ (u i ) : ε(v i ) dΩ - (Γ 1 C ∪Γ 1 D ∪Γ 2 D ) θ γ 1 R ρ (u 1 ) • R ρ (v 1 ) dΓ . The new term R ρ (v h,i
) is an approximation of σ i (v i )n i built in order to recover an optimal order of convergence.

In [START_REF] Fabre | A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method[END_REF] this operator is defined as the polynomial extension of the displacement field of neighbour elements having a sufficiently large intersection with the domain Ω i , as proposed initially in [START_REF] Haslinger | A new fictitious domain approach inspired by the extended finite element method[END_REF]. Note that an alternative is the use of the so-called ghost penalty introduced in [START_REF] Burman | Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche's method[END_REF]. Remark finally that the presentation here differs slightly from [START_REF] Fabre | A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method[END_REF] in which the second Newton's law is used to reformulate differently the contact conditions [START_REF] Fabre | A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method[END_REF]. This second Newton's law reveals as well to be a key ingredient for the derivation of an unbiased formulation, as detailed in next section.

Unbiased formulation for self-and multi-body contact

If the master/slave formulation consists in a natural extension of the contact treatment between a deformable body and a rigid ground, it has no complete theoretical justification, and induces detection difficulties in the case of self-contact and multi-body contact. We provide in this section a short description of the unbiased formulation presented in [START_REF] Chouly | An unbiased Nitsche's approximation of the frictional contact between two elastic structures[END_REF] that circumvents these difficulties. We do not distinguish between a master surface and a slave one since we impose the non-penetration and the friction conditions on both of them. Unbiased contact and friction formulations have been considered before in [START_REF] Sauer | An unbiased computational contact formulation for 3D friction[END_REF] and references therein.

In order to obtain an unbiased method we prescribe the contact condition on the two surfaces in a symmetric way. Thus, the conditions describing contact on Γ i C (i = 1, 2) are:

u i n ≤ 0, σ i n (u i ) ≤ 0, σ i n (u i ) u i n = 0. ( 45
)
Let

s i ∈ L 2 (Γ i C ), s i ≥ 0, u 1 t = u 1 t -u 2 t • Π 1 and u 2 t = u 2 t -u 1 t • Π 2 = -u 1 t • Π 2 . The Tresca friction condition on Γ 1 C and Γ 2 C reads:        σ i t (u i ) ≤ s i if u i t = 0, σ i t (u i ) = -s i u i t u i t otherwise. (46)
Finally, we need to consider the second Newton law between the two bodies:

       γ 1 C σ 1 n (u 1 )ds - γ 2 C σ 2 n (u 2 )ds = 0, γ 1 C σ 1 t (u 1 )ds + γ 2 C σ 2 t (u 2 )ds = 0,
where

γ 1 C is any subset of Γ 1 C and γ 2 C = Π 1 (γ 1 C
). Mapping all terms on γ 1 C allows writing:

σ 1 n (u 1 ) -J 1 σ 2 n (u 2 • Π 1 ) = 0, σ 1 t (u 1 ) + J 1 σ 2 t (u 2 • Π 1 ) = 0, on Γ 1 C . (47)
Remark 13. : A similar condition holds on Γ 2 C :

σ 2 n (u 2 ) -J 2 σ 1 n (u 1 • Π 2 ) = 0, σ 2 t (u 2 ) + J 2 σ 1 t (u 1 • Π 2 ) = 0.
Let us mention that, due to second Newton law, we need to fix s 1 and s 2 such that:

-s 1 u 1 t u 1 t = σ 1 t (u 1 ) = -J 1 σ 2 t (u 2 • Π 1 ) = J 1 s 2 u 2 t • Π 1 u 2 t • Π 1 = -J 1 s 2 u 1 t | u 1 t | .
It results that the following compatibility condition on s 1 and s 2 needs to be satisfied:

s 1 = J 1 s 2 . ( 48
)
The reformulation of the contact and friction conditions ( 45)-( 46) reads:

σ i n (u i ) = σ i n (u i ) -γ i u i n R -, (49)
σ i t (u i ) = σ i t (u i ) -γ i u i t s i , (50)
where γ i plays the same role as γ on each contact surface Γ i C . Still using the Green formula and the different equations considered as well as the Nitsche's writing of the contact and friction we obtain the following unbiased finite element approximation:

                 Find u h ∈ V h such that, A 1,2 θ γ (u h , v h ) + 1 2 ∑ i=1,2 Γ i C 1 γ i P i,n 1,γ i (u h ) R - P i,n θ ,γ i (v h ) dΓ + 1 2 ∑ i=1,2 Γ i C 1 γ i P i,t 1,γ i (u h ) s i • P i,t θ ,γ i (v h ) dΓ = L(v h ), ∀v h ∈ V h , (51)
where this time:

A 1,2 θ γ (u, v) := 2 ∑ i=1 Ω i σ (u i ) : ε(v i ) dΩ - 1 2 Γ i C θ γ i σ i (u i )n • σ i (v i )n dΓ .
Note that in the above method the two contact surfaces are treated similarly. This small strain formulation allows an extended mathematical study which is impossible to perform in a large transformations framework. Indeed, all the mathematical analysis carried out for the unilateral contact problem can be adapted for this method, and yields the same theoretical properties of well-posedness and convergence. This analysis as well as a complete numerical study in the small strain framework can be found in [START_REF] Chouly | An unbiased Nitsche's approximation of the frictional contact between two elastic structures[END_REF]. However, the true potential of this unbiased paradigm is to able the design of a similar method in large transformations. We quickly describe it in the rest of this section (see [START_REF] Mlika | An unbiased Nitsche's formulation of large deformation frictional contact and self-contact[END_REF] for more details). Let Ω ⊂ R d be a potential multi-body domain representing the reference configuration of one or several hyper-elastic bodies. Let Γ C be the part of the boundary of Ω where contact may occur. Let J H (•) be the potential of the hyper-elastic law (more details can be found in [START_REF] Poulios | An unconstrained integral approximation of large sliding frictional contact between deformable solids[END_REF]). Then, to each (material) point x ∈ Γ C and to the displacement u : Ω → R d , we associate g(x, u) the (nonnegative) gap distance (in the deformed configuration) between the surface Γ C at point x and the nearest other surface in potential contact. The main difference with the small strain case is that the point in potential contact at x ∈ Γ C is a priori unknown. Moreover it changes during the deformation and may be quite difficult to define. This gap calculation can be performed in different ways, see [START_REF] Poulios | An unconstrained integral approximation of large sliding frictional contact between deformable solids[END_REF] and the references therein for more details.

Then, an extension of unbiased Nitsche's method to large transformation contact in the frictionless case can be written

DJ H (u h ), v h + 1 2 Γ C σ n (u h ) + γg(u h ) R - Dg(u h ), v h dΓ = 0, (52)
for all admissible increment of displacement v h . The notation σ n (u) stands for the contact stress in large transformations.

Note that, similarly to the method presented in the small strain case, the integral on Γ C in (52) being calculated on the whole boundary Γ C , it is notably calculated on both the two corresponding surfaces of a potential contact whether a multi-body contact is considered or a self-contact. Nitsche's method allows thus to approximate both multi-body-and self-contact in the same simple formalism. An example of application of this method to the contact and self-contact of two deformable tubes is presented in Fig. 7.

Contact in elastodynamics

This section shows an adaptation of Nitsche's method to elastodynamic frictionless contact problems. In this context Nitsche-based approximation produces well-posed space semi-discretizations contrary to standard finite element discretizations. We recall the main results of the analysis for the space semi-discretization in terms of well-posedness and energy conservation. We consider some time-marching schemes (theta-scheme, Newmark and a new hybrid scheme) for which the well-posedness and the stability results are given. The details can be found in the papers [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and timemarching schemes[END_REF] and [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 2. Stability of the schemes and numerical experiments[END_REF]. This sections ends with an extension of the method to the frictional case which is actually under investigation in [START_REF] Chouly | Nitsche method for dynamic contact with Coulomb friction[END_REF].

The setting is the same as in section 2, but in the evolution case: an elastic body Ω in R d with d = 2, 3 of density ρ (constant to simplify). In the reference configuration, the body is in contact on Γ C with a rigid foundation and we suppose that the unknown contact zone during deformation is included into Γ C . We consider the unilateral contact problem in linear elastodynamics during a time interval [0, T ) where T > 0 is the final time. We denote by Ω T := (0, T ) × Ω the time-space domain, and similarly Γ DT := (0, T ) × Γ D , Γ NT := (0, T ) × Γ N and Γ CT := (0, T ) × Γ C . The problem consists in finding the displacement field u : [0, T ) × Ω → R d verifying the equations and conditions ( 53)-( 54):

div σ (u) + f = ρ ü, σ (u) = C ε(u) in Ω T , u = 0 on Γ DT , σ (u)n = f N on Γ NT , u(0, •) = u 0 , u(0, •) = u0 in Ω , ( 53
)
where u is the velocity of the elastic body and ü its acceleration; u 0 is the initial displacement and u0 is the initial velocity.

The conditions describing unilateral contact without friction on Γ CT are:

u n ≤ 0, σ n (u) ≤ 0, σ n (u) u n = 0, σ t (u) = 0. ( 54
)      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 . (57)
The following theorem together with the boundedness of (M h ) -1 γ (see [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and timemarching schemes[END_REF]) show that Problem (56) (or equivalently Problem ( 57)) is well-posed.

Theorem 7. 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 γ ≤ C(1 + γ -1 0 )(1 + |θ |) v h 1 -v h 2 γ . ( 58
)
As a consequence, for every value of θ ∈ R and γ 0 > 0, Problem (56) admits one unique solution

u h ∈ C 2 ([0, T ], V h ).
Remark 15. Note that, conversely to the static case (see [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF][START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF][START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF]) and the fully-discrete case there is no condition on γ 0 for the space (semi-)discretization, which remains well-posed even if γ 0 is arbitrarily small. Now we consider the energy estimates which are counterparts of the equation ( 55), in the semi-discretized case.

Let us define 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 ].
which is associated to the solution u h (t) to Problem [START_REF] Hild | Residual error estimators for Coulomb friction[END_REF]. Note that this is the direct transposition of the mechanical energy E(t) for the continuous system. Set also

E h θ (t) := E h (t) - θ 2γ 0 σ n (u h (t)) 2 -1 2 ,h,Γ C -[P n 1,γ (u h (t))] R - 2 -1 2 ,h,Γ C := E h (t) -θ R h (t),
that corresponds to a modified energy in which a consistent term is added. This term denoted R h (t) represents, roughly speaking, the nonfulfillment of the contact condition (12) by u h .

Theorem 8. Suppose that the system associated to ( 53)-( 54) is conservative, i.e., that L(t) ≡ 0 for all t ∈ [0, T ].

The solution u h to (56) then satisfies the following identity:

d dt E h θ (t) = (1 -θ ) Γ C 1 γ [P n 1,γ (u h (t))] R - uh n (t) dΓ .
Notably, when θ = 1, we get for any t ∈ [0, T ]: E h 1 (t) = E h 1 (0).

Fully discrete formulations

Now we fully discretize the dynamic contact problem by combining Nitsche's method with some classical schemes (theta-scheme, Newmark) as well as a new hybrid scheme. We focus on the well-posedness and the stability of the schemes.

Let τ > 0 be the time-step, and consider a uniform discretization of the time interval [0, T ]: (t 0 , . . . ,t N ), with t n = nτ, n = 0, . . . , N. Let θ ∈ [0, 1], we use the notation:

x h,n+ θ = (1 -θ )x h,n + θ x h,n+1
for arbitrary quantities x h,n , x h,n+1 ∈ V h . Hereafter we denote by u h,n (resp. uh,n and üh,n ) the resulting discretized displacement (resp. velocity and acceleration) at time-step t n . We next define the following energy:

E h,n := 1 2 ρ uh,n 2 0,Ω + 1 2 a(u h,n , u h,n ),
which is associated with the solution u h,n to Problems ( 59), ( 61) or [START_REF] Juntunen | Nitsche's method for general boundary conditions[END_REF]. Set also

E h,n θ := E h,n - θ 2γ 0 σ n (u h,n ) 2 -1 2 ,h,Γ C -[P n 1,γ (u h,n )] R - 2 -1 2 ,h,Γ C := E h,n -θ R h,n .
Note that the energies E h,n and E h,n θ are the fully discrete counterparts of the semi-discrete energies E h (t) and E h θ (t). • Theta-scheme. We discretize in time Problem (56) using a theta-scheme, of parameter θ ∈ [0, 1]. For n ≥ 0, the fully discretized problem reads:

               Find u h,n+1 , uh,n+1 , üh,n+1 ∈ V h such that: u h,n+1 = u h,n + τ uh,n+ θ , uh,n+1 = uh,n + τ üh,n+ θ , (ρ üh,n+1 , v h ) 0,Ω + A n θ γ (u h,n+1 , v h ) + Γ C 1 γ [P n 1,γ (u h,n+1 )] R -P n θ ,γ (v h ) dΓ = L n+1 (v h ), ∀ v h ∈ V h , (59)
with initial conditions u h,0 = u h 0 , uh,0 = uh 0 , üh,0 = üh 0 and where L n+1 (•) = L(t n+1 )(•). The following proposition concerns well-posedness and stability of the theta-scheme: Proposition 3. 1. If θ = 0, existence and uniqueness of (59) always holds since the scheme is fully explicit.

2. Let θ > 0. If (1 + θ ) 2 ≤ Cγ 0 1 + ρh 2 τ 2 θ 2
where C is a positive constant, then at each time-step n, Problem (59) admits one unique solution. 3. Suppose that L n ≡ 0 for all n ≥ 0. Then, for γ 0 sufficiently large, θ = 1 and θ = 1 (backward Euler scheme), the following stability estimate holds for the solution to Problem (59), for all n ≥ 0:

E h,n+1 1 ≤ E h,n 1 . (60)
Consequently, the scheme (59) is unconditionally stable when θ = 1, θ = 1 (i.e., stable for all h > 0 and all τ > 0).

• Newmark scheme. We discretize in time Problem (56) using a Newmark scheme of parameters β ∈ [0, 1/2] and γ ∈ [0, 1]. For n ≥ 0, the fully discretized problem reads:

                 Find u h,n+1 , uh,n+1 , üh,n+1 ∈ V h such that: u h,n+1 = u h,n + τ uh,n + τ 2 2 üh,n+2 β , uh,n+1 = uh,n + τ üh,n+ γ , (ρ üh,n+1 , v h ) 0,Ω + A n θ γ (u h,n+1 , v h ) + Γ C 1 γ [P n 1,γ (u h,n+1 )] R -P n θ ,γ (v h ) dΓ = L n+1 (v h ), ∀ v h ∈ V h , (61)
with initial conditions u h,0 = u h 0 , uh,0 = uh 0 , üh,0 = üh 0 . The following proposition holds: Proposition 4. 1. If β = 0, existence and uniqueness of (61) always holds since the scheme is fully explicit.

2. Let β > 0. If (1 + θ ) 2 ≤ Cγ 0 1 + ρh 2 τ 2 β
where C is a positive constant, then at each time-step n, Problem (61) admits one unique solution. 3. Suppose that L n ≡ 0 for all n ≥ 0. Then, for γ 0 sufficiently large, θ = 1, γ = 1, and β = 1/2, the following stability estimate holds for the solution to Problem (61) for all n ≥ 0:

E h,n+1 1 ≤ E h,n 1 .
So the scheme (61) is unconditionally stable when θ = 1, γ = 1 and β = 1/2 (i.e., stable for all h > 0 and all τ > 0).

• A new time-marching scheme. We followed the ideas from [START_REF] Gonzalez | Exact energy and momentum conserving algorithms for general models in nonlinear elasticity[END_REF] in the context of incompressible nonlinear elasticity and from [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF] for penalized contact, where the authors design modified time-marching schemes to enforce energy conservation. The idea is to propose an hybrid discretization of the Nitsche-based contact term: the linear part of Problem ( 56) is treated with a conservative Crank-Nicolson scheme, whereas the non-linear part arising from contact is discretized with a linear combination of Crank-Nicolson and Midpoint schemes. This strategy is of interest since the resulting scheme is unconditionally stable in the symmetric case (for more properties, see [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and timemarching schemes[END_REF][START_REF] Chouly | A Nitsche finite element method for dynamic contact: 2. Stability of the schemes and numerical experiments[END_REF]).

For n ≥ 0, the fully discretized problem reads:

                 Find u h,n+1 , uh,n+1 , üh,n+1 ∈ V h such that: u h,n+1 = u h,n + τ uh,n+ 1 2 , uh,n+1 = uh,n + τ üh,n+ 1 2 , (ρ üh,n+ 1 2 , v h ) 0,Ω + A n θ γ (u h,n+ 1 2 , v h ) + Γ C 1 γ Φ(u h,n , u h,n+1 )P n θ ,γ (v h ) dΓ = L n+ 1 2 (v h ), ∀ v h ∈ V h , (62)
with the initial conditions u h,0 = u h 0 , uh,0 = uh 0 , üh,0 = üh 0 and with the following expression for Φ(u h,n , u h,n+1 ):

Φ(u h,n , u h,n+1 ) := H(-P n 1,γ (u h,n ))[P n 1,γ (u h,n+ 1 2 )] R -+ H(P n 1,γ (u h,n ))[P n 1,γ (u h )] n+ 1 2 R -. Remark that [P n 1,γ (u h )] n+ 1 2 R -= 1/2([P n 1,γ (u h,n )] R -+ [P n 1,γ (u h,n+1
)] R -) represents the Crank-Nicolson part, whereas [P n 1,γ (u h,n+ 1 2 )] R -= [1/2(P n 1,γ (u h,n ) + P n 1,γ (u h,n+1 ))] R -stands for the Midpoint part. So, when P n 1,γ (u h,n ) < 0, the Midpoint scheme is applied, and when P n 1,γ (u h,n ) > 0, the Crank-Nicolson scheme is applied instead. When P n 1,γ (u h,n ) = 0 both schemes coincide.

The well-posedness and stability properties of the fully discrete scheme (62) are stated below.

Proposition 5. 1. If the condition below is satisfied

(1 + θ ) 2 ≤ Cγ 0 1 + ρh 2 τ 2
where C is a positive constant, then at each time-step n, Problem (62) admits one unique solution.

2. Suppose that L n ≡ 0 for all n ≥ 0 and that Problem (62) is well-posed. Suppose also that θ = 1. The following stability estimate holds for the solution to Problem (62) for all n ≥ 0:

E h,n+1 1 ≤ E h,n 1 .
So the scheme (62) is unconditionally stable when θ = 1 (i.e., stable for all h > 0 and all τ > 0).

Remark 16. From Proposition 5 we observe that the scheme ( 62) is not exactly conservative. Nevertheless the energy dissipation caused by the scheme is almost negligible in practice, and reduces when the discretization parameters h, τ are taken smaller.

We next show an extension of the method to the Coulomb frictional case where F ≥ 0 still denotes the friction coefficient (the case F = 0 corresponds to the former frictionless case). Now the condition σ t (u) = 0 on Γ CT in [START_REF] Heintz | Stabilized Lagrange multiplier methods for bilateral elastic contact with friction[END_REF] has to be changed with

ut = 0 =⇒ |σ t (u)| ≤ -F σ n (u) ut = 0 =⇒ σ t (u) = F σ n (u) ut | ut | . (63)
It is easy to check that the total mechanical energy E(t) satisfies

d dt E(t) ≤ L(t) u(t). ( 64
)
In particular, when L vanishes, we get energy dissipation contrary to the frictionless case: E(t 2 ) ≤ E(t 1 ), for all 0 ≤ t 1 ≤ t 2 ≤ T . Next, the space semi-discretized Nitsche-based method for frictional unilateral contact problems in elastodynamics reads:

                  
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 1 γ [P n 1,γ (u h (t))] R -P n θ ,γ (v h ) dΓ + Γ C 1 γ [σ t (u h (t)) -γ uh (t)] -F [P n 1,γ (u h (t))] R - • P t θ ,γ (v h ) dΓ = L(t)(v h ), ∀ v h ∈ V h , u h (0, •) = u h 0 , uh (0, •) = uh 0 , (65)
where u h 0 (resp. uh 0 ) is a finite element approximation in V h of the initial displacement u 0 (resp. the initial velocity u0 ). Using similar tools as in the frictionless case allows us to prove that for any F ≥ 0, θ ∈ R, γ 0 > 0, Problem (65) admits one unique solution u h ∈ C 2 ([0, T ], V h ). Further details and schemes can be found in [START_REF] Chouly | Nitsche method for dynamic contact with Coulomb friction[END_REF]. We end this section with a numerical result that comes from [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 2. Stability of the schemes and numerical experiments[END_REF] (see this reference for more details and other numerical tests). The corresponding test-case concerns multiple impacts of an elastic bar: in the absence of external volume forces, the bar is initially compressed. It is then released without initial velocity. It first impacts the rigid ground, located at x = 0, then gets compressed once again and moves back to its initial position. This problem admits a closed-form solution which derivation and expression are detailed in [START_REF] Dabaghi | A robust finite element redistribution approach for elastodynamic contact problems[END_REF]. Especially it has a periodic motion of period 3. At each period, the bar remains in contact with the rigid ground during one time unit (see Fig. 8). The chosen simulation time is T = 12, so that we can observe 4 successive impacts.

x = 0 L t = 0 t 2 = 2 t 3 = 3 t 1 = 1 Fig. 8 Multiple impacts of an elastic bar. The bar is clamped at x = L and the contact node is located at the bottom. The solution is periodic of period 3, with one impact during each period (here between t = 1 and t = 2).

We make comparisons to the modified mass method (see, e.g., [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF]). The chosen method to compute the modified mass matrix is the simplest possible, since we set the entries associated with the contact node to 0 and no mass redistribution is considered (see also [START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF]). We combine the modified mass method either to the standard (mixed) method or to Nitsche's method for the treatment of contact conditions. The chosen time-marching scheme is Crank-Nicolson (Newmark with γ = 1/2, β = 1/4). These two methods are compared to the Nitsche-Hybrid scheme, with the same parameters and without modified mass. The results are depicted in Fig. 9. The three methods compare well and there is no significant difference: the displacement is free of spurious oscillations, the energy is quite well preserved, with only 1 or 2 % of dissipation and some small spurious oscillations are still present on the contact pressure, that are of similar magnitude. This behaviour is well-known for modified mass combined to mixed discretization of the contact (see, e.g., [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF][START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF]). These results show that, also for Nitsche's discretization of the contact condition, the modified mass improves the quality of the solution in terms of spurious oscillations and energy conservation. In this test-case, the treatment through modified mass produces almost the same effects as the Hybrid time-marching scheme.

Conclusion and perspectives

Nitsche's method allows a numerical treatment of contact and friction in a simple manner since it remains a primal method, and is more robust than penalty, since it is consistent. For the Signorini problem, a rather complete numerical analysis can be carried out, to establish well-posedness and optimal convergence under assumptions on Nitsche's parameters similar to those commonly encountered for Dirichlet boundary conditions. The method can be extended to take into account various situations such as multi-body contact, large transformations and contact in elastodynamics. Most common friction's law such as Coulomb or Tresca can be formulated as well within this framework. Non-symmetric methods corresponding to θ = 0 and θ = -1 reveal to be more robust numerically, though they imply a loss of symmetry in the tangent system. Forthcoming studies may deal with numerical analysis in some situations in which results are lacking, such as contact in elastodynamics, or Coulomb's friction. The same method can be considered as well to discretize other categories of contact / friction problems, or other types of non-linear boundary conditions associated to variational inequalities.

Fig. 1

 1 Fig. 1 Elastic body that occupies the domain Ω . The boundary ∂ Ω is divided into three non-overlapping parts: Γ D (the body is clamped), Γ N (tractions are imposed) and Γ C (contact boundary).

Theorem 2 . 3 2

 23 Suppose that the solution u to Problem (6) belongs to (H +ν (Ω )) d with ν > 0 and d = 2 or d = 3. 1. Let θ ∈ R. Suppose that the parameter γ 0 > 0 is sufficiently large. Then the solution u h to Problem (16) satisfies the following abstract error estimate:

Fig. 2

 2 Fig. 2 Example of two-dimensional mesh and reference solution (with color plot of the Von-Mises stress).

24 γ

 24 0 /E H 1 (Ω ) L 2 (Γ C )

Fig. 3 Fig. 4

 34 Fig. 3 Rate of convergence for uniform and adaptive refinement methods. Parameters γ 0 = E, θ = -1 and Lagrange P 2 elements.

Fig. 7

 7 Fig. 7 Example of use of the unbiased Nitsche's method for the contact with large transformations of two tubes. Only a quarter of each tube is taken into account. Left panel : geometry and mesh in undeformed configuration. Right panel : deformed configuration with Von Mises contour plot.

Fig. 9

 9 Fig. 9 Multiple impacts. Comparison of mixed with modified mass, Nitsche with modified mass and Nitsche-Hybrid methods. Displacement u (left), zoom on displacement u during the fourth impact (center) and discrete energy E h (right).

Setting and Nitsche-based methodIn this section we first present a frictional contact problem in 2.1. Its Nitsche-based finite element approximation is detailed in 2.2. We show in 2.3 that, for the problem under consideration, symmetric Nitsche formulation can be derived from an energy functional. Finally we discuss the relationship with other methods in 2.4: Nitsche for Dirichlet boundary conditions, Barbosa and Hughes stabilized method and the augmented lagrangian of Alart and Curnier.2.1 Unilateral contact with Tresca friction.Ω fondation Γ Γ N Γ C Γ D n Γ N .

see www.getfem.org

Acknowledgements. The authors thank Erik Burman, the editorial board and Springer for the invitation to write this paper for their special volume. Moreover they thank the two anonymous referees for their constructive comments that helped to improve the paper. They thank also Thomas Boiveau and Susanne Claus for the GUFEM meeting and discussions. The first author thanks Région Bourgogne Franche-Comté for partial funding ("Convention Région 2015C-4991. Modèles mathématiques et méthodes numériques pour l'élasticité non-linéaire"), as well as Erik Burman and Miguel A. Fernández for some inspiring discussions on Nitsche's method.

Note additionally that the initial displacement u 0 should satisfy the compatibility condition u 0n ≤ 0 on Γ C . To our knowledge, the well-posedness of Problem ( 53)- [START_REF] Heintz | Stabilized Lagrange multiplier methods for bilateral elastic contact with friction[END_REF] is still an open issue. For partial results, see e.g., [70,[START_REF] Kim | 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] Eck | Unilateral contact problems[END_REF]. Remark 14. The (total) mechanical energy associated with the solution u to the dynamic contact problem ( 53)- [START_REF] Heintz | Stabilized Lagrange multiplier methods for bilateral elastic contact with friction[END_REF] is:

Formally, we get from [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF], after multiplication by u(t), integration by parts, with the boundary conditions on Γ DT , Γ NT , the absence of friction and the persistency condition σ n (u(t)) u(t) = 0 (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]) :

In particular, when L vanishes, we get energy conservation:

Semi-discretization in space with a Nitsche-based finite element method

We consider the same family of finite element spaces V h as defined in Section 2.2. The space semi-discretized Nitsche-based method for unilateral contact problems in elastodynamics then reads:

where u h 0 (resp. uh 0 ) is a finite element approximation in V h of the initial displacement u 0 (resp. the initial velocity u0 ). Nitsche's formulation leads to a well-posed (Lipschitz) system of differential equations, as it will be shown below. This feature is shared with the standard penalty method, the difference being that Nitsche's method remains consistent (for details concerning consistency, we refer the reader to [START_REF] Chouly | A Nitsche finite element method for dynamic contact: 1. Space semi-discretization and timemarching schemes[END_REF]). Since we consider the frictionless case, we modify slightly the definition of (•, •) γ provided initially in Definition 2:

In order to prove well-posedness we reformulate (56) as a system of (non-linear) second-order differential equations. To this purpose, using Riesz's representation theorem in (V h , (•, •) γ ) we first introduce the mass operator M h : V h → V h , which is defined for all v h , w h ∈ V h by (M h v h , w h ) γ = (ρv h , w h ) 0,Ω . Still using Riesz's representation theorem, we define the (non-linear) operator B h : V h → V h , by means of the formula

for all v h , w h ∈ V h . 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 ) γ = L(t)(w h ). With the above notation, Problem (56) reads: