Franz Chouly
email: franz.chouly@univ-fcomte.fr

Mathieu Fabre
email: mathieu.fabre@insa-lyon.fr

Patrick Hild
email: patrick.hild@math.univ-toulouse.fr

Jérôme Pousin
email: jerome.pousin@insa-lyon.fr

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

Residual-based a posteriori error estimation for contact problems approximated by Nitsche's method

Keywords: unilateral contact, finite elements, Nitsche's method, a posteriori error estimates, residuals AMS Subject Classification: 65N12, 65N30, 74M15

We introduce a residual-based a posteriori error estimator for contact problems in two and three dimensional linear elasticity, discretized with linear and quadratic finite elements and Nitsche's method. Efficiency and reliability of the estimator are proved under a saturation assumption. Numerical experiments illustrate the theoretical properties and the good performance of the estimator.

Introduction

The computations of contact problems between deformable bodies are usually obtained with the finite element method [START_REF] Laursen | Computational contact and impact mechanics[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF]. An important aspect for the user is to quantify the quality of the simulations by evaluating the discretization errors coming from the finite element approximation. This quantification requires the definition of a posteriori error estimators which can be of different types (residual based, equilibrated fluxes, smoothing of the stress fields . . . ). The main aim of the estimators is to furnish some information on the local error in order to adapt or refine the mesh and to reduce the computational costs.

Among the finite element discretizations for contact problems, a recent effort was devoted to Nitsche's method which can be seen as a consistent penalty formulation with only one primal unknown (like the penalty or the variational inequality formulation): the displacement field. In contrast to Nitsche's method, the Lagrange (stabilized or standard) methods admit the contact pressure as a supplementary unknown. Nitsche's method was introduced for the frictionless unilateral contact problem in a simple (symmetric) form in [START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF] then generalized and numerically investigated in [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF]. In the latter references the theoretical results deal with well-posedness of the discrete problems and a priori error estimates in two and three space dimensions with linear and quadratic finite elements. A generalization to frictional contact problems is carried out in [START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF] (numerical analysis for Tresca friction) and [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF] (numerical study for Coulomb friction). To our knowledge the a posteriori quantification of the discretization errors committed by the Nitsche finite element approximation has not been considered for unilateral contact problems up to now.

Nevertheless there are several studies concerning a posteriori error analyses for frictionless or frictional contact problems in [START_REF] Carstensen | Adaptive finite elements for elastic bodies in contact[END_REF][START_REF] Fernández | A posteriori error analysis for the normal compliance problem[END_REF][START_REF] Lee | A posteriori error estimation of h-p finite element approximations of frictional contact problems[END_REF][START_REF] Wriggers | Different a posteriori error estimators and indicators for contact problems[END_REF] (residual approach using a penalization of the contact condition or the normal compliance law), in [START_REF] Coorevits | A posteriori error control of finite element approximations for Coulomb's frictional contact[END_REF][START_REF] Coorevits | A posteriori error estimation for unilateral contact with matching and nonmatching meshes[END_REF][START_REF] Louf | Constitutive error estimator for the control of contact problems involving friction[END_REF][START_REF] Weiss | A posteriori error estimator and error control for contact problems[END_REF][START_REF] Wohlmuth | An a posteriori error estimator for two body contact problems on nonmatching meshes[END_REF] (equilibrated residual method), in [START_REF] Eck | A residual-based error estimator for BEM-discretizations of contact problems[END_REF][START_REF] Maischak | Adaptive hp-versions of BEM for Signorini problems[END_REF] (residual approach for BEM-discretizations), in [START_REF] Blum | An adaptive finite element discretization for a simplified Signorini problem[END_REF][START_REF] Schröder | Goal-oriented error control in adaptive mixed FEM for Signorini's problem[END_REF] (error technique measure developed in [START_REF] Becker | A feed-back approach to error control in finite element methods: basic analysis and examples[END_REF]). Moreover a residual type estimator for the Signorini problem in its common formulation (variational inequality or mixed method) can be found in [START_REF] Hild | A posteriori error estimations of residual type for Signorini's problem[END_REF][START_REF]Residual a posteriori error estimators for contact problems in elasticity[END_REF][START_REF] Schröder | A posteriori error estimates of higher-order finite elements for frictional contact problems[END_REF] and in the recent work [START_REF] Krause | An efficient and reliable residual-type a posteriori error estimator for the Signorini problem[END_REF].

Finally we mention that only few works are devoted to a posteriori error estimates for Nitsche's method, and all concern linear boundary / interface conditions. For interface conditions and elliptic problems, Hansbo and Hansbo [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF] introduce a residual type estimator for a Nitsche's unfitted treatment of the interface condition. They prove an upper bound on a linear functional of the error, in the spirit of Becker and Rannacher [START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF]. Note as well an early work of Becker [START_REF] Becker | Mesh adaptation for Dirichlet flow control via Nitsche's method[END_REF] in the context of optimal control for Navier-Stokes equations, with a Nitsche treatment of Dirichlet boundary condition and an a posteriori error estimate for the functional to minimize. Residual error estimates are introduced as well by Becker, Hansbo and Stenberg in [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF] for a Nitsche-based domain decomposition with non-matching meshes. Upper bounds both in H 1 and L 2 -norms are established, with help of a saturation assumption (as in [START_REF] Wohlmuth | A residual based error estimator for mortar finite element discretizations[END_REF]) for the H 1 -norm. In the context of composite grids two variants of residual-based error estimates are proposed by Hansbo, Hansbo and Larson in [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF]. Upper bounds in H 1 -norm without any saturation assumption are proposed for both of them. Later on Juntunen and Stenberg provide in [START_REF] Juntunen | On a mixed discontinuous Galerkin method[END_REF] a residual-based error estimator for the stabilized Bassi-Rebay discontinuous Galerkin method, that relies on Nitsche's treatment of continuity. Upper and lower bounds are proven for this method. The same authors in [START_REF]Nitsche's method for general boundary conditions[END_REF] introduce a Nitsche's method for a general boundary condition, and an associated residual error estimator. They prove an upper bound in H 1 -norm under a saturation assumption (as in [START_REF] Braess | A posteriori error estimators for the Raviart-Thomas element[END_REF]), and they establish a lower bound too. Finally let us mention two recent papers on the Brinkman problem by Juntunen, Könnö and Stenberg [START_REF]Analysis of finite element methods for the Brinkman problem[END_REF][START_REF] Könnö | H(div)-conforming finite elements for the Brinkman problem[END_REF].

The paper is outlined as follows: in Section 2 the Nitsche's finite element discretization for contact problems in linear elasticity is described and the results dealing with well-posedness are recalled from [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF]. In Section 3 a residual a posteriori error estimator is introduced, and we prove its reliability and efficiency. In Section 4 numerical experiments in 2D and 3D illustrate the theoretical results and allow to assess the quality of the estimator, for different values of the numerical parameters.

