Stéphane Abide
email: stephane.abide@univ-perp.fr

Mikaël Barboteu
email: barboteu@univ-perp.fr

David Danan
email: david.danan@univ-perp.fr

Analysis of two active set type methods to solve unilateral contact problems

Keywords: AMS Subject Classification : 74M15, 74G15, 74B05, 74B20, 74S30, 49M15, 90C53, 90C25, 65K15 Unilateral constraint, Hertzian contact, Hyperelasticity, Fixed point equation, Weak solution, Finite element approximation, Projection iterative method, Primal dual Active set, Augmented Lagrangian, Numerical simulations, Error estimates

In this work two active set type methods are considered in order to solve a mathematical problem which describes the frictionless contact between a deformable body and a perfectly rigid obstacle, the so-called Signorini Problem. These methods are the primal dual active set method and the projection iterative method. Our aim, here, is to analyse these two active set type methods and to carry out a comparison with the well-known augmented Lagrangian method by considering two representative contact problems in the case of large and small deformation. After presenting the mechanical formulation in the hyperelasticity framework, we establish weak formulations of the problem and the existence result of the weak solution is recalled. Then, we give the finite element approximation of the problem and a description of the numerical methods is presented. The main result of this work is to provide a convergence result for the projection iterative method. Finally, we present numerical simulations which illustrate the behavior of the solution and allow the comparison of the numerical methods.

Introduction

Various contact boundary conditions have been used to model contact phenomena, both in engineering and mathematics in the literature; see for instance [START_REF] Eck | Unilateral Contact Problems: Variational Methods and Existence Theorems[END_REF][START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF][START_REF] Panagiotopoulos | Inequality Problems in Mechanics and Applications[END_REF][START_REF] Shillor | Models and Variational Analysis of Quasistatic Contact[END_REF][START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF] and the references therein for the mathematical analysis and [START_REF] Haslinger | Numerical Methods for Unilateral Problems in Solid Mechanics[END_REF][START_REF] Hlaváček | Solution of Variational Inequalities in Mechanics[END_REF][START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF][START_REF] Laursen | Computational Contact and Impact Mechanics[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF] for the numerical analysis. One of the most popular remains the so-called Signorini condition, introduced in [START_REF] Signorini | Sopra alcune questioni di elastostatica[END_REF], which describes the contact with a perfectly rigid foundation. Expressed in terms of unilateral constraints for the displacement field, this condition leads to a nonlinear and nonsmooth mathematical problem. To overcome these difficulties, several methods were used. For instance, Alart and Curnier presented in [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF] an augmented Lagrangian formulation combined with a Generalized Newton method to solve non differentiable but continuous equations arising from frictional contact problems. Over the last few years, several other methods emerged. Amongst them, the active set strategies have been relatively popular for the last decade or so. The aim of this kind of method is to find the correct subset A of all nodes that are currently in contact with the perfectly rigid obstacle; those nodes are called active while the others are inactive, see [START_REF] Hintermuller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF][START_REF] Hueber | A primal dual active set strategy for non-linear multibody contact problems[END_REF] and references therein for more details. Also, one of the most interesting aspects of such a method is that it does not require the use of the Lagrange multipliers and, therefore, could facilitate the implementation of the algorithm and improve the condition number of the system. Indeed, since the boundary conditions on the contact boundary for the active nodes are directly enforced, we only get to solve a series of problems with simple boundary conditions, such as Dirichlet, Neumann or Robin boundary condition, based on the active set type method used. The purpose of the present work is to study and compare two active set type methods, the primal dual active set method and the projection iterative method, in the case of large and small deformation theories. The first one is the primal dual active set method used in various works such as [START_REF] Hintermuller | The primal dual active set strategy as a semismooth Newton Method[END_REF][START_REF] Hintermuller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF] and [START_REF] Hueber | A primal dual active set strategy for non-linear multibody contact problems[END_REF] and the reference therein. In [START_REF] Hintermuller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF], the Authors study a class of semismooth Newton methods for quadratic minimization problems to non negativity constraints in which global and local non-linear convergence results of the resulting primal dual active set strategy were established under strong assumptions on the matrix of the linearized systems. The proof is based on the M-matrix properties of the discrete operator; while such a statement is true for the Laplace operator it is not for the discrete elasticity operator, as the authors admitted. Furthermore, in that work, only linear elasticity problems with unilateral boundary constraints are considered. The work in [START_REF] Hueber | A primal dual active set strategy for non-linear multibody contact problems[END_REF] is devoted to provide an inexact primal dual active set approach to solve nonlinear multibody contact problems for linear elasticity. The non-linearity of such problems arise only from the non penetration conditions for the involved bodies in contact. The second active set type method was used for instance in [START_REF] Zhang | A projection iterative algorithm boundary element method for the Signorini Problem[END_REF] to solve the so-called Signorini problem for the Laplace equations and to obtain a convergence analysis of the method. This active set method is based on a fixed point equation and can be formulated as a projection iterative algorithm. The projection iterative method consists to solve a sequence of Dirichlet or Robin boundary conditions according to a contact criteria in order to find the correct sets of active and inactive contact nodes. Our aim in this work is twofold. First, we present and analyse the active set type methods compared to the well-known augmented Lagrangian approach to solve a Signorini contact problem for hyperelasticity. Furthermore, we extend the convergence result obtained in [START_REF] Zhang | A projection iterative algorithm boundary element method for the Signorini Problem[END_REF] for the projection iterative method to a Signorini problem in the case of linear elasticity. Next we provide a numerical comparison of the different methods by considering simulations on two-dimensional test problems: one in the small deformation framework, the Hertz contact problem, and one in the large deformation framework, the contact between a hyperelastic ring and a rigid foundation under strong compressions. In particular, we analyze and compare the numerical convergence of the different methods with respect to several parameters such as the number of degrees of freedom or the number of iterations. By doing so, we also illustrate the behavior of the solution related to the contact condition.

The rest of the paper is structured as follows. In Section 2 we describe the contact conditions and introduce the mechanical problem in the large deformation framework. Then we introduce the notation and some preliminary material, list the assumptions on the data and state the variational formulation of the problem in the hyperelasticity framework. In Section 3, we provide the finite element approximation of the variational formulation. After reformulating the problem into a minimization one, an augmented Lagrangian method is recalled to treat the unilateral constraints. In Section 4, we describe the two active set type methods used to solve unilateral contact problems in hyperelasticity: the primal dual active set method and the projection iterative method. The Section 5 is devoted to the convergence analysis of the projection iterative method to the solution of the Signorini Problem, in the small deformation hypothesis. After that, in Section 6 we present several numerical simulations to illustrate and compare the behavior of the two active set type methods. Finally, in Section 7, we conclude by recalling the results obtained and discussing about the conceivable works in the continuation of this one.

Formulations of the contact problem

The purpose of this section is to present both the mechanical problem and the variational formulation in the framework of hyperelasticity.

Hyperelastic contact model

The physical setting of the mechanical problem is as follows. A hyperelastic body occupies a bounded domain Ω ⊂ R d (d = 1, 2, 3) with a Lipschitz continuous boundary Γ, divided into three disjoint measurable parts Γ 1 , Γ 2 and Γ 3 such that meas (Γ 1 ) > 0. The notation x = (x i ) for a typical point in Ω ∪ Γ is used and we denote by ν = (ν i ) the outward unit normal at Γ. Here and below the indices i, j, k, l run between 1 and d and, unless stated otherwise, the summation convention over repeated indices is used. We denote by M d the space of second order tensors on R d or, equivalently, the space of square matrices of order d. The inner product and norm on R d and M d are defined by

u • v = u i v i , v = (v • v) 1 2 ∀ u, v ∈ R d , Π • τ = Π ij τ ij , τ = (τ • τ ) 1 2 ∀ Π, τ ∈ M d .
We use the notation u and Π for the displacement field and the first Piola-Kirchoff stress tensor, respectively. Also, we denote by u ν and u τ the normal stress and tangential components of u on Γ given by

v ν = v • ν, v τ = v -v ν ν.
Finally, Π ν and Π τ will represent the normal and the tangential stress on Γ, defined by Π ν = (Πν) • ν and Π τ = Πν -Π ν ν. Furthermore, an index that follows a comma represents the partial derivative with respect to the corresponding component of the spatial variable x, e.g. u i,j = ∂u i /∂x j . Moreover, we recall that the divergence operator is defined by the equality Div Π = (Π ij,j ).

In the problems studied below, the material's behavior is described with a hyperelastic constitutive law. We recall that hyperelastic constitutive laws are characterized by the first Piola-Kirchhoff tensor Π which derives from an internal hyperelastic energy density

W : Ω × M d + → R, Π = ∂ ∂F W (x, F) = ∂ F W (x, F), for all x ∈ Ω and F ∈ M d .
Here F is the deformation gradient defined by F = I + ∇u and ∂ F represents the differential with respect to the variable F, see [START_REF] Ciarlet | I : Three-Dimensional Elasticity[END_REF][START_REF] Tallec | Numerical Methods for Nonlinear Three-dimensional Elasticity[END_REF] for details on hyperelasticity. In what follows, we consider the contact without friction of a hyperelastic body with a perfectly rigid obstacle, the so-called foundation. The hyperelastic body is in equilibrium under the action of body forces of density f 0 and surface tractions of density f 2 which act on Γ 2 . We also assume that the body is fixed on Γ 1 and may come in contact over Γ 3 with the foundation. The frictionless contact conditions are based on the unilateral contact conditions without friction and can be written following the Karush-Kuhn-Tucker conditions on Γ 3 :

u ν ≤ 0, (2.1)
Π ν ≤ 0, (2.2)

u ν Π ν = 0, (2.3)
Π τ = 0. (2.4)
With these preliminaries, the formulation of the contact problem we consider in this work is as follows.

Problem P. Find a displacement field u : Ω → R d and a stress field Π : Ω → M d such that

Π = ∂ F W (F) in Ω, (2.5
)

Div Π + f 0 = 0 in Ω, (2.6
)

u = 0 on Γ 1 , (2.7
)

Πν = f 2 on Γ 2 , (2.8
)

u ν ≤ 0, Π ν ≤ 0, u ν Π ν = 0 on Γ 3 , (2.9
)

Π τ = 0 on Γ 3 . (2.10)
Here and below, in order to simplify the notation, we do not indicate explicitly the dependence of various functions on the spatial variable x. Equation (2.5) represents the hyperelastic constitutive law of the material. Equation (2.6) is the equation of equilibrium in the static case. Conditions (2.7) and (2.8) represent the displacement and traction boundary conditions, respectively. Finally, conditions (2.9) and (2.10) represent the frictionless contact conditions with unilateral constraint previously described in this section. Note that the conditions (2.9) are equivalent to the following subdifferential inclusion (cf [START_REF] Moreau | Application of convex analysis to some problems of dry friction[END_REF])

-Π ν ∈ ∂Ψ IR -(u ν ) on Γ 3 , (2.11)
where ∂ represents the subdifferential operator in the sense of the convex analysis and Ψ A denotes the indicator function of the set A ⊂ R.

Variational formulation

In order to derive the variational formulation of Problem P, also called in the literature the Signorini Problem, further notation and some preliminary material are needed. Everywhere in this paper, the standard notations for Sobolev and Lebesgue spaces associated to Ω and Γ are used. The following spaces are considered

