Z Belhachmi

F Hecht

ADAPTIVE FINITE ELEMENTS FOR VARIATIONAL INEQUALITIES WITH NON-SMOOTH COEFFICIENTS

Keywords: 1991 Mathematics Subject Classification. 35J85, 65N30, 74M15 Adaptive finite elements method, unilateral cracks in contact mechanics, linear elasticity, a posteriori analysis

We consider an elliptic variational inequality with discontinuous coefficients arising in unilateral contact mechanics in linearized elasticity. The contact zone is an internal boundary separating sub-domains with difference elastic properties. We study some discrete formulations with mixed finite element methods and we give optimal error estimates in appropriate norms, independent of the variation of the elasticity coefficients. The focus of the article is the a posteriori analysis with residual error indicators. It is achieved both for the conforming and nonconforming discretization, in a unified framework. The residual error indicators are well suited to handle non-matching meshes and the contact conditions, and they allow us to obtain sharp and robust a posteriori estimates. An adaptive solution algorithm is proposed and few numerical experiments confirming the theory are presented.

Introduction

The numerical implementation of contact and impact problems in solid mechanics is usually based on finite element tools [START_REF] Zhong | Finite Element Procedures for Contact-Impact Problems[END_REF][START_REF] Haslinger | Nečas-Numerical Methods for Unilateral Problems in Solid Mechanics[END_REF][START_REF] Kikuchi | Oden--Contact problems in elasticity: A study of variational inequalities and finite element methods[END_REF]. The choice of the finite element methods which are both easy to implement, accurate from the theoretical point of view and of low cost is crucial for these simulations. Such choice requires the use of a priori and a posteriori analysis tools to design efficient discretization strategies.

There has been a lot of progress in the numerical solution of such problems since the pioneering works in the 70's, see [START_REF] Glowinski | Lectures on numerical methods for nonlinear variational problems[END_REF][START_REF] Haslinger | Nečas-Numerical Methods for Unilateral Problems in Solid Mechanics[END_REF][START_REF] Kikuchi | Oden--Contact problems in elasticity: A study of variational inequalities and finite element methods[END_REF][START_REF] Zhong | Finite Element Procedures for Contact-Impact Problems[END_REF] for an exhaustive bibliography. Working with linear (and quadratic) finite elements, various discrete formulations, depending on the modelling of the contact conditions at discrete level, were addressed in many studies [START_REF] Haslinger | Nečas-Numerical Methods for Unilateral Problems in Solid Mechanics[END_REF][START_REF] Brezzi | Raviart-Error estimates for the finite element solution of variational inequalities[END_REF][START_REF] Brezzi | Raviart-Error estimates for the finite element solution of variational inequalities, Part 2: Mixed methods[END_REF][START_REF] Kikuchi | Oden--Contact problems in elasticity: A study of variational inequalities and finite element methods[END_REF][START_REF] Zhong | Finite Element Procedures for Contact-Impact Problems[END_REF][START_REF] Falk | Error estimates for the approximation of a class of variational inequalities[END_REF][START_REF] Lhalouani | Sassi-Nonconforming mixed variational formulation and domain decomposition for unilateral problems[END_REF][START_REF] Coorevits | Sassi-Mixed finite element methods for unilateral problems: convergence analysis and numerical studies[END_REF][START_REF] Ben Belgacem | Renard-Hybrid finite element methods for the Signorini problem[END_REF][START_REF] Belhachmi | Ben Belgacem-Quadratic finite element for Signorini problem[END_REF][START_REF] Slimane | Laborde-Mixed formulations for a class of variational inequalities[END_REF][START_REF] Hüeber | Wohlmuth-An optimal error estimate for nonlinear contact problems[END_REF]. In particular, a priori error estimates and numerical algorithms for solving such variational inequalities have been extensively studied. On the contrary, the a posteriori analysis and adaptive strategies have not been sufficiently developed, particularly for problems of Signorini type. In fact, most of the existing studies are devoted to obstacle problems where various error estimators are studied. We refer the reader to [START_REF] Hoppe | Adaptive multilevel methods for obstacle problems[END_REF][START_REF] Ainsworth | Lee-Local a posteriori error estimators for variational inequalities[END_REF][START_REF] Johnson | Adaptive finite element methods for obstacle problems[END_REF][START_REF] Chen | Nochetto-Residual type a posteriori error estimates for elliptic obstacle problems[END_REF] and the references therein.

For frictionless unilateral contact problems, the residual based method using a penalized approach is considered in [START_REF] Carstensen | Wriggers-Adaptive finite elements for elastic bodies in contact[END_REF] (see also the references therein). The study of the error in the constitutive law is performed in [START_REF] Coorevits | Pelle-A posteriori error estimators for unilateral contact with matching and nonmatching meshes[END_REF], an error indicator based on equilibrated fluxes in [START_REF] Wohlmuth | An a posteriori error estimator for two-body contact problems on nonmatching meshes[END_REF], and a residual error estimator for the Signorini problem is considered in [START_REF] Hild | Nicaise-A posteriori error estimations of residual type for Signorini problems[END_REF][START_REF] Hild | Nicaise-Residual a posteriori error estimators for contact problems in elasticity[END_REF].

Most previous works dedicated to unilateral contact problems consider a contact between a deformable and a rigid body or the contact between two deformable bodies, in homogeneous linearized elasticity. Our aim in this paper is to consider contact problems for nonhomogeneous linearized elasticity that rely on the study of a variational inequality with discontinuous coefficients: the elasticity coefficients are smooth on a finite number of subdomains. These models describe, for example, the unilateral crack propagation in dissimilar media or the delamination in composite materials (see [2]).

It is our goal to prove a priori and a posteriori estimates that are independent of the size of the jumps in the coefficients across the interfaces between the subdomains. We consider both the conforming and nonconforming discretization, and we perform their analysis in a unified framework. In particular we introduce new residual error estimators that take into account the contact conditions and the non-matching grids at the interfaces of the sub-domains. The residual error indicators are easily computable and we obtain sharp estimates without any saturation assumption or any extra regularity assumption. The scaling factors of the error estimators allow to handle correctly the local ratio between adjacent coefficients. As the meshes should be aligned with the discontinuities, we build a new quasi-interpolation operator of Clément type [START_REF] Clément | Approximation by finite element functions using local regularization[END_REF]. This operator and the use of appropriate norms allow us to obtain estimates for the interpolation error that are independent of the size of the jumps of the coefficients. The a posteriori analysis presented in the article covers, as largely as possible, the most realistic model of frictionless unilateral contact in the framework of nonhomogeneous linear elasticity and could be extended straightforwardly to the anisotropic elasticity. The numerical experiments show the convergence of the adaptive strategy and are in accordance with the theoretical results.

The outline of the paper is as follows. We consider the variational formulations and prove well posedness results in Section 2. In section 3, we describe the discrete problems, we prove a priori estimates and we perform the convergence analysis. In Section 4, we perform the a posteriori analysis, we introduce the residual error indicators and prove upper and lower bounds for the error in an appropriate norm. The details of the implementation and some numerical results are given in Section 5. The appendix, is devoted to the construction of the appropriate quasi-interpolant operator.

Variational formulation and discrete problems

Let Ω be a bounded domain of R 2 with smooth boundary Γ, and (γ i ) i , 1 ≤ i ≤ I, a given number of Lipschitz continuous curves in Ω without self-intersections, such that Ω \ I i=1 γ i is connected. We assume given a family Ω j , 1 ≤ j ≤ J of pairwise disjoint Lipschitzian open subsets which constitutes a non overlapping decomposition of Ω, and the family (γ i ) i , 1 ≤ i ≤ I is made of (parts) of interfaces of such decomposition.

(1) Ω = J j=1

Ω j ∪ ( I i=1 γ i )
Let α be a function which is equal to the constant α j on each Ω j , 1 ≤ j ≤ J (see figure 1). We define the two parameters

α min = min 1≤j≤J α j α max = max 1≤j≤J α j ,
and we assume α min positive.

1 Ω j2

Ω j1 αj 1 j 1 α j2
-div σ = f in Ω C , (2)
σ = A α ε(u) in Ω C , (3)
u = 0 on Γ , (4)
[u] n ≤ 0, [σ n ] = 0, σ n [u] n = 0 on γ i , 1 ≤ i ≤ I (5) σ n ≤ 0, σ τ = 0 on γ ± i , 1 ≤ i ≤ I. (6)
Here [u] = u+u -denotes the jump in the displacement field across Γ C , and the signs ± indicate the positive and negative directions with respect to the external normal n. The unknown is the displacement field u = (u 1 , u 2 ). The symmetric stress tensor σ = (σ ij ), i, j = 1, 2 is linked to the displacement by Hooke's law [START_REF]Belhachmi-A posteriori error estimates for the 3D stabilized mortar finite element method applied to the Laplace equation[END_REF] where ε ij (u) = 1 2 (u i,j + u j,i ) denotes the symmetric strain tensor. We set

A α = α j A in Ω j , 1 ≤ j ≤ J,
where A is the elasticity tensor. The body is subjected to volume forces f. We have used the following standard notation

σn = {σ ij n j } 2 i=1 , σ n = σ ij n j n i , σ τ = σn -σ n n = σ i τ 2 i=1 , {A α ε(u)} ij = A α ijkℓ ε kℓ , A α ijkℓ = A α jikℓ = A α kℓij , A α ijkℓ ∈ L ∞ (Ω)
. The fourth order tensor A α satisfies the ellipticity condition (7)

A α ijkℓ ξ ji ξ kℓ ≥ a 0 |ξ| 2 , ∀ξ ji = ξ ij , a 0 > 0. We use the summation convention over repeated indices.

2.1.

Variational formulation and well-posedness. We assume that the data f are in L 2 (Ω C ) and we introduce the Sobolev space

H 1 Γ (Ω C ) 2 = v = (v 1 , v 2 ) ∈ H 1 (Ω C ) 2 ; v i = 0 on Γ, i = 1, 2
, and the closed convex set (of admissible displacements) ( 8)

K C = v = (v 1 , v 2 ) ∈ H 1 Γ (Ω C ) 2 ; [v] n ≤ 0 a.e on Γ C .
Problem ( 2)-( 6) admits the following equivalent variational form:

(9) find u ∈ K C such that, Ω C σ(u) : ε(v -u) dx ≥ Ω C f(v -u) dx, ∀v ∈ K C .
Due to the boundedness of α, the bilinear form on the right-hand side is continuous, and the positivity of α min ensures its coercivity thanks to the Korn inequality which holds in Ω C . The following result is a consequence of Stampacchia theorem [START_REF] Haslinger | Nečas-Numerical Methods for Unilateral Problems in Solid Mechanics[END_REF] Proposition 2.1.

For any data f ∈ L 2 (Ω C ) 2 , problem (9) admits a unique solution u ∈ K C .
The equivalence between problem (2)-( 6) and the variational formulation [START_REF] Ben Belgacem | Renard-Hybrid finite element methods for the Signorini problem[END_REF] as well as the precise mathematical meaning of the boundary, respectively the contact condition on Γ C , requires some care. Let Σ i be one of ∂Ω j1,i or ∂Ω j2,i (the subdomains that share the interface γ i , and that we denote for brevity Ω i ). The space H 1 2 (Σ i ) for Lipschitz closed curves Σ i is defined in [START_REF] Girault | Raviart-Finite element methods for the Navier-Stokes equations, Theory and algorithms[END_REF]Chapter1,relation (1.16)]. For the convenience of the reader, we recall that its norm is

v 2 H 1 2 (Σi) = v 2 L 2 (Σi) + Σi Σi |v(x) -v(y)| 2 |x -y| 2 ds(x) ds(y).
We denote by

H -1 2 (Σ i ) its dual. For σ ∈ L 2 (Ω i ) 4 , div σ ∈ L 2 (Ω i ) 2
, the traces (σ.n i ) ± can be defined as elements of H -1 2 (Σ i ), and the trace operator with values in H -1 2 (Σ i ) is continuous [START_REF] Girault | Raviart-Finite element methods for the Navier-Stokes equations, Theory and algorithms[END_REF]. We denote by H 1 2 00 (γ i ) the subspace of H 1 2 (Σ i ) consisting of functions with support in γi . Note that the definition of H 1 2 00 (γ i ) does not depend on the choice of Σ i being the boundary of Ω j1 or Ω j2 . This space is defined for smooth γ i in [42, Chapter 1]. When γ i are smooth enough (say C 1,1 ), we obtain, by standard arguments, the precise interpretation of conditions ( 5)-( 6):