Let us introduce some useful notations. In what follows, bold letters like u, v, indicate vector or tensor valued quantities, while the capital ones (e.g., V, K . . .) represent functional sets involving vector fields. As usual, we denote by (H s (.)) d , s ∈ R, d = 1, 2, 3, the Sobolev spaces in one, two or three space dimensions (see [START_REF] Adams | Sobolev spaces[END_REF]) with the convention H 0 = L 2 . The usual norm (resp. semi-norm) of (H s (D)) d is denoted by For two scalar quantities a and b, the notation a b means there exists a constant C, independent of the mesh size parameters and of the Nitsche parameter γ 0 (see Section 2.2), such that a ≤ Cb. Moreover, a ∼ b means that a b and b a.

Setting

The unilateral contact problem

We consider an elastic body whose reference configuration is represented by the domain Ω in R d with d = 2 or d = 3. Small strain assumptions are made, as well as plane strain when d = 2. The boundary ∂Ω of Ω is polygonal or polyhedral and we suppose that ∂Ω consists in three nonoverlapping parts Γ D , Γ N and the (candidate) contact boundary Γ C , with meas(Γ D ) > 0 and meas(Γ C ) > 0. The (candidate) contact boundary is supposed to be a straight line segment when d = 2 or a polygon when d = 3 to simplify. The unit outward normal vector on ∂Ω is denoted n. In its initial stage, the body is in contact on Γ C with a rigid foundation (the extension to two elastic bodies in contact can be easily made, at least for small strain models) and we suppose that the unknown final contact zone after deformation will be included into Γ C . The body is clamped on Γ D for the sake of simplicity. It is subjected to volume forces f ∈ (L 2 (Ω)) d and to surface loads

g ∈ (L 2 (Γ N )) d .
The unilateral contact problem in linear elasticity consists in finding the displacement field u : Ω → R d verifying the equations and conditions (1)-( 2):

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

v = v n n + v t and σ(v)n = σ n (v)n + σ t (v),
where v t (resp. σ t (v)) are the tangential components of v (resp. σ(v)n). The conditions describing unilateral contact without friction on Γ C are:

u n ≤ 0, (i) σ n (u) ≤ 0, (ii) σ n (u) u n = 0, (iii) σ t (u) = 0. (iv) (2)
We introduce the Hilbert space V and the convex cone K of admissible displacements which satisfy the noninterpenetration on the contact zone Γ C :

V := v ∈ (H 1 (Ω)) d : v = 0 on Γ D , K := {v ∈ V : v n = v • n ≤ 0 on Γ C } .
We define as well

a(u, v) := Ω σ(u) : ε(v) dΩ, L(v) := Ω f • v dΩ + Γ N g • v dΓ,
for any u and v in V. From the previous assumptions, we deduce that a(•, •) is bilinear, symmetric, V-elliptic and continuous on V × V. Likewise we observe that L(•) is a continuous linear form on V. The weak formulation of Problem ( 1)-( 2), as a variational inequality (see [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF][START_REF] Haslinger | Numerical methods for unilateral problems in solid mechanics[END_REF][START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]), reads as:

Find u ∈ K such that: a(u, v -u) ≥ L(v -u), ∀ v ∈ K. (3)
Stampacchia's Theorem ensures that Problem (3) admits a unique solution.

Finite element setting and Nitsche-based method

To approximate problem (3) we fix a family of meshes (T h ) h>0 , regular in Ciarlet's sense (see [START_REF] Ciarlet | Handbook of Numerical Analysis[END_REF]), made of closed elements and assumed to be subordinated to the decomposition of the boundary ∂Ω into Γ D , Γ N and Γ C . For K ∈ T h we recall that h K is the diameter of K and h := max

K∈T h h K .
The regularity of the mesh implies notably that for any edge (or face when

d = 3) E of K one has h E := |E| ∼ h K .
Let us define E h as the set of edges (or faces when d = 3) of the triangulation and define For an element K, we will denote by E K the set of edges/faces of K and according to the above notation, we set

E int h = {E ∈ E h : E ⊂ Ω}
E int K = E K ∩ E int h , E N K = E K ∩ E N h , E C K = E K ∩ E C h .
For an edge/face E of an element K, introduce ν K,E the unit outward normal vector to K 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 ).
Note that the sign of v E depends on the orientation of ν E . Finally we will need local subdomains (also called patches). As usual, let ω K be the union of all elements having a nonempty intersection with K. Similarly for a node x and an edge/face E, let

ω x = ∪ K:x∈K K and ω E = ∪ x∈E ω x .
The chosen finite element space V h ⊂ V involves standard Lagrange finite elements of degree k with k = 1 or k = 2 (see [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF][START_REF] Ciarlet | Handbook of Numerical Analysis[END_REF][START_REF] Ern | Theory and practice of finite elements[END_REF]) i.e.:

V h := v h ∈ (C 0 (Ω)) d : v h |K ∈ (P k (K)) d , ∀K ∈ T h , v h = 0 on Γ D .
Let us introduce the notation [•] + for the positive part of a scalar quantity a ∈ R: [a] + = a if a > 0 and [a] + = 0 otherwise. The monotonicity property below holds:

([a] + -[b] + )(a -b) ≥ ([a] + -[b] + ) 2 ≥ 0. ( 4
)
Note that the condition (4) can be straightforwardly extended to real valued functions.

Let γ be a positive piecewise constant function on the contact interface Γ C which satisfies

γ K∩Γ C = γ 0 h K ,
for every K that has a non-empty intersection of dimension d -1 with Γ C , and where γ 0 is a positive given constant which we call the "Nitsche parameter". Note that the value of γ on element intersections has no influence. We introduce the discrete linear operator

P γ : V h → L 2 (Γ C ) v h → v h n -γ σ n (v h )
, and the bilinear form where θ ∈ R is a fixed parameter:

A θγ (u h , v h ) := a(u h , v h ) - Γ C θγ σ n (u h )σ n (v h ) dΓ.
Our Nitsche-based method then reads:

     Find u h ∈ V h such that: A θγ (u h , v h ) + Γ C 1 γ [P γ (u h )] + P θγ (v h ) dΓ = L(v h ), ∀ v h ∈ V h . (5)
We consider the quasi-interpolation (regularization) operator introduced in e.g. [8, formula (4.11)] and its straightforward extension to the vectorial case, that we denote R h : V → V h . This operator has the following approximation and stability properties:

Lemma 2.1. For any v ∈ V ∩ (H l (Ω)) d , 1 ≤ l ≤ k + 1, the following estimates hold v -R h v 0,K h l K |v| l,ω K , ∀K ∈ T h , (6)
v -R h v 0,E h l-1/2 E |v| l,ω E , ∀E ∈ E h . ( 7
)
Moreover R h is stable in the H 1 -norm, i.e. R h v 1,Ω v 1,Ω , ∀v ∈ V. (8)
Proof. Estimates ( 6) and ( 7) are provided in [8, Theorem 4.8, Remark 8]. The stability of R h in H 1 -norm is proven in [START_REF] Bernardi | A local regularization operator for triangular and quadrilateral finite elements[END_REF]Theorem 4.4] (in all cases, it suffices to apply the results of [START_REF] Bernardi | A local regularization operator for triangular and quadrilateral finite elements[END_REF] component-wise).