V = {v ∈ H 1 (Ω; R d ) : v = 0 on Γ 1 }, H = L 2 (Ω; R d ).
These are real Hilbert spaces endowed with their standard inner products (u, v) V and (Π, τ ) H and their associated norms • V and • H , respectively. Note that V ⊂ H ⊂ V * is an evolution triple, with all embeddings being continuous, compact and dense. The duality pairing between V * and V will be denoted by u, v V * ×V . Also, recall that H 1 = H 1 (Ω) d algebraically and topologically. Completeness of the space (V, • V ) follows from the assumption meas(Γ 1 ) > 0, which allows the use of Korn's inequality. We also recall that there exists c 0 > 0 which depends on Ω, Γ 1 and Γ

3 such that v L 2 (Γ3) d ≤ c 0 v V for all v ∈ V.
(2.12) Inequality (2.12) represents a consequence of the Sobolev trace theorem. It is well known that the weak solution of problem P can be defined from a variational inequality [START_REF] Hlaváček | Solution of Variational Inequalities in Mechanics[END_REF][START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF], or can be obtained equivalently from a minimization problem on a convex set, see [START_REF] Boieri | Existence, Uniqueness, and Regularity Results for the Two-Body Contact Problem[END_REF][START_REF] Haslinger | Contact between two elastic bodies -I. continuous problems[END_REF]. In this section, a saddle point formulation within the form of a hybrid variational formulation is considered. In this case, the Lagrange multiplier λ ν related to the contact stress verifies an extended subdifferential inclusion derived from the pointwise subdifferential inclusion defined in (2.11). Note that, for convenience, the Lagrange multiplier λ ν is taken as equal to -Π ν . To this end, a trace space is defined by

X ν = { v ν | Γ3 : v ∈ V }
equipped with its usual norms. We denote by X ′ ν the duals of the spaces X ν . Moreover, we denote by •, • X ′ ν ,Xν the corresponding duality pairing mapping. Then a function ϕ ν : X ν → (-∞, +∞] is introduced, defined by

ϕ ν (v ν ) = Γ3 Ψ IR -(v ν ) dΓ ∀ v ν ∈ X ν .
Therefore, the pointwise inclusion condition (2.11) leads to the extended subdifferential inclusion

λ ν ∈ ∂ϕ ν (u ν ) in X ′ ν . (2.13)
For an element v ∈ H 1 , v is still written for the trace of v. Also, for a regular stress function Π the following Green's formula holds:

Ω Π : ∇v dx + Ω Div Π • v dx = Γ Πν • v da for all v ∈ H 1 . (2.14)
Moreover, in the study of the mechanical problem (2.5)-(2.10) we assume that the body forces and tractions densities have the regularity

f 0 ∈ L 2 (Ω) d , f 2 ∈ L 2 (Γ 2 ) d . (2.15)
Now we turn to the hybrid variational formulation of Problem P. To this end, we introduce the non linear hyperelastic operator A : V → V * , and the element f defined by

A(u), v V * ×V = Ω Π : ∇v dx ∀ u, v ∈ V, (2.16)
(f , v) V = (f 0 , v) H + (f 2 , v) L 2 (Γ2) , ∀v ∈ V.
(2.17)

We use Green's formula (2.14) and definition (2.17) to see that

A(u), v V * ×V = (f , v) V + Γ3 Π ν v ν da + Γ3 Π τ • v τ da.
Then by using the Lagrange multiplier λ ν , related to the normal contact stress Π ν , and the condition (2.10), we obtain the variational formulation of the frictionless contact problem P in terms of two unknown fields.

Problem P V . Find a displacement field u and a normal stress field λ = λ ν ν such that u ∈ V, A(u), v V * ×V + λ ν , v ν X ′ ν ,Xν = (f , v) V ∀ v ∈ V, (2.18)
λ ν ∈ ∂ϕ ν (u ν ) in X ′ ν . (2.19)
A pair (u, λ) which satisfies (2.18), (2.19) and the hyperelastic constitutive law in the form Π = ∂ F W (F) is called a weak solution to the frictionless contact problem P.

The solvability of the Problem P V can be obtained by using an existence result for a three dimensional hyperelastic Signorini contact problem, provided in [8, p. 381]. The result recalled below is based on the polyconvexity and the coerciveness of the stored density energy W .

Theorem 2.1. (existence result for a Signorini contact problem) Let Ω be a domain in R 3 , and let W : Ω × M 3 + → R + be a stored energy function with the following properties: (a) Polyconvexity: For almost all x ∈ Ω, there exists a convex function W(x, •) :

M 3 ×M 3 ×R + * → R such that W(x, F, Cof F, detF) = W (x, F) for all F ∈ M 3 + ; (2.20) the function W(•, F, H, δ) : Ω → R is measurable for all (F, H, δ) ∈ M 3 × M 3 × R + * . (b) lim detF→0 + W (x, F) = +∞, for almost all x ∈ Ω.
(c) Coerciveness: There exist constants α, β, p, q, r such that

α > 0, p ≥ 2, q ≥ p p -1 , r > 1, (2.21) W (x, F) ≥ α( F p + Cof F q + detF r ) + β (2.22)
for almost all x ∈ Ω and for all F ∈ M 3 .

Let Γ 1 , Γ 2 and Γ 3 be disjoint relatively open subset of Γ with meas(Γ 1 ) > 0 and meas({Γ

-(Γ 1 ∪ Γ 2 ∪ Γ 3 )}) = 0
, and let u 0 : Γ 1 → R 3 be a measurable function such that the set

Φ = {v ∈ W 1,p (Ω); Cof ∇v ∈ L q (Ω), det ∇v ∈ L r (Ω), v = u 0 on Γ 1 , det ∇v > 0 a.e. in Ω} is nonempty. Let the linear form L defined by L(v) = (f , v) V be continuous on W 1,p (Ω). Finally, let I(v) = W (x, ∇v(x))dx + Γ3 Ψ IR -(v ν ) dΓ -L(v)
, and assume that inf v∈Φ I(v) < +∞. Then there exists at least one function u such that v ∈ Φ and

I(v) = inf v∈Φ I(v). (2.23)
Remark 1. The solvability of problem P V can be obtained directly by applying the Theorem 2.1 in the case where p = 2, q = 2, r = 2 and u 0 = 0. It is easy to see that problem P V with d = 3 can be derived from the minimization problem (2.23), see [START_REF] Ayyad | Formulation and analysis of two energy-consistent methods for nonlinear elastodynamic frictional contact problems[END_REF][START_REF] Ciarlet | I : Three-Dimensional Elasticity[END_REF][START_REF] Tallec | Numerical Methods for Nonlinear Three-dimensional Elasticity[END_REF].

Remark 2. In the case of linear elasticity, if we assume that the elasticity operator is a strongly monotone Lipschitz continuous operator on the space V (which is not true for hyperelasticity) then, as proved in [27, p. 126], for each f ∈ V there is a unique weak solution to the contact problem P V .

Variational approximation and Lagrangian method

This section is devoted to provide a discrete minimization approach in order to solve the variational problem P V . To this end, first the discretization of the variational Problem P V is presented and next, the augmented Lagrangian approach based on an optimization formulation is recalled.

Discretization step

Let Ω be a polyhedral domain. Moreover, let {T h } be a regular family of triangular finite element partitions of Ω that are compatible with the boundary decomposition Γ = Γ 1 ∪ Γ 2 ∪ Γ 3 , i.e., if one side of an element T r ∈ T h has more than one point on Γ, then the side lies entirely on Γ 1 , Γ 2 or Γ 3 . The space V is approximated by the finite dimensional space V h ⊂ V of continuous and piecewise affine functions, that is,

V h = { v h ∈ [C(Ω)] d : v h | T r ∈ [P 1 (T r)] d ∀ T r ∈ T h , v h = 0 at the nodes on Γ 1 },
where P 1 (T r) represents the space of polynomials of degree less or equal to one in T r and h > 0 denotes the spatial discretization parameter. For the discretization of the Lagrange multiplier space X ′ ν , we use discontinuous piecewise constant functions as it is done in [START_REF] Barboteu | A Hyperelastic Dynamic Frictional Contact Model with Energy-Consistent Properties[END_REF][START_REF] Khenous | Hybrid discretization of the Signorini problem with Coulomb friction. Theoretical aspects and comparison of some numerical solvers[END_REF]. The discrete Lagrange multiplier space is denoted by Y h ν . Thus, let us consider the following discrete variational problem. Problem P h V . Find a discrete displacement field u h and a discrete normal stress field λ h = λ h ν ν such that

u h ∈ V h , A(u h ), v h V * ×V + λ h ν , v h ν X ′ ν ,Xν = (f , v h ) V ∀ v h ∈ V h , (3.1)
λ h ν ∈ ∂ϕ ν (u h ν ) in Y h ν . (3.2)
For details about the discrete formulation, we refer, e.g., to [START_REF] Ben Belgacem | Extension of the mortar finite element method to a variational inequality modeling unilateral contact[END_REF][START_REF] Belgacem | Hybrid finite element methods for the Signorini problem[END_REF][START_REF] Laursen | Computational Contact and Impact Mechanics[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF].

Minimization formulation and augmented Lagrangian method

The discrete weak solutions of Problem P h V can also be obtained from a minimization problem on a convex set which leads to an unconstrained augmented Lagrangian method. This method was described in detail in [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF][START_REF] Laursen | Computational Contact and Impact Mechanics[END_REF][START_REF] Wriggers | Computational Contact Mechanics[END_REF]. The starting point of this formulation is to consider the following optimization problem.

u h = argmin v h ∈V h J(v h ) + ϕ ν (v h ν ) , (3.3)
where J(v h ) is the hyperelastic strain energy potential defined by

J(v h ) = 1 2 A(v h ), v h V * ×V -(f , v h ) V . (3.4)
Therefore, with the formalism of the standard Lagrangian, the problem (3.4) is equivalent to the minmax problem

(u h , λ h ν ) ∈ arg(min v h max γ h ν ∈IR + {L(v h , γ h ν )}), (3.5)
where λ h ν ∈ Y h ∩ R + stands for the Lagrange multiplier while γ h ν is the virtual Lagrange variable, which represents the normal contact stress. The Lagrangian function L is defined by

L(v h , γ h ν ) = J(v h ) + v h ν , γ h ν Xν ×X ′ ν . (3.6)
We use now an augmented Lagrangian approach which permits to derive a smooth minimization problem without constraints. At this point, the regularized augmented Lagrangian expression is given by

L r (v h , γ h ν ) = J(v h ) - 1 2r γ h ν 2 + 1 2r [γ h ν + rv h ν ] 2 + , (3.7)
with r > 0, a penalty factor and [.] + , the positive part of an element of R i.e. [a] + = max(a, 0) for a ∈ R . Then the following augmented Lagrangian minimization problem is obtained

(u h , λ h ν ) ∈ arg(min{L r (v h , γ h ν )}), where v h ∈ V h and γ h ν ∈ Y h ν . (3.8)
Assuming that J is continuously differentiable with respect to v and strictly convex, L r is continuously differentiable with respect to v h and γ h ν and therefore, a saddle-point (u h , λ h ν ) verifies the following stationary condition

∇ u h ,λ h ν L r (u h , λ h ν ) = 0. (3.9)
It leads to a system of non linear equations of the form

           ∇ u h L r (u h , λ h ν ) = ∇ u h J(u h ) + [λ h ν + ru h ν ] + ν = 0, ∇ λ h ν L r (u h , λ h ν ) = - 1 r λ h ν -[λ h ν + ru h ν ] + = 0.