for each i, 1 ≤ i ≤ I, (10) < σ n (u) --σ n (u) + , ϕ > γi = 0, ∀ϕ ∈ H 1 2 00 (γ i ), ( 11
) < σ n (u), ϕ > γi ≤ 0, ∀ϕ ∈ H 1 2 00 (γ i ), ϕ ≥ 0, ( 12
) < σ τ (u), ϕ > γi = 0, ∀ϕ = (ϕ 1 , ϕ 2 ) ∈ H 1 2 00 (γ i ), ϕ • n = 0, and (13) < σ n (u) [u] n, ϕ > γi = 0, ∀ϕ ∈ H 1 2 00 (γ i ),
where < ., . > γi denotes the duality product between H 1 2 00 (γ i ) and its dual H -1 2 (γ i ) (we will not distinguish between scalar and vector valued cases).

2.2. Hybrid variational formulation. We introduce the spaces

V(Ω j ) = v ∈ H 1 (Ω j ), v = 0 on Γ ∩ ∂Ω j , 1 ≤ j ≤ J,
and the space

V = v ∈ L 2 (Ω); v |Ωj ∈ V(Ω j ) ∩ H 1 (Ω C ).
We define the Lagrange multiplier convex set

M = µ = (µ i ) ∈ I i=1 H -1 2 (γ i ), I i=1 < µ i , ψ i > γi ≥ 0, ∀ψ = (ψ i ) ∈ I i=1 H 1 2 00 (γ i ), ψ i ≥ 0, 1 ≤ i ≤ I .
We denote

µ -1 2 , * = ( I i=1 µ i 2 H -1 2 (γi) ) 1 2
.

In order to obtain error estimates which are independent of α, we will work with the norm ( 14)

v α = J j=1 |α 1 2 j v| 2 H 1 (Ωj ) 2 + v 2 L 2 (Ωj ) 2 1 2 .
We introduce the following notations :

for u = (u j ) 1≤j≤J , v = (v j ) 1≤j≤J in V and µ in I i=1 H -1 2 (γ i ) a α (u, v) = J j=1 Ωj (A α ) mnkh ε j mn (u j ) ε j kh (v j ) dx, (f , v) = J j=1 Ωj f j • v j dx, and
b(µ, v) = I i=1 < µ i , v j1,i • n j1,i + v j2,i • n j2,i > γi = I i=1 < µ, [v] n i > γi .
The hybrid variational formulation of problem (9) consists of finding u ∈ V and λ ∈ M , such that

(15) a α (u, v) + b(λ, v) = (f , v), ∀v ∈ V, b(µ -λ, u) ≤ 0, ∀µ ∈ M.
The existence and uniqueness of the solutions of problem [START_REF] Brezzi | Raviart-Error estimates for the finite element solution of variational inequalities, Part 2: Mixed methods[END_REF] follow, in a standard way, from the ellipticity of a α (., .) and the usual Brezzi-Babuska inf-sup condition on b(., .) (see [START_REF] Haslinger | Nečas-Numerical Methods for Unilateral Problems in Solid Mechanics[END_REF], III. Theorem 9.4) Proposition 2.2. Problem (15) admits a unique solution (u, λ) ∈ V × M . Moreover, we have

(16) λ = (λ i ) 1≤i≤I = -(σ ℓm ) j1,i n j1,i ℓ n j1,i m = -(σ ℓm ) j2,i n j2,i ℓ n j2,i m .

The discrete problems

We assume that each of the subdomains Ω j , 1 ≤ j ≤ J is polygonal and the interfaces γ i , 1 ≤ i ≤ I are straight lines . More general cases require additional and non-essential technicalities that we omit for simplicity. We consider regular families (T j h ) h of partitions of Ω j , 1 ≤ j ≤ J, into a finite number of triangles which satisfy the usual admissibility conditions [START_REF]Ciarlet-Basic Error Estimates for Elliptic Problems[END_REF].

We denote by h j = max K∈T j h h K the discretization parameter on Ω j , and h = max(h j ) 1≤j≤J . We also assume that the endpoints of each γ i , c i 1 and c i 2 are common vertices of the triangulations T j1,i h and T j2,i h , and the traces of the triangulations T j1,i h and T j2,i h on γ i are one-dimensional triangulations that we denote by T ℓ h,γi , ℓ = 1, 2, 1 ≤ i ≤ I. The set of vertices of T ℓ h,γi , ℓ = 1, 2 is denoted by

ζ ℓ h,i = c i 1 = x ℓ 0,i , x ℓ 1,i , . . . , x ℓ m ℓ ,i = c i 2 ,
and their elements are

t ℓ k,i = x ℓ k,i , x ℓ k+1,i , 0 ≤ k ≤ m ℓ -1.
We consider the (affine) finite element spaces

V h (Ω j ) = v j h ∈ C(Ω j ), ∀K ∈ T j h v j h|K ∈ P 1 (K) 2
, v j h = 0 on Γ , and the space

V h = v h ∈ L 2 (Ω); v h|Ωj ∈ V h (Ω j ) ∩ H 1 (Ω C ). Note that V h ⊂ V.
For the approximation of the Lagrange multipliers, for each 1 ≤ i ≤ I, and ℓ = 1, 2 we introduce the affine, respectively piecewise, finite element spaces

W 1,ℓ h (γ i ) = µ h ∈ C(γ i ), ∃v h ∈ V h such that v h|Ωj ℓ(i) • n j ℓ(i) = µ h on γ i , = µ h ∈ C(γ i ), µ h|t k,i ∈ P 1 (t k,i ), 1 ≤ k ≤ m ℓ -2 µ h|t0,i ∈ P 0 (t 0,i ), µ h|tm ℓ -1,i ∈ P 0 (t m ℓ -1,i ) ,
and the discrete convex cones

M 1,ℓ h (γ i ) = µ h ∈ W 1,ℓ h (γ i ), µ h ≥ 0 on γ i , M 1,ℓ, * h (γ i ) = µ h ∈ W 1,ℓ h (γ i ), γi µ h ψ h ≥ 0 ∀ψ h ∈ M 1,ℓ h (γ i ) .
These two convex cones are commonly used to express the jump conditions on the contact zone. The choice of one is linked to the way one selects to enforce the nonnegativity conditions (either on the displacement field or on the normal component of the stress tensor). We set

M h (γ i ) = M 1,ℓ h (γ i ) or M 1,ℓ, * h (γ i ), ℓ = 1,

and we define

W 1 h = I i=1 W 1,ℓ h (γ i ), M h = I i=1 M h (γ i ). Note that if M h (γ i ) = M 1,ℓ h (γ i ) then M h ⊂ M , M 1,ℓ h (γ i ) ⊂ M 1,ℓ, * h (γ i ) but if M h (γ i ) = M 1,ℓ, * h then M h ⊂ M . The discrete problem reads : find (u h , λ h ) ∈ V h × M h such that (17) a α (u h , v h ) + b(λ h , v h ) = (f , v h ), ∀v h ∈ V h , b(µ h -λ h , u h ) ≤ 0, ∀µ h ∈ M h .
The V h -ellipticity of the bilinear form holds thanks to the Korn inequality, still valid in the case of unilateral contact cracks. It is also readily checked that

{µ h ∈ M h , b(µ h , v h ) = 0, ∀v h ∈ V h } = {0} .
Therefore, the next proposition follows from standard saddle-point theory in the finite dimensional setting.

Proposition 3.1. With M h based on the local choices M h (γ i ) = M 1,ℓ h (γ i ), or M 1,ℓ, * h (γ i ), ℓ = 1, 2, 1 ≤ i ≤ I,
there exists a unique solution (u h , λ h ) of the discrete problem [START_REF] Chen | Nochetto-Residual type a posteriori error estimates for elliptic obstacle problems[END_REF].

In all what follows we denote by C all the positives constants independent of α and h.

In order to obtain optimal approximation properties, we need the bilinear form b(., .) to satisfy a uniform inf-sup condition (with respect to h). This condition requires some assumptions on the triangulations on each γ i , 1 ≤ i ≤ I. The simplest and sufficient one is the quasi-uniformity, however, we do not make this assumption which is too stringent especially for adaptive mesh refinement which is one of the objectives of this work. Instead, we will make the following assumption: the 1D triangulations T ℓ h,γi , ℓ = 1, 2, satisfy the Crouzeix-Thomée criterion [START_REF] Crouzeix | Thomée-Resolvent estimates in lp for discrete Laplacians on irregular meshes and maximum-norm stability of parabolic finite difference schemes[END_REF] (18)

|t ℓ k,i | |t ℓ k ′ ,j | ≤ Cβ |k-k ′ | , ∀k, k ′ (0 ≤ k, k ′ ≤ m ℓ -1),
where 1 ≤ β ≤ 4.

Remark 3.2. Meshes satisfying Assumption 1 allows the adaptive mesh refinement unlike the quasi-uniform meshes [START_REF] Ben Belgacem | Renard-Hybrid finite element methods for the Signorini problem[END_REF] We recall some approximation tools which will be used in the following analysis. Let R j h , and r i h , be the Lagrange interpolation operators with values in

V h (Ω j ) (1 ≤ j ≤ J) and W 1,ℓ h (γ i ) (ℓ = 1, 2), 1 ≤ i ≤ I, respectively. There exists a constant C > 0, such that ∀v ∈ (H 2 (Ω j )) 2 and v ∈ H 3 2 (γ i ) ([18]) the following estimates hold (19) v -R j h v (H 1 (Ωj )) 2 ≤ Ch v (H 2 (Ωj )) 2 and v -r i h v L 2 (γi) ≤ Ch 3 2 v H 3 2 (γi)
.

Let us denote by γ any γ i , 1 ≤ i ≤ I. We define the projection operator π 1 h :

L 2 (γ) → W 1,ℓ
h (γ), with respect to the scalar product in L 2 (γ), which satisfies the following properties (see [START_REF] Bernardi | Patera-A new nonconforming approach to domain decomposition: the mortar element method[END_REF], [START_REF] Ben Belgacem | Renard-Hybrid finite element methods for the Signorini problem[END_REF]). Given µ ∈ [0, 1] and ν ∈ 1 2 , 2 , there exists a constant c > 0 which is independent of h, such that for all functions ϕ ∈ H ν (γ), [START_REF] Coorevits | Sassi-Mixed finite element methods for unilateral problems: convergence analysis and numerical studies[END_REF] ϕ -

π 1 h ϕ H -µ (γ) + h µ+ 1 2 ϕ -π 1 h ϕ H 1 2 (γ) ≤ ch µ+ν ϕ H ν (γ) . Proposition 3.3.
Under assumption (18) on the triangulations T ℓ h,γi , ℓ = 1 or 2, 1 ≤ i ≤ I, the following inf-sup condition holds [START_REF] Coorevits | Pelle-A posteriori error estimators for unilateral contact with matching and nonmatching meshes[END_REF] inf

µ h ∈W 1 h sup v h ∈V h b(µ h , v h ) µ h H -1 2 v h α ≥ δ ( α m α M ) 1 2 .
The constant δ is independent of h and α.

Proof. Let µ h ∈ W 1,ℓ h . We want to construct v h ∈ V h satisfying (22) b(µ h , v h ) ≥ c 1 (α) µ h 2 -1 2 , * and c 2 (α) v h α ≤ µ h -1 2 , * , where c ℓ (α), ℓ = 1, 2 are constants depending, on α but not on h. Let us consider v, the solution of the problem    div α(x)∇ v = 0, in Ω C α(x)∇ v • n = 0, on Γ α(x)∇ v • n ± = µ h,i n ± , on γ i , 1 ≤ i ≤ I.
The existence of such a v comes from the direct method of the calculus of variations. It is clear that v is also the solution of the equivalent variational problem

Ω C α(x) ∇ v ∇ w dx = b(µ h , w), ∀w ∈ H 1 (Ω C ) 2 ,
and satisfies, by the usual stability inequality and the continuous inf-sup condition on b(., .),

(23) c -µ h -1 2 , * ≤ v α ≤ c + µ h -1 2 , *
, where c + and c -are constants independent of h and α.

We set