We next define a convenient mesh dependent norm which is in fact a weighted

L 2 (Γ C )-norm (since γ/γ 0 = h K ). Definition 2.2. For any v ∈ L 2 (Γ C ), we set v -1/2,h,Γ C := γ γ 0 1 2 v 0,Γ C .
We end this subsection with a discrete trace inequality, that will be useful for the analysis (for the proof, see, e.g., [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF]):

Lemma 2.3. For any v h ∈ V h , we have σ n (v h ) -1/2,h,Γ C v h 1,Ω . (9)

Consistency and well-posedness of the Nitsche-based method

We recall two theoretical properties for the Nitsche-based method (5): consistency and wellposedness. These properties, together with optimal a priori error estimates in H 1 (Ω)-norm, are proven in [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF]. Like Nitsche's method for second order elliptic problems with Dirichlet boundary conditions or domain decomposition [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF], our Nitsche-based formulation (5) for unilateral contact is consistent: Lemma 2.4. The Nitsche-based method for contact is consistent: suppose that the solution u of (1)-( 2) lies in (H

3 2 +ν (Ω)) d with ν > 0 and d = 2, 3. Then u is also solution of A θγ (u, v h ) + Γ C 1 γ [P γ (u)] + P θγ (v h ) dΓ = L(v h ), ∀ v h ∈ V h .
Problem ( 5) is well-posed in the following sense and under the assumptions below:

Theorem 2.5. Suppose that either θ = -1 and γ 0 > 0 is sufficiently small, or θ = -1 and γ 0 > 0. Then Problem (5) admits one unique solution u h in V h .

Remark 2.6. When γ 0 is large and θ = -1 we can neither conclude to uniqueness, nor to existence of a solution. In reference [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF] there are some simple explicit examples of nonexistence and nonuniqueness of solutions.

3 A posteriori error analysis

Definition of the residual error estimator

The element residual of the equilibrium equation in ( 1) is defined by div σ(u h ) + f in K.

Remark 3.1. For linear elements (k = 1), the term div σ(u h ) vanishes.

As usual this element residual can be replaced by some finite dimensional approximation, called approximate element residual (see, e.g., [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF])

div σ(u h ) + f K , f K ∈ (P l (K)) d , l ≥ 0. A current choice is to take f K = K f (x) /|K| dx since for f ∈ (H 1 (Ω)) d , scaling arguments yield f -f K 0,K h K f 1,K
and it is then negligible with respect to the estimator η defined hereafter. In the same way g is approximated by a computable quantity denoted g E on any E ∈ E N h .

Definition 3.2. The local error estimators η K and the the global estimator η are defined by

η K = 4 i=1 η 2 iK 1/2 , η 1K = h K div σ(u h ) + f K 0,K , η 2K = h 1/2 K   E∈E int K ∪E N K J E,n (u h ) 2 0,E   1/2 , η 3K = h 1/2 K   E∈E C K σ t (u h ) 2 0,E   1/2 , η 4K = h 1/2 K   E∈E C K 1 γ [P γ (u h )] + + σ n (u h ) 2 0,E   1/2 , η =   K∈T h η 2 K   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 -g E , ∀E ∈ E N h . (10)
The local and global approximation terms are given by

ζ K =   h 2 K K ⊂ω K f -f K 2 0,K + h E E⊂E N K g -g E 2 0,E   1/2 , ζ =   K∈T h ζ 2 K   1/2 .

Upper error bound

First we state a "saturation" assumption that we need in order to prove the estimate (see also [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF] in the case of Nitsche for domain decomposition, and [START_REF] Wohlmuth | A residual based error estimator for mortar finite element discretizations[END_REF] for mortar methods).

Assumption 3.3. The solution u of (3) and the discrete solution u h of ( 5) are such that:

σ n (u -u h ) -1/2,h,Γ C u -u h 1,Ω . (11)
Remark 3.4. Note that for a Nitsche's treatment of (linear) interface conditions, an upper bound for a residual-based estimator has been derived without such an assumption in [START_REF] Hansbo | A finite element method on composite grids based on Nitsche's method[END_REF]. Similarly for some classes of mixed non-conforming finite element approximations, an assumption such as 3.3 has revealed to be superfluous, see, e.g., references [START_REF] Carstensen | A posteriori error estimate for the mixed finite element method[END_REF][START_REF] Kim | A posteriori error analysis for locally conservative mixed methods[END_REF]. However for the method (5) the derivation of an upper bound without such kind of saturation assumption remains an open issue.

The following statement guarantees the reliability of the a posteriori error estimator given in Definition 3.2: Theorem 3.5. Let u be the solution to the variational inequality (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] Becker | A finite element method for domain decomposition with non-matching grids[END_REF]. Assume that, for θ = -1, γ 0 is sufficiently small, and otherwise that γ 0 > 0 for θ = -1. Assume that the saturation assumption [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF] holds as well. Then we have

u -u h 1,Ω + σ n (u) + 1 γ [P γ (u h )] + -1/2,h,Γ C + σ n (u) -σ n (u h ) -1/2,h,Γ C (1 + γ 0 )η + ζ. Proof: Let v h ∈ V h .
To lighten the notation, we denote e := uu h . We first use the V-ellipticity of a(•, •), together with the Green formula, equations ( 1) and ( 5) to obtain:

α e 2 1,Ω ≤ a(u -u h , u -u h ) = a(u, u -u h ) -a(u h , u -v h ) -a(u h , v h -u h ) = L(u -u h ) + Γ C σ n (u)(u n -u h n ) dΓ -a(u h , u -v h ) -L(v h -u h ) + Γ C 1 γ [P γ (u h )] + P θγ (v h -u h ) dΓ -θ Γ C γ σ n (u h )σ n (v h -u h ) dΓ = T 1 + T 2 ( 12
)
where α is the V-ellipticity constant of a(•, •) and

T 1 := L(u -v h ) -a(u h , u -v h ) + Γ C 1 γ [P γ (u h )] + (v h n -u n ) dΓ T 2 := Γ C σ n (u)(u n -u h n ) dΓ + Γ C 1 γ [P γ (u h )] + P θγ (u -u h ) dΓ -θ Γ C 1 γ [P γ (u h )] + γσ n (v h -u) dΓ -θ Γ C γ σ n (u h )σ n (v h -u h ) dΓ.
The quantity T 1 is an expression which is handled hereafter in a classical way. Namely, by integrating by parts on each triangle K, using the definition of J E,n (u h ) in [START_REF] Braess | A posteriori error estimators for the Raviart-Thomas element[END_REF] and splitting up the integrals on Γ C into normal and tangential components we get:

T 1 = K∈T h K (div σ(u h ) + f ) • (u -v h ) dΓ + E∈E C h E 1 γ [P γ (u h )] + + σ n (u h ) (v h n -u n ) dΓ + E∈E C h E σ t (u h ) • (v h t -u t ) dΓ - E∈E int h ∪E N h E J E,n (u h ) • (u -v h ) dΓ + E∈E N h E (g -g E ) • (u -v h ) dΓ. (13)
We now need to estimate each term of this right-hand side. For that purpose, we take

v h = u h + R h (u -u h ), ( 14
)
where R h is the quasi-interpolation operator defined in Section 2.2. We start with the integral term on elements K. Cauchy-Schwarz's inequality implies

K∈T h K (div σ(u h ) + f ) • (u -v h ) dΓ ≤ K∈T h div σ(u h ) + f 0,K u -v h 0,K ,
and it suffices to estimate uv h 0,K for any triangle K. From the definition of v h and ( 6) we get:

u -v h 0,K = e -R h e 0,K h K e 1,ω K . As a consequence Ω (div σ(u h ) + f ) • (u -v h ) dΓ (η + ζ) e 1,Ω .
We now consider the interior and Neumann boundary terms in [START_REF] Carstensen | Adaptive finite elements for elastic bodies in contact[END_REF]. As we previously noticed, the application of Cauchy-Schwarz's inequality leads to

E∈E int h ∪E N h E J E,n (u h ) • (u -v h ) dΓ ≤ E∈E int h ∪E N h J E,n (u h ) 0,E u -v h 0,E .
Therefore using expression [START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF] and estimate [START_REF] Becker | An optimal control approach to a posteriori error estimation in finite element methods[END_REF], we obtain

u -v h 0,E = e -R h e 0,E h 1/2 E e 1,ω E .
Inserting this estimate in the previous one we deduce that

E∈E int h ∪E N h E J E,n (u h ) • (u -v h ) dΓ η e 1,Ω .
Moreover,

E∈E N h E (g -g E ) • (u -v h ) dΓ ζ e 1,Ω .
The two following terms are handled in a similar way as the previous ones so that

E∈E C h E 1 γ [P γ (u h )] + + σ n (u h ) (v h n -u n ) dΓ η e 1,Ω ,
and

E∈E C h E σ t (u h ) • (v h t -u t ) dΓ η e 1,Ω .
Collecting the previous results, we deduce:

T 1 (η + ζ) e 1,Ω . (15)
The first two terms in T 2 are split using the definition of P γ (•) and P θγ (•), and the last one is split using relationship

σ n (v h -u h ) = σ n ((v h -u) + (u -u h ))
. This leads to:

T 2 = Γ C σ n (u)P γ (u -u h ) dΓ + Γ C σ n (u)γσ n (u -u h ) dΓ + Γ C 1 γ [P γ (u h )] + P γ (u -u h ) dΓ + (1 -θ) Γ C 1 γ [P γ (u h )] + γσ n (u -u h ) dΓ -θ Γ C 1 γ [P γ (u h )] + + σ n (u h ) γσ n (v h -u) dΓ -θ Γ C γ σ n (u h )σ n (u -u h ) dΓ.
Then we split the second term in the above expression using 1 = θ + (1 -θ) and we gather the resulting terms:

T 2 = Γ C 1 γ [P γ (u h )] + + σ n (u) P γ (u -u h ) dΓ + (1 -θ) Γ C 1 γ σ n (u) + [P γ (u h )] + γσ n (u -u h ) dΓ -θ Γ C 1 γ [P γ (u h )] + + σ n (u h ) γσ n (v h -u) dΓ + θ γ 1/2 σ n (u -u h ) 2 0,Γ C .
Now we substitute σ n (u) using the reformulation of contact conditions (2 [START_REF] Alart | A generalized Newton method for contact problems with friction[END_REF][START_REF] Chouly | A Nitsche-based method for unilateral contact problems: numerical analysis[END_REF]). This reformulation makes sense in L 2 (Γ C ) due to the regularity assumption u ∈ (H 3 2 +ν (Ω)) d . Afterwards we apply the bound (4) in the first term as well as Cauchy-Schwarz inequality in the second one:

) (i)-(iii) as σ n (u) = -1 γ [P γ (u)] + (see for instance
T 2 ≤ -γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 2 0,Γ C + |θ -1| γ 1 2 (σ n (u) + 1 γ [P γ (u h )] + ) 0,Γ C γ 1/2 σ n (u -u h ) 0,Γ C -θ Γ C 1 γ [P γ (u h )] + + σ n (u h ) γσ n (v h -u) dΓ + θ γ 1/2 σ n (u -u h ) 2 0,Γ C .
The expression ab ≤ a 2 + b 2 /4 yields, for any β > 0:

T 2 ≤ |θ -1| 2 4 γ 1/2 σ n (u -u h ) 2 0,Γ C -θ Γ C 1 γ [P γ (u h )] + + σ n (u h ) γσ n (v h -u) dΓ + θ γ 1/2 σ n (u -u h ) 2 0,Γ C = (θ + 1) 2 4 γ 1/2 σ n (u -u h ) 2 0,Γ C -θ Γ C 1 γ [P γ (u h )] + + σ n (u h ) γσ n (v h -u) dΓ ≤ (θ + 1) 2 4 γ 1/2 σ n (u -u h ) 2 0,Γ C + |θ|γ 1/2 0 η γ 1/2 σ n (v h -u) 0,Γ C ≤ (θ + 1) 2 4 γ 1/2 σ n (u -u h ) 2 0,Γ C + βθ 2 γ 0 η 2 + 1 2β γ 1/2 σ n (v h -u h ) 2 0,Γ C + 1 2β γ 1/2 σ n (u h -u) 2 0,Γ C = 1 2β + (θ + 1) 2 4 γ 0 σ n (u -u h ) 2 -1/2,h,Γ C + βθ 2 γ 0 η 2 + γ 0 2β σ n (v h -u h ) 2 -1/2,h,Γ C .
Using [START_REF] Blum | An adaptive finite element discretization for a simplified Signorini problem[END_REF] and the H 1 -stability of R h (see [START_REF] Bernardi | A local regularization operator for triangular and quadrilateral finite elements[END_REF] in Lemma 2.1) we bound:

σ n (v h -u h ) -1/2,h,Γ C ≤ C v h -u h 1,Ω = C R h (u -u h ) 1,Ω ≤ C u -u h 1,Ω .
We combine this last bound with the saturation assumption (11) and get

T 2 ≤ Cγ 0 (θ + 1) 2 4 + 1 β u -u h 2 1,Ω + βθ 2 γ 0 η 2 . ( 16
)
Now we combine estimates ( 12), ( 15) and ( 16):

α e 2 1,Ω ≤ C(η + ζ) e 1,Ω + Cγ 0 (θ + 1) 2 4 + 1 β e 2 1,Ω + βθ 2 γ 0 η 2 .
We treat the first term on the right-hand side with Young's inequality and obtain:

α 2 -Cγ 0 (θ + 1) 2 4 + 1 β e 2 1,Ω ≤ C α (η 2 + ζ 2 ) + βθ 2 γ 0 η 2 .
When θ = -1, we choose γ 0 sufficiently small, and for θ = -1, we can choose for instance β = 4Cγ 0 α (for a fixed value of γ 0 > 0, that does not need to be small in this case). We obtain the upper bound on the error in natural norm:

e 1,Ω (1 + γ 0 )η + ζ.
The saturation assumption [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF] provides directly a bound on the contact stress error:

σ n (u -u h ) -1/2,h,Γ C e 1,Ω .
For the contact error we make use of the triangle inequality, and of the above inequality:

σ n (u) + 1 γ [P γ (u h )] + -1/2,h,Γ C ≤ σ n (u -u h ) -1/2,h,Γ C + σ n (u h ) + 1 γ [P γ (u h )] + -1/2,h,Γ C e 1,Ω + η.
Collecting the three previous results allows to prove the theorem.

Lower error bound

We now consider the local lower error bounds of the discretization error terms.

Theorem 3.6. For all elements K ∈ T h , the following local lower error bounds hold:

η 1K u -u h 1,K + ζ K , (17)
η 2K u -u h 1,ω K + ζ K . ( 18
)
For all elements K such that K ∩ E C h = ∅, the following local lower error bounds hold:

η 3K u -u h 1,K + ζ K , (19)
η 4K E∈E C K h 1/2 K σ n (u) + 1 γ [P γ (u h )] + 0,E + σ n (u -u h ) 0,E . (20)
Proof: The estimates of η 1K , η 2K in ( 17)-( 18) are standard (see, e.g., [START_REF] Verfürth | A review of a posteriori error estimation techniques for elasticity problems[END_REF]). The estimate η 3K in [START_REF] Coorevits | A posteriori error estimation for unilateral contact with matching and nonmatching meshes[END_REF] is handled in a standard way as in [START_REF]Residual a posteriori error estimators for contact problems in elasticity[END_REF].

The estimate of η 4K in ( 20) is obtained from the Definition 3.2 and triangular inequality.

Remark 3.7. Remark that, from Theorem 3.6, optimal convergence rates of order O(h min(k, 1 2 +ν) ) are expected for the estimator of Definition 3.2. Remark 3.8. An extension of the above analysis for the Tresca friction case is sketched in Appendix 1.

Numerical experiments

We illustrate numerically the theoretical properties of the error estimator η given in Definition 3.2 and compute its convergence order when h vanishes. To study separately global contributions of each component of η we introduce the notation

η i =   K∈T h η 2 iK   1/2 1 ≤ i ≤ 4,
where the expressions of η iK are provided in Definition 3.2. In all examples below Hooke's law is considered: E and ν P will denote respectively Young's modulus and Poisson's ratio. Moreover a dimensional analysis allows to deduce that γ 0 is the inverse of a stiffness parameter. Consequently we choose in our discussion γ 0 = C/E where C is a constant which does not depend on E. The finite element method (5) as well as the residual estimator η are implemented under the open source finite element library GetFEM++ 1 . For details on numerical solving we refer to [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] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF].

To measure the quality of the estimator η we introduce the effectivity index:

ef f E = η E u -u h 1,Ω .
As in [START_REF]Residual a posteriori error estimators for contact problems in elasticity[END_REF] this index has been normalized with respect to Young's modulus E. Indeed, we remark that if u(E) denotes the solution of a (linear) Lamé system with a Young modulus E then u(mE) = u(E)/m whereas σ(u(mE)) = σ(u(E)). Thus the error estimator η is independent of E (for η 4 this property comes from the scaling γ 0 = C/E). In contrast there holds u(mE) -u h (mE)

1,Ω = u(E) -u h (E) 1,Ω
/m that becomes independent of E for the choice m = 1/E.

First example: a square with slip and separation

Description

We first consider a test case taken from [START_REF]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 = 1/E and θ = -1.

First of all we illustrate in Figure 1 the difference between uniform and adaptive refinement. For the latter we refine only the mesh elements K in which the local estimator η K 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.

1 see http://download.gna.org/getfem/html/homepage/ The solution obtained with adaptive refinement and θ = -1 is depicted in Figure 2. 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. As expected we observe that all the nodes on Γ C have a negative tangential displacement and that Γ C is divided into two parts: the upper part where the body remains in contact (slipping nodes) and the lower part where it is separated, with a transition point near (1, 0.685). The value is close to the transition point (1, 0.69 ± 0.01) found in [START_REF]Residual a posteriori error estimators for contact problems in elasticity[END_REF] and (1, 0.65) found in [START_REF] Lleras | Modélisation, analyse et simulation de problèmes de contact en mécanique des solides et des fluides[END_REF]. The slight difference with [START_REF] Lleras | Modélisation, analyse et simulation de problèmes de contact en mécanique des solides et des fluides[END_REF] should be due to Coulomb friction. Remark 4.1. Note that the solution in the case θ = 1 (see Figure 3) has an error estimator on the contact zone which is larger than in the case θ = -1. In the case θ = -1, the discrete solution is less dependent on the parameter γ 0 than for the other methods (see [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF]) and we obtain a better approximation of the problem on the contact boundary Γ C .

Numerical convergence

We perform a numerical convergence study for three variants of method ( 5) corresponding to θ = 1, θ = 0 and θ = -1. The Nitsche parameter γ 0 is fixed to 1/E, which should ensure wellposedness and optimal convergence in each case. Lagrange P 1 finite elements are chosen. The reference solution for error computations corresponds to the one described in Section 4.1.1 and depicted in Figure 2 (P 2 finite elements, θ = -1 and adaptive finest mesh). No mesh adaptation is carried out anymore and only uniform refinement is imposed, with a sequence of decreasing mesh sizes h.

First the estimator η, the L 2 and the H 1 -norms of the error uu h are depicted in Figure 4. One can note a sub-optimality of the convergence rate in L 2 and H 1 -norms of the error. They are caused by the Neumann-Dirichlet transition on the left corners of Ω (the same observation has been reported in [START_REF] Fabre | A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method[END_REF]). Then the different contributions of η are reported on Tables 1, 2 and 3. The convergence rate of η 1 is strictly equal to 1 since, for piecewise linear finite elements, the expression of this estimator reduces to η 1K = h K f K 0,K . More generally, all the estimators η i converge towards zero as h vanishes, and they behave identically whatever the value of θ is (this is due to the low value of γ 0 ). Moreover, the convergence rate of η 2 is slightly less than the one of the H 1 -norm of the error whereas the convergence rates of η 3 and η 4 are far greater, and higher than 1 (we do not have a clear interpretation of this). In all cases, we obtain an effectivity index between 1.2 and 1.6 (the average is close to 1.45 and the standard deviation is close to 0.12). These overall results are quite similar to those presented in [START_REF]Residual a posteriori error estimators for contact problems in elasticity[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]. On Figure 5 the numerical experiment is performed for a larger parameter γ 0 = 1000/E. In the case θ = 1 and in the case θ = 0 the convergence rate is degraded compared to the case γ 0 = 1/E. Conversely, in the case θ = -1, the convergence is not deteriorated which confirms the theoretical results obtained in both the a priori analysis in [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF] and the a posteriori analysis in Section 3 (see Theorem 3.5). Additionally we present a numerical convergence study for θ = 1, 0, -1 and for a very large value of the parameter γ 0 = 10 6 /E, far from its reference value of 1/E. In this case for θ = 1 and θ = 0 there is no more guarantee of well-posedness and optimal convergence (see [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF]). The error estimator η, the L 2 and H 1 -norms of the error uu h are plotted in Figure 6, while Tables 4,5 and 6 present the different contributions of η. For the method θ = 0 the solution does not converge while the effectivity index Ef f E tends to 0. This is consistent with our theoretical results since Theorem 3.5 is no more applicable and no upper bound is guaranteed. The estimator η converges despite of the term η 4 , slightly increasing but still very small. For the method θ = 1, even though γ 0 is large, the method converges in L 2 and H 1 -norms of the error with an acceptable effectivity index, however with a deteriorated convergence rate. Conversely, for the method θ = -1, both convergence and effectivity index are optimal and are not deteriorated compared to the case γ 0 = 1/E. This supports its theoretical property of robustness with respect to γ 0 . The previous experiment for θ = 0 reveals the bad behavior of η for a very large γ 0 . An heuristics to recover a meaningful estimator is to decouple the value of γ 0 for the problem [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF] and for the estimator η 4 . A last experiment shows in Figure 7 the convergence curves in the same case γ 0 = 10 6 /E yet with an error estimator that makes use of a Nitsche's parameter γ0 = 1/E. For the methods θ = 1 and θ = -1 this has no visible influence on the effectivity index Ef f E . For the method θ = 0, a better effectivity index is obtained: at least the estimator does not tend to zero for a non-convergent solution, in contrast with what happens in Figure 6. To summarize, this study for a large γ 0 confirms the analysis provided in Section 3 which requires a sufficiently small parameter γ 0 to obtain a reliable and efficient a posteriori estimator when θ = -1.

Second example: Hertz's contact

Description

We consider Hertz's contact problems of a disk/a sphere with a plane rigid foundation (see, e.g., the numerical examples in [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]Residual a posteriori error estimators for contact problems in elasticity[END_REF]). The parameters have been fixed as θ = -1 and γ 0 = 10 -3 /E. The disc (resp. the sphere) is of center (0, 20) (resp. of center (0, 0, 20)) and radius 20. The lower part of the boundary Γ C is potentially in contact with the rigid support ; the remaining (upper part) of the boundary Γ N is subjected to a homogenous Neumann condition. To overcome the non-definiteness coming from free rigid motions, the horizontal displacement is prescribed to be zero on the two points of coordinates (0, 20) and (0, 25) (resp. the horizontal displacement components u 1 and u 2 on the point (0,0,20), the component u 1 on the point (0,5,20) and the component u 2 on the point (5,0,20)): this blocks the horizontal translation and the rigid rotation. Young's modulus is fixed at E = 25 and Poisson's ratio is ν P = 0.25. A vertical density of volume forces of intensity 20 is applied in Ω. The reference solutions are depicted in Figure 8. There are refined solutions with an average mesh size h = 0.10 for the disc (resp. h = 1.27 for the sphere), Lagrange P 2 elements, θ = -1 and γ = 10 -3 /E.

The initial gap between Γ C and the obstacle is computed as gap(x) := x • n obs , where x ∈ Γ C and with n obs the unit outward normal vector on the boundary of the plane obstacle. The error estimator associated to the contact condition on Γ C is then slightly modified:

η 4K = h 1/2 K   E∈E C K 1 γ [P γ (u h ) -gap] + + σ n (u h ) 2 0,E   1/2
.

Numerical convergence in 2D

The error curves in the 2D case are depicted in Figure 9, both for linear and quadratic finite elements. In the case of P 1 finite elements, and as in [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF], a slight super-convergence is observed in the H 1 -norm of the error (1.5 instead of 1). This behavior is not recovered by the error estimator η, that converges with a rate close to 1. The origin of this difference is unknown. For P 2 finite elements the agreement between η and the error in H 1 -norm is better: for the H 1 -norm, the convergence rate is close to 1.7, while of approximately 1.5 for η. We observe the same results for the variants θ = 0, 1. In Table 7 the contribution of each component η i of η is detailed. Each term of the error estimator converges towards zero when h becomes smaller. Note however the increasing values of the effectivity index, due to the super-convergence in H 1 -norm and the convergence rate of the contribution η 4 , that is close to 1.5.

Numerical convergence in 3D

The error curves in the 3D case are depicted in Figure 10, both for linear and quadratic finite elements. For P 1 finite elements the convergence rates for η and for the error in H 1 -norm are close (around 1.3) and slightly above the expected rate of 1. For P 2 finite elements we observe a sub-optimality of the error estimator η, which converges but with a rate of 1, while the error in H 1 -norm remains optimal, with a convergence rate around 1.5. Lemma 4.3. For any v h ∈ V h , we have

σ n (v h ) -1/2,h,Γ C + σ t (v h ) -1/2,h,Γ C v h 1,Ω . (25)
For contact with Tresca friction, the following statement guarantees the reliability of the a posteriori error estimator: Theorem 4.4. Let u be the solution to (1)-(2 (i,ii,iii))-( 21), 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] Fabre | A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method[END_REF]. Assume that, for θ = -1, γ 0 is sufficiently small, and otherwise that γ 0 > 0 for θ = -1. Assume that the saturation assumption [START_REF] Fernández | A posteriori error analysis for the normal compliance problem[END_REF] Proof: The proof is a direct adaptation of Theorem 3.5. Let v h ∈ V h . To lighten the notation, we denote e := uu h . We start as in Theorem 3.5 and get