Note that from now on, in order to simplify the notation, we do not indicate explicitly the spatial discretization parameter on the unknown fields. We denote by u and λ the vectors representing the generalized vectors of displacement and Lagrange multiplier, defined by u

= {u i } N h T ot i=1 and λ = {λ i } N h Γ 3
i=1 , respectively, with N h T ot , the number of nodes of the discretization of Ω and N h Γ3 , the number of nodes of the discretization of Γ 3 . Here, u i represents the value of the corresponding function u h at the i th nodes of T h while λ i denotes the value of the corresponding function λ h = λ h ν ν at the i th nodes arising from the discretization of the contact interface Γ 3 denoted by Γ h 3 . Therefore, by adopting a vector form, we have the following system of non linear equations:

Problem P L . Find a displacement field u ∈ R d×N h T ot and a normal stress field λ = λ ν ν ∈ R d×N h Γ 3 such that G(u) + F(u, λ ν ) = 0 d×N h T ot +N h Γ 3 , (3.10)
where

G(u) ∈ R d×N h T ot × R N h Γ 3 is the hyperelastic part and F(u, λ) ∈ R d×N h T ot × R N h Γ 3
is the non differentiable contact part defined as follows

G(u) = A(u) -f 0 N h Γ 3
and

F(u, λ ν ) =    [λ ν + ru ν ] + ν - 1 r λ ν -[λ ν + ru ν ] +    , here 0 N h Γ 3 is the zero of R N h Γ 3 and 0 d×N h T ot +N h Γ 3 is the zero of R d×N h T ot +N h Γ 3
. The solution of the non linear Problem P L is based on a generalized Newton method which permits to treat both variables (u, λ) simultaneously. Now, by considering the pair x = (u, λ), the iterative scheme of index k can be summarized as follows. (i) Choose x 0 , set k = 0.

(ii) Set x (k+1) , such that

x (k+1) = x (k) -(K k + L k ) -1 (G(x (k) ) + F(x (k) )), L k ∈ ∂F(x (k) ), (3.11)
where ∂F(x (k) ) is the generalized Jacobian of F at x (k) and K k is the Jacobian of G at x (k) .

(iii) If x (k+1)x (k) ≤ ǫ and R L (x (k) ) ≤ ǫ then stop, else goto (ii).

Here, ǫ represents a small value and R L (x (k) ) is the nonlinear operator belonging to R d×N h T ot ×R N h Γ 3 , defined by

R L (x (k) ) = G(u (k) ) + F(u (k) , λ (k) ).
Remark 3. This method requires the consideration of additional immaterial nodes associated to the Lagrange multipliers, for the treatment of the contact operator F(u, λ). The construction of these fictitious nodes depends on the contact element used for the geometrical discretization of the interface Γ 3 . In the case of the numerical examples presented in Section 6, the discretization is based on "node-to-rigid" contact element, which is composed by one node of Γ h 3 and one Lagrange multiplier node.

Active set type methods

The aim of this section is to present the main traits of two active set type methods used to solve unilateral contact problems without friction between a hyperelastic body and a rigid obstacle. These two methods are the primal dual active set method and the projection iterative method. See [START_REF] Hueber | A primal dual active set strategy for non-linear multibody contact problems[END_REF] and [START_REF] Zhang | A projection iterative algorithm boundary element method for the Signorini Problem[END_REF] for a presentation of the methods in the case of linear elasticity and Laplace equations, respectively. Then, in each case, we present the method and detail the associated iterative algorithm for the hyperelastic Signorini problem.

Primal dual active set

The primal dual active set type method can be seen as a semismooth Newton method based on the reformulation of conditions (2.9), see [START_REF] Hintermuller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF]. As a matter of fact, the unilateral contact conditions (2.9) can be formulated in terms of a fixed point problem, as mentioned in [START_REF] Chouly | A Nitsche finite element method for dynamic contact : 1. Semi-discrete problem analysis and time-marching schemes[END_REF]. Let us recall and prove this result. Proposition 4.1. Let γ > 0, the unilateral contact conditions (2.9) are equivalent to:

Π ν = -[-Π ν + γu ν ] + on Γ 3 . (4.1)
Proof. Let's assume that (2.9) hold. Then, either u ν < 0 or u ν = 0. First, if u ν < 0, (2.3) implies that Π ν = 0. So

-[-Π ν + γu ν ] + = -[γu ν ] + = 0 and Π ν = 0 =⇒ Π ν = -[-Π ν + γu ν ] + .
Suppose now that u ν = 0 and Π ν < 0

-[-Π ν + γu ν ] + = -[-Π ν ] + = Π ν .
Conversely, if (4.1) holds, it implies that (2.2) holds. Next, if Π ν = 0 we have

-[γu ν ] + = 0, which means that u ν < 0, since γ > 0. At last, if Π ν < 0, -Π ν + γu ν > 0 and Π ν = -[-Π ν + γu ν ] + = Π ν -γu ν ,
so γu ν = 0 and, since γ > 0, u ν = 0. Therefore (2.1) and (2.3) hold, which concludes the proof. Now, let (u, λ) be the solution of the discrete variational problem P h V , with λ ν the Lagrange multiplier introduced in the previous part and defined by λ ν = -Π ν . The Problem P h V can be rewritten in a discrete algebraic form with the condition (4.1) for the Lagrange multiplier. Then, it leads to consider the system of non linear equations of the following form

R(u, λ) =   R u (u, λ) = A(u) + λ ν ν -f R λ (u, λ) = λ ν -[λ ν + γu ν ] +   =   0 0   , (4.2)
where, R(u, λ) is the generalized non linear operator belonging to R d×N h T ot × R N h Γ 3 . Therefore, by adopting the semismooth Newton formalism, the solution (u, λ) can be formulated as follows:

(i) Choose (u (0) , λ (0) ), set k = 0. (ii) Set (u (k+1) , λ (k+1) ) such that G u (u (k) , λ (k) )δu (k) = -R u (u (k) , λ (k) ), (4.3) ∆R λ (u (k) , λ (k) ) = -R λ (u (k) , λ (k) ), (4.4)
u (k+1) = u (k) + δu (k) , (4.5)
λ (k+1) = λ (k) + δλ (k) , (4.6)
where G u (u (k) , λ (k) ) stands for the differential of R u with respect to the variable u (k) as follows

G u (u (k) , λ (k) ) = ∂R u (u (k) , λ (k) ) ∂u (k) ,
and ∆R λ (u (k) , λ (k) ) represents the total differential of the complementarity function R λ with respect to the variables u (k) and λ (k) as follows

∆R λ (u (k) , λ (k) ) = ∂R λ (u (k) , λ (k) ) ∂(u (k) , λ (k) ) (δu (k) , δλ (k) ). (iii) If (u (k+1) , λ (k+1) ) -(u (k) , λ (k) ) ≤ ǫ and R(u (k) , λ (k)
) ≤ ǫ then stop, else goto (ii). The main traits of this method is to consider separately the solution of R u (u, λ) = 0 and the fixed point λ ν = [λ ν + γu ν ] + which is determined by equation (4.4). Let us denote by S the set of all nodes of the finite element mesh belonging to Γ h 3 and p a node of S. After some calculus of variations, it is easy to see that the equation (4.4) leads directly to enforce u ν,p = 0 if the nodes p ∈ S are in status of contact by checking an active set condition that is λ ν,p + γu ν,p ≥ 0 for all p ∈ S. Furthermore, equation (4.4) leads also to the condition λ ν,p = 0 in the case of non contact, see [START_REF] Hintermuller | The primal dual active set strategy as a semismooth Newton Method[END_REF][START_REF] Hintermuller | Semismooth Newton methods for a class of unilaterally constrained variational problems[END_REF] for more details. Thereby, let us consider the active subset A of contact nodes defined by A = {p ∈ S : λ ν,p + γu ν,p ≥ 0}, and the associated inactive subset: I = S \ A. Now, we turn to the description of the iterative active set algorithm of index k.

(i) Choose (u (0) , λ (0) ), set k = 0. (ii) Set A k+1 = {p ∈ S : λ (k) ν,p + γu (k) ν,p ≥ 0}, I k+1 = S \ A k+1 . (iii) Find (u (k+1) , λ (k+1) ) such that u (k+1) = u (k) -[G u (u (k) , λ (k) )] -1 R u (u (k) , λ (k) ) u (k+1) ν,p = 0 for all p ∈ A k+1 , (4.7)
λ (k+1) ν,p = 0 for all p ∈ I k+1 . (4.8) (iv) If (u (k+1) , λ (k+1) ) -(u (k) , λ (k) ) ≤ ǫ, R(u (k+1) , λ (k+1) ) ≤ ǫ and A k+1 = A k then stop, else goto (ii).
We recall that the aim of the active set type strategy is to find the correct subset A of all nodes that are currently in contact with the perfectly rigid obstacle.

Remark 4. Unlike the augmented Lagrangian approach, another traits of the active set type method is the non-use of additional nodes for the determination of the normal contact stress. The computation of λ (k+1) ν,p for all p ∈ A k+1 is obtained a posteriori by using the solution u (k+1) in the equation

A(u (k+1) ) + λ (k+1) ν ν = f
. Nevertheless, we can note a similarity with the augmented Lagrangian method. Indeed, the second equation of (4.2) corresponds exactly to the second equation of (3.10).

Projection iterative method

In this part, a projection iterative algorithm based on a fixed point equation is presented. Note that in (4.1) we used a reformulation of the unilateral contact conditions in the term of a nonlinear complementary problem function. Here, in order to obtain the adequate fixed point equation, this function can be used slightly differently to obtain another formulation of (2.9). Indeed, let c > 0, the Signorini boundary conditions (2.9) are equivalent to the following fixed point problem:

u ν = -[-u ν + cΠ ν ] + on Γ 3 . (4.9)
Remark that the equivalence between (4.9) and (4.1) can be easily proved by similar arguments that those used in Section 4.1. Then, the discrete variational problem P h V can be rewritten in another equivalent discrete algebraic form, as follows

A(u) + λ ν ν = f , (4.10)
u ν,p + [-u ν,p -cλ ν,p ] + = 0 for all p ∈ S. (4.11)
Using (4.11), an implicit scheme of index k on S is used and it can be formulated as follows

u (k+1) ν,p = -[-u (k) ν,p -cλ (k+1) ν,p
] + , for all p ∈ S. (4.12)

Note that an explicit scheme for the numerical treatment of (4.11) could also be used, but with all the drawbacks inherent in this type of method. Now, we turn to the description of the projection iterative algorithm (see [START_REF] Zhang | A projection iterative algorithm boundary element method for the Signorini Problem[END_REF] for a presentation in the case of Laplace equations with a unilateral constraint).

(i) Choose (u (0) , λ (0) ), set k = 0.

(ii) Set

A k+1 = {p ∈ S : u (k) ν,p + cλ (k) ν,p ≥ 0}, I k+1 = S \ A k+1 . (iii) Find (u (k+1) , λ (k+1) ) such that A(u (k+1) ) + λ (k+1) ν ν = f , (4.13)
u (k+1) ν,p
= 0 for all p ∈ A k+1 , (4.14)

u (k+1) ν,p = u (k) ν,p + cλ (k+1) ν,p for all p ∈ I k+1 . (4.15) (iv) If (u (k+1) , λ (k+1) ) -(u (k) , λ (k) ) ≤ ǫ, R(u (k+1) , λ (k+1) ) ≤ ǫ and A k+1 = A k then stop, else goto (ii).
Even if it was not referred to as an active set type method in [START_REF] Zhang | A projection iterative algorithm boundary element method for the Signorini Problem[END_REF], the algorithm fulfills the same objective than the active set strategy presented in Section 4.2, that is to say to find the correct subset A of all nodes in contact with the perfectly rigid obstacle. However, the establishment of the subset I = S \ A changes according to the algorithm computed. In the primal dual active set method, the inactive set is obtained by using homogeneous Neumann boundary condition while, in the projection iterative method, the Robin condition u

(k+1) ν,p = u (k) ν,p + cλ (k+1) ν,p
for all p ∈ I k+1 is used instead. As a consequence, the non linear equation (4.13) can be solved by a linearized Newton method in which the Robin condition (4.15) has to be included beforehand. It leads to a differential of R u (u (k) , λ (k) ) with respect to u (k) different from the case of the primal dual active set method.

Remark 5. We also insist on the fact that for the projection iterative method, it is not necessary to consider additional fictitious nodes in the initial mesh as for the augmented Lagrangian method. This particular trait represents one of the main advantages of the active set type methods.

Convergence analysis of the projection iterative method

Since the convergence of active set type methods remains a challenging and open problem in the framework of hyperelasticity, in this section we investigate the convergence of the projection iterative algorithm in the infinitesimal strain theory. Therefore, under the hypothesis of small displacementgradient theory, the first Piola-Kirchhoff tensor Π can be replaced by the Cauchy stress tensor σ, since u << 1 and ∇u << 1 (cf [8, p. 71]). In this part only, additional notations compatible with the small strain theory are introduced. We denote by S d the space of second order symmetric tensor on R d and we consider the space Q defined by

Q = σ = (σ ij ) : σ ij = σ ji ∈ L 2 (Ω) .
Note that ε represents the linearized strain tensor in the small displacement-gradient theory

ε(v) = (ε ij (v)), ε ij (v) = 1 2 (v i,j + v j,i ).
Then, we can consider an elastic constitutive law given by the following relation σ = Eε(u).

(5.1)

The domain Ω represents the cross section of a three-dimensional deformable body subjected to the action of traction in such a way that a plane deformation hypothesis is assumed. The horizontal component of the displacement field vanishes on Γ 1 . Vertical traction of density f 2 acts on the boundary Γ 2 . No body forces are assumed to act on the body during the process. The body is in contact without friction with an obstacle on the part Γ 3 of its boundary. For the discretization, we have control over a parameter N h Γ3 that describes the number of nodes on Γ h 3 (the number of nodes on Γ h 2 is, therefore, chosen accordingly in order to obtain regular-enough elements). In the configuration N h Γ3 = 256 depicted in Figure 1, 30356 elastic finite elements were used and the total number of degrees of freedom is equal to 30810. We model the material's behavior with an elastic linear constitutive law in which the elasticity tensor E satisfies

(Eε) αβ = Eκ (1 + κ)(1 -2κ) (ε 11 + ε 22 )δ αβ + E 1 + κ ε αβ , 1 ≤ α, β ≤ 2.
Here E is the Young modulus, κ the Poisson ratio of the material and δ αβ denotes the Kronecker symbol.

For the computation below the following data were used:

E = 150N/m 2 , κ = 0.3, f 0 = (0, 0)N/m 2 , f 2 = (0, -2) N/m on Γ 2 , R = 8m, stopping criterion : ǫ = 10 -5 .
In Figure 2, the deformed configuration as well as the normal contact stresses on the boundary Γ 3 are plotted.

Accuracy of the methods comparing to the analytical solution. According to the Hertz theory, it is possible to calculate the analytical solution given by the following normal pressure distribution on Γ 3 :

     if |x 1 | < b, p ν = - 4Rf 2y πb 2 b 2 -x 2 1 (contact), else p ν = 0 (no contact), (6.1
)

with b = 2 2R 2 f 2y (1 -κ 2 )
Eπ . We present now the comparison between the numerical solutions, obtained by the projection iterative method, the primal dual active set method, the augmented Lagrangian method, and the analytical solution (6.1), given by the Hertz theory.

In Figure 3, we investigate the accuracy of the methods on Γ 3 by plotting the normal contact stress σ ν with respect to the abscissa for each method. We note that the analytical solution and the numerical solutions are already almost indistinguishable. Therefore, as expected, the projection iterative method and the primal dual active set method are reliable on this test case. In Figure 4, the difference (∆ Analytic ) between the analytical solution and the normal contact stress σ ν obtained with the 3 numerical methods are also plotted. Note that, while there is a very slight difference between the projection iterative method and the augmented Lagrangian method, no difference is visible at first sight between the primal dual active set method and the augmented Lagrangian method.

Performances of the algorithms. In Table 1, we give the number of iterations for the convergence of the projection iterative method and the augmented Lagrangian method with respect to the arbitrary positive parameter c used in Section 4.2 and with respect to the number of nodes of Γ h 3 , N h Γ3 . In Table 2, we conduct the same study for the primal dual active set method with respect to the arbitrary positive parameter γ used in Section 4.1.

According to Tables 1-2, as expected, the arbitrary positive parameters c and γ do not have a significant effect on the number of iterations for the projection iterative method and the primal dual active set method. However, such statements are not exactly verified for the projection iterative method in the cases where c is less than or equal to 10. This is due to the fact that the condition (4.15) tends to the condition (4.8) when c tends to infinity. Furthermore, the number of iterations to reach the convergence is identical for the primal dual active set method and the augmented Lagrangian method, and a such trait comes from that the second equation of (4.2) corresponds exactly to the second equation of (3.10). Besides, note that the number of iterations for the augmented Lagrangian method and the two active set type methods is not degraded with the increasing of N h Γ3 . Next, in Figure 5-7 we provide the evolution of the error E k = x (k) -x (k-1) 2 with respect to the iteration index k for each method, with c = γ = 10 4 , and for 3 values of N h Γ3 (N h Γ3 = 128, 256, 512). There is almost no difference between the curve obtained with the projection iterative method and the curve obtained with primal dual active set. Also it seems that both methods have an error evolution similar to that of the augmented Lagrangian method.

In Table 3, we provide the CPU time of the simulation for each method with respect to N h Γ3 . The two active set type methods are faster by far, in term of CPU time, than the augmented Lagrangian method. This behavior can be explained by the use of the Lagrange multipliers for the augmented Lagrangian method which worsen the condition number of the linearized systems coming from the Newton method and therefore this leads to an increase of the solution's CPU time of the linearized systems. Such a result must be highlighted by the actual size of the tangent linear system, i.e. the number of degrees of freedom, in both cases with respect to the discretization considered:

• N h Γ3 = 128: 7921 dof for the augmented Lagrangian method compared to 7792 dof for the active set type methods,

• N h Γ3 = 256: 31067 dof for the augmented Lagrangian method compared to 30810 dof for the active set type methods,

• N h Γ3 = 512: 123065 dof for the augmented Lagrangian method compared to 122552 dof for the active set type methods.

Note that, as already implied in the presentation of the algorithms, the tangent linear system arising from the active set type methods is smaller since we do not use the Lagrange multiplier in those cases. In fact, given the small differences between the primal dual active set method and the projection iterative method, we can consider that they are equally fast, both in term of iteration and in term of CPU time.

Evolution of the deformed meshes during the iterative process. We display now in Figures 8-9 the deformed meshes at different steps of convergence of the active set methods and the augmented Lagrangian method. More precisely, the deformed configurations are plotted for the iteration number 1, 3, 5 and 8. It is relevant to compare them since the three methods converge with the same number of iterations for c = γ = 10 4 . Unexpectedly, the deformed meshes during the iterative process of the augmented Lagrangian method and the active set type methods are identical to the point that they are stackable. This behavior can be explained by the fact that the contact conditions (4.1) and (4.9) used for the two active set methods correspond exactly to the second equation of the non linear system (3.10) arising from the augmented Lagrangian method.

Error estimates. At last, we provide and check the numerical convergence order of the active set methods. In order to do that, a sequence of numerical solutions is computed by using triangulations of the body according to the spatial discretization parameter h. Since the solutions are virtually identical, we provide this numerical study for the projection iterative method. The numerical estimated error values uu h V , denoted by E h in Figure 10, are computed for several discretization parameters of h. Here, the boundary Γ 3 of Ω is divided into 1/h equal parts. The numerical solution corresponding to h = 1/512 is taken as the "exact" solution, which is used to compute the errors of the numerical solutions with 6 higher values of h. Therefore, the results obtained for the following values of h : h = 1/64, h = 1/96, h = 1/128, h = 1/192, h = 1/256 and h = 1/384 are plotted. The numerical results are presented in Figure 10 where the dependence of the error estimate uu h V with respect to h is plotted. According to Figure 10, the curve of the numerical error estimates is asymptotically linear with respect to h, which is consistent with the theoretical error estimates established in [START_REF] Belgacem | Hybrid finite element methods for the Signorini problem[END_REF][START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF] for the case of the small deformation theory.

Compression of a hyperelastic ring against a foundation

In this section, we consider a second representative application in order to assess the performances of the active set type methods in the case of large deformation framework. It concerns an academic frictionless contact problem based on the compression of a hyperelastic ring against a rigid foundation. Details on the physical setting of the problem, depicted in Figure 11, are given below:

Ω = (x 1 , x 2 ) ∈ R 2 : R 2 i ≤ x 2 1 + (x 2 -R e ) 2 ≤ R 2 e , Γ 1 = AB, Γ 2 = ∅, Γ 3 = (x 1 , x 2 ) ∈ R 2 : x 2 1 + (x 2 -R e ) 2 = R 2 e , x 2 ≤ √ 2R e .
The domain Ω represents the cross section of a three-dimensional hyperelastic body subjected to the action of displacement in such a way that the plane stress hypothesis is assumed. The foundation is given by (x 1 , x 2 ) ∈ R 2 : x 2 ≤ 0 . On the part Γ 1 , a vertical displacement u 0 is imposed. No forces are assumed to act on the hyperelastic body during the process. The body is in frictionless contact with an obstacle on the part Γ 3 of the boundary. For the discretization, we have control over a parameter N h Γ that describes the number of nodes on Γ (the number of nodes on the thickness of the ring is, therefore, chosen accordingly in order to obtain regular-enough elements). In the configuration N h Γ = 256 depicted in Figure 11, the boundary Γ 3 is divided in 192 equal parts, 6144 hyperelastic finite elements were used for a total number of degrees of freedom equal to 6656. The compressible material response, considered here, is governed by a variant of the Ogden constitutive law (see [START_REF] Ciarlet | Sur les lois de comportement en élasticité non-linéaire compressible[END_REF]) defined by the following energy density

W (C) = c 1 (I 1 -3) + c 2 (I 2 -3) + a(I 3 -1) -(c 1 + 2c 2 + a) ln I 3 .
Here I 1 , I 2 and I 3 represent the three invariants of the tensor C, with C = F T F. For the numerical experiments, the data are:

u 0 = (0, -14)m, f 0 = (0, 0)N/m 2 , f 2 =
(0, 0) N/m , (0, 0) N/m , c 1 = 0.5M P a, c 2 = 0.5 × 10 -2 M P a, a = 0.35M P a,

R i = 9m, R e = 10m.
In Figure 12, the deformed configuration as well as the normal contact stresses on the contact boundary Γ 3 are plotted.

As in the previous Section 6.1, we present the comparison between the numerical solutions obtained by the active set type methods and the augmented Lagrangian method with respect to the numerical accuracy, convergence and performances of the methods.

First in Figure 13, we plot the normal contact stress Π ν with respect to the abscissa on Γ 3 for the three mentioned methods. It is clear that the methods do not present any significant differences, which leads us to say that the active set type methods are accurate even in the large deformation framework on this numerical test compared to the augmented Lagrangian methods. Now, in Tables 4-5, we study the evolution of the number of iterations needed for the convergence of both the projection iterative method and the primal dual active set method compared to the augmented Lagrangian method. In the case of this hyperelastic contact example, the number of iterations is more dependent to the value of c than in the case of the Hertz contact problem of Section 6.1. We remark that, when c ≤ 10, the algorithm even fails to converge (Ncv). However, according to Tables 4 and5, for a large value of c, greater or equal to 10 4 , the number of iterations needed for the convergence of the projection iterative method becomes constant and, is substantially equal to that obtained for the convergence of the primal dual active set method. Furthermore, we remark that the values of γ hardly impact the number of iterations for the primal dual active set method. Moreover, unlike the numerical example studied in Section 6.1, both the primal dual active set method and the projection iterative method take extra iterations to converge, compared to the augmented Lagrangian method.

Next, we provide in Figures 14-16 the evolution of E k = x (k)x (k-1) 2 with respect to the iteration index k for the augmented Lagrangian method and the active set methods for several values of discretization (N h Γ = 128, 256, 512). Although the projection iterative method requires some extra iterations to converge, it can be seen that the convergence of the projection iterative method is smoother than the convergence of the primal dual active set method.

At last, we provide in Table 6 a comparison of the methods in term of CPU time to reach the convergence with respect to the size of the problem. It leads us to say that the two active set type methods have a faster convergence in term of CPU time than the augmented Lagrangian method. This behavior is due to the fact that the two active set methods do not use any Lagrange multipliers. As for the first numerical example (section 6.1), such a result must be highlighted by the actual size of the tangent linear system, i.e. the number of degrees of freedom, in both cases with respect to the discretization considered:

• N h Γ = 128: 1889 dof for the augmented Lagrangian method compared to 1792 dof for the active set type methods,

• N h Γ = 256: 6849 dof for the augmented Lagrangian method compared to 6656 dof for the active set type methods,

• N h Γ = 512: 25985 dof for the augmented Lagrangian method compared to 25600 dof for the active set type methods. This is why the avoidance of the Lagrange multipliers remains a solid argument to rather choose the active set type methods over the augmented Lagrangian method.

Conclusion

This paper provides the analysis of two active set type methods through a classical problem from the contact mechanics literature, that is to say the Signorini problem both in the case of the large and small deformation. First, a variational formulation is derived from the mechanical problem and an existence result based on stored energy function properties is provided in the framework of hyperelasticity. Then, we derived a numerical approximation of the problem and provided a minimization formulation through the augmented Lagrangian formalism. After that, we presented two active set type methods with their algorithms. Next, we investigated the convergence of the projective iterative algorithm in the infinitesimal strain theory; the proof is based on considerations about projection operators and the properties of the elasticity operator. Note that the same properties were used for the uniqueness of the weak solution in the small deformation theory. Then, we presented several numerical simulations in order to compare the behavior of the two active set type methods. We carried them out on two test problems, with the augmented Lagrangian method taken as a reference: one in the small deformation framework with the Hertzian contact and one in the large deformation framework with the compression of a hyperelastic ring against a rigid foundation. While the results showed quite clearly that the active set type methods are effective on the first test case, both in term of convergence rate and in term of CPU time, the second test case must be studied with caution. Indeed, these active set methods took more iterations than the augmented Lagrangian to converge. However, both of them were faster in term of CPU time than the augmented Lagrangian. Nevertheless, as already said, the fact that the active set type methods does not require the use of the Lagrange multipliers cannot be neglected from an implementation point of view. Therefore, for all these reasons, it seems consistent to consider an active set type method over the augmented Lagrangian method.

At this point, we can already highlight some prospects which could be considered in the continuation of this paper. From a theoretical point of view, we showed the convergence of the projection iterative algorithm in the small deformation framework, using only the hypotheses needed for the uniqueness of the weak solution in the same framework. Since we presented an existence result in the large deformation framework based on several different hypotheses, a challenging prospect would be to prove the convergence of this algorithm in this difficult case and this point remains for now an open question, at the best of our knowledge. Besides, we only studied frictionless static contact problem. It could be interesting to consider the challenging case of a frictional contact problem or the dynamic case.

Figure 1 :

 1 Figure 1: Physical setting of the Hertz contact problem.

Figure 2 :Figure 3 : 1 Figure 4 :

 2314 Figure 2: Deformed mesh and normal contact stresses on Γ 3 .

Figure 5 :Figure 6 :Figure 7 :Figure 8 :Figure 9 :

 56789 Figure 5: Evolution of E k with respect to k for N h Γ 3 = 128 (c = γ = 10 4 ).

Figure 11 :

 11 Figure 11: Physical setting of the hyperelastic contact problem.

Figure 12 :

 12 Figure 12: Deformed mesh with normal contact stresses on Γ 3 .

Figure 14 :Figure 15 :Figure 16 :

 141516 Figure 14: Evolution of E k with respect to k for N h Γ = 128 (c = γ = 10 4 ).

Table 1 :

 1 Number of iterations for the convergence of the projection iterative method with respect to c and N h Γ 3 .

	4

Table 2 :

 2 Number of iterations for the convergence of the primal dual active set method with respect to γ and N h Γ 3 .

	E k
	Iteration index k

Table 4 :

 4 Number of iterations for the convergence of the projection iterative method with respect to c and N h Γ . = 10 2 γ = 10 3 γ = 10 4 γ = 10 5 γ = 10 6

	6

Table 5 :

 5 Number of iterations for the convergence of the primal dual active set method with respect to γ and N h Γ .

	E k
	Iteration index k

We also assume that the elasticity tensor E satisfies the following conditions.

∀ ε 1 , ε 2 ∈ S d , a.e. x ∈ Ω. (d) The mapping x → E(x, ε) is measurable on Ω, for any ε ∈ S d . (e) The mapping x → E(x, 0 S d ) belongs to Q.

(5.2) Remark 6. Note that the previous hypotheses were used in [27, p. 126] to prove the uniqueness of the weak solution.

Let (., .) denote the extension of the usual L 2 (Γ) scalar product to H 1/2 (Γ) × H -1/2 (Γ). In order to establish the convergence analysis of the projection iterative method for the case of small strain theory, we have to consider and to prove preliminary results presented in the form of two lemmas.

Lemma 5.1. Denote by (u, σ) the solution of the Signorini problem P, and (u (k+1) , σ (k+1) ) the solution obtained at the iteration k + 1 of the projection iterative method we have

Proof. With the complementary condition (2.3), we have

Also, since

(5.4)

Applying now the Green formula, we have:

By using (2.6), (2.7), (2.8) and (2.10) we obtain

besides, with (5.2)(c) and (5.1)

and, by combining (5.5) and (5.6),

Then, according to (5.4) and (5.7)

From there, we deduce Lemma 5.1

Lemma 5.2. Denote by (u, σ) the solution of the Signorini problem P, and (u (k+1) , σ (k+1) ) the solution obtained at the iteration k + 1 of the projection iterative method; then

(5.8)

Proof. Let us consider K, a nonempty closed convex subset, x, a vector of K, and P K (x), the projector of x on K. Therefore, we have the following proposition proved in [27, p. 13]

Here, we consider

Using the above Lemmas, we can now prove the following theorem Theorem 5.3. The sequence u

, defined in Section 4.2 by the projection iterative method, converges to the unique solution u ν of the Signorini problem as k → ∞.

Proof. First, by using the Lemma 5.2, we have

(5.10)

then, by the Lemma 5.1,

(5.12)

It means that the infinite sum of the terms of the sequence u

ν . Such a couple verifies the fixed point problem (4.9) that is the solution of the Signorini problem; therefore u * = u ν and σ * = σ ν , (5.14) which concludes the proof.

Note that the Theorem 2.1 does not provide a coercivity assumption, or at least a strong monotonicity assumption, for the Piola-Kirchhoff tensor Π, needed in (5.7). Therefore, whether Theorem 5.3 still holds in the case of large deformation framework remains, for now, an open question which deserves to be investigated in the future.

Numerical simulations

The aim of this section is to present several numerical simulations which illustrate, amongst other, the behavior of the numerical solutions of the two active set type methods described in Section 4 compared to the augmented Lagrangian method. Moreover, in order to conduct an objective comparison of methods, the linearized subproblems resulting from each method are solved by the same numerical algorithm, namely a conjugate gradient method with incomplete LU factorization preconditioning. We carried out the numerical simulations based on two static frictionless contact problems: the Hertz half-sphere contact problem and the compression of a hyperelastic ring against a rigid foundation. Note that the different numerical methods have been implemented in a Fortran computer code which is based on a MODULar Finite Element library (MODULEF) developed by INRIA (Institut National de Recherche en Informatique et en Automatique, Rocquencourt, France). For more details, we refer to https://www.rocq.inria.fr/modulef/english.html.

Hertz half-sphere contact problem

Generally, it is quite difficult to find analytical solutions for contact problems, because of their complexity, however there are several exceptions studied in the Hertz theory [START_REF] Hertz | Über die Berührung fester elastischer Körper[END_REF]. For this reason, we describe in what follows the compression of a linear elastic half ball against a foundation that is included in the Hertz theory. The physical setting is depicted in Figure 1. There, the following notations are used:


Mikaël Barboteu

David Danan

Mircea Sofonea

A Hyperelastic Dynamic Frictional Contact Model with Energy-Consistent Properties

published or not. The documents may come


Mircea Sofonea

David Danan

Cong

Cong Zheng

Primal and dual variational formulation of a frictional contact problem

Keywords: AMS Subject Classification : 74M15, 74G25, 74G30, 49J40 elastic material, frictional contact, normal compliance, unilateral constraint, primal variational formulation, dual variational formulation, variational inequality, weak solution

HAL is a

Introduction

Modelling and analysis of contact process is an important topic which, currently, is still under investigation, as illustred in the recent literature. The reason is that contact phenomena abound in industry and everyday life and lead to interesting and non standard mathematical models. Here and below by a mathematical model of contact we mean a system of partial differential equations, associated to boundary value conditions and, eventually, to initial conditions, which describe a specific contact process.

The mathematical models of contact are constructed by considering various constitutive laws (i.e. materials) associated to different contact and frictional conditions. One of the most popular contact condition is the so-called the Signorini condition, introduced in [START_REF] Signorini | Sopra alcune questioni di elastostatica[END_REF] to describe the contact with a rigid foundation. The normal compliance contact condition represents a regularization of the Signorini contact condition and is used to describe the contact with an elastic foundation. It was introduced in [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF], and used in a large number of papers, see [START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF][START_REF] Klarbring | Frictional contact problems with normal compliance[END_REF][START_REF] Klarbring | On friction problems with normal compliance[END_REF][START_REF] Martins | Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws[END_REF][START_REF] Shillor | Models and Variational Analysis of Quasistatic Contact[END_REF] and the references therein. A more general contact condition, called the normal compliance condition with unilateral constraint, was introduced in [START_REF] Jarušek | On the solvability of dynamic elastic-viscoplastic contact problems[END_REF]. It contains as particular cases both the Signorini contact condition and the normal compliance condition, and models the contact with an elastic-rigid foundation. On the other hand, friction is usually modelled with the Coulomb law of dry friction and its versions or regularizations.

The analysis of mathematical models of contact is based on the weak variational formulation of the corresponding models. In most of the cases it is given by a variational or hemivariational inequality in which the unknown is the displacement or the velocity field, as illustrated in [START_REF] Eck | Unilateral Contact Problems: Variational Methods and Existence Theorems[END_REF][START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF][START_REF] Migórski | Nonlinear Inclusions and Hemivariational Inequalities. Models and Analysis of Contact Problems[END_REF][START_REF] Panagiotopoulos | Inequality Problems in Mechanics and Applications[END_REF][START_REF] Shillor | Models and Variational Analysis of Quasistatic Contact[END_REF][START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF], for intance. There, the unique solvability of the models was obtained by using arguments of monotonicity, convexity and fixed point and, in the nonconvex case, by using the properties of the Clarke subdifferential. The progress made in the variational analysis of contact models in the last period was impressive, as it results from the references above. Nevertheless, many open problems still remain to be investigated and resolved. The reason is that, owing to inherent complexity, contact phenomena are modelled by nonlinear problems that are difficult to analyze and, therefore, new and nonstandard mathematical tools are required.

Two mathematical models of contact with normal compliance and unilateral constraint were recently considered in [START_REF] Barboteu | Analysis of a contact problem with normal compliance, finite penetration and nonmonotone slip dependent friction[END_REF][START_REF] Barboteu | Analysis of a contact problem with unilateral constraint and slip-dependent friction[END_REF]. In [START_REF] Barboteu | Analysis of a contact problem with normal compliance, finite penetration and nonmonotone slip dependent friction[END_REF] the material's bahavior was described with a linear elastic constitutive law and friction was modeled with a slip-dependent version of Coulomb's law; the weak solvability of the model was proved by using arguments on pseudomonotone operators followed by a passage to the limit procedure; a convergence result was proved and its numerical validation was also provided; the question of the uniqueness of the solution was left open. In contrast, in [START_REF] Barboteu | Analysis of a contact problem with unilateral constraint and slip-dependent friction[END_REF] the material's bahavior was described with a nonlinear elastic constitutive law and friction was modeled with a new condition, able to describe the transition from the Coulomb'law to Tresca's friction law. In addition, besides the existence of the weak solution to the model, its uniqueness was proved, a finite element approximation of the problem was considered and error estimates have been derived. In both [START_REF] Barboteu | Analysis of a contact problem with normal compliance, finite penetration and nonmonotone slip dependent friction[END_REF] and [START_REF] Barboteu | Analysis of a contact problem with unilateral constraint and slip-dependent friction[END_REF] the analysis of the models was carried out based on variational formulations in terms of displacement, the so-called the primal variational formulation of the corresponding contact problems.

This paper represents a continuation of [START_REF] Barboteu | Analysis of a contact problem with unilateral constraint and slip-dependent friction[END_REF]. Its aim is to study the mathematical model in [START_REF] Barboteu | Analysis of a contact problem with unilateral constraint and slip-dependent friction[END_REF] by using a weak formulation in terms of the stress, the so-called dual variational formulation. Such formulations are used in the literature at least for three reasons. First, in many contact problems the main interest lies to the contact stress, since the behavior of the system and especially the surface integrity and wear depend on it. Second, such formulations can be handled numerically, save errorprone postprocessing and, therefore, lead to a better prediction to the real contact situation. And, finally, in most of the cases, they lead to interesting and non standard mathematical models. References concerning dual variational formulations in the study of various contact problems include [START_REF] Anders | Dual-dual Formulations for Frictional Contact Problems in Mechanics[END_REF][START_REF] Awbi | Dual formulation of a quasistatic viscoelastic contact problem with Tresca's friction law[END_REF][START_REF] Hlaváček | Solution of Variational Inequalities in Mechanics[END_REF][START_REF] Kikuchi | Contact Problems in Elasticity: A Study of Variational Inequalities and Finite Element Methods[END_REF][START_REF] Matei | Dual Formulation of a Viscoplastic Contact Problem with Unilateral Constraint[END_REF][START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF][START_REF] Telega | Topics on unilateral contact problems of elasticity and inelasticity[END_REF].

The rest of the paper is structured as follows. In Section 2 we introduce the contact problem and describe the frictional contact conditions. Next, in Section 3 we introduce the notation and some preliminary material, list the assumptions on the data and state the primal and the dual variational formulation of the problem. In Section 4 we present an equivalence result, Theorem 4.1. It states the equivalence between the primal and dual variational formulation of the contact problem. Its proof is based on arguments of monotonicity and convexity. Then, in Section 5, we provide an existence and uniqueness results, Theorem 5.1. Its proof is based on arguments on quasivariational inequalities.

The model

Everywhere below we denote by S d the space of second order symmetric tensors on R d or, equivalently, the space of symmetric matrices of order d. The inner product and norm on R d and S d are defined by

u • v = u i v i , ∥v∥ = (v • v) 1 2 ∀ u, v ∈ R d , σ • τ = σ ij τ ij , ∥τ ∥ = (τ • τ ) 1 2 ∀ σ, τ ∈ S d .
Let Ω ⊂ R d (d = 2, 3) be a bounded domain. We denote by Γ its boundary, assummed to be Lipschitz continuous. We use the notation x = (x i ) for a typical point in Ω ∪ Γ and we denote by ν = (ν i ) the outward unit normal at Γ. Here and below the indices i and j run between 1 and d and, unless stated otherwise, the summation convention over repeated indices is used. An index that follows a comma represents the partial derivative with respect to the corresponding component of the spatial variable x, e.g. u i,j = ∂u i /∂x j .

The physical setting we consider in this paper is as follows. An elastic body occupies the domain Ω ⊂ R d with a Lipschitz continuous boundary divided into three measurable parts Γ 1 , Γ 2 and Γ 3 such that meas (Γ 1 ) > 0. The body is in equilibrium under the action of body forces of density f 0 and surface tractions of density f 2 which act on Γ 2 . We also assume that the body is fixed on Γ 1 and, on Γ 3 , it is in frictional contact with an obstacle, the so-called foundation. The classical formulation of the contact problem we consider in this paper is as follows.

Problem P. Find a displacement field u : Ω → R d and a stress field σ : Ω → S d such that σ = Fε(u)
in Ω, (2.1)

Div σ + f 0 = 0 in Ω, (2.2
)

u = 0 on Γ 1 , (2.3) σν = f 2 on Γ 2 , ( 2.4
)

u ν ≤ g, σ ν + p(u ν ) ≤ 0, (u ν -g)(σ ν + p(u ν )) = 0 } on Γ 3 , (2.5
)

∥σ τ ∥ ≤ µ p(u ν ), -σ τ = µ p(u ν ) uτ ∥uτ ∥ if u τ ̸ = 0 } on Γ 3 . (2.6)
We now provide a short description of the equations and conditions (2.1)-(2.6) in which, in order to simplify the notation, we do not indicate explicitly the dependence of various functions on the spatial variable x. For more details concerning the mechanical assumptions which lead to the construction of this contact model we send the reader to [START_REF] Barboteu | Analysis of a contact problem with unilateral constraint and slip-dependent friction[END_REF].

First, Equation (2.1) represents the elastic constitutive law of the material in which F is a (possible nonlinear) constitutive function while equation (2.2) is the equation of equilibrium; we use it here since we assume that process is static. Here and below ε and Div represent the deformation and the divergence operators, respectively, i.e.

ε(v) = (ε ij (v)), ε ij (v) = 1 2 (v i,j + v j,i ), Div σ = (σ ij,j ).
Conditions (2.3), (2.4) represent the displacement and traction boundary conditions, respectively. Finally, conditions (2.5) and (2.6) represent the frictional contact condition with normal compliance and unilateral constraint in which u ν , u τ are the normal and tangential components of u on Γ given by u ν = u • ν, u τ = u -u ν ν and σ ν , σ τ represent the normal and the tangential stress on Γ, that is

σ ν = (σν) • ν and σ τ = σν -σ ν ν.
Note that the contact condition (2.5) models the contact with a foundation is made by a rigid body covered by a layer made by of elastic material, say asperities. It shows that the penetration is restricted, since u ν ≤ g where g > 0 represents the thickness of the elastic layer. Also, combined with the friction law (2.6), condition (2.5) shows that when there is penetration, as far as the normal displacement does not reach the bound g, the contact is described with a normal compliance condition associated to the Coulomb's law of dry friction. Here p represents a positive normal compliance function such that p(0) = 0 and µ denotes the coefficient of friction. When the normal displacement reaches the bound g, then the normal stress is larger than the given positive value µp(g) and, moreover, friction follows the Tresca law with the friction bound µp(g). To conclude, this model describes a natural transition from the Coulomb law (which is valid as far as 0 < u ν < g) to the Tresca frcition law (which is valid when u ν = g).

Variational Formulation

In the study of Problem P we need further notation and preliminary material that we introduce in what follows. Everywhere in this paper we use the standard notation for Sobolev and Lebesgue spaces associated to Ω and Γ. In addition, we consider the spaces

Q = { σ = (σ ij ) : σ ij = σ ji ∈ L 2 (Ω) } , H 1 = { u = (u i ) : ε(u) ∈ Q } .
The spaces Q and H 1 are real Hilbert spaces endowed with the canonical inner products given by

(σ, τ ) Q = ∫ Ω σ ij τ ij dx, (u, v) H 1 = (u, v) L 2 (Ω) d + (ε(u), ε(v)) Q .
The associated norms are denoted by ∥ • ∥ Q and ∥ • ∥ H 1 , respectively. Also, recall that H 1 = H 1 (Ω) d algebraically and topologically.

For an element v ∈ H 1 we still write v for the trace of v and we denote by v ν and v τ the normal and tangential components of v on Γ given by

v ν = v • ν, v τ = v -v ν ν.
Also, for a regular stress function σ the following Green's formula holds:

∫ Ω σ • ε(v) dx + ∫ Ω Div σ • v dx = ∫ Γ σν • v da for all v ∈ H 1 . (3.1)
Consider the space

V = { v ∈ H 1 (Ω) d : v = 0 on Γ 1 } .
It is well known that V is a real Hilbert space endowed with the inner product

(u, v) V = (ε(u), ε(v)) Q = ∫ Ω ε(u) • ε(v) dx,
and the associated norm ∥ • ∥ V , respectively. Completeness of the space (V, ∥ • ∥ V ) follows from the assumption meas (Γ 1 ) > 0, which allows the use of Korn's inequality. We also recall that there exists c 0 > 0 which depends on Ω, Γ 1 and Γ 3 such that

∥v∥ L 2 (Γ 3 ) d ≤ c 0 ∥v∥ V for all v ∈ V. (3.2)
Inequality (3.2) represents a consequence of the Sobolev trace theorem.

In the study of the mechanical problem (2.1)-(2.6) we assume that the elasticity operator F satisfies the following conditions.

                         (a) F : Ω × S d → S d . (b) There exists L F > 0 such that ∥F(x, ε 1 ) -F(x, ε 2 )∥ ≤ L F ∥ε 1 -ε 2 ∥ ∀ ε 1 , ε 2 ∈ S d , a.e. x ∈ Ω. (c) The mapping x → F (x, ε) is measurable on Ω, for any ε ∈ S d . (d) There exists m F > 0 such that (F(x, ε 1 ) -F(x, ε 2 )) • (ε 1 -ε 2 ) ≥ m F ∥ε 1 -ε 2 ∥ 2 ∀ ε 1 , ε 2 ∈ S d , a.e. x ∈ Ω. (e) The mapping x → F (x, 0 S d ) belongs to Q. (3.3)
We also assume that the body forces and tractions densities have the regularity

f 0 ∈ L 2 (Ω) d , f 2 ∈ L 2 (Γ 2 ) d . ( 3.4)
The normal compliance function and the coefficient of friction satisfy

                 (a) p : Γ 3 × R → R + . (b) There exists L p > 0 such that |p(x, r 1 ) -p(x, r 2 )| ≤ L p |r 1 -r 2 | ∀ r 1 , r 2 ∈ R, a.e. x ∈ Γ 3 . (c) (p(x, r 1 ) -p(x, r 2 )) (r 1 -r 2 ) ≥ 0 ∀ r 1 , r 2 ∈ R, a.e. x ∈ Γ 3 . (d)
The mapping x → p(x, r) is measurable on Γ 3 , for any r ∈ R. (e) p(x, r) = 0 for all r ≤ 0, a.e. x ∈ Γ 3 .

(3.5) µ ∈ L ∞ (Γ 3 ), µ(x) ≥ 0 a.e. x ∈ Γ 3 . (3.6)
Note that assumptions (3.3)-(3.6) are used in the analysis of the primal variational formulation of the contact model P. Examples of constitutive functions F and normal compliance functions p which satisfy assumptions (3.3) and (3.5), respectively, can be found in [START_REF] Han | Quasistatic Contact Problems in Viscoelasticity and Viscoplasticity[END_REF][START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF], for instance. Now, to derive a dual variational of the contact model P we need an additional assumption on the geometry of the problem. Thus, we assume that there exists θ ∈ V such that θ = ν on Γ 3 .

(3.7)

Examples of domains Ω together with subsets Γ 1 , Γ 3 ⊂ Γ for which assumption (3.7) holds are presented below, both in the two and three dimensional case.

Example 1.

Let

Ω = { x = (x 1 , x 2 ) ∈ R 2 : x 2 1 + x 2 2 < 1 }, Γ 1 = { x = (x 1 , x 2 ) ∈ R 2 : x 2 1 + x 2 2 = 1, x 1 ≤ 0, x 2 ≤ 0 }, Γ 3 = { x = (x 1 , x 2 ) ∈ R 2 : x 2 1 + x 2 2 = 1, x 1 ≥ 0, x 2 ≥ 0 } and consider the function θ : R 2 → R defined by θ(x 1 , x 2 ) =              (x 1 , x 2 ) if x 1 ≥ 0, x 2 ≥ 0, (0, x 2 ) if x 1 ≤ 0, x 2 ≥ 0, (0, 0) if x 1 ≤ 0, x 2 ≤ 0, (x 1 , 0) if x 1 ≥ 0, x 2 ≤ 0.
Then, it is easy to see that θ ∈ H 1 (Ω) 2 , θ = 0 on Γ 1 and θ = ν on Γ 3 , i.e. θ satisfies assumption (3.7).

Example 2. Let a, b, c be positive constants and let

Ω = { x = (x 1 , x 2 , x 3 ) ∈ R 3 : 0 < x 1 < a, 0 < x 2 < b, 0 < x 3 < c }, Γ 1 = { x = (x 1 , x 2 , x 3 ) ∈ R 3 : 0 ≤ x 1 ≤ a, 0 ≤ x 2 ≤ b, x 3 = c }, Γ 3 = { x = (x 1 , x 2 , x 3 ) ∈ R 3 : 0 ≤ x 1 ≤ a, 0 ≤ x 2 ≤ b, x 3 = 0 }.
Consider the function θ : IR 3 → IR defined by

θ(x 1 , x 2 , x 3 ) = ( 0, 0, x 3 -c c
) .

Then, it is easy to see that θ ∈ H 1 (Ω) 2 , θ = 0 on Γ 1 and θ = ν on Γ 3 , i.e. θ satisfies assumption (3.7).

We now turn to the variational formulation of Problem P. To this end we introduce the set of admissible displacements U , the function j : V × V → IR and the element f ∈ V defined by

U = { v ∈ V : v ν ≤ g a.e. on Γ 3 }, (3.8) j(u, v) = ∫ Γ 3 p(u ν )v ν da + ∫ Γ 3 µ p(u ν )∥v τ ∥ da ∀ u ∈ V, v ∈ V, (3.9) (f , v) V = ∫ Ω f 0 • v dx, + ∫ Γ 2 f 2 • v da ∀ v ∈ V. (3.10)
Assume that (u, σ) are sufficiently regular functions which satisfy (2.1)-(2.6) and let v ∈ U . We use Green's formula (3.1) and definition (3.10) to see that

(σ, ε(v) -ε(u)) Q = (f , v -u) V + ∫ Γ 3 σ ν (v ν -u ν ) da (3.11) + ∫ Γ 3 σ τ • (v τ -u τ ) da.
Next, we use the frictional contact conditions (2.5), (2.6) and the definition (3.8) to obtain ∫

Γ 3 σ ν (v ν -u ν ) da ≥ - ∫ Γ 3 p(u ν )(v ν -u ν ) da, (3.12) ∫ Γ 3 σ τ • (v τ -u τ ) da ≥ ∫ Γ 3 µ p(u ν )∥u τ ∥ da (3.13) - ∫ Γ 3 µ p(u ν )∥v τ ∥ da.
Finally, we combine equality (3.11) with inequalities (3.12) and (3.13), then we use the constitutive law (2.1) and the definitions (3.8)-(3.9). As a result we obtain the following variational formulation of the frictional contact problem P.

Problem P V . Find a displacement field u such that u ∈ U, (Fε(u), ε(v) -ε(u)) Q + j(u, v) -j(u, u) ≥ (f , v -u) V ∀ v ∈ U. (3.14)
Note that Problem P V is formulated in terms of the displacement field and, for this reason, we refer to it as the primal variational formulation of Problem P. Once the displacement field is known, the stress field can be easily obtained by using the constitutive law (2.1).

Our aim in what follows is to derive a second variational formulation of the contact problem P, in terms of stress, the so-called dual variational formulation. To this end we note that, since meas (Γ 1 ) > 0, the range of the deformation operator ε : V → Q, denoted ε(V ), is a closed subspace of Q. A proof of this preliminary result can be found in [START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF] page 87. It is a direct consequence of the equality

∥v∥ V = ∥ε(v)∥ Q ∀ v ∈ V. (3.15)
Denote by P : Q → ε(V ) the projection operator on ε(V ) ⊂ Q and note that equality (3.15) shows that ε : V → ε(V ) is an invertible operator. We denote in what follows by ε -1 : ε(V ) → V the inverse of ε. Also, we recall that assumption (3.3) implies that F : Q → Q is a strongly monotone Lipschitz continuous operator and, therefore, using Proposition 1.25 in [START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF] it follows that it is invertible. Moreover, the inverse of F, denoted F -1 , is still a strongly monotone Lipschitz continuous operator. The ingredients above allow us to define the operator Λ :

Q → V defined by Λσ = ε -1 P F -1 σ ∀ σ ∈ Q. (3.16)
The importance of this operator arises in the fact that it helps to inverse the elastic constitutive law, as it results from the following elementary result. Lemma 3.1. Assume (3.3) and let σ ∈ Q, v ∈ V be such that σ = Fε(u). Then u = Λσ.

Proof. Equality σ = Fε(u) shows that ε(u) = F -1 σ and, therefore, since ε(u) ∈ ε(V ) we have ε(u) = P ε(u) = P F -1 σ which implies that u = ε -1 P F -1 σ. We combine this equality with the definition of (3.16) to see that u = Λσ, which concludes the proof.

Besides the operator Λ we consider the element g ∈ V given by

g = gθ ∈ V (3.17)
where, recall, g is a positive constant. In addition, for every η ∈ V we consider the set Σ(η) ⊂ Q defined by

Σ(η) = { τ ∈ Q : (τ , ε(v) -ε( g)) Q + j(η, v -g) ≥ (f, v -g) V ∀ v ∈ U } . (3.18)
Then we consider the following variational problem.

Problem P D V . Find a stress field σ such that σ ∈ Σ(Λσ), (F -1 σ, τ -σ) Q ≥ (ε( g), τ -σ) Q ∀ τ ∈ Σ(Λσ). (3.19)
We refer in what follows to P D V as the dual formulation of problem P. Its link with the primal variational formulation, P V , together with its unique solvability, will be studied in the next two sections. Here we restrict ourselves to note that both problem P V and P D V are expressed in term of elliptic quasivariational inequalities in which the unknown is the dispalcement and the stress field, rerspectively. In the case of Problem P V the set of constraints is given but the problem is governed by the functionl j which depends on the solution. In contrast, Problem P D V there is no such functional but, there, the set of constraints depends of the solution.

We end this section with two additional preliminary results. Lemma 3.2. Assume (3.5), (3.6) and (3.7). Then, the element g defined by (3.17) satisfies the following properties.

                 (a) g ∈ U. (b) 2v -g ∈ U ∀ v ∈ U. (c) λ(v -g) ∈ U ∀ λ ≥ 0, ∀ v ∈ U. (d) j(v, 2v -g) + j(v, g) = 2j(v, v) ∀ v ∈ U. (e) j(u, v) -j(u, g) = j(u, v -g) ∀ u, v ∈ U.
(3.20)

Equivalence results

We now study the link between the variational problems Problem P V and P D V . Our main result in this section is the following equivalence result.

(σ, ε(v) -ε(u)) Q + j(u, v) -j(u, u) ≥ (f , v -u) V ∀ v ∈ U. (4.2)
We take v = 2u -g and v = g in (4.2) and use the equalities (3.20)(d), (e) to obtain

(σ, ε(u) -ε( g)) Q + j(u, u -g) = (f , u -g) V . ( 4.3)
Then, we add (4.2) and ( 4.3) and use (3.20)(e) to see that

(σ, ε(v) -ε( g)) Q + j(u, v -g) ≥ (f , v -g) V ∀ v ∈ U. ( 4
(τ -σ, ε(u) -ε( g)) Q ≥ 0 ∀ τ ∈ Σ(u)
and, since ε(u) = F -1 σ, we deduce that

(F -1 σ, τ -σ) Q ≥ (ε( g), τ -σ) Q ≥ 0 ∀ τ ∈ Σ(u). (4.6)
We now combine (4.1), (4.5) and (4.6) to deduce that σ is a solution to P D V , which concludes the first part of the proof.

2) Assume now that σ satisfies (3.19) and consider an element z ∈ Q such that

(z, ε(v)) Q = 0 ∀ v ∈ V (4.7)
Then, using (3.18) we deduce that τ = σ ± z ∈ Σ(Λσ) and, testing in (3.19) with τ = σ ± z ∈ Σ(Λσ), we deduce that

(F -1 σ -ε( g), z) Q = 0. (4.8)
On the other hand, since ε(V ) is a closed subspace of Q we deduce that

ε(V ) ⊥⊥ = ε(V ) (4.9)
where the superscript ⊥ indicates the orthogonal complement in Q. Using (4.7) and (4.8) we see that F -1 σε( g) ⊥ z for all z ∈ ε(V ) ⊥ and, by (4.9) it follows that Then u ̸ = P U u, i.e.

F -1 σ -ε( g) ∈ ε(V ). Therefore, there exists u ∈ V such that ε( u) = F -1 σ -ε( g). Let u = u + g. Then ε(u) = F -
(P U uu, P U uu) V > 0. (4.12)

Next, using (4.12) and the variational characterization of the projection we have

(P U u -u, v) V ≥ (P U u -u, P U u) V > (P U u -u, u) V ∀ v ∈ U.
These inequalities imply that there exists α ∈ IR such that

(P U u -u, v) V > α > (P U u -u, u) V ∀ v ∈ U and, using (3.20)(c) with λ = 1 it follows that (P U u -u, v -g) V > α > (P U u -u, u) V ∀ v ∈ U. (4.13) Let τ = ε(P U u -u) ∈ Q. Then (4.13) yields ( τ , ε(v) -ε( g)) Q > α > ( τ , ε(u)) Q ∀ v ∈ U (4.14)
and, taking v = g in the previous inequality, we obtain that

α < 0. ( 4
λ( τ , ε( v) -ε( g)) Q > α + ( τ , ε( g)) Q ∀ λ > 0.
Therefore, passing to the limit as λ → ∞, and using (4.16) we deduce that

-∞ ≥ α + ( τ , ε( g)) Q which contradicts α ∈ R. We conclude from above that ( τ , ε(v) -ε( g)) Q ≥ 0 ∀ v ∈ U. (4.17)
Next, equality (4.10) combined with the regularity σ ∈ Σ(Λσ) and the definition (3.18) imply that

(σ, ε(v) -ε( g)) Q + j(u, v -g) ≥ (f, v -g) V ∀ v ∈ U.
Therefore using (4.17) we deduce that τ + σ ∈ Σ(u). This allows to test in (3.19) with τ = τ + σ and, as a result, we find that

(F -1 σ, τ ) Q ≥ (ε( g), τ ) Q .
We now use equality σ = Fε(u) to obtain

(ε(u), τ ) Q ≥ (ε( g), τ ) Q and, since τ = ε(P U u -u) ∈ Q, we find that (P U u -u, u) V ≥ (P U u -u, g) V . ( 4.18)
On the other hand, taking v = g in (4.13), we obtain

(P U u -u, g) V > (P U u -u, u) V (4.19)
The inequalities (4.19) and (4.18) lead to a contradiction. Therefore, we deduce that assumption (4.11) is not valid, and we conclude that

u ∈ U. (4.20)
Next, we denote by

τ (u) = ε(f -ξ(u)) (4.21)
where ξ(u) is the element defined by the Lemma 3.3. Then, using (3.22), we have We combine now inequality (4.30) with (4.20) to see that u is a solution to Problem P V , which completes the proof.

(τ (u), ε(v) -ε(u)) Q ≥ (f , v -u) V + j(u, u) -j(u, v) ∀ v ∈ V. ( 4
constitutive law σ = Fε(u). For this reason we refear to the pair (u, σ) as a weak solution to the frictional contact problem P and we conclude that this problem has a unique weak solution. Note that condition (5.1) represents a smallness inequality for the coefficient of friction. Whether this condition represents an intrinsic feature of this frictional contact problem P or it is only a limitation of our mathematical tools remains an open question. Clearly, this question deserves more investigation in the future.

Theorem 4 . 1 . 1 ) 2 )

 4112 Assume (3.3) -(3.7). Then, the following statements hold: If u is a solution to Problem P V and σ = Fε(u), then σ is a solution to Problem P D V . Conversely, if σ is a solution to Problem P D V , then there exists a unique u ∈ V such that σ = Fε(u) and, moreover, u is a solution to Problem P V . Proof. 1) Let u be a solution to Problem P V and let σ = Fε(u). Then Lemma 3.1 implies that u = Λσ (4.1) and, moreover, (3.14) yields

. 4 )

 4 We now combine (4.4) and the definition (3.18) to deduce that σ ∈ Σ(u).(4.5)Moreover,(3.18) and (4.3) imply that

  .22)We note that (3.20)(a), (b) allow to take v = 2u -g and v = g in (4.22) . Therefore, using equality (3.20)(d) we obtain(τ (u), ε(u) -ε( g)) Q + j(u, ug) = (f , ug) V .(4.23)Then, we combine (4.22) and (4.23) and use(3.20)(e) to see that(τ (u), ε(v) -ε( g)) Q + j(u, vg) ≥ (f , vg) V ∀ v ∈ U. (4.24)This inequality implies that τ (u) ∈ Σ(u) and, using (4.10), it shows that τ (u) ∈ Σ(Λσ). Therefore, testing in(3.19) with τ = τ (u) we have(F -1 σ, τ (u) -σ) Q ≥ (ε( g), τ (u) -σ) Q and, since F -1 σ = ε(u), it follows that (τ (u), ε(u) -ε( g)) Q ≥ (σ, ε(u) -ε( g)) Q . (4.25)Next, we combine (4.23) and (4.25) to deduce that(f , ug) V ≥ (σ, ε(u) -ε( g)) Q + j(u, ug). (4.26)Finally, since σ ∈ Σ(Λσ), (4.10) implies that σ ∈ Σ(u) and, therefore, since u ∈ U , (3.18) shows that(f , ug) V ≤ (σ, ε(u) -ε( g)) Q + j(u, ug). (4.27) Inequalities (4.26) and (4.27) yield(f , ug) V = (σ, ε(u) -ε( g)) Q + j(u, ug), i.e. (σ, ε( g) -ε(u)) Q -j(u, ug) = (f , gu) V . (4.28)On the other hand, since σ ∈ Σ(u) we have(σ, ε(v) -ε( g) Q + j(u, vg) ≥ (f , vg) V ∀ v ∈ U, (4.29)We now add (4.28) and (4.29) to obtain(σ, ε(v) -ε(u) Q + j(u, vg) -j(u, ug) ≥ (f , vu) V ∀ v ∈ U. and, using (3.20)(e) yields (σ, ε(v) -ε(u)) Q + j(u, v) -j(u, u) ≥ (f , vu) V ∀ v ∈ U.Finally, since σ = Fε(u), we find that that (Fε(u), ε(v) -ε(u)) Q + j(u, v) -j(u, u) ≥ (f , vu) V ∀ v ∈ U . (4.30)

  1 σ, i.e. σ = Fε(u). This proves the existence part. The uniqueness of u follows from equality(3.15). In addition, Lemma 3.1 shows that Next, we prove that u is a solution to Problem P V and, to this end, we start proving that u ∈ U . Let P U : V → U be the projector on the nonempty closed convex subset U ⊂ V . Arguing by contradiction, we assume in what follows that

	Λσ = u.	(4.10)
	u / ∈ U.	(4.11)

Acknowledgment. This work was supported by a the Marie Curie International Research Staff Exchange Scheme Fellowship within the 7th European Community Framework Programme under Grant Agreement No. 295118.

Proof. We use (3.7) to see that θ ν = 1 on Γ 3 and, therefore, (3.17) implies that g ν = g on Γ 3 , which shows that g ∈ U . Let v ∈ V and λ ≥ 0. Then 2v ν -g ν ≤ g and λ(v ν -g ν ) ≤ 0 on Γ 3 which imply that 2v -g ∈ U and λ(v -g) ∈ U , respectively. This proves the inclusions (3.20)(a)-(c).

Next, we use (3.17) to see that g τ = 0 on Γ 3 which implies that ∥2v τ -g τ ∥ = 2 ∥v τ ∥ -∥ g τ ∥, ∥v τ -g τ ∥ = ∥v τ ∥ -∥ g τ ∥ on Γ 3 , ∀ v ∈ V. (3.21) We now combine the definition (3.9) with equalities (3.21) to see that (3.20) (d)-(e) hold, which concludes the proof. Lemma 3.3. Assume (3.5) and (3.6) and let j be the functional given by (3.9). Then, for each u ∈ V there exists ξ(u) ∈ V such that

Proof. Let u ∈ V and denote by ξ(u) the function defined by

Then, it is easy to see that for all v ∈ V we have

Therefore, we deduce that ∫

and, using (3.9) yields

Next we apply the Riesz representation theorem to obtain that there exists a unique element ξ(u) ∈ V such that

We combine now (3.23) and (3.24) to obtain (3.22) which concludes the proof.

Note that Lemma 3.2 states the subdifferentiability of the function v → j(u, v) in the point u, for each u ∈ V . Its statement follows from a classical result concerning the subdifferentiability of convex lower semicontinuous functions. Nevertheless, for the convenience of the reader, we decided to present above a constructive proof of this result.

Existence and Uniqueness results

We now proceed with the study of the existence and uniqueness of the solution for problems P V and P D V . Our main results in this section is the following.

Theorem 5.1. Assume (3.3)- (3.6). Then the following statements hold.

1) Problem P V has at least a solution.

2) There exist a constant α 0 which depends only on Ω, Γ 1 , Γ 3 and F such that the solution of Problem P V is unique, if V , obtained at points 2) and 4), respectively, are related by the elastic constitutive law σ = Fε(u).

Proof. 1) The proof of this point follows from arguments of elliptic quasivariational inequalities, as shown in [START_REF] Barboteu | Analysis of a contact problem with unilateral constraint and slip-dependent friction[END_REF]. Nevertheless, for the convenience of the reader, we resume below its main steps.

First, we consider the operator

We use assumption (3.3) to see that A is a strongly monotone Lipschitz continuous operator on the space V . Moreover, we use assumptions (3.5), (3.6) and the property of the trace map to see that the functional j defined by (3.9) satisfies the following conditions:

For all sequences {η n } ⊂ U and {u n } ⊂ U such that

the inequality below holds:

(5.4)

Finally, we recall that U is a closed convex subset of V such that 0 X ∈ U . The existence of a solution of Problem P V follows now from a classical argument of quasivariational inequalities, see for instance [19, p. 51].

2) In order to prove the uniqueness part, let u 1 , u 2 ∈ V be two solution to Problem P V . Then, using (3.14) and (3.3), after a standard computation we find that

Moreover, using the definition (3.9) of the function j, the properties (3.5) of the function p and the trace inequality (3.2) it turns that

and note that, clearly, α 0 depends only on Ω, Γ 1 , Γ 3 and F. Assume that (5.1) holds.

Then it follows from (5.5)-(5.7) that u 1 = u 2 , which concludes the proof.

3) Assume (3.7) holds. Denote by u a solution to Problem P V whose existence is guaranteed by 1) and let σ = Fε(u). Then, using Theorem 4.1 1) it follows that σ is a solution to Problem P D V which ends the proof. 4) Assume now that (3.7) and (5.1) hold and let σ 1 and σ 2 be two solutions to Problem P D V . Then, it follow from Theorem 4.1 1) that there exists two elements

Moreover, u 1 and u 2 are solutions to Problem P V . Using now the smallness assumption (5.1) we deduce by 2) that u 1 = u 2 .

(5.9)

We combine (5.8) and (5.9) to see that σ 1 = σ 2 , which concludes the proof. Assume in what follows that (3.3)-(3.7) and (5.1) hold. Then Theorem 4.1 provides the unique solvability of both the primal variational formulation P and the dual variational formulation P D V . Moreover, it shows that if the displacement field u is the solution of the primal variational formulation P and the stress field σ in the solution of the variational formulation P D V , then u and σ are connected by the elastic