v h ∈ V h such that v h|γi n = π 1 h ( v |γi n ), 1 ≤ i ≤ I and (24) v h α ≤ c α 1 2 M v h H 1 (Ω C ) ≤ c α 1 2 M I i=1 π 1 h v |γi n H -1 2 (γi) ≤ c α 1 2 M I i=1 v |γi n H -1 2 (γi) ≤ c α 1 2 M α -1 2 m v α .
Such a v h is built using a stable finite element extension operator similar to the standard local regularization operator studied in [START_REF] Bernardi | Girault-A local regularization operator for triangular and quadrilateral finite elements[END_REF]. Next, we note that

b(µ h , v h ) = b(µ h , v) = v 2 α ≥ (c -) 2 µ h 2 - 1
2 , * , and thanks to ( 23)-( 24) and the trace theorem we have the second statement of (22).

A Remark on the use of piecewise constant Lagrange multipliers.

The discrete Lagrange multiplier space can also be defined with piecewise constants instead of the one dimensional affine finite elements: for 1

≤ i ≤ I, W 0,ℓ h (γ i ) = µ h , µ h|t ℓ k,i ∈ P 0 (t ℓ k,i ), 0 ≤ k ≤ m ℓ,i -1 , which yields the discrete convex cone M 0,ℓ h (γ i ) = µ h ∈ W 0,ℓ h (γ i ), µ h,i ≥ 0 . Denoting by M 0 h = I i=1 M 0,ℓ h (γ i ), and W 0 h = I i=1 W 0,ℓ h (γ i ),
we have M 0 h ⊂ M and the resulting discrete problem also admits a unique solution (u h , λ h ). However, this choice in the present case leads to non optimal approximation properties because of the presence of spurious modes. Following [START_REF] Ben Belgacem | Renard-Hybrid finite element methods for the Signorini problem[END_REF], we use, in this case, a stabilization technique by adding two bubble functions on each γ i , 1

≤ i ≤ I ϕ t ℓ k,i (x) = 6 |t ℓ k,i | λ 1 (x)λ 2 (x), ∀x ∈ K k , i = 0, m ℓ,i -1,
where K k , is the triangle having t ℓ k,i , k = 0 or m ℓ,i -1 as an edge and λ 1 , λ 2 the barycentric coordinates associated to the vertices of t ℓ k,i . We then replace V h by

V h = V h ⊕ (⊕ I i=1 ⊕ t ℓ k,i ∈T ℓ h,γ i Rϕ t ℓ k,i ). We denote by π 0 h the L 2 -projection operator L 2 (γ) -→ W 0,ℓ h (γ)
, defined as follows:

(25)

γ vψ h dσ = γ π 0 h (v)ψ h dσ, ∀ψ h ∈ W 0,ℓ h (γ),
where π 0 h satisfies the following estimates (see [START_REF] Ben Belgacem | Renard-Hybrid finite element methods for the Signorini problem[END_REF]). Namely, for the functions ϕ ∈ H ν (γ), with ν = 1 2 , or with ν = 1, there exists a constant c > 0 independent of h such that [START_REF] Girault | Raviart-Finite element methods for the Navier-Stokes equations, Theory and algorithms[END_REF] ϕ

-π 0 h ϕ L 2 (γ) ≤ ch ν ϕ H ν (γ) . Moreover, if ϕ ∈ L 2 (γ), then (27) ϕ -π 0 h ϕ H -1 2 (γ) ≤ ch 1 2 ϕ -π 0 h ϕ L 2 (γ) .
The uniform inf-sup condition is obtained similarly to proposition 3.3

Proposition 3.4. Under assumption (18) on the triangulations T ℓ h,γi , ℓ = 1 or 2, 1 ≤ i ≤ I, the following inf-sup condition holds [START_REF] Grisvard | Elliptic Problems in Nonsmooth Domains[END_REF] inf

µ h ∈W 0 h sup v h ∈ V h b(µ h , v h ) µ h -1 2 , * v h α ≥ δ ( α m α M ) 1 2 > 0,
where the constant δ is independent of h or α.

Convergence Analysis.

The convergence analysis and the a priori error estimates come from the standard approximation theory for unilateral contact problems [START_REF] Falk | Error estimates for the approximation of a class of variational inequalities[END_REF][START_REF] Haslinger | Nečas-Numerical Methods for Unilateral Problems in Solid Mechanics[END_REF][START_REF] Coorevits | Sassi-Mixed finite element methods for unilateral problems: convergence analysis and numerical studies[END_REF][START_REF] Lhalouani | Sassi-Nonconforming mixed variational formulation and domain decomposition for unilateral problems[END_REF][START_REF] Ben Belgacem | Renard-Hybrid finite element methods for the Signorini problem[END_REF][START_REF] Ben | Belgacem-Numerical simulation of some variational inequalities arisen from unilateral contact problems by finite element method[END_REF][START_REF] Slimane | Laborde-Mixed formulations for a class of variational inequalities[END_REF][START_REF] Hüeber | Wohlmuth-An optimal error estimate for nonlinear contact problems[END_REF]. For the specific case of unilateral cracks we refer the reader to [START_REF] Tahir | Problèmes de contact unilatéral et maillages incompatibles[END_REF][START_REF] Belhachmi | Sokolowski-Mixed finite element methods for a smooth domain formulation of a crack problem[END_REF][START_REF] Tahir | Belhachmi-Mixed finite elements discretizations of some variational inequalities arising in elasticity problems in domains with cracks[END_REF]. We give a brief summary for the reader convenience, we skip the proofs which are rather long but standard adaptations of arguments in [START_REF] Coorevits | Sassi-Mixed finite element methods for unilateral problems: convergence analysis and numerical studies[END_REF][START_REF] Ben Belgacem | Renard-Hybrid finite element methods for the Signorini problem[END_REF][START_REF] Hüeber | Wohlmuth-An optimal error estimate for nonlinear contact problems[END_REF] for the homogeneous case. We start with the following lemmas, Lemma 3.5. Let (u, λ) be the solution of problem [START_REF] Brezzi | Raviart-Error estimates for the finite element solution of variational inequalities, Part 2: Mixed methods[END_REF] and (u h , λ h ) the solution of problem [START_REF] Chen | Nochetto-Residual type a posteriori error estimates for elliptic obstacle problems[END_REF]. Then the following estimate holds: for any

(v h , µ h ) ∈ V h × M h , (29) a
(u -u h , u -u h ) ≤ a(u -u h , u -v h ) + b(λ -µ h , u h -u) + b(λ -λ h , u -v h ) + b(λ -µ h , u) + b(λ h , u) .
From [START_REF] Haslinger | Nečas-Numerical Methods for Unilateral Problems in Solid Mechanics[END_REF] and standard computations, we deduce Lemma 3.6. Let (u, λ) be the solution of problem [START_REF] Brezzi | Raviart-Error estimates for the finite element solution of variational inequalities, Part 2: Mixed methods[END_REF]. Assume that

u |Ωj ∈ H 2 (Ω j ), 1 ≤ j ≤ J. Let (u h , λ h ) be the solution of problem (17) with M h (γ i ) = M 1,ℓ, * h (γ i ) or M h (γ i ) = M ℓ,0 h (γ i ), 1 ≤ i ≤ I.
Then the following estimate holds

(30) u -u h α ≤ C    J j=1 α 1 2 j h j u H 2 (Ωj ) 2    λ -λ h -1 2 , * +    J j=1 α 1 2 j h 3 4 j u H 2 (Ωj ) 2    .
Remark 3.7. It is well known that a priori error estimates are not optimal if the number of points where the constraint changes from binding to non binding [START_REF] Haslinger | Nečas-Numerical Methods for Unilateral Problems in Solid Mechanics[END_REF][START_REF] Brezzi | Raviart-Error estimates for the finite element solution of variational inequalities[END_REF][START_REF] Brezzi | Raviart-Error estimates for the finite element solution of variational inequalities, Part 2: Mixed methods[END_REF]. In 2D (which is the case in this article), if

M h (γ i ) = M 1,ℓ, * h (γ i )
, and the number of points on each γ i where the constraint changes from binding to non binding, is finite, 1 ≤ i ≤ I, the following (optimal) estimate holds

(31) u -u h α ≤ C    J j=1 α 1 2 j h j u H 2 (Ωj ) 2    λ -λ h -1 2 , * +    J j=1 α 1 2 j h j u H 2 (Ωj ) 2    .
Remark 3.8. Note that the regularity assumption is local and thus generally is satisfied. The cases for which such an assumption is not valid arises only when some subdomains Ω j contain corners or if some changes in the boundary conditions (from Neumann to Dirichlet) occur. We exclude such situations in order to focus on the contact zone.

The proof of the following result is standard (see [START_REF] Coorevits | Sassi-Mixed finite element methods for unilateral problems: convergence analysis and numerical studies[END_REF] for example) Lemma 3.9. Under the same assumptions as in the previous lemma, the following estimate holds

(32) ( α m α M ) 1 2 λ -λ h -1 2 , * ≤ C( u -u h α +    J j=1 α j h 2 j u 2 H 2 (Ωj ) 2    1 2
), with a constant C independent of α.

Remark 3.10. The choice of the multiplier spaces yields the same result in each case. However, in the case of smooth solutions, namely λ ∈ H s (γ i ), 1 ≤ i ≤ I, s > 1 2 , the approximation order for piecewise constant multipliers does not change, while it increases for the other choices.

Finally, assembling the estimates of the two previous lemmas, we deduce Theorem 3.11. Let (u, λ) be the solution of problem [START_REF] Brezzi | Raviart-Error estimates for the finite element solution of variational inequalities, Part 2: Mixed methods[END_REF]. Assume that

u |Ωj ∈ H 2 (Ω j ), 1 ≤ j ≤ J. Let (u h , λ h ) be the solution of problem (17) with M h = M ℓ, * h or M h = M 0 h , the following estimate holds (33) u -u h α + ( α m α M ) 1 2 λ -λ h -1 2 , * ≤ C    J j=1 α 1 2 j h 3 4 j u H 2 (Ωj ) 2    . Remark 3.12. ii) In 2D (our case), if M h (γ i ) = M 1,ℓ, * h (γ i ),
and the number of points on each γ i where the constraint changes from binding to non binding, is finite,1 ≤ i ≤ I, the following estimate holds

(34) u -u h α + ( α m α M ) 1 2 λ -λ h -1 2 , * ≤ C    J j=1 α 1 2 j h j u H 2 (Ωj ) 2    .
Remark 3.13.

(1) In the case of homogeneous media, we retrieve the expected rate of convergence of O(h).

(2

) The case M h (γ i ) = M 1,ℓ h (γ i ) yields the rate of convergence O(α 1 2 h 1 2
) by adapting the argument of [START_REF] Ben Belgacem | Renard-Hybrid finite element methods for the Signorini problem[END_REF].

A posteriori analysis

The a posteriori analysis by residual error indicators that we develop in this section aims to present, in a relatively general and unified framework, the most used approximation approaches (conforming and nonconforming) for contact problems. Since these two approximation methods resort to different techniques to perform the analysis and to derive the a posteriori estimates we will presents the nonconforming case in details and we refere to the literature for the conforming one (e.g. [START_REF] Chen | Nochetto-Residual type a posteriori error estimates for elliptic obstacle problems[END_REF][START_REF] Hild | Nicaise-Residual a posteriori error estimators for contact problems in elasticity[END_REF]).

There are many difficulties to deal with in this analysis: discontinuous coefficients with high contrast, incompatible meshes at the interfaces of the subdomains, handling the contact conditions, and possible low regularity in these contact zones. In particular, a key point to conduct such an analysis is to build a Clément type operator: an interpolation operator for singular functions, that includes all these constraints and that satisfies optimal approximation properties. As usual, the analysis consists to define the well suited error indicators and to obtain upper and lower bounds of the error with them. In addition, we will give, as far as possible, an interpretation of the theoretical error estimates, in order to highlight how each error indicator acts in the adaptive process.

For the a posteriori error estimates we assume that f ∈ L 2 (Ω) 2 and we fix f h to be a finite element approximation of it associated with T h (usually this approximation is the same as for the displacement).

Notation . Given K ∈ T j h , we denote by E K the set of its edges not contained on the boundary ∂Ω j , 1 ≤ j ≤ J. The union of all E K , K ∈ T j h is denoted by E h,j and the union of E h,j , 1 ≤ j ≤ J will be denoted E h . We denote by E h,j,γi the set of edges of T j h which are contained in γ i and we set E h,γi = E h,jℓ(i),γi , 1 ≤ i ≤ I, ℓ = 1 or 2 with the same choice of ℓ as for M h . With each edge e ∈ E h or e ∈ E h,γi , 1 ≤ i ≤ I, we associate a unit vector n e normal to e and we denote by [ϕ] e the jump of the piecewise continuous (vector valued) function ϕ across e in the direction n e . For each K ∈ T h we denote by h K the diameter of K and we denote by h e the diameter of e, e ∈ E K . Remark 4.1. In all that follows when we write "contact zone" it means the zone where the contact is active, that is where the jump (u n) + -(u n) -= 0. In practice, we know dynamically that zone during the computations(see Section. (68)).

We define two kinds of residual error indicators

• An error indicator for the elements of the mesh. For each element K ∈ T h , we set

(35) η K = α -1 2 K h K f h + div ασ(u h ) L 2 (K) 2 + 1 2 e∈E K α -1 2 e h 1 2 e [ασ(u h ) n e ] e L 2 (e) 2 .
• Error indicators for the edges on the contact zone. For each e ∈ E h,γi , 1 ≤ i ≤ I, we set ( 36)

η e = h -1 2 e [u h n e ] e L 2 (e) 2 ,

and

(37)

η C,e = α -1 2 e h 1 2
e β e λ h + n e • ασ(u h ) n e L 2 (e) ,

with

β e = 1 if M h (γ i ) = M 0 h (γ i ) or M 1,ℓ h (γ i ) and β e = 0 if M h = M 1,ℓ, * h (γ i ).
We denote indifferently the scalar product x • ασ(u h ) y or ασ(u h ) x y for any vectors x, y.

Remark 4.2. The residual error indicators appear from the computation of the residu of the equation in the appropriate norm (the dual of the energy norm see [START_REF]Verfürth-A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques[END_REF]). The first indicator η K , is the standard one for a nonhomgeneous material in linear elasticity [START_REF]Verfürth-A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques[END_REF]. The second one, η e is due to the nonconformity of the meshes at the contact zone. Basically, it is the jump in the H 1 2 -norm of the normal component of the displacement. Note that it is zero, except at locations where the contact is active and where the meshes or not compatible. The last error indicator is specific to the contact condition it could be expressed in η K -in this case it looks like a Neumann condition on each γ i -but for the clarity and ease of implementation we define it separately.

Remark 4.3. The error indicator η C,e measures the norm H -1 2 of the multiplier (the pressure on the contact zone). Its interpretation is the following : in the conforming case, it expresses the residual λ h = n e .(ασ(u h ) n e ) . In the nonconforming case, it expresses the residual n e .(ασ(u h ) n e ) (which should vanishes for the exact solution if the contact is active (the jump is zero) and if the contact is not active, the non penetration condition imply that the pressure (on each side) should be zero (for the exact solution) and so is the "jump".

4.1.

Abstract upper bound for the error : the nonconforming case. We set M h (γ i ) = M 1,ℓ, * h (γ i ) and we will derive an abstract bound of the error by the error indicators.

The ellipticity of a α (., .) on V gives [START_REF] Brezzi | Raviart-Error estimates for the finite element solution of variational inequalities, Part 2: Mixed methods[END_REF] and subtracting the first equation of [START_REF] Chen | Nochetto-Residual type a posteriori error estimates for elliptic obstacle problems[END_REF], we obtain [START_REF] Kikuchi | Oden--Contact problems in elasticity: A study of variational inequalities and finite element methods[END_REF] a α (u

(38) u -u h 2 α ≤ C a α (u -u h , u -u h ). Since V h ⊂ V, choosing v = v h ∈ V h in the first equation of
-u h , v h ) + b(λ -λ h , v h ) = 0, ∀v h ∈ V h .
We set w = uu h , and we fix an approximation w h of w in V h , then we deduce from ( 38) and ( 39) that

u -u h 2 α ≤ C (a α (u -u h , w -w h ) -b(λ -λ h , w h )).
Integrating by parts, inserting f h and considering separately the cases where e ∈ E K belongs to ∪ I i=1 γ i or not, yields (40)

u-u h 2 α ≤ C K∈T h ( K (f h +div α σ(u h ))•(w-w h ) dx+ K (f -f h )•(w-w h ) dx - 1 2 e∈E K e [n e • ασ(u h )] (w -w h ) dτ ) - I i=1 e∈E h,γ i e n e • ασ(u h ) [w -w h ] dτ -b(λ, w -w h ) -b(λ -λ h , w h ) .
Choosing w h as a conforming approximation of w (The construction of such w h is not obvious, it is given in proposition 4.7) yields

(41) u -u h 2 α ≤ C K∈T h ( K (f h + div α σ(u h )) • (w -w h ) dx + K (f -f h ) • (w -w h ) dx - 1 2 e∈E K e [n e • ασ(u h )] (w -w h ) dτ ) - I i=1 e∈E h,γ i e n e • ασ(u h ) n e [w • n e ] dτ -b(λ, w) . Note that b(λ, w) = b(λ, u h ) = b(λ -λ h , u h ).
Moreover, for each e ∈ E h,γi , 1

≤ i ≤ I, e (λ -λ h ) [u h • n] dτ ≤ λ -λ h H -1 2 (e) [u h • n] H 1 2 (e) ≤ C λ -λ h H -1 2 (e) (h -1 2 e [u h • n] L 2 (e) ),
summing up for all e ∈ E h,γi and taking the square root, we obtain (42) (

I i=1 e∈E h,γ i e (λ -λ h ) [u h • n] dτ ) 1 2 ≤ C λ -λ h -1 2 , * ( I i=1 e∈E h,γ i h -1 e [u h • n] 2 L 2 (e) ) 1 2
.

Next, invoking the inf-sup condition and (39), we have

( α m α M ) 1 2 λ h -µ h -1 2 , * ≤ sup w h ∈V h b(λ h -µ h , w h ) w h α ≤ sup w h ∈V h b(λ h -λ, w h ) w h α + sup w h ∈V h b(λ -µ h , w h ) w h α ≤ sup w h ∈V h a α (u -u h , w h ) w h α + sup w h ∈V h b(λ -µ h , w h ) w h α ≤ C u -u h α + C ′ inf µ h ∈M h λ -µ h -1 2 , * .
Using the triangle inequality

λ -λ h -1 2 , * ≤ λ h -µ h -1 2 , * + λ -µ h -1
and combining it with [START_REF] Kinderlehrer | Stamppachia-An introduction to variational inequalities and their applications[END_REF], we obtain ( 43)

u-u h α ( u-u h α +( α m α M ) 1 2 λ-λ h -1 2 , * ) ≤ C K∈T h ( K (f h +div α σ(u h ))•(w-w h ) dx + K (f -f h ) • (w -w h ) dx - 1 2 e∈E K e [n e • ασ(u h )] (w -w h ) dτ ) - I i=1 e∈E h,γ i e n e • ασ(u h ) [w] dτ -b(λ -λ h , u h ) .
Using [START_REF] Lions | Magenes-Problèmes aux limites non homogènes[END_REF], Cauchy-Schwartz inequality, and dividing by w α , except the last term that we divide by (

α m α M ) 1 2 λ -λ h -1 2 , * , this gives (44) u -u h α ≤ C K∈T h ( f h + div α σ(u h )) L 2 (K) 2 w -w h L 2 (K) w α + f -f h L 2 (K) 2 w -w h L 2 (K) 2 w α + 1 2 e∈E K [ασ(u h )n e ] L 2 (e) 2 w -w h L 2 (e) 2 w α ) - I i=1 e∈E h,γ i e ασ(u h )n e [w] dτ w α + C ′ ( α m α M ) - 1 2 ( I i=1 e∈E h,γ i h -1 e [u h • n] 2 L 2 (e) ) 1 2
.

Finally, we obtain the abstract upper bound in the nonconforming case Lemma 4.4. The upper bound on the error in the nonconforming case is given by

(45) u-u h α +( α m α M ) 1 2 λ-λ h -1 2 , * ≤ C K∈T h ( f h +div α σ(u h )) L 2 (K) 2 w -w h L 2 (K) w α + f -f h L 2 (K) 2 w -w h L 2 (K) 2 w α + 1 2 e∈E K [ασ(u h )n e ] L 2 (e) 2 w -w h L 2 (e) 2 w α ) - I i=1 e∈E h,γ i e ασ(u h )n e [w] dτ w α + C ′ ( α m α M ) - 1 2 ( I i=1 e∈E h,γ i h -1 e [u h • n] 2 L 2 (e) ) 1 2 + inf µ h ∈M h λ -µ h -1 2 , * .
For the convenience of the reader, we only give the upper bound in the conforming case Lemma 4.5. The upper bound on the error in the conforming case is given by

(46) u -u h α + λ -λ h -1 2 , * ≤ C K∈T h f h + div α σ(u h ) L 2 (K) 2 w -w h L 2 (K) 2 w α + f -f h L 2 (K) 2 w -w h L 2 (K) 2 w α + 1 2 e∈E K [α σ(u h )n e ] L 2 (e) 2 w -w h L 2 (e) 2 w α - e∈E hΓ C ( e (λ h + n e • α σ(u h ) n e ) [w • n e ] dτ w α + 1 µ -1 2 , * b(µ -µ h , u h ).
Remark 4.6. Comparing these two upper bounds [START_REF] Tahir | Problèmes de contact unilatéral et maillages incompatibles[END_REF] and (45), we note that the fourth term in (45) (the nonconforming case) do not appears in the second inequality. In fact, it is the "variational crime", similar to the consistency error in the a priori analysis. The two last terms in both inequalities are different since they involve the contact conditions and thus are inherent to the discretization in this region and particularly to the choice of the Lagrange multipliers discrete cone.

4.2.

Final upper bound for the error. In this section, we will derive the upperbound of the discretization error by the error indicators. The common factors in the error estimates ( 45) and ( 46) will be handled by introducing a well suited Clément type operator, which requires some technical assumptions that we introduce below and which imposes some constraints on the topology of the meshes. The specific factors will be estimated directly in Lemma 4.8 and Lemma 4.9.

Notations and definitions. For any K ∈ T h , we denote by ∆ K , resp. ∆ e , the union of all elements that share at least one vertex, resp. one edge, with K. We denote by N h , N K , and N e , the set of all vertices of elements of T j h , 1 ≤ j ≤ J, of a given element K, and of a given edge e, respectively. With each vertex z we associate the corresponding unique continuous, piecewise affine function that takes the value 1 at z and vanishes at all other vertices. We denote by E - h,γi the set E h,γi (recall that it coincides with E h,ℓ,γi when ℓ is such that

M h (γ i ) = M 1,ℓ h (γ i ) or M h (γ i ) = M 1,ℓ, * h (γ i )) and E + h,γi = E h,3-ℓ,γi .
We denote also by N - C and N + C the set of the vertices which are endpoints of e ∈ E - h,γi and E + h,γi respectively. We need the following assumption (see Figure 2), that allows to handle the adaptivity with nonconforming meshes Assumption 1 . Each element of E - h,γi is the union of some entire elements of E + h,γi and the number of such elements is independent of h.

This assumption is easy to satisfy in the practical implementation of the adaptive strategy.

Assumption 2 . For any two subdomains Ω i and Ω j sharing at least one point, there is a connected path passing from Ω i to Ω j trough adjacent subdomains such that α is monotone along this path (adjacent means they share a common edge).

The following proposition is proved in the Appendix Proposition 4.7. Under assumption 2, there exists an operator R h from V into V h ∩ H 1 0 (Ω) 2 , and a constant C > 0 depending only on the shape parameter of T h , such that: for any v ∈ V, every element K and every edge e of K, the following estimates hold .

(47) v -R h v L 2 (K) 2 ≤ Ch K α -1 2 K v 1,α,∆ K , v -R h v L 2 (e) 2 ≤ ch
Using the inverse inequality [START_REF]Ciarlet-Basic Error Estimates for Elliptic Problems[END_REF], we get

e (β e λ h + n e • ασ(u h ) n e ) [w • n e ] dτ ≤ C h 1 2 e (β e λ h + n e • ασ(u h ) n e L 2 (e) K∈Ke [w • n e ] H 1 2 (∂K) ≤ C h 1 2 e α -1 2 e (β e λ h + n e • ασ(u h ) n e L 2 (e) K∈Ke u -u h 1,α,K .
Lemma 4.9. We have

(49) b(µ -µ h , u h ) µ -1 2 , * ≤ C( I i=1 e∈E h,γ i η 2 e ) 1 2
.

Proof. For each e ∈ E h,γi , 1 ≤ i ≤ I, we have

e (µ -µ h ) [u h • n e ] dτ ≤ µ -µ h H -1 2 (e) [u h • n e ] H 1 2 (e)
, and using the inverse inequality [START_REF]Ciarlet-Basic Error Estimates for Elliptic Problems[END_REF], we note that

[u h • n e ] H 1 2 (e) ≤ h -1 2 e [u h • n e ] L 2 (e) .
Extending the operator π 1 h , respectively π 0 h , to H -1 2 (γ i ), we have

µ -µ h H -1 2 (e) ≤ c µ -µ h (H -1 2 (γi)) ≤ c µ (H -1 2 (γi))
Summing up for all e ∈ E h,γi , 1 ≤ i ≤ I yields the result.

Finally, choosing 46) and ( 45), using Lemma 4.7, Lemma 4.8 and Lemma 4.9 yield the upper bound on the error. Theorem 4.10. If Assumption 1 is satisfied, there exists a constant C independent of h and α, such that if

w h = R h (w) ∈ V ∩ H 1 (Ω) 2 in (
M h (γ i ) = M 1, * ,ℓ h (γ i ), (50) u -u h α + ( α m α M ) 1 2 λ -λ h -1 2 , * ≤ C K∈T h (η 2 K + α -1 K h 2 K f -f h 2 L 2 (Ω C ) 2 ) + I i=1 e∈E h,γ i (η 2 e + η 2 C,e ) 1 2
.

Corollary 4.11. If Assumption 1 is satisfied, there exists a constant C independent of h and α, such that if

M h (γ i ) = M 0,ℓ h (γ i ) or = M 1,ℓ h (γ i ) (51) u -u h α + λ -λ h -1 2 , * ≤ C K∈T h (η 2 K + α -1 K h 2 K f -f h 2 L 2 (Ω C ) 2 ) + I i=1 e∈E h,γ i (η 2 e + η 2 C,e )
+ I i=1 e∈E h,γ i e e λ [w • n e ] dτ = (f , w) -a α (u h , w) -b(λ h , w) = K∈T h K (f h + div α σ(u h )) • w dx + K (f -f h ) • w dx - e∈E K n e • (ασ(u h )) w dx -b(λ h , w) = K∈T h K (f h + div α σ(u h )) • w dx + K (f -f h ) • w dx - 1 2 e∈E h e [n e • (ασ(u h )] w dτ - I i=1 e∈E h,γ i e (ασ n (u h ) + λ h ) [w • n e ] dτ.
We denote by r h w = (r i h w) 1≤i≤I . If M h = M 1,ℓ, * h (γ i ), following the same lines than (52), and taking into account b(λ h , r h w) ≤ 0 (by the definition of M h ), we have (53)

a α (e u , w) + b(e λ , w) = K∈T h K (f h + div α σ(u h )) • w dx + K (f -f h ) • w dx - e∈E K n e • (α σ(u h )) w dx -b(λ h , w) ≥ K∈T h K (f h + div α σ(u h )) • w dx + K (f -f h ) • w dx - e∈E K n e • (ασ(u h )) w dx -b(λ h , w -r h w) = K∈T h K (f h + div α σ(u h )) • w dx + K (f -f h ) • w dx - 1 2 e∈E h e
[n e • (ασ(u h )] w dτ -

I i=1 e∈E h,γ i e α σ n (u h ) [w • n e ] dτ - I i=1 e∈E h,γ i e λ h ( w • n e -r i h (w • n e ) dτ .
We will obtain the desired estimates by appropriate choices of w. With each K ∈ T h , we associate the bubble function ψ K equal to the product of the three barycentric coordinates on K. For each e ∈ E Ω , we associate the bubble function ψ e equal to the product of the two barycentric coordinates on e. We introduce a lifting operator defined as follows: On the reference element K, we fix a lifting operator P from polynomial traces on ê on K that vanish at the endpoints of ê, into polynomials on K that vanish on ∂ K \ ê. A similar operator is obtained on each K by an affine transformation.

Proposition 4.12. There exists a constant C independent of h and α such that:

∀K ∈ T h (54) η K ≤ C u -u h 1,α,G K + ( K ′ ∈G K α -1 K ′ h 2 K ′ f -f h 2 L 2 (K ′ ) 2 ) 1 2 ,
where G K is the union of K and all triangles containing an edge of K.

Proof. We will bound each term of η K which, as standard in a posteriori analysis ( [START_REF]Verfürth-A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques[END_REF]), relies on appropriate choices of the test function w in (52) [START_REF] Ainsworth | Lee-Local a posteriori error estimators for variational inequalities[END_REF] We take w in (52) equal to

w = (f h + div α σ(u h ))ψ K in K, 0 elsewhere.
Since ψ K vanishes on ∂K, this yields

(f h + div α σ(u h ))ψ 1 2 K 2 L 2 (K) 2 ≤ K α σ(e u ) : ε((f h + div α σ(u h ))ψ K ) dx - K (f -f h ) • (f h + div α σ(u h ))ψ K ) dx.
It follows that

(f h + div α σ(u h ))ψ 1 2 K 2 L 2 (K) 2 ≤ C( u -u h 1,α,K |(f h + div α σ(u h ))ψ 1 2 K | 1,α,K + f -f h L 2 (K) 2 (f h + div α σ(u h ))ψ 1 2 K L 2 (K) 2 .
It can be checked by going to the reference element that for any polynomial ϕ of degree at most k, the following inequalities hold [START_REF]Ciarlet-Basic Error Estimates for Elliptic Problems[END_REF] :

ϕ L 2 (K) ≤ C ϕψ 1 2 K L 2 (K) , |ϕψ K | H 1 (K) ≤ Ch -1 K ϕ L 2 (K)
, with constants depending only on k and the shape parameter of K. Noting that ψ K is ≤ 1 (and with obvious extension to vector valued functions), we obtain

(55) α -1 2 K h K f h + div α σ(u h ) L 2 (K) 2 ≤ c( u -u h 1,α,K + α -1 2 K h K f -f h L 2 (K) 2 ).
(2) We denote by e an edge in E K . We distinguish between two cases.

• First, if e is not contained in ∪ I i= γ i , it is a common edge of the two adjacent elements K and K ′ . Now, we choose w in (52) to be equal to

w =    P K,e ([α σ(u h )n e ] ψ e ) in K, P K ′ ,e ([α σ(u h )n e •] ψ e ) in K ′ , 0 elsewhere, This yields
[ασ(u h )n e ] ψ 1 2 e 2 L 2 (e) 2 ≤ κ∈(K,K ′ ) |u -u h | 1,α,κ |P κ,e ([ασ(u h )n e ] ψ e )| 1,α,κ + ( f h + div α σ(u h ) L 2 (κ) 2 + f -f h L 2 (κ) 2 ) P κ,e ([ασ(u h )n e ] ψ e ) L 2 (κ) 2 .
The following inequalities, obtained by going to the reference element, hold

[α σ(u h ) n e ] L 2 (e) 2 ≤ c [α σ(u h )n e ] ψ 1 2 e L 2 (e) 2 , |P κ,e ([α σ(u h )n e ] ψ e )| H 1 (κ) 2 +h -1 e P κ,e ([α σ(u h )n e ] ψ e ) L 2 (κ) 2 ≤ C h -1 2 e [α σ(u h ) n e ] L 2 (e) 2 .
Noting that ch κ ≤ h e ≤ h κ , we obtain

(56) α -1 2 h 1 2 e [α σ(u h ) n e ] L 2 (e) 2 ≤ C κ∈(K,K ′ ) (|u -u h | 1,α,κ + α -1 2 K h κ f -f h L 2 (κ) 2 + α -1 2 K h κ f h + div α σ(u h ) L 2 (κ) 2 ). • Let e ∈ ∪ I i=1 γ i . If e ∈ E +
γi , for a given i, we denote by K ′ the element such that e is contained in e ′ ∈ E K ′ . We extend [ασ(u h )n e ] ψ e to the entire e ′ by zero and we make the same choice of w as previously and we obtain (56). If e ∈ E - γi , then according to assumption 1, we replace K ′ by the finite number of elements K i , which share the edge e, we define w as previously with respect to this change and we proceed similarly to the previous case to obtain the estimate (56) Corollary 4.13. If Assumption 1 is satisfied, there exists a constant C independent of h and α such that

(57) ( K∈T h η 2 K ) 1 2 ≤ C ( u -u h α + ( K∈T h α -1 K h 2 K f -f h 2 L 2 (K) 2 ) 1
2 ), Proposition 4.14. If Assumption 1 is satisfied, we take

β e = 0 if M h (γ i ) = M 1,ℓ, * h (γ i )
and β e = 1 otherwise, then, there exists a constant C independent of h and α such that

(58) η C,e ≤ C u -u h 1,α,Ge + α -1 2 e h 1 2 e λ -λ h L 2 (e) + (1 -β e ) α -1 2 e h 5 4 e λ h L 2 (e) + ( κ∈Ge α -1 κ h 2 κ f -f h 2 L 2 (κ) 2 ) 1 2 ,
where G e is the union of all triangles having a non null intersection with e.

Remark 4.15. Note that either β e = 1 and the third term in the right-hand side vanishes, or β e = 0 and it is of high order. In all cases it could be neglected.

Proof. Let us denote by σ n (u) = (n • σ(u))n. For e ∈ E h,γi , 1 ≤ i ≤ I, we take w in (52) equal to

w = P K,e (ψ e (β e λ h n e + α σ(u h )n e )) in K, 0 elsewhere.
with

β e = 1 if M h (γ i ) = M 0 h (γ i ) or M 1,ℓ h (γ i ) and β e = 0 M h (γ i ) = M 1,ℓ, * h (γ i ). • If M h (γ i ) = M 0 h (γ i ) or M 1,ℓ h (γ i
) then arguing as in the previous proposition, we deduce that

(λ h + ασ n (u h ) ψ 1 2 e 2 L 2 (e) 2 ≤ u -u h 1,α,K |P κ,e ([λ h n e + ασ(u h )n e ] ψ e )| 1,α,K + ( f h + div α σ(u h ) L 2 (K) 2 + f -f h L 2 (K) 2 ) P κ,e ([λ h n e + ασ(u h )n e ] ψ e ) L 2 (K) 2 + e (λ -λ h ) (λ h + ασ n (u h )) ψ e dτ ).
Thus, we obtain (59) α 52), we deduce as in the previous case

-1 2 e h 1 2 e λ h + ασ n (u h ) L 2 (e) ≤ C( u -u h | 1,α,K + α -1 2 K h K f -f h L 2 (K) 2 + α -1 2 K h K f h + div ασ(u h ) L 2 (K) 2 ) + α -1 2 e h 1 2 e λ -λ h L 2 (e) . • If M h (γ i ) = M 1,ℓ, * h (γ i ): using (
ασ n (u h ) ψ 1 2 e 2 L 2 (e) 2 ≤ u -u h 1,α,K |P κ,e ([ασ n (u h )] ψ e )| 1,α,K + ( f h + div ασ(u h ) L 2 (K) 2 + f -f h L 2 (K) 2 ) P κ,e ([ασ n (u h )] ψ e ) L 2 (K) 2 + e (λ -λ h ) (ασ n (u h )) ψ e dτ ) + e λ h (ασ n (u h ) -r i h (ασ n (u h ))) ψ e dτ ).
The last term and ( 19) yield the additional extra term α

-1 2 e h 5 4
e λ h L 2 (e) on the right-hand side of (59). Corollary 4.16. If Assumption 1 is satisfied, there exists a constant C independent of h and α such that, (60) (

I i=1 e∈E h,γ i η 2 C,e ) 1 2 ≤ C ( u -u h α + ( I i=1 e∈E h,γ i α -1 e h e λ -λ h 2 L 2 (e) ) 1 2
+ (

K∈T h α -1 K h 2 K f -f h 2 L 2 (K) 2 ) 1 2
).

Remark 4.17. The second term in the right-hand side is bounded by a constant times the norm H -1 2 (Γ C ). However, the constant depends on the ratio ( αm α M ). Considering the a posteriori estimate (60) (with regard to the error for the Lagrange multipliers), shows that when the ratio between adjacent bodies (sharing the same contact zone) is large, nonconforming approximation is more suited.

It now remains to bound η e . Proposition 4.18. If Assumption 1 is satisfied, there exists a constant C independent of h and α such that the following estimate holds

(61) η e ≤ C γ e u -u h H 1 (Ge) 2
where G e denotes, as before, the union of all elements having a non null intersection with e and γ e = 1 + max(α -1 K , α -1 K ′ ). Proof. We define a function µ h on ∪ I i=1 E h,γi by

µ h = [u h • n e ]
ψ on e, 0, elsewhere.

(62) b(µ h , u h ) = b(µ h , u hu)

• If E h,γi = E + h,γi then we consider the two elements K, K ′ from both sides of γ i such that e is an entire edge of K and is contained on an entire edge e ′ of K ′ . We extend the product µ h ψ e by zero to e ′ and then we solve the problem: for κ ∈ (K, K ′ ), find ϕ ∈ H 1 0,∂κ\e (κ) (63)

   ∆ϕ = 0, in κ, ϕ = µ h , on e, ϕ = 0, elsewhere,
where

H 1 0,∂κ\e (κ) = v ∈ H 1 (κ); v = 0 on ∂κ \ e .
From the definition of µ h , b(., .), (62) and integration by parts, it follows that

[u h • n e ] ψ 1 2 e 2 L 2 (e) = b(µ h , u h ) = b(µ h , u h -u) = e µ h [(u h -u) • n e ] dτ = κ∈(K,K ′ ) ( κ grad ϕ • (u h -u) dx - κ ϕ div (u h -u) dx)
and we deduce that

[u h • n e ] 2 L 2 (e) ≤ [u h • n e ] ψ 1 2 e 2 L 2 (e) ≤ C κ∈(K,K ′ ) (|ϕ| H 1 (κ) u h -u L 2 (κ) 2 + h -1 κ ϕ L 2 (κ) h κ div (u h -u) L 2 (κ) |.
By going to the reference element and using direct estimates on the boundary value problem (63), we have

h -1 κ ϕ L 2 (κ) + |ϕ| H 1 (κ) ≤ C [u h • n e ] ψ 1 2 e L 2 (e) ≤ C h -1 2 e [u h • n e ] L 2 (e) .
Thus, we have

[u h • n e ] 2 L 2 (e) ≤ C h -1 2 e [u h • n e ] L 2 (e) κ∈(K,K ′ ) ( u h -u L 2 (κ) 2 + h κ div (u h -u) L 2 (κ) ) + [u h • n e ] L 2 (e)
Noting that ch κ ≤ h e ≤ c ′ h κ we have

[u h • n e ] L 2 (e) ≤ Ch -1 2 e κ∈(K,K ′ ) h κ u h -u H 1 (κ) 2 Ch 1 2 e κ∈(K,K ′ ) u h -u H 1 (κ) 2 , thus, (64) h
-1 2 e [u h • n e ] L 2 (e) ≤ C δ e κ∈(K,K ′ ) u h -u 1,α,κ ,
where δ e = 1 + max(α -1 K , α -1 K ′ ). which is the disred estimate.

• If E h,γi = E - h,γi
, then it follows from Assumption 1, that e is an entire edge of an element K on one side and is a union of edges e ′ i , i = 1, . . . , i * . Since u h belongs to V h it is continuous at the endpoints of e ′ i , therefore we can still solve (63) with K ′ replaced by ∆ e ∩ Ω + C . Arguing as in the previous case we obtain the desired estimate.

I i=1 e∈E h,γ i η 2 e ) 1 2 ≤ C u -u h α .

Numerical experiments

First, we describe briefly the implementation of discrete problem [START_REF] Chen | Nochetto-Residual type a posteriori error estimates for elliptic obstacle problems[END_REF]. All of the numerical experiments are achieved with FreeFem++ [START_REF] Hecht | Pironneau-FreeFem++[END_REF]. The discrete solution (u h , λ h ) is a saddle-point of the Lagrangian functional defined over K h by

L(v h , µ h ) = 1 2 a α (v h , v h ) -L(v h ) -b(µ h , v h ),
which means that it satisfies the min-max principle (66)

L(u h , µ h ) ≤ L(u h , λ h ) ≤ L(v h , λ h ), ∀(v h , µ h ) ∈ K h .
Let V, U α denote the vectors with the entries given by the nodal values of the functions (v h , µ h ) and (u h , λ h ), respectively. Let M and Λ be the vectors with the entries given by the nodal values of µ h and λ h , respectively, for the three different choices of the space M h , namely

M h (γ i ) = M 1 h (γ i ), M h (γ i ) = M 1, * h (γ i ) or M h (γ i ) = M 0 h (γ i ), 1 ≤ i ≤ I.
The saddle-point problem for the Lagrangian (66) can be rewritten in the finite dimensional setting : Find U α = (u h , λ h ) and Λ, solution of the following max-min problem (67) max

P M ≥0 min V 1 2 t VKV -t VF + ( t V L)P M ,
where the matrix K denotes the stiffness matrix, L the coupling matrix and P expresses the sign conditions for the multipliers. F denotes the vector corresponding to the external loads.

Given the triangulations T j h of Ω j , 1 ≤ j ≤ J, let N j , 1 ≤ j ≤ J denote the number of nodes in Ω j . We introduce the finite element basis of V h (Ω j ):

(η 1 , . . . , η 2Nj ) = ((w 1 e 1 , w 1 e 2 ), . . . , (w Nj e 1 , w Nj e 2 )), j = 1, . . . , J, where (w i ) i denotes the (scalar) affine Lagrange finite element basis and (e 1 , e 2 ) the canonical basis in R 2 . Then the matrix K is defined by

K =     K 1 0 . . . 0 K 2 0 . . . . . . . . . K J     , (K j ) ns = a α (η n , η s ) = Ωj ε(η n ) : A α ε(η s ) dx, n, s = 1, . . . , 2N j ,
and the right-hand side takes the form F = (F j ), with F j = (f j n ) n , D = ( Ωj η n η s dx) ns , n, s = 1, . . . , 2N j .

Since each contact zone γ i occurs between two bodies Ω j1(i) and Ω j2(i) , we will drop the index i, and we denote by Ω 1 and Ω 2 , the subdomains in contact at γ i . Let m ℓ denote the number of nodes of Ω ℓ on γ. We fix ℓ = 1 as a Lagrange multiplier side and we define (ψ k ) k , 1 ≤ k ≤ m 1 -1 to be the finite element basis associated with W

1,1 h (γ) and (ϕ k ) k , 1 ≤ k ≤ m 1 -1, to be the finite element basis associated with W 0,1 h (γ). If M h (γ) is = M 0,1 h (γ) or = M 1,1 h (γ)
, then P is given by the identity matrix, else M h (γ) = M 1,1, * h (γ), and

P ij = γ ψ i ψ j dτ , 1 ≤ i, j ≤ m 1 -1.
Finally, the coupling matrix L = L 1 L 2 is defined in the following way

• If M h (γ) = M 1,1 h (γ) or M h (γ) = M 1,1, * h (γ), then (L ℓ ) ij = γ ψ j (η i • n ℓ ) dτ, 1 ≤ i ≤ N ℓ , 1 ≤ j ≤ m 1 -1. • If M h (γ) = M 0,1 h (γ), then (L ℓ ) ij = γ ϕ j (η i • n ℓ ) dτ, 1 ≤ i ≤ N ℓ , 1 ≤ j ≤ m 1 -1.
In order to solve the discrete problem we will use an adapted semi-smooth Newton method [START_REF] Hintermüller | Kunisch-The dual-primal active set strategy as a semi-smooth Newton method[END_REF][START_REF] Ito | Semi-smooth Newton methods for variatonal inequalities of the first kind[END_REF] which consists in finding the contact zone (68) S = λ h ∈ W k,1 h (γ), k = 0, 1 P λ h < 0 . We use the algorithm of primal-dual active set strategy [ [START_REF] Ito | Semi-smooth Newton methods for variatonal inequalities of the first kind[END_REF], algorithm (2.9)]. The stopping criterion consists in checking if the residual vanishes at the endpoints of each connected component of a stable set S. In practice the algorithm is very fast and very efficient.

Remark 5.1. One has to pay attention to the initial guess in the algorithm which has to contain the full reached contact zone. In practice, such a choice is always possible because the candidate area to the contact is always known. If the initial guess does not fulfill this condition, it may occur that the missed part of the contact zone is not obtained when convergence is reached.

5.1. Adaptive strategy. We start always with a fixed uniform or quasi-uniform triangulation T h,n , n = 1. Next, we perform iteratively the following adaptivity step:

On the triangulation T h,n , we compute the solution (u h , λ h ) of problem [START_REF] Chen | Nochetto-Residual type a posteriori error estimates for elliptic obstacle problems[END_REF], the corresponding error indicators as defined in ( 35), ( 36) and ( 37) and the mean value (69)

η n h = 1 N n h K∈T h,n η K ,
where N n h is the number of triangles in T h,n . Next, each triangle K such that η K ≥ η n h , is divided in such a way so that the diameters of the new triangles inside it are very close to h K times the ratio η n h η K .

When η e + η C is of the same order as η K we perform a second step of the adaptivity on the edges of Γ C following the same strategy as above.

The adaptivity is performed either a fixed (small) number of times or until the quantity η n h becomes smaller than a given tolerance. 5.2. Example 1. In this case the data are as follows: we consider the L-shaped domain Ω

1 = ]-10, 0[ × ]0, 1[ ∪ ]-1, 0[ × ]-3, 1[, Ω 2 = ]-15, -5[ × ]-2, 0[, Γ c = ]-10, -5[ × {0}, ν = 0.3, E = 215. The surface force is f (x) = -1 and is defined on the boundary Γ N = ]-1, 0[ × {-3}. Also, α(x, y) = 1 in Ω 1 , 10 -3 in Ω 2 ,
We show the initial and the final deformed mesh after 10 cycles of mesh refinement, as well as the error indicators on Figures 34and on Figures 5-6

Example 2. The T-shaped domain is a union of two rectangles Ω

1 = ]-5, 5[× ]0, 1[ ∪ ]-1, 1[ × ]-4, 1[ and the support part is Ω 2 = ]-8, -1.5[ × ]-2, 0[ and Ω 3 = ]1.5, 8[ × ]-2, 0[, Γ c = ]-5, -1.5[ × {0} ∪ ]1.5, 5[ × {0}, ν = 0.3, E = 2150.
The gravity force is g = -0.1. To remove free sliding in the x-direction we impose zero displacement on {0} × ]-5, 1[. We give the initial deformed mesh in Figure 7 and the final one in Figure 8 for α = 1 in the entire domain.

We give the same plots in Figures 9-10 for the case The convergence of the adaptive strategy is given on Figure 11. We have computed a reference solution on a very fine uniform mesh with a constant mesh size of 0.05. The mesh size in the coarse initial mesh is 0.3 in all examples. We plot, as a function of iteration numbers, the global error indicator η, the error in the energy norm and the error indicator η-mortar, the Hilbertian sum of the error indicators on the contact zone. We observe in this example that both η and η-mortar decrease with the energy norm of the error. We also remark that the η-mortar is negligible for the adaptive process in this example and this is not surprising since solutions of the Signorini unilateral contact problems with straight contact zones are smooth, see [START_REF] Moussaoui | Régularité des solutions d'un problème mêlé Dirichlet-Signorini dans un domaine polygonal plan[END_REF] for more details. in Ω 1 , 10 -1 in Ω 2 , 10 -2 in Ω 3 .

α(x, y) = 1 in Ω 1 , 10 -2 in Ω 2 ∪ Ω 3 .
We plot the final deformed mesh in Figure 12 where we can see the influence of the parameter α in the adaptive strategy Final mesh after 10 cycle of adaptivity α 2 = 10 -1 , α 3 = 10 -2

• Consider now z which is not in ∂Ω but is in ∂Ω ℓ(K) where ℓ(K) is such that K ∈ Ω ℓ(K) . If ℓ(K) = ℓ(z), then the previous argument with ω z replaced by ω z ∩ Ω ℓ(K) still applies. • If ℓ(K) = ℓ(z), we have

ϕ z (v i -R h v i ) L 2 (K) ≤ ϕ z (v i - ωz∩Ω ℓ(z) v i dx) L 2 (K) ≤ ϕ z (v i - ωz∩Ω ℓ(K) v i dx) L 2 (K) + ϕ z ( ωz∩Ω ℓ(K) v i - ωz∩Ω ℓ(z) v i dx) L 2 (K) .
The first term is estimated as previously. The second term is estimated as follows

ϕ z ( ωz∩Ω ℓ(K) v i - ωz∩Ω ℓ(z) v i dx) L 2 (K) ≤ ϕ z L 2 (K) | ωz∩Ω ℓ(K) v i - ωz∩Ω ℓ(z) v i dx)| ≤ ch K | ωz∩Ω ℓ(K) v i - ωz∩Ω ℓ(z) v i dx)|.
Consider first the case where the two subdomains Ω ℓ(z) and Ω ℓ(K) are adjacent, i.e. they share a common edge that we denote by e • As a subcase, we still have ℓ(K) = ℓ(z) and we consider z such that z ∈ γ, γ ∈ ∪ I i=1 γ i being a crack between the two subdomains Ω ℓ(z) and Ω ℓ(K) . Then it is an endpoint of e which is the entire edge of two elements, each in one subdomain. Using the regularity of the triangulations T ℓ h , ℓ = ℓ(K), ℓ(z), we have e v i -

h K | ωz∩Ω ℓ(K) v i - ωz∩Ω ℓ(z) v i dx)| ≤ c h
ωz∩Ω ℓ(κ) v i L 2 (e) ≤ c v i - ωz∩Ω ℓ(κ) v i dx L 2 (κ) + h e |v i | H 1 (κ) ≤ c h K α -1 2 K v i 1,α,∆ K .
• If z ∈ γ and z ∈ N - C \ N + C , it follows from the Assumption 1, that z is an endpoint of e which is an entire edge of K and only a part of ẽ an edge of K. In addition v = (v i ), i = 1, 2 is not continuous through γ, therefore, we have to modify the previous argument. Inserting and subtracting h -1 e e v i|Ω ℓ(κ) , κ = (K, K), we obtain

h K | ωz∩Ω ℓ(K) v i - ωz∩Ω ℓ(z) v i dx| ≤ h K | ωz∩Ω ℓ(K) v i -h -1
v i dx|) dτ + h -1 e | e [v i ] dτ | , ≤ c h K h -1 2 e v i|Ω ℓ(K) - ωz∩Ω ℓ(K) v i dx L 2 (e) + h -1 2 e v i|Ω ℓ(z) - ωz∩Ω ℓ(z) v i dx L 2 (e) + h -1 2 e [v i ] L 2 (e) .
Invoking the trace theorem (73) once, we obtain

h K | ωz∩Ω ℓ(K) v i - ωz∩Ω ℓ(z) v i dx| ≤ c h K (h -1 e v i - ωz∩Ω ℓ(K) v i dx L 2 (κ) + |v i | H 1 (κ) ) + (h -1 e v i - ωz∩Ω ℓ(z) v i dx L 2 (κ) + |v i | H 1 (κ) ) + h -1 2 e [v i ] L 2 (e) .
The two first terms yield as previously

h -1 e v i - ωz∩Ω ℓ(κ) v i dx L 2 (κ) + |v i | H 1 (κ) ) ≤ c h K α -1 2 κ v i 1,α,∆κ .
Thanks to assumption 1, the last term is bounded as follows

h K (h -1 2 e [v i ] L 2 (e) ) ≤ h K (h -1 2 e v i|K -R h v i L 2 (e) + ( h e he ) 1 2 ( h-1 2 e v i| K -R h v i L 2 (ẽ) )) ≤ c (h 1 2 e v i|K -R h v i L 2 (e) + h 1 2 e v i| K -R h v i L 2 (ẽ) )
The second term on the right hand side of this inequality is bounded for small δ > 0 as follows

(74) h 1 2 e v i| K -R h v i L 2 (ẽ) ≤ δ v i| K -R h v i 2 L 2 ( K) + 1 4δ h ẽ,
and we have obtained for K

v i| K -R h v i L 2 ( K) ≤ c h K α -1 2 K v i 1,α,∆ K .
Therefore the worst term is O(h K ). The first term is simply added to v i -R h v i L 2 (K) on the left hand side of the estimate thanks to the decomposition (74) which will still add to the right hand side a term of O(h K ).

• If z ∈ N + C , z is an endpoint of e which is an entire edge of K and only a part of ẽ an edge of K. We denote by a j , j = 1, 2, 3 the vertices of K, and by λ i , the associated barycentric functions,

v i|Ω ℓ(K) -R h v i L 2 (K) ≤ 3 j=1 |(v i -R h v i )(a j )| λ j L 2 ( K) ≤ c h K 3 j=1 |(v i|Ω ℓ(K) -R h v i )(a j )|,
Using Assumption 1, and the regularity of the triangulations, we obtain

v i|Ω ℓ(K) -R h v i L 2 (K) ≤ h 1 2 ẽ v i|Ω ℓ(z) -R h v i L 2 (ẽ) + h 1 2 e [v i ] L 2 (e) .
When the subdomains Ω ℓ(K) and Ω ℓ(z) are not adjacent, by using Assumption 2, we introduce the subdomains (Ω ℓ ) ℓ which are on the path between Ω ℓ(K) and Ω ℓ(z) . We apply the previous argument each pair of adjacent subdomains. This establishes the first estimate of the Lemma. The second one is proven in exactly the same way by noting that ϕ z L 2 (e) ≤ ch

Figure 1 .

 1 Figure 1. Zoom on the partition of the domain Ω. The plain curve is the contact zone γ

Figure 2 .

 2 Figure 2. Graphic representation of assumption 1

Lemma 4 . 8 . 1 2

 481 The following inequality holds(48) e (β e λ h + (n e • ασ(u h ) n e )) [w • n e ] dτ ≤ λ h + (n e • ασ(u h ) n e ) L 2 (e) K∈Ke uu α h 1,α,K ,where K e denotes the union of elements K sharing the edge e and β e = 0 if M h = M 1,ℓ, * h and = 1 otherwise. Proof. We have e (β e λ h +n e •ασ(u h ) n e ) [w • n e ] dτ ≤ β e λ h +n e •ασ(u h ) n e H -

1 2 . 4 . 3 .

 1243 An upper bound for the indicators. To complete the a posteriori analysis, we have to show the equivalence of the error and the error indicators which means that we should obtain upper bounds to each error indicator in term of the local discretization error. In order to bound the indicators η K , η e and η C,e , we set E = (e u , e λ ) = (uu h , λλ h ), then we take a test function w ∈ V and compute (52) a α (e u , w) + b(e λ , w) = Ω C α σ(e u ) : ε(w) dx

Corollary 4 . 19 .

 419 If Assumption 1 is satisfied, there exists a constant C independent of h and α such that (65) (

Figure 3 .

 3 Figure 3. Initial mesh

Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 8 .

 45678 Figure 4. Final mesh after 10 cycle of adaptivity

Figure 9 .

 9 Figure 9. Initial mesh Figure 10. Final mesh α = 10 -2 in Ω 2 ∪ Ω 3 .

Figure 11 .

 11 Figure 11. Convergence of the adaptivity

Figure 12 .

 12 Figure 12.Final mesh after 10 cycle of adaptivity α 2 = 10 -1 , α 3 = 10 -2

1 2 eϕ L 2 (κ) + h 1 2 e |ϕ| H 1 (

 12221 v iv i L 2 (e) + v i -ωz∩Ω ℓ(z) v i L 2 (e) ).Denoting by κ any of the two elements K, K ′ which share the edge e, and ℓ(κ) = ℓ(z) or ℓ(K), and thanks to the trace theorem [[START_REF]Verfürth-A Review of A Posteriori Error Estimation and Adaptive Mesh-Refinement Techniques[END_REF], Lemma 3.2] we have (73) ϕ L 2 (e) ≤ c(h -

v

  i|Ω ℓ(z) dτ -ωz∩Ω ℓ(z) v i | + h -1 e | e [v i ] dτ | , ≤ h K h -1 e e (|v i|Ω ℓ(K) -ωz∩Ω ℓ(K) v i dx|) dτ + h -1 e e (|v i|Ω ℓ(z) -ωz∩Ω ℓ(z)

, *

Appendix

In this appendix, we give the proof of proposition 4.7. The proof consists in constructing a regularized interpolation operator with the appropriate approximation properties. Such operators are based on local modified quasi-interpolation operators of Clément type [START_REF] Bernardi | Girault-A local regularization operator for triangular and quadrilateral finite elements[END_REF]. In our case the construction has to take into account the nonhomogeneous materials in the spirit of [START_REF] Bernardi | Verfürth-adaptive finite element methods for elliptic equations with nonsmooth coefficients[END_REF], the non-compatibility of the meshes at the cracks ∪ I i=1 γ i , and the new (and major) difficulty which is the discontinuity (even) of the solution of the continuous problem at the contact zone

, let ω z denote the support of the nodal basis function ϕ z . It is the union of all elements that have z as a vertex. With each z we associate ℓ(z) in {1, . . . , J}, such that

• z belongs to Ω ℓ(z) .

• α ℓ(z) is maximal among α j , j = 1, . . . , J such that Ω j contain z.

We denote by

the mean value of the function v on the set ω. Then, we set (70)

. In order to enforce the continuity at the vertices of N + C , we define the affine piecewise continuous function on γ

We have to distinguish between several cases.

• We consider a vertex z which is not contained in the boundary of any of two subdomains. Then from the Bramble-Hilbert inequality, we deduce

When z ∈ ∂Ω, similar computations with the Poincaré-Friedrichs inequality lead to the same estimate.


O Bareille
email: olivier.bareille@ec-lyon.fr

M Chabchoub

S Besset

M N Ichchou

AN INVERSE MID-HIGH FREQUENCY ENERGY METHOD : FORMULATION AND APPLICATIONS

An energy method called Simplified Energy Method (MES) has already been proposed to predict energy densities repartition for structural-acoustic problems in the mid-high frequency range. In order to illustrate but also to present one of the applications of this method, this latter is used here to solve inverse structural problems. The main application, in this paper, is the structural forces identification thanks to energy densities calculation and analysis. The injected forces estimation and localization are hence obtained in the mid-high frequency range. Internal as well as boundary sources are here detected through 2D Kirchhoff-plate numerical tests.

Introduction

Identification of forces features exciting structural components is an important concern in structural engineering. Indeed, localization and quantification of input loads is of great interest and often the direct measurements of sources reveals complicated to perform. Indirect means to define the excitation is often an alternative. A number of work have been done and published around that question [1, [START_REF] Ma | Input force estimation of beam structures by an inverse method[END_REF][START_REF] Ji | A study on an estimation method for applied force on the rod[END_REF][START_REF] Pezerat | Identification of vibration sources[END_REF][START_REF] Djamaa | Reconstruction of a distributed force applied on a thin cylindrical shell by an inverse method and spatial filtering[END_REF][START_REF] Pezerat | Identification of vibration excitations from acoustic measurements using near field acoustic holography and the force analysis technique[END_REF][START_REF] Perotin | An inverse method for the identification of a distributed random excitation acting on a vibrating structure: Part 2, flow induced vibration application[END_REF][START_REF] Pavic | Measurement of structure borne wave intenisty: Part I, formulation of the methods[END_REF][START_REF] Bai | Application of BEM-based acoustic holography to radiation analysis of sound sources with arbitrarily shaped geometries[END_REF]. The problem considered here can have also some applications in structural health monitoring. The damage, whose features are often unknown, leads to modifications of the vibration energy distribution of structural members. The damage size and position detection is then a key issue.

By and large, two ways of thinking are available in the open literature concerning input forces identifications issues. The first way of thinking is roughly inspired from identification issues of blackbox like behavior of dynamical systems and the use of estimation theories [1,[START_REF] Ma | Input force estimation of beam structures by an inverse method[END_REF]. The second way of thinking consists in defining the input force that allows the measured set of inputs to verify the equilibrium equations given either in a continuous or discretised format [START_REF] Pezerat | Identification of vibration sources[END_REF][START_REF] Djamaa | Reconstruction of a distributed force applied on a thin cylindrical shell by an inverse method and spatial filtering[END_REF][START_REF] Pezerat | Identification of vibration excitations from acoustic measurements using near field acoustic holography and the force analysis technique[END_REF]. The problem of vibration input force localization was also dealt with in structural intensity community. Moreover, instead of using analytical model [START_REF] Pezerat | Identification of vibration sources[END_REF] or modal modal [START_REF] Perotin | An inverse method for the identification of a distributed random excitation acting on a vibrating structure: Part 2, flow induced vibration application[END_REF] in the inversion process, discretized finite element (FE) or boundary element (BE) can be introduced in the identification [START_REF] Bai | Application of BEM-based acoustic holography to radiation analysis of sound sources with arbitrarily shaped geometries[END_REF]. Ultimately, in [START_REF] Pavic | Measurement of structure borne wave intenisty: Part I, formulation of the methods[END_REF] intensity vectors indicating the energy flow and the localization of the sources were used.

It should be however noticed that most of what was proposed in the literature concerns a limited frequency band of interest; precisely, the so-called low frequency region where modal overlap is weak and the resonant behavior of the structure still dominant. Energy methods are often used as alternative methods in high and medium frequency range when numerical methods like finite element or boundary equation based formulations run to important limitations. Among these energy methods, the most recognized remains Statistical Energy Analysis (SEA) [START_REF] Lyon | Statistical Energy Analysis of Dynamical Systems: Theory and Application[END_REF] Many energy methods have been based on SEA, attempting to enhance the Statistical Energy Analysis robustness and predictivity , leading to the Simplified Energy Method (MES) [START_REF] Nefske | Power Flow Finite Element Analysis of Dynamic Systems: Basic Theory and Application to Beams[END_REF][START_REF] Wohlever | Mechanical Energy Flow Models of Rods and Beams[END_REF][START_REF] Lase | Energy Analysis of bars an beams: Theoretical Formulations[END_REF], as a "direct" prediction tools [START_REF] Lase | Energy Analysis of bars an beams: Theoretical Formulations[END_REF][START_REF] Besset | A coupled BEM and energy flow method for mid-high frequency internal acoustic[END_REF].

In this paper, the main concern is the use of MES to predict energy sources [START_REF] Hardy | Absorption coefficient and energy flow path identification by means of inverse local energy method[END_REF]. Such an issue needs the inversion of the direct formulation. The paper is structured as follows. For a sake of clarity, the direct MES formulation is reminded in section 2. Then the formulation of the inverse MES is given. A two dimensional case is considered and the algebraic matrix format of the inversion problem is offered in section 3. In order to validate the proposed method, numerical simulations are processed section 4, using estimated direct energies or a finite-element simulated input energy.

Overview of the direct energy flow method MES

Assumptions

Simplified Energy Method is based on the description of two local energy quantities. The first one is the total energy density W defined as the sum of the potential energy density and the kinetic energy density. The second energy variable I is the energy flow. The energy balance of a system can be written as follows:

∇ • I + π diss = π inp (1)
where ∇ is the gradient operator, π diss is the dissipated power density and π inp is the input power density. MES adopts the same damping model as SEA. It has already been discussed in the literature [START_REF] Lyon | Statistical Energy Analysis of Dynamical Systems: Theory and Application[END_REF] and can be expressed as follow:

π diss = ηωW ( 2
)
where η is the damping loss factor and ω the circular frequency. The propagating waves considered in the MES formulation are made of partial energy quantities corresponding to direct and reverberated fields. As the considered fields are quadratic variables, the superposition principle can be applied:

£ = i £ i (3)
where (£, £ i ) are quadratic variables corresponding to global energy or a partial energy associated to wave fields (W, W i ) or I, I i . In the following, an intrinsic energy law is used, which is often introduced to define the wave velocity c:

I i = c • W i n (4)
This expression simply defines the energy velocity, which is in some cases the group velocity, as the ratio of the energy flow and convected energy density. This expression was demonstrated in a number of textbooks and is respected for any kind of linear elastic propagative wave.

Direct formulation

In the following, the derivation of the MES formulation is offered considering a nDimensional space (nD) space. Considering a symmetrical waves field and equations ( 4) and (1), the MES approach [START_REF] Ichchou | Comments on simple models of the energy flow in vibrating membranes and transversely vibrating plates[END_REF] leads to the following relation:

1 r n-1 ∂ ∂r r n-1 I + ηωW n = 0 (5)
which can also be written using equation ( 4):

-c 2 1 r n-1 ∂ ∂r r n-1 W n = ηω I (6)
Finally, the equation can be rewritten considering only the W field and provides a local energy equation for symmetrical waves energies.

1 r n-1 ∂ ∂r r n-1 W + ηω c W = 0 (7)
Energy variables are defined as the superposition of a direct field and a reverberated field which can be added to obtain the total field.

The elementary solutions in terms of energy density and active intensity are denoted G and H.

G(S, M) = 1 SM e -ηω c SM H(S, M) = c SM e -ηω c SM u SM (8)
The kernels are determined in each point M of the considered space and created by a source located in S. Hence, using a superposition principle the energy solution is then expressed using the contribution of the primary source ρ (direct field) and fictitious sources σ. The fictitious sources are energy unknown parameters representing the contribution of the reverberant field. Precisely:

W (M) = Ω ρ(S)G(S, M)dS + ∂Ω σ(P ) u SM n P G(P, M)dP (9)
Ω is the system surface and ∂Ω is the system boundaries.

The active energy flow can be described by a similar expression:

I(M) = Ω ρ(S)H(S, M)dS + ∂Ω σ(P ) u SM n P H(P, M)dP (10)
where ρ corresponds to internal sources to the considered domain (here a two dimensional one) and σ corresponds to sources at the boundaries of the two dimensional space.

The inverse energy flow method IMES formulation

In what follows, a discretized format of the energy integral equations are employed. To that end, we define firstly T as the (N, N) geometric interaction matrix between N two dimensional edge elements. This operator is given below referring to figure 1. Referring to the same figure 1, another needed (N, k) matrix is defined and named here Q. This operator is representing the geometric interaction between N boundary elements and k sources. It is expressed as follows:

T ij = ∂Ω ∂Ω 1 2 1 P i P j e -ηω c P i P j u P i P j n P i dP i dP j (11)
Q ij = ∂Ω 1 2 1 S j P i e -ηω c S j P i u S j P i n P i dP i (12)
Then, we define two matrices named Y W and Y I (having a dimension (k, m)) representing a geometric interaction between sources and measurement points. These matrices are determined respectively for the energy density W and the energy flow I as follows:

Y W ij = ∂Ω 1 S j P i e -ηω c S j P i dP i and Y I ij = 1 c ∂Ω 1 S j P i e -ηω c S j P i u S j P i dP i (13)
Ultimately, we define A and B as follows (Id being the identity matrix and ρ injected power vector):

A = Id -T and B = Qρ (14)
The power vector σ can thus be readily expressed as a function of injected power vector ρ as:

σ = A -1 B = A -1 Qρ (15)
Combining equations ( 9), ( 10) and ( 15), the energy density vector W and the energy flow vector I can be obtained from:

W = S W ρ = R W ρ + Y W σ and I = S I ρ = R I ρ + Y I σ (16)
where R W and R I are defined as follows:

R W ij = 1 S i M j e --ηω c S i M j and R I ij = 1 S i M j e --ηω c S i M j u S i M j (17)
S W and S I given by:

S W = R W + Y W A -1 Q and S I = R I + Y I A -1 Q (18)
Finally, both W and I fields can be expressed thanks to a linear operator denoted F :

W, I = F (ρ) (19)
The IMES formulation aims at inverting operator F proposed in equation ( 19). Sources have to be detected thanks to measurements made on the plate. As the number of measurements is not equal to the number of sources, the problem is not directly invertible. Thus, we look for the G operator that gives σ and ρ as a function of W and I.

ρ = G W, I (20)
Let now introduce the quantity R (X ) defined as follows:

R (X ) = W, I -F oX W, I (21)
We now have to minimize the quantity R (X ) among operator X . In the following, the IMES will be applied to a 2D plate. Thus, Ω corresponds to the surface of the plate, whereas ∂Ω corresponds to the boundaries. In this case, only surface sources are treated. In what follows, The access to structural intensities experimentally being a hard task to achieve, simulations have been made considering only energy densities W measurements in some locations of the structure. Thus, we will replace S W by S for the sake of clarity. MES equation ( 9) is discretized to provide the following matrix formulation:

           W 1 W 2 . . . W m            =     
                ρ 1 ρ 2 . . . ρ k            ( 22
)
where m is the number of the measurements and k is the number of the sources. To invert the formulation, as explained in equation ( 21), we have to find a matrix X that minimize the quantity R(X):

R(X) = (I m -SX)            W 1 W 2 . . . W m            (23)
In fact, we choose matrix X as a pseudo-inverse of matrix S:

X = S T S -1 S T (24)
The inverse formulation can then be performed. From a set of input energy parameters, the inversion can be achieved in order to characterize input sources.

Numerical simulations and Results

In what follows, a number of numerical tests are performed in order to validate the formulation. The first considered test case is very simple. It consists in achieving the inverse energy flow approach when a set of energy prediction using the direct energy flow method are used. This is an MES/IMES simulation. A more realistic test case is then considered. The input variables are in this situation computed from a full finite element model.

MES/IMES simulation

A first simulation named MES/IMES simulation is made. IMES calculations are processed using energy densities W coming from the direct MES formulation. The model studied here is a thin plate whose characteristics are mentioned in table (1). 2. Error on the input power ξ Error 0.06 30 % 0.14 25 % 0.5 18 % It appears that the source are correctly identified in this simulation (localization and level). So, the inversion of S T S is validated. To illustrate the effect of the number of measurements m, we repeated the MES/IMES simulation while decreasing m. Let's then introduce the factor ξ = m N to quantify the efficiency of the method (ξ = 1 corresponds to a square matrix to be inverted).

It appears that by increasing the number of measurements m spread over the plate surface and corresponding to ξ = 0.06, ξ = 0.14 and ξ = 0.5 respectively, the source is well located and the error on the input power summarized in Table [START_REF] Ma | Input force estimation of beam structures by an inverse method[END_REF] decrease. Nevertheless, the input power of the source is equal to 0.7 W/m 2 for ξ = 0.06, 0.75 W/m 2 for ξ = 0.14 and 0.82 W/m 2 for ξ = 0.5 whereas it should be 1 W/m 2 . The right value of the input power is obtained for higher values of ξ. The IMES seems to correctly localize the input structural force but underestimates the injected power.

FEM/IMES simulation

The IMES/MES subsection allowed mainly the validation S T S inversion. In this subsection the validation of approach from a FEM simulated energies is considered. In fact, IMES calculations are processed using energy densities W coming from a finite element software (Comsol in the present study). An FEM model of the plate (characteristics provided below) are processed with 1W/m 2 excitation and the same number of elements N. Given the input power ρ, the force F applied to a point of the plate is recalculated by:

ρ = 1 2 F 2 8 √ Dρ s h ( 25
)
Where D is the plate stiffness defined by D = Eh 3 12(1-ν 2 ) . FEM calculation is processed on the plate, and provides energy density W field used as an input for IMES algorithm. FEM calculations are first averaged first space and frequency averaged. Octave bands are considered and various center frequencies f c are considered. The values of W recovered by FEM calculation and realized around f c feed the IMES algorithm. Results found for ξ = 0.06, ξ = 0.14, ξ = 0.5 and ξ = 1 are presented in figure 3.

Moreover, the influence of f c is also studied. Results found for ξ = 1 and for different values of f c are given in figure 4. This figure shows a good localization of the source whatever the value of f c belonging to the range of medium and high frequencies.

Conclusion

The aim of this paper is to analyze the effectiveness of the Simplified Energy Method (MES) to identify sources from a set of given measured energies in the plate. The MES formulation and assumptions were first detailed. Then, the Inverse formulation identifying sources through data energy recovered from the plate is given. After that, the MES is validated with FEM/IMES simulations taking numerical data from FEM software.

Validation of the strategy for 3D acoustic cavities is under investigation.

Figure 1 .

 1 Figure 1. element/element interaction

Figure 2 .

 2 Figure 2. MES/IMES simulation for the plate

Figure 3 .

 3 Figure 3. FEM/IMES simulations for ξ = 0.06 (a), ξ = 0.14 (b), ξ = 0.5 (c), ξ = 1 (d)

Figure 4 .

 4 Figure 4. FEM/IMES simulations: f c = 950Hz (a) and f c = 1500Hz (b)



  S 11 S 12 . . . S 1k S 21 S 22 . . . S 2k

				
	. . .	. . .	. . .	. . .
	S m1 S m2 . . . S mk

Table 1 .

 1 plate caracteristicsThis plate is excited by an input power of 1W/m 2 . Considering N to be the number of edge plate elements, IMES result found for N = m = 64 is presented in figure2.

	parameter	unit	value
	Length	m	1
	Width	m	1
	Height Young's modulus of elasticity E	m P a	2e -3 2.1e11
	Poisson's ratio ν		0.3
	Solid density ρ s	Kg/m 3 7800