α e 2 1,Ω ≤ T 1 + T 2
where α is the V-ellipticity constant of a(•, •) and

T 1 := L(u -v h ) -a(u h , u -v h ) + Γ C 1 γ [P γ (u h )] + (v h n -u n ) dΓ + Γ C 1 γ P t γ (u h ) γg • (v h t -u t ) dΓ T 2 := Γ C σ(u)n • (u -u h ) dΓ + Γ C 1 γ [P γ (u h )] + P θγ (u -u h ) dΓ + Γ C 1 γ P t γ (u h ) γg • P t θγ (u -u h ) dΓ -θ Γ C 1 γ [P γ (u h )] + γσ n (v h -u) dΓ -θ Γ C 1 γ P t γ (u h ) γg • γσ t (v h -u) dΓ -θ Γ C γ σ(u h )n • σ(v h -u h )n dΓ.
The quantity T 1 is bounded almost exactly as in Theorem 3.5, except for the new Tresca friction term, that is bounded as follows:

E∈E C h E 1 γ P t γ (u h ) γg + σ t (u h ) • (v h t -u t ) dΓ η e 1,Ω .
Note that the remaining terms in T 2 can be split as

T 2 =T C 2 + T T 2 ,
where T C 2 represents the contact terms and T T 2 contains the Tresca friction terms. The contact terms T C 2 are handled as in Theorem 3.5. Moreover, we can bound the friction terms T T 2 in a similar fashion, following step by step the proof of Theorem 3.5 and using the bound [START_REF] Ern | Theory and practice of finite elements[END_REF]. We get finally for any β > 0:

T T 2 ≤ 1 2β + (θ + 1) 2 4 γ 0 σ t (u -u h ) 2 -1/2,h,Γ C + βγ 0 θ 2 η 2 + γ 0 2β σ t (v h -u h ) 2 -1/2,h,Γ C .
Using [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: Il problema di Signorini con ambigue condizioni al contorno[END_REF] and the H 1 -stability of R h (see [START_REF] Bernardi | A local regularization operator for triangular and quadrilateral finite elements[END_REF] in Lemma 2.1) we bound:

σ t (v h -u h ) -1/2,h,Γ C ≤ C v h -u h 1,Ω = C R h (u -u h ) 1,Ω ≤ C u -u h 1,Ω .
We combine this last bound with the saturation assumption [START_REF] Fernández | A posteriori error analysis for the normal compliance problem[END_REF] and get (remembering that the same result holds for the contact terms T C 2 ):

T 2 ≤ Cγ 0 (θ + 1) 2 4 + 1 β u -u h 2 1,Ω + βγ 0 θ 2 η 2 .
From now on the proof is exactly the same as in Theorem 3.5.

Remark 4.5. An extension of Theorem 3.6 holds as well for Problem [START_REF] Fabre | A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method[END_REF] and similar local lower error bounds can be derived following the same method. The only difference is that the term η 3K is bounded as follows:

η 3K E∈E C K h 1/2 K σ t (u) + 1 γ P t γ (u h ) γg 0,E + σ t (u -u h ) 0,E
.

  • s,D (resp. | • | s,D ) and we keep the same notation for any d = 1, 2, 3. In the sequel the symbol | • | will either denote the Euclidean norm in R d , or the measure of a domain in R d . The letter C stands for a generic constant, independent of the discretization parameters.

Figure 1 :Figure 2 :Figure 3 :

 123 Figure 1: Rate of convergence for uniform and adaptive refinement methods. Parameters γ 0 = 1/E, θ = -1 and Lagrange P 2 elements.

θ = - 1 Figure 4 :

 14 Figure 4: First example. Convergence curves of the error estimator η, the L 2 and H 1 -norms of the error uu h , for γ 0 = 1/E.

θ = - 1 Figure 5 :

 15 Figure 5: First example. Convergence curves of the error estimator η, the L 2 and H 1 -norms of the error uu h , for γ 0 = 1000/E.

θ = - 1 Figure 6 :

 16 Figure 6: First example. Convergence curves of the error estimator η, the L 2 and H 1 -norms of the error uu h , for γ 0 = 10 6 /E.

θ = - 1 Figure 7 :

 17 Figure 7: First example. Convergence curves of the error estimator η, the L 2 and H 1 -norms of the error uu h , with γ 0 = 10 6 /E in the Nitsche's method and γ0 = 1/E in the error estimator η 4 .

Figure 8 :

 8 Figure 8: Reference solutions with von Mises stresses, in 2D (left) and 3D (right).

Figure 9 :

 9 Figure 9: Hertz's contact in 2D. Error estimator η and the H 1 -norm of the error uu h , for Lagrange P 1 (left) and P 2 (right) finite elements.

Figure 10 :

 10 Figure 10: Hertz's contact in 3D. Error estimator η and H 1 -norm of the error u-u h , for Lagrange P 1 (left) and P 2 (right) finite elements.

  as the set of interior edges/faces of T h (the edges/faces are supposed to be relatively open). We denote by E N h = {E ∈ E h : E ⊂ Γ N } the set of boundary edges/faces that correspond to Neumann conditions, and similarly E C h = {E ∈ E h : E ⊂ Γ C } is the set of boundary edges/faces included into the contact boundary.

Table 1 :

 1 First example, θ = 1 and γ 0 = 1/E.

	Mesh size h	1/4	1/8	1/16	1/32	1/64	1/80	slope
	Degrees of freedom	32	128	512	2048	8192	12800
	u -u h	0,Ω (×10 -4 )	104.7551 48.2436 17.3689	5.9666	2.0366	1.4262	1.4589
	u -u h	1,Ω (×10 -3 )	51.3896	28.8563 16.1335	9.0627	4.9777	4.1489	0.8412
		η 1	16719.8	8359.9	4179.95 2089.97 1044.99	835.99	1.0000
		η 2	60779.5	38076.7	22698	13222.3 7724.01 6507.89	0.7522
		η 3	7626.32	3209.18 1207.19 427.694 157.242 118.467	1.4107
		η 4	13501	4604.89 1395.58 370.912	100.73	77.2	1.7646
		η	64916.4	39385.6 23153.3 13398.4 7796.61 6562.89	0.7779
	Effectivity index Ef f E	1.2632	1.3649	1.4351	1.4784	1.5661	1.5816

Table 2 :

 2 First example, θ = 0 and γ 0 = 1/E.

	Mesh size h	1/4	1/8	1/16	1/32	1/64	1/80	slope
	Degrees of freedom	32	128	512	2048	8192	12800
	u -u h	0,Ω (×10 -4 )	113.6807 47.1350 17.0780	5.9262	2.0312	1.4229	1.4757
	u -u h	1,Ω (×10 -3 )	48.8181	28.0213 15.9877	9.0359	4.9716	4.1459	0.8251
		η 1	16719.8	8359.9	4179.95 2089.97 1044.99	835.99	1.0000
		η 2	57305.3	37374.7 22547.2 13200.7 7720.86 6505.24	0.7356
		η 3	3938.22	1852.35 720.951 256.135 95.0474	71.047	1.3686
		η 4	11946.5	4002.56 1154.11 324.915 89.6552	61.026	1.7809
		η	61005.6	38551.4 22971.7 13371.5 7792.35	6559.4	0.7779
	Effectivity index Ef f E	1.2496	1.3758	1.4368	1.4798	1.5672	1.5819

Table 3 :

 3 First example, θ = -1 and γ 0 = 1/E.

	Mesh size h	1/4	1/8	1/16	1/32	1/64	1/80	slope
	Degrees of freedom	32	128	512	2048	8192	12800
	u -u h	0,Ω (×10 -4 )	120.9371 48.9718 17.3613	5.9619	2.0360	1.4255	1.4952
	u -u h	1,Ω (×10 -3 )	49.3705	28.1269 16.0087	9.0385	4.9714	4.1467	0.8283
		η 1	16719.8	8359.9	4179.95 2089.97 1044.99	835.99	1.0000
		η 2	58846.3	37649.9 22607.7 13213.2 7723.58 6506.99	0.7428
		η 3	2690.5	1464.81 558.637 192.194 70.7559 53.7733	1.3544
		η 4	9202.06	2854.93 832.228 229.683	62.842	44.0949	1.8004
		η	61922.2	38700.1 23012.7 13380.8 7794.52 6560.84	0.7779
	Effectivity index Ef f E	1.2542	1.3759	1.4375	1.4804	1.5677	1.5820

Table 4 :

 4 First example, θ = 1 and γ 0 = 10 6 /E.

	Mesh size h	1/4	1/8	1/16	1/32	1/64	1/80	slope
	Degrees of freedom	32	128	512	2048	8192	12800
	u -u h	0,Ω (×10 -4 )	122.6500	58.4959	46.4511	16.6143	3.0112	2.0808	1.3628
	u -u h	1,Ω (×10 -3 )	57.7770	30.5558	30.8275	15.1381	7.5190	7.2669	0.6963
		η 1	16719.8	8359.9	4179.95	2089.97	1044.99	835.99	1.0000
		η 2	62073.4	38335.4	41033.3	22552.4	10916.7	9635.48	0.6172
	η 3 (×10 -2 )	2.13709	1.01961	0.63268 0.768462 0.414401	0.42499	0.4878
	η 4 (×10 -2 )	2.52415	0.571842 2.55605	1.38494	0.525521 0.446467	0.4177
		η	64285.7	39236.3	41245.6	22649	10966.6	9671.68	0.6272
	Effectivity index Ef f E	1.1127	1.2841	1.3380	1.4962	1.4564	1.3266

Table 5 :

 5 First example, θ = 0 and γ 0 = 10 6 /E.

	Mesh size h	1/4	1/8	1/16	1/32	1/64	1/80	slope
	Degrees of freedom	32	128	512	2048	8192	12800
	u -u h	0,Ω (×10 -4 )	321.9223 518.2042 592.2560 635.8190 615.1016 577.5225	-0.1689
	u -u h	1,Ω (×10 -3 )	130.6775 175.6706 192.4574 203.4624 197.2235 186.4225	-0.1058
		η 1	16719.8	8359.9	4179.95	2089.97	1044.99	835.99	1.0000
		η 2	75562.4	49342.9	29582.9	17606.7	10284.5	8546.95	0.7339
		η 3	2468.7	908.957	316.586	107.411	38.2804	28.2145	1.5045
	η 4 (×10 -1 )	2.00548	3.65123	5.58251	8.52101	1.38006	1.37438	-0.6465
		η	77429.5	50054.3	29878.5	17730.7	10337.5	8587.78	0.7399
	Effectivity index Ef f E	5.9252	2.8493	1.5525	0.8714	0.5242	0.4607
	(×10 -1 )

Table 6 :

 6 First example, θ = -1 and γ 0 = 10 6 /E.

	Mesh size h	1/4	1/8	1/16	1/32	1/64	1/80	slope
	Degrees of freedom	32	128	512	2048	8192	12800
	u -u h	0,Ω (×10 -4 )	110.5852 47.6266 16.9809	5.9093	2.0290	1.4216	1.4709
	u -u h	1,Ω (×10 -3 )	50.6403	29.1195 16.2386	9.0861	4.9803	4.1565	0.8386
		η 1	16719.8	8359.9	4179.95 2089.97 1044.99	835.99	1.0000
		η 2	62292.2	38204	22809.5 13249.4 7732.01	6512.98	0.7582
	η 3 (×10 -4 )	143.671	83.7405 26.7592 8.72031 3.05775 0.0236947	1.4405
	η 4 (×10 -4 )	168.808	67.5774 14.2866 3.96942 1.17445	0.948677	1.8030
		η	64497	39108	23189.4 13413.2 7802.31	6566.42	0.7666
	Effectivity index Ef f E	1.2736	1.3430	1.4280	1.4762	1.5663	1.5796

Table 7 :

 7 Hertz's contact in 2D, θ = -1, γ 0 = 10 -3 /E and Lagrange P 1 elements

	Mesh size h	6.04766 5.23002	2.7327	1.64637 0.482414 0.246359	slope
	u -u h	1,Ω	50.5984 12.3996	4.1083	1.5561	0.3399	0.1534	1.6392
	η 1		7781.84 7376.92 4066.11 2379.44	728.236	359.817	0.9715
	η 2		18000.7 12350.7 9279.79 5866.83	2009.88	1029.86	0.8525
	η 3		2523.15 1055.64 852.542 458.121	90.2956	38.1934	1.2132
	η 4		21999.5 10276.6 2537.53 1735.77	321.871	152.501	1.4597
	η		29579.2 17711.1 10479.2 6580.59	2163.73	1102.18	0.9643
	Effectivity index Ef f E	0.2338	0.5713	1.0203	1.6916	2.5467	2.8735

  holds as well. Then we have+ σ n (u) -σ n (u h ) -1/2,h,Γ C + σ t (u) -σ t (u h ) -1/2,h,Γ C (1 + γ 0 )η + ζ.

	u -u h	1,Ω + σ n (u) +	1 γ	[P γ (u h )] +	-1/2,h,Γ C	+ σ t (u) +	1 γ	P t γ (u h )	γg -1/2,h,Γ C

Acknowledgements

We thank Roland Becker, Daniela Capatina and Robert Luce for helpful discussions and comments.

Appendix 1. Extension to the Tresca friction case

We extend in this appendix the analysis of Section 3 to the case of unilateral contact with Tresca friction. About a posteriori error estimates for the Tresca friction problem one may refer, e.g., to

Setting and Nitsche-based finite element method for Tresca

Let g ∈ L 2 (Γ C ), g ≥ 0 be a given threshold. The Tresca friction problem with unilateral contact consists in finding the displacement field u : Ω → R d verifying the equations and conditions (1)-(2 (i,ii,iii))- [START_REF] Eck | A residual-based error estimator for BEM-discretizations of contact problems[END_REF], with [START_REF] Eck | A residual-based error estimator for BEM-discretizations of contact problems[END_REF] given by:

where

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 α. The following property holds for all x, y ∈ R d-1 :

where • is the euclidean scalar product in R d-1 .

Let us introduce the discrete linear operator

) and the bilinear form:

The extension of our Nitsche-based method for unilateral contact with Tresca friction then reads:

Consistency, well-posedness and a priori error estimates for the method (23) are established in [START_REF] Chouly | An adaptation of Nitsche's method to the Tresca friction problem[END_REF].

Residual error estimator, upper and lower bound

The Definition 3.2 still holds for problem [START_REF] Fabre | A fictitious domain method for frictionless contact problems in elasticity using Nitsche's method[END_REF], except for η 3K whose expression is now:

First we provide counterparts of the Assumption 3.3 and of the discrete trace inequality of Lemma 2.3. Assumption 4.2. The solution u of (1)-(2 (i,ii,iii))-( 21) and the discrete solution u h of ( 23) are such that: