Michel Raous

E Cocu

Pratt

Quasistatic Signorini problem with Coulomb friction and coupling to adhesion

the study of the mathematical aspects of the previous problems with Professor M

Introduction

In this course, we propose to outline some of the classical results obtained in contact mechanics and to present in addition some contributions resulting from the research carried out in our group "Mécanique et Modélisation du Contact" in the "Laboratoire de Mécanique et d'Acoustique" in Marseille. These research topics concern : -the study of the dynamic instabilities associated with friction in small or finite elastic deformations, in collaboration with Professor J.A.C. Martins from the IST in Lisbon, with applications to the modelling of stress waves occurring in some sliding contact or squeal phenomena involving rubber-glass contact ( [START_REF] Barbarin | Instabilité et frottement en élasticité -Application à un problème d'ondes de contrainte[END_REF] [START_REF] Vola | Frottement et instabilités en dynamique : bruits de crissement[END_REF]), -the development of models coupling adhesion, unilateral contact and friction with applications to the modelling of the fiber/matrix interface of composite materials ( [START_REF] Cangémi | Frottement et adhérence : modèle, traitement numérique et application à l'interface fibre/matrice[END_REF] [88]), -the modelling of frictional unilateral contact in finite plastic deformations conducted by P. Chabrand with applications to metal forming ( [START_REF] Pinto | Contact et frottement en grandes déformations plastiques[END_REF] [START_REF] Dubois | Contact, frottement, grandes déformations élastoplastiques[END_REF] [24] [START_REF] Martinet | Etude et mise en oeuvre de formulation ALE -Application au contact lubrifié en mise en forme[END_REF]), -the development of accelerating numerical approaches in contact mechanics, using methods such as multigrid or subdomain decomposition methods (Fast Adaptive Composite grids) or Arbitrary Lagrangian Eulerian formulations ( [START_REF] Grégo | Méthodes multiniveaux pour des problèmes de contact unilatéral avec frottement[END_REF] [START_REF] Martinet | Etude et mise en oeuvre de formulation ALE -Application au contact lubrifié en mise en forme[END_REF]), [START_REF] Aitken | On the iterative solution of a system of linear equations[END_REF] Constitutive models for frictional contact

Introduction

In this section, a few comments on the real physical complexity of a contact surface are first presented before giving an overview of the various basic formulations of frictional unilateral contact. In this short presentation of some of the tribological aspects of the contact, I would like to make it clear that on the one hand the models are only rough approximations of reality, and that on other hand, the regularization parameters introduced into some models should be based on physical considerations. Various formulations for unilateral contact and different friction laws are given. Among them, the Signorini problem (strictly unilateral conditions) and the Coulomb friction law which will be treated in the following chapters. Most of the analytical models in the past were based on microscopic considerations. References to some of them will be given. Nowadays, some finite element models are being developed with a view to correlate friction with the plastic deformation of the asperities.

Some tribological aspects 1.2.1 Surface asperities

To let the reader have a look at the real appearance of a contact surface, an experimental profile of a metal sheet is given in Fig. 1. These surfaces are classically characterized by statistical parameters (Ra, Rp, etc).

Figure 1: Roughness of a metal surface (metal forming) [START_REF] Felder | La tribologie de l'emboutissage[END_REF]

Real contact surface

A real physical contact is not a full contact but a contact between asperities. During loading, the asperities will be deformed and the real contact area changes as schematized in Fig. 2, where the obstacle is assumed to be flat which is not the case in reality (see [START_REF] Kikuchi | Contact problems in Elasticity: a study of variational inequalities and finite element methods[END_REF]).

Complexity of the surface material

Regarding the physico-chemical aspects, a contact surface is very complex and some of the components are not well known. A scheme of a metal surface is given in Fig. 3 (see [START_REF] Rabinowicz | Friction and wear of materials, 2nd editionJ[END_REF]).

Unilateral contact

First, various formulations for unilateral contact are given. They take into account the non penetration (or the weak penetration) of the solid into the obstacle, which is assumed to be rigid and fixed for the sake of simplicity.

The Signorini problem

Let Ω be an open bounded set of R d (with d = 2, 3), which is the interior of an elastic body, with a sufficiently smooth boundary Γ = Γ F ∪ Γ D ∪ Γ C , in contact with a rigid body. The boundary Γ D (respectively Γ F ) is the part of Γ on which the displacements (respectively the forces) are prescribed. Γ C is the part (of the boundary) initially in contact with the rigid body (see Fig. 4). Let F , the unknown contact force density, and u, the displacement defined on Γ C , be decomposed into normal and tangential components, where n is the outward unit vector to Γ.

u = u N n + u T (1) F = F N n + F T (2)
The strictly unilateral contact problem, known as the Signorini problem, can then be written:

u N ≤ 0 F N ≤ 0 u N F N = 0              (3)
This satisfactorily describes the unilateral contact with the two possible conditions : -no contact, then u N < 0 and F N = 0, -contact, then u N = 0 and F N ≤ 0. This is a non smooth problem, and the complementarity problem [START_REF] Antes | The boundary integral approach to static and dynamic contact problems -Equality and inequality methods[END_REF] implies that F N is a multivalued application of u N , see the graph Fig. 5. In the previous problem, the size of the contact area cannot exceed the size of Γ C . In the case where an extension of the contact area is possible during the loading process, the condition u N ≤ 0 will be replaced by the condition u N ≤ d, where d is the initial gap between the solid and the obstacle. In small deformations, this initial gap can be defined using the normal related either to the solid or to the obstacle.

b -Contact between two deformable solids

In small deformations, the contact between two deformable bodies can be easily dealt with by using a point to point procedure, making a simple change of variables ( [START_REF] Raous | Numerical methods in mechanics of contact involving friction[END_REF]).

When dealing with finite deformations or large displacements, the problem is much more complex.

c -Orientation of the surfaces Let us mention that, from the computational point of view, instead of making the change of variables [START_REF] Aitken | On the iterative solution of a system of linear equations[END_REF], a local change of referential will be made once and for all in the finite element matrix for each contact node. [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF] This is another class of models which can be said to be a regularization of the Signorini conditions. The "unilateral" contact is then accounted for by a nonlinear behaviour law relating the normal force F N to the normal displacement u N as follows (where (u N ) + is zero when u N < 0 (no contact), and (u N ) + = u N when u N ≥ 0) :

Compliance model

-F N = C n (u N ) mn + ( 4
)
This problem is much more regular, because we now have F N as a function of u N , the graph is given in Fig. 6. This model is simpler in terms of both the computational and the mathematical. Nevertheless, some penetration of the solid into the obstacle occurs. This penetration depends on the choise of the constitutive parameters C n and m n . When they are chosen for the sake of computational comfort, the penetration can be large. Whereas when they are obtained from mechanical estimates (see [START_REF] Raous | Numerical simulation of the behaviour of surface asperities for metal forming[END_REF]) where these parameters are deduced by calculating the flattening of asperities with a finite plasticity model, the penetration is slight but computation may become very uncomfortable because F N is then a very stiff function of u N . This compliance law, which as many similarities with penalty methods, has the advantage of having a mechanical meaning. This model and these numerical solvers are very convenient for dealing with problems where the solution is only slightly sensitive to the penetration. When subtle phenomena such as instabilities or the adhesion coupling are studied, the amount of penetration one accepts can strongly affect the results ; that is different values of the penalty or constitutive parameters induce slightly different contact solutions which may have large consequences in the mentioned phenomena. For these reasons, we will deal here mainly with the Signorini formulation, i.e. with the strictly unilateral conditions.

"Third body" model

Some models treat the contact by assuming that a third body is connected to the solid and the obstacle (or to the two deformable bodies) by boundary conditions. This third body approach is very rarely used now. In this case, the contact is not modelled in terms of inequalities, but in terms of the behaviour of this third body. One advantage of this approach is that the numerical treatment can be attempted using a classical finite element code with various behaviour laws. Another attractive point is that it takes into account some specific physical aspects of the phenomena involved in a contact (fluid or solid lubricants, wear fragments, etc). Some difficulties arise however : the microbehaviour of these third bodies (except for fluids) are very complex and little known and the numerical treatment has to deal with evanescent thicknesses which leads to ill-posed problems.

Friction

Friction in the analysis and the numerical methods presented in the following chapters is modelled by the Coulomb law or a variant with a given sliding threshold, Tresca friction. In this section, after giving these two models, a compliance law developed by Martins-Oden [START_REF] Martins | Existence and uniqueness results for dynamic contact problem with nonlinear normal and friction interface laws[END_REF] is recalled and some variants of the Coulomb law are given. Non local friction laws have been introduced by Oden-Pires [START_REF] Oden | Non local and non linear friction laws and variational principles for contact problem in elasticity[END_REF] ; the normal force is then defined at each point through a convolution on a small area surronding the point. These formulations are helpful for establishing mathematical results (see section 2). Another interesting formulation in terms of non associate plasticity has been given by Curnier [START_REF] Curnier | A theory of friction[END_REF].

Tresca friction

Let us first introduce an elementary friction law where the function giving the sliding threshold g is given. The term "Tresca friction" is often used because of the similarity with the corresponding plasticity law. It is written as follows, where uT denotes the time derivative of u T and F T denotes the modulus of vector F T (the 3D formulation is given here ; in 2D it becomes an absolute value).

F T ≤ g with if F T < g then uT = 0 if F T = g then ∃λ ≥ 0 such that uT = -λF T              (5)
This law will be realistic only in those cases where the normal pressure is known and constant or weakly changing (polymer injection, for example). It cannot be directly written with unilateral contact because in this case, the normal force is not known a priori (although it has to be zero under non-contact conditions).

On the other hand, this law is interesting however, because of its simplicity as regards both the mathematical and the numerical aspects. As we will see in chapter 2, the convexity of this law makes it possible to associate it with a minimum principle, and will give a result of existence and uniqueness of the solution. As the problem can also be written as a minimization one, the class of numerical methods which can be used will be much larger than in the case of the Coulomb law. In some of the algorithms presented in chapter 3, the solution of this problem will be used as an intermediate step in treating the Coulomb friction.

Coulomb friction

This law was introduced by Amontons in 1699 [START_REF] Amontons | Sur l'origine de la résistance dans les machines[END_REF] and developed by Coulomb in 1785 [START_REF] Coulomb | Théorie de machines simples[END_REF]. The sliding limit is proportional to the unknown normal contact force and the velocity of the tangential displacement is colinear to the tangential force.

F T ≤ µ | F N | with if F T < µ | F N | then uT = 0 if F T = µ | F N | then ∃λ ≥ 0 such that uT = -λF T              (6)
This is a non associate law because the sliding direction is not normal to the Coulomb cone but colinear to the friction force. As given in [START_REF] Auatt | An interior point optimization algorithm for contact problem in linear elasticity[END_REF], a friction law has to be written in terms of velocities because the solution depends not only on the current loading condition but also on the previous evolution of the system. Nevertheless a "static" friction problem will be also introduced in the next section. It will be considered as either an intermediate problem in solving problem (6) (using an incremental formulation), or a convenient formulation in the case of monotonic and proportional loading. The Coulomb law is a multivalued application as can be seen from the graph given in R 2 on Fig. 7.

Other equivalent forms of the Coulomb law

Here we briefly give some other formulations of the Coulomb law (details can be found in [Jean,1993], [START_REF] Moreau | Unilateral contact and dry friction in finite freedom dynamics[END_REF]).

-maximum dissipation principle

F T ∈ C ∀S T ∈ C (S T -F T ) uT ≥ 0 (7)
where, in 2D :

C = [-µF N , +µF N ] (8)
-subdifferential formulation

-uT ∈ ∂I C T (F N ) (F T ) (9)
where ∂I C T (F N ) (F T ) is the subdifferential of the indicator function of C T (F N ) with:

C T (F N ) = {P / |P| ≤ -µF N } (10)
-dual subdifferential formulation

F T ∈ ∂φ F N (-uT ) (11)
with :

φ F N (-uT ) = µF N uT (12)
1.4.3 Compliance model [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF] As in the case of unilateral contact, a compliance law has been introduced (see [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF]) for the friction law.

F T ≤ C T (u N ) m T + with if F T < C T (u N ) m T + then uT = 0 if F T = C T (u N ) m T + then ∃λ ≥ 0 such that uT = -λF T              (13)
The case where C T = C N and m T = m N corresponds to the classical Coulomb law associated with the compliance model for the unilateral contact.

Regularization of the Coulomb law

To obtain a smoother formulation, some classical regularizations of (6) may be used, (F T becomming a function of uT instead of a multivalued application) :

F T = -µ ϕ ε ( uT ) | F N | (14)
with for 2D problems :

• Square root : ϕ 1 ε ( uT ) = uT u2 T + ε 2 , • Hyperbolic tangent : ϕ 2 ε ( uT ) = tanh( uT ε ) ,
• Piecewise polynomial :

ϕ 3 ε ( uT ) =              -1 , if uT < -ε uT 2ε , if -ε ≤ uT ≤ ε +1 , if uT > ε
Graphs of the various regularizations are given on Fig. 8. To take account of high pressure levels (as in metal forming [START_REF] Oudin | Approches expérimentales et numériques des conditions de contact et de frottement[END_REF]) or changes in the nominal contact surface, the following form with various choices for the function f has been proposed (graphs are given on Fig. 9) :

F T ≤ f with if F T < f then uT = 0 if F T = f then ∃λ ≥ 0 such that uT = -λF T              (15) i -Coulomb-Orowan law f = M in(µ | F N |, k) (16)
where k is usually the elastic limit of the material.

ii -Shaw law

f = αk with α = Γ R /Γ C ( 17
)
where Γ R is the nominal (real) contact surface involved in the flattening of the asperities (depending on the normal compression).

iii -other forms Other forms have been introduced by Cristescu, Wanheim, Bay [START_REF] Bay | Real area of contact and friction stress at high pressure sliding contact[END_REF] The Coulomb law is sometimes used with a variable coefficient. As observed experimentally the force necessary to initiate the movement is often larger than the friction force occurring during steady sliding. More generally, a continuous evolution of the friction coefficient relatively to the velocity is sometimes adopted. A classical law in metal forming problems (in the presence of lubricants) is the Stribeck law given in Figure 10: Stribeck law Fig. 10 where p is the pressure and β the lubricant viscosity. In spite of the attractiveness of this kind of model, we have to be aware of the loss of uniqueness (even for small friction coefficients) and the lack of existence results. This can lead to serious problems as far as the reliability of the numerical results is concerned.

1.5 Some models based on microscopic considerations

Analytical modelling

Apart from the numerical analysis presented in this course, a lot of studies have been carried out with analytical models based on micro-analysis conducted at the asperity scale : see Archard [START_REF] Archard | Elastic deformation and the laws of friction[END_REF] [5], Bay-Wanheim-Avitzur [START_REF] Bay | Models for friction in metal forming[END_REF] [11], Challen-McLean-Oxley [START_REF] Challen | An explanation of the different regimes of friction and wear -Effect of lateral tension and compression on plan srain flattening processes of surface asperities lying over a plastically deformable bulk[END_REF] [23], Greenwood-Tabor-Williamson [START_REF] Greenwood | Deformation properties of friction junctions[END_REF] [53] [START_REF] Tabor | Friction-The present state of our understanding[END_REF], Wilson-Sheu [START_REF] Wilson | Real area of contact and boundary friction in metal forming[END_REF], etc. Numerous references can be found in [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF]].

In some of these studies, the asperities are assumed to have elastic behaviours, as in the case of the Archard model summarized in Fig. 11. The asperities are assumed to be spherical with a degressive radius and Hertz theory is used for the contact. Still using micro-analysis approach, Greenwood and Williamson have proposed a model based on the hypothesis that the asperity heights have a Gaussian distribution.

In other models, the asperities are taken to have plastic behaviour [START_REF] Challen | Plastic deformation of a metal surface in sliding contact with a hard wedge: its relation to friction and wear[END_REF]] [ [START_REF] Challen | An explanation of the different regimes of friction and wear -Effect of lateral tension and compression on plan srain flattening processes of surface asperities lying over a plastically deformable bulk[END_REF]). The plastic wave occurring when an asperity indents and ploughs into a deformable obstacle has been described in terms of a sliding line theory by Wanheim, Bay, Rault, Entringer, etc. More recently, some numerical approaches have been used to account for the finite plastic deformation of asperities, thus introducing some micro-mechanical information into the macroscopic unilateral contact and friction models [START_REF] Makinouchi | A finite element analysis of flattening of surface asperities by perfectly lubricated rigid dies in metal working processes[END_REF], [START_REF] Raous | Numerical simulation of the behaviour of surface asperities for metal forming[END_REF]], [START_REF] Ike | Plastic deformation of surface asperities associated with bulk deformation of metal workpiece in contact with rigid tool[END_REF], [Chabrand et al, 1996]). This makes it possible to treat any shape of asperity. In [START_REF] Raous | Numerical simulation of the behaviour of surface asperities for metal forming[END_REF], by calculating the flattening of an asperity, it was possible to estimate the normal compliance parameters C N and m N of the law 6. The identification is conducted on the computed evolution of the normal contact force relative to the decrease in the height of the asperity given in Fig. 13. In a study on the shear off of two asperities (see Fig. 13 and14) in finite plasticity, the evolution of the micro tangential force during the sliding was estimated (see Chabrand et al, 1996]).

Homogenization

To correlate the micro-behaviour of the asperities with macroscopic unilateral contact and friction has been an attractive topic both to tribologists and contact mechanics researchers. This has been done as far as material behaviour is concerned, using the very powerful theory of homogenization. In contact mechanics, the problem is more complex because the laws are not smooth. Some mathematical results have been established in the case of elastic asperities when the height goes to zero by Sanchez-Palencia and Suquet [109], of plastic asperities by Bouchitté, Lidouh and Suquet [START_REF] Bouchitté | Homogénéisation de frontière pour la modélisation du contact entre un corps déformable non linèaire et un corps rigide[END_REF], and with a thin dissipative layer by Licht [START_REF] Licht | Comportement asymptotique d'une bande dissipative mince de faible rigidité[END_REF]. With plastic asperities, when both the amplitude and the period are approching zero an associated law which is similar to the Tresca one is obtained.

Variational formulations for static and quasistatic problems

This section is devoted to the variational formulation of the Signorini problem (strict unilateral conditions) with Coulomb friction and the main mathematical results will be briefly recalled. These formulations are written in terms of implicit variational inequalities. The academic "static" friction problem will first be analyzed, and recent formulations for quasi-static problems will then be presented. Dynamic formulations for discrete problems can be found in [START_REF] Jean | Dynamics in the presence of unilateral contact and dry friction: a numerical approach[END_REF] [85] [90] [START_REF] Oden | Models and computational methods for dynamic friction phenomena[END_REF] [115] [START_REF] Wronski | Some somputational aspects of structural dynamics problems with frictional contact[END_REF].

Even though chapter 3 and 4 will focus on 2D problems, we give in this chapter 3D formulations. The presentation remains here very formal in order to give a short survey of the mathematical aspects fo contact formulations. We do not give all the details on the functional framework since the proofs are not given : they can be found in the papers given in the references.

The Signorini problem with Coulomb friction

As previously mentioned in section 1.3.1 and in figure 4, the part of the boundary initially in contact with the obstacle is Γ C . Boundary conditions are given on Γ D . A volumic (respectively surface) force, the density of which is Φ 1 (respectively Φ 2 ), is applied in Ω (respectively Γ F ). The behaviour is assumed to be linear elastic and the deformations to be small. The problem can be written as follows.

Problem (P c P c P c ) Let Φ 1 , Φ 2 be the given forces, find the displacement field u, the strain ε, the stress σ and the contact force F such that :

ε = grad s u σ = Kε divσ = -Φ 1              on Ω u = 0 on Γ D , σ.n = Φ 2 on Γ F σ.n = F u N ≤ 0, F N ≤ 0, u N F N = 0 F T ≤ µ | F N | with if F T < µ | F N | then uT = 0 if F T = µ | F N | then ∃λ ≥ 0 such that uT = -λF T                            on Γ C 2.2
The "static" Coulomb problem

Formulation

Let us first consider a "static" Coulomb problem where, in Problem (P c P c P c ), the friction is written on the tangential displacement instead of the velocity. The relations (6) are replaced by :

F T ≤ µ | F N | with if F T < µ | F N | then u T = 0 if F T = µ | F N | then ∃λ ≥ 0 such that u T = -λF T              (18)
This law is not realistic when dealing with general cases, because the friction phenomena clearly depend on the previous history of the loading and have to be formulated in terms of velocities. Nevertheless, this static formulation will be acceptable in some specific situations (for example, if the loadings are monotonic and proportional). The main interest of this formulation is to constitute an intermediate problem when an incremental formulation of the Coulomb problem will be written. On the other hand, it should be noted that most of the mathematical results available so far have been established using this static formulation.

Variational formulation : an implicit variational inequality

As announced before, the details of the rigorous mathematical analysis are not given, but they can be found in [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF] for example and in other references given in the sections below. The functional framework is that of the Sobolev spaces, which ensures that the various fields are sufficiently regular.

Let K be the convex of the admissible displacements :

K = {v ∈ U/ v N ≤ 0 on Γ C } (19)
with :

U = v ∈ [H 1 (Ω)] 3 / v = 0 on Γ D .
The variational form of the Signorini problem with static friction is then written as follows.

Problem (P v P v P v ) : Let Φ 1 ∈ [L 2 (Ω)] 3 , Φ 2 ∈ [L 2 (Γ F )] 3 be given, find u ∈ K such that : a(u, v -u) + j 1 (u, v) -j 1 (u, u) ≥ L(v -u) ∀v ∈ K (20)
with :

a(u, v) = Ω σ(u)ε(u) dx (21)
=

Ω K ijkl ε ij (u)ε kl (v) dx ∀u, v ∈ U (22)
L(v) = Ω Φ 1 v dx + Γ F Φ 2 v ds ∀v ∈ U (23) j 1 (v, w) = Γ C µ|F N (v)| w T ds in the sense of the traces on Γ C (24)
The bilinear form a(u, v) is associated with the elasticity mapping, the linear form L(v) with the loadings and the non differentiable functional j 1 (v, w) with the friction. This is an implicit variational inequality. A dual formulation, set on the stresses, would lead to a quasi-variational inequality because of the dependence of the convex K ′ on the solutions in that case.

Indications about the equivalence in a weak sense with the initial problem (details can be found in [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF])

(i) Implication
The following generalized Green formula is applied as usual to the integration on Ω of the product of the volumic equilibrium equation by the test function (v -u).

∀ v ∈ W (Γ) < σ.n, v > = Ω σε dx + Ω (div σ)v dx (25)
It is easy to eliminate the terms on Γ D by using the boundary conditions and to introduce the form L(v -u) by using the surfacic equilibrium equation. Eliminating the unknown normal contact force leads to the variational inequality by using the Signorini conditions. To deal with the friction term on Γ C , the friction law is written in the following equivalent form (see [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF]) :

Γ C F T (v T -w T ) ds + Γ C µ|F N (v)| ( v T -w T ) ds ≥ 0 ∀v ∈ U (26)
(ii) Reverse implication This is established by making various convenient choices of test functions vand using again [START_REF] Cocu | Formulation and approximation of quasistatic frictional contact[END_REF].

Fixed point and sequence of Tresca problems

Because of the non convexity of the Coulomb problem, which is due to the dependence of the sliding limit on the normal contact force, no minimum principle can be associated with the previous problem. A fixed point method on the sliding limit is now introduced in order to treat the problem as a sequence of Tresca problems, the convexity of which is ensured. A condition on the friction coefficient µ to ensure that the following application S(.) is a contraction, is given in [START_REF] Licht | Remarks on a numerical method for unilateral contact including friction[END_REF]. In [START_REF] Licht | Remarks on a numerical method for unilateral contact including friction[END_REF], the Coulomb friction problem is shown to be equivalent to problem (P f P f P f ).

Problem (P f P f P f ) : Find the fixed point of the application S :

S(g) = -µF N (u g ) (27)
with u g solution of the following problem (P T P T P T ) :

Problem (P T P T P T ) : For a given g, find u g ∈ K such that :

a(u g , v -u g ) + j g (v) -j g (u g ) ≥ L(v -u g ) ∀v ∈ K (28)
with :

j g (v) = Γ C g v T ds (29)
Problem (P T P T P T ) is the variational formulation of the Tresca problem 18. It is a classical variational inequality, where the implicit character disappears because the sliding limit is given in this case. Because of the convexity of this problem, it can be shown that Problem (P T P T P T ) is equivalent to the following minimization problem (P m P m P m ) :

Problem (P m P m P m ) : For a given g, find u g ∈ K such that

J(u g ) ≤ J(v) ∀ v ∈ K (30)
with:

J(v) = 1 2 a(v, v) + j g (v) -L(v) (31)
In the next chapter, a numerical algorithm will be deduced from this formulation. At each step, the normal contact force is computed by solving the Tresca problem with a previous value of the function g, and a new value of this sliding threshold g is obtained for the next step from [START_REF] Demkowicz | On some existence and uniqueness results in contact problem with non local friction[END_REF].

2.3

The quasistatic problem The problem has been defined in section 2.1. With evolutive problems, the main difficulty that arises in establishing the variational formulation of this quasistatic problem, is due to the simultaneous presence of the displacement field, which appears in the unilateral conditions, and the velocity field, which appears in the friction law. For this reason, the following variational form combines two inequalities : a variational inequality, which represents the unilateral conditions (with the displacement field used as test function) and an implicit variational inequality, which describes the friction conditions (with the velocity field used as test function). Details about this formulation are given in Cocu-Pratt-Raous [START_REF] Cocu | Existence d'une solution du problème quasi statique de contact unilatéral avec frottement non local[END_REF] [26] and Telega [START_REF] Telega | Quasi-static Signorini's contact problem with friction and duality[END_REF].

Problem (P time P time P time ) : Let Φ 1 ∈ W 1,2 (0, T ; [L 2 (Ω)] 2 ), Φ 2 ∈ W 1,2 (0, T ; [L 2 (Γ F )]
2 ) be given, find u ∈ W 1,2 (0, T ; U ) such that u(0) = u 0 given, and, for almost all t ∈ [0, T ], u(t) ∈ K such that :

                   a(u(t), v -u(t)) + j 1 (u(t), v) -j 1 (u(t), u(t)) ≥ L(v -u(t), t) + < F N (u(t)), v N -uN (t)) > ∀v ∈ U < F N (u(t)), z N -u N (t)) > ≥ 0 ∀z ∈ K (32)
where u 0 belongs to K and satisfies the following compatibility conditions :

a(u 0 , w -u 0 ) + j 1 (u 0 , w -u 0 ) ≥ L(w -u 0 , 0) ∀w ∈ K (33)
The first inequality is obtained from the integration on Ω of the product of the equilibrium equation by the function test (w -u), where w is homogenous to a velocity field. This time, the normal contact force is kept as an unknown variable, and the following velocity form of ( 26) is used to obtain the implicit variational inequality.

Γ C F T (v T -uT ) ds + Γ C µ|F N (v)| ( v T -uT ) ds ≥ 0 ∀v ∈ U (34)
A weak formulation of the unilateral conditions (3), where the test function z is homogenous to a displacement, gives the second variational inequality .

If the solution u of Problem (P time P time P time ) is sufficiently regular, u will also be a solution of Problem (P c P c P c ) because : -the equilibrium equation can easily be obtained by choosing [START_REF] Delamare | Tribologie et mise en forme des matériaux-Actes des Journées Scientifiques de la Société de Tribologie de France-9/10 avril[END_REF] then holds, which yields the friction relations (6), -the choice of u ∈ K ensures the unilateral condition u N ≤ 0 and by choosing z N = 0 and z N = 2 u N in the second variational inequality we obtain :

v = u + φ and v = u -φ on Ω with φ ∈ (C ∞ 0 (Ω)) 2 and

the surface equilibrium equation is obtained by choosing

v = u + ξ and v = u -ξ on Γ C with ξ ∈ U and v = u on Γ D , -inequality
Γ C F N u N = 0 ( 35
)
It follows that Γ C F N z N ≥ 0 ∀z ∈ K which implies that F N ≤ 0 and finally, using equation [START_REF] Piero | Unilateral Problems in Structural Analysis[END_REF], that the complementarity relation u N F N = 0 holds.

Incremental formulation

By using a time discretization t i+1 = t i + ∆t, the following incremental formulation is obtained.

Problem (P inc1 P inc1 P inc1 ) : Find u i+1 ∈ K such that :                    a(u i+1 , v -u i+1 -u i ∆t ) + j 1 (u i+1 , v) -j 1 (u i+1 , u i+1 -u i ∆t ) ≥ L i+1 (v -u i+1 -u i ∆t ) + < F N (u i+1 ), v N -u N i+1 -u N i ∆t ) > ∀v ∈ U < F N (u i+1 ), z N -u N i+1 > ≥ 0 ∀z ∈ K (36)
where u i approximates u(t i ) and

L i (v) denotes L(v, t i ).
Equivalence with the following problem is shown in [START_REF] Cocu | Formulation and approximation of quasistatic frictional contact[END_REF] :

Problem (P inc2 P inc2 P inc2 ) : Find u i+1 ∈ K such that : a(u i+1 , w -u i+1 ) + j 1 (u i+1 , w -u i ) -j 1 (u i+1 , u i+1 -u i ) ≥ L i+1 (w -u i+1 ) ∀w ∈ K (37)
This is an important point because it establishes that the incremental formulation of the quasistatic problem can be written in the form of unique implicit variational inequality, which was not the case for the continuous problem (P time P time P time ). This implicit variational inequality is very similar to the one [START_REF] Chabrand | Modelling drawbeads in sheet metal forming[END_REF] resulting from the static problem, but the test function of the friction term is now (w -u i ) where u i is the solution obtained for the previous time step t i . The previous history is included in this term.

As stated in the case of the static problem, a fixed point method is introduced on the sliding limit and Problem (P inc2 P inc2 P inc2 ) is shown to be equivalent to the following one.

Problem (P inc3 P inc3 P inc3 ) : At each time step t i+1 , find the fixed point of the application S:

S(g) = -µF N (u g i+1 ) (38)
where u g i+1 is the solution of the following problem (P inc4 P inc4 P inc4 ) :

Problem (P inc4 P inc4 P inc4 ) : For a given g, find u g i+1 ∈ K such that : a(u g i+1 , v -u g i+1 ) + j(v -u i ) -j(u g i+1 -u i ) ≥ L(v -u g i+1 ) ∀v ∈ K (39)
with :

j(v -u i ) = Γ C g v T -u i ds (40)
As previously, Problem (P inc4 P inc4 P inc4 ) is shown to be equivalent to the following minimization problems under constraint, either (P opt1 P opt1 P opt1 ) (set on the displacements) or problem (P opt2 P opt2 P opt2 ).

Problem (P opt1 P opt1 P opt1 ) : For a given g, find u g i+1 ∈ K such that

J(u g i+1 ) ≤ J(v) ∀ v ∈ K (41)
with:

J(v) = 1 2 a(v, v) + j(v -u g i+1 ) -L i+1 (v) (42)
This is a displacement formulation. For each value of the function g, at each step in the fixed point iteration, the minimization Problem (P opt1 P opt1 P opt1 ) under the constraint v ∈ K of the non differentiable functional [START_REF] Dubois | Contact, frottement, grandes déformations élastoplastiques[END_REF] has to be solved. This problem is very similar to Problem (P m P m P m ) except for the fact that the known solution u i of the previous step is present in the friction term. The loading history, due to the velocity formulation of the friction, is characterized by this extra term. The convex K does not change from one step to the next.

An alternative formulation of Problem (P opt1

P opt1 P opt1 ) can also be obtained. Instead of choosing the displacement u g i+1 as the unknown, we can take ∆u g i+1 where :

∆u g i+1 = u i+1 g -u i g ∆t (43)
This leads to an incremental (or velocity) formulation. It leads to the minimization problem (P opt2 P opt2 P opt2 ). In this case, the convex K changes at each time step, but it does not depend on the current solution because K i is defined at the previous step by [START_REF] Felder | La tribologie de l'emboutissage[END_REF], and an extra loading term R i (v) (given by ( 47)) has to be considered : it consists the contact forces at the previous step. This time, the memory of the previous history of the solution is taken into account through two contributions. The problem set on the increments of displacement is written as follows.

Problem (P opt2 P opt2 P opt2 ) : For a given g, find ∆u g i+1 ∈ K i such that

J(∆u g i+1 ) ≤ J(v) ∀ v ∈ K i (44)
with :

J(v) = 1 2 a(v, v) + j(v) -∆L i+1 (v) -R i (v) (45)
where :

K i = v ∈ U/ v N + u N i ≤ 0 on Γ C ( 46
) (R i , v) = L i (v) -a(u i , v) (47)
The problems (P opt1 P opt1 P opt1 ) and (P opt2 P opt2 P opt2 ), associated with the fixed point problem (P inc3 P inc3 P inc3 ) are therefore two alternative ways of solving the quasistatic frictional problem at each step t i+1 of the loading. Numerical methods on these lines will be presented in the next chapter.

Mathematical results

Properties and hypotheses

The following results are based on the basic properties of the bilinear form a(u, v):

-symmetry a(u, v) = a(v, u) ∀u, v ∈ U (48) -coerciveness ∃α > 0 such that a(u, u) ≥ α u 2 H 1 v ∈ U (49)
-continuity

∃M > 0 such that | a(u, v) |≤ M u H 1 v H 1 u, v ∈ U (50)
and of the loadings assumed to be regular

Φ 1 ∈ W 1,2 (0, T ; L 2 (Ω) 2 ) Φ 2 ∈ W 1,2 (0, T ; L 2 (Γ F ) 2 ) (51)

The Signorini problem

Theorem :

The Signorini problem has a unique solution.

Proof: see [START_REF] Lions | Variational inequalities[END_REF][START_REF] Fichera | Boundary value problems in elasticity with unilateral constraints[END_REF]

The static Tresca problem

Theorem :

The Tresca problem has a unique solution.

Proof: see [START_REF] Duvaut | Les inéquations en mécanique et en physique[END_REF] where the normal contact force is given.

b -The regularized problem (non local friction)

To overcome this difficulty, a regularized problem has been proposed by Duvaut [START_REF] Duvaut | Equilibre d'un solide élastique avec contact unilatéral et frottement de Coulomb[END_REF].

As introduced in chapter 1, the contact force is then defined at each point by a convolution acting on a small set surrounding this point. This introduces the notion of non local friction. Apart from of its mathematical value, this hypothesis also has tribological relevance, since it introduces local interactions between the asperities.

Using this regularization procedure, various authors have established an existence result for this non local friction problem and an uniqueness result if the friction coefficient µ is small (see [Cocu, 1984] [ [START_REF] Demkowicz | On some existence and uniqueness results in contact problem with non local friction[END_REF]). For systems with a finite number of degrees of freedom (springs and mass systems), Klarbring [START_REF] Klarbring | Examples of non-uniqueness and non existence of solutions to quasistatic contact problem with friction[END_REF] has shown the existence of several simultaneous solutions when the friction coefficient is larger than a given limit.

Results of the same kind were obtained by Klarbring-Mikelic-Shillor [START_REF] Klarbring | Frictional contact problems with normal cpmpliance[END_REF] using the compliance model.

Theorem : For regularized problems (either non local friction or compliance models), the existence of the solutions is proved but uniqueness is obtained only under the assumption that the friction coefficient µ is small enough (this limit depends on mathematical constants which make it difficult to estimate this limit mechanicaly).

This weakness of the mathematical results for the Coulomb problem has to be kept in mind because it is also of interest to engineers undertaking computations. Large friction coefficients may lead to bad conditioning and numerical problems may arise. Nevertheless, within the classical range of mechanical estimations of the friction coefficient, the problem is generally well posed.

On the contrary, if we take a variable friction coefficient (one depending on the sliding velocity, for example), the lack of mathematical results means that great care should be taken with the computations and with the physical significance of the results. When variable friction is used, it is easy to obtain different solutions for a same problem simply by changing some of the parameters in the computations.

The loss of uniqueness, when the friction coefficient becomes large, shows that : -if µ tends towards 0, the problem will tend towards a frictionless contact (free boundary conditions), as mechanical intuition suggests, -but if µ tends towards infinity, the problem will tend towards an ill-posed problem, and not towards a problem with clamped boundary conditions, as mechanical intuition might have suggested.

On the other hand, it can be shown that the possibility of instabilities occurring with friction is enhanced by large values of the friction coefficient (see [START_REF] Martins | to appear, Dynamic stability of finite dimensional linearly elastic systems with unilateral contact and Coulomb friction[END_REF] [105] [106] [START_REF] Barbarin | Friction and instabilities : stress waves in a sliding contact[END_REF]).

c -The discretized problem with local Coulomb friction

Theorem : For the finite element problem P v P v P v (triangles P 1 with 3 nodes), existence is proved and uniqueness obtained if the friction coefficient is small enough (µ ≤ µ 0 (h), where h characterizes the meshsize of the finite element discretization). In addition, the fixed point process defined by [START_REF] Demkowicz | On some existence and uniqueness results in contact problem with non local friction[END_REF] converges if µ ≤ µ 0 (h).

Proof: see [START_REF] Licht | Remarks on a numerical method for unilateral contact including friction[END_REF] No result on the continuous problem can be deduced from this because µ 0 (h) tends to zero when h tends to zero.

The quasistatic problem with a non local friction or compliance model

Theorem : If the friction coefficient is small enough, the existence of a solution of the quasistatic problem is proved.

Proof for non local friction: see [START_REF] Cocu | Formulation and approximation of quasistatic frictional contact[END_REF] With the incremental solution (for which existence, and uniqueness can be established when the friction is small), we can construct a solution to the quasistatic problem by proving the weak convergence of a subsequence of mappings interpolating the incremental solution.

Proof for compliance model : see [Klarbring-Mikelić-Shillor, 1991].

3 Numerical methods for the Signorini problem with Coulomb friction

Introduction

A large panel of methods has been developed for solving friction contact problems. Some of them involve a regularization of the contact conditions. Without attempting to give an exhaustive list, the following methods can be mentioned.

a -Penalty methods [55] [65] [21] [51]

These are widely used in the standard finite element codes dealing with contact problems. They tolerate a penetration of one solid into the other, which can be convenient when an accurate contact solution is not needed. When penetration occurs, a high normal contact force is introduced to push the node out of the obstacle. This penetration can be greatly reduced with a good choice of the penalty parameters, but in that case, the computational time may increase drastically.

b -Lagrangian formulations and especially the augmented Lagrangian method [START_REF] Simo | An augmented Lagrangian treatment of contact problems involving friction[END_REF] [55] [START_REF] Wriggers | Penalty an augmented Lagrangian formulations for contact problems[END_REF] [120] [START_REF] Chabrand | Various numerical methods for solving unilateral contact problems with friction[END_REF] After introducing the contact force as a Lagrange multiplier (mixed formulation), a saddle point problem is obtained. It can be solved with the Uzawa algorithm. By combinating this with a penalty method, an augmented Lagrangian formulation is obtained. This is a powerful method, the accuracy of which depends on the number of augmentations. The case with only one augmentation corresponds to a penalty method.

c -Complementarity formulations and mathematical programming method

[12] [19] [21] [66] [68] [69] [77] [100] [76]
This kind of method is widely presented by A. Klarbring in the present volume. The problem has first to be set under a complementarity form. This is a natural way of formulating the unilateral conditions but new variables have to be introduced to write the friction law under this form. In addition, the problem has to be condensed, i.e. reduced to the contact variable alone by inversing the linear part. Mathematical programming methods, which are direct pivoting methods, can then be used. These methods are very powerful (fast and accurate). They can be used on 3D problems by polygonalizing the Coulomb cone [START_REF] Klarbring | A mathematical programming approach for threedimensional contact problem with friction[END_REF].

d -Projection techniques [51] [100] [76] [74] [21]

These methods are associated with the constrained minimization problem introduced in the previous section. They do not require any regularization of the contact condi-tions (except for one of them, the conjugate gradient method). They constitute very robust iterative methods : no specific difficulties arise when the friction coefficient becomes large and they yield an accurate determination of the contact forces and contact displacements.

e -Newton methods [START_REF] Curnier | A generalized Newton method for contact problems with friction[END_REF] As in finite plasticity, a Newton method combined with radial return mapping can be used. This would be a good choise when treating contact problems in finite plasticity.

[20] [START_REF] Chabrand | Various numerical methods for solving unilateral contact problems with friction[END_REF] f -interior point methods [START_REF] Auatt | An interior point optimization algorithm for contact problem in linear elasticity[END_REF] [30] These are efficient optimization methods which were recently introduced for solving contact problems.

g-Successive approximation methods [START_REF] Francavilla | A note on numerical computation of elastic contact problems[END_REF] The solution is obtained through iterates on boundary conditions on the contact zone. At each iterate, the next boundary conditions are defined by writing the unilateral and the frictional conditions deduced from the previous solution. This algorithm is fast and deals accurately with linear problems. For non linear problems, the algorithm may not converge.

In this paper, we have focused on projection and mathematical programming methods. Both methods involve Signorini conditions (strictly unilateral conditions) and the Coulomb law. They are both presented for 2D problems but the second one can be extended to 3D formulations.

The projection technique is used in three methods : Successive Overrelaxation (SSOR), Gauss-Seidel with Aitken acceleration and Preconditioned Conjugate Gradient (which needs a regularization in order to obtain differentiability). With the second of these, there is no need to determine an optimum acceleration parameter as it is the case with SSOR. A characteristic of the two last methods is that they do not require determining any computational or convergence parameters. This is a significant advantage. Among the various mathematical programming methods avalaible, the Lemke algorithm will be presented here.

Comparisons will be made between these methods by giving several examples. The problem of the identification of the friction coefficient is discussed and illustrated in the case of an industrial problem.

In this section, the various numerical methods are implemented on the Signorini problem with static Coulomb friction. In the previous section, we saw that this can be easily extended to the real Coulomb friction, by using the incremental formulation.

Projection methods

Sequence of minimization problems under constraint

As stated previously, by using a fixed point iterate on the sliding threshold, we have to solve a sequence of minimization problems under constraints (v ∈ K) of a non differentiable functional (the absolute value in the friction term j(v)). Let us first recall this problem.

Problem (P f P f P f ) : Find the fixed point of the mapping S :

S(g) = -µF N (u g ) with g ∈ G ( 52
)
where u g is the solution of the following Problem (P m P m P m ) :

Problem (P m P m P m ) : For a given g, find u g ∈ K such that

J(u g ) ≤ J(v) ∀ v ∈ K h ( 53
)
with

J(v) = 1 2 a(v, v) + j(v) -L(v) (54)

Resolution of the fixed point problem (P

f P f P f )
a -Local coordinates for the contact variables Rather than evaluating the normal and tangential components from the values of these variables in the Cartesian coordinates, we perform an initial rotation of coordinates in the finite element matrix at the beginning. The normal direction is given for each contact node (to deal with curve boundaries) and the finite element matrix is adjusted by multiplying it with the elementary rotation matrix.

b -Evaluation of the contact forces

The contact forces are not computed from the stress field by using σ.n = F , but directly from the equilibrium residues :

(F (u g ), v) = L(v) -a(u g , v). (55)
This increases the accuracy, because they are not computed by using the (rough) approximation of the stresses.

c -Diagonal algorithm

In order to reduce the computational time, the first iterates of g are obtained by solving problem P m P m P m roughly, i.e. little accuracy is demanded. The accuracy of the solution of P m P m P m is increased progressively for the following iterates. The accuracy is tested on the relative variation of g according to a norm L ∞ (ε = 10 -3 is used for the final iterates). This diagonal process reduces the computational time by two.

d -Result of convergence

In [START_REF] Licht | Remarks on a numerical method for unilateral contact including friction[END_REF] for the discretized problem, we proved that S is a contraction on condition that µ is small, therefore the iterates g converge to a solution of P f P f P f . This condition is the same as the one which ensures uniqueness of the continuous problem. In practice about 6 to 8 iterations on g are necessary.

Discretization of Problem (P m

P m P m ) (minimization under constraints)

a -The approximate problem

In order to give an approximate form of the problem (P m P m P m ), we introduce the classical finite dimensional space U h generated by the shape functions associated with a finite element approximation. The presentation is restricted here to linear elements (P 1 or Q 1 ) because the approximation is internal in that case. The following algorithms can of course be extended to high order finite elements. The space U h is the space of linear piecewise functions and h characterizes the meshsize of the discretization. We have proposed the following approximate problem in [Licht et al, 1991] [79].

Problem (P h g P h g P h g ) : Find u h ∈ K h such that : a(u h , v h -u h ) + g, Π h (|v h T | -|u h T |) ≥ L(v h -u h ) ∀v h ∈ K h (56)
The mapping Π h is the projector on the space U h . Taking the extra projection of the term |v h T | -|u h T | on U h , we can write the friction term as follows : i g i |v i | (see [START_REF] Licht | Remarks on a numerical method for unilateral contact including friction[END_REF]). It greatly simplifies the discrete problem.

The convex K h is defined by :

K h = v h ∈ U h / v h N ≤ 0 (57)
b -The discrete problem Problem (P discret P discret P discret ) : Find u ∈ Kh such that :

J(u) ≤ J(v) ∀ v ∈ Kh ( 58
)
with :

J(v) = 1 2 2N h i=1 2N h j=1 a ij v i v j - 2N h i=1 f i v i + i∈I T g m i |v i | (59)
with :

a ij = a(ω i , ω j ) (60)
f i = (f, ω i ) (61)
g m i = Γ C g m ω i dl (62)
Kh = v h ∈ R 2N h / (v h N ) i ≤ 0 if i ∈ I N = 2N h i=1 K i (63)
with,

K i = R if i ∈ I N = R -if not ( 64
)
where I N and I T are the sets of numberss defined by 2i -

1 ∈ I N (normal component) and 2i ∈ I T (tangential component) if i is a contact node number.
Three algorithms are presented in the next subsection for solving Problem (P discret P discret P discret ).

Successive Over Relaxation with Projection (SSORP)

a -The algorithm This method belongs to the class of punctual relaxation methods. The minimization is conducted on one component after the other. As it is a sequence of one dimensional minimizations, it is very simple to cope with the absolute value (as can be seen below).

A punctual relaxation method is an iterative method where the step n + 1 is written :

for i = 1, . . . , 2N h ,      Find u n+1 i ∈ K i such that ∀ v ∈ K i J(u n+1 1 , . . . , u n+1 i , u n i+1 , . . . , u n 2N h ) ≤ J(u n+1 1 , . . . , v, u n i+1 , . . . , u n 2N h ) (65)
The relaxation algorithm (Gauss-Seidel method accelerated by the relaxation coefficient ω ∈]0, 2[) is adapted to the constraint v ∈ Kh and to the treatment of the absolute value in i∈I T g m i |v i | in the following way :

                         u n+ 1 2 i = 1 a ii   f i - i-1 j=1 a ij u n+1 j - 2N h j=i+1 a ij u n j -sgn(u n+ 1 2 i ) g m i   if u n+ 1 2 i ∈ R * u n+ 1 2 i = 0 if not u n+1 i = P K i (1 -ω)u n i + ωu n+ 1 2 i (66)
where :

• If i ∈ I N and i ∈ I T (not a component of a contact node), we have the classical relaxation algorithm :

         u n+ 1 2 i = 1 a ii   f i - i-1 j=1 a ij u n+1 j - 2N h j=i+1 a ij u n j   u n+1 i = (1 -ω)u n i + ωu n+ 1 2 i (67)
• if i ∈ I N (normal component of a contact node), a projection P K i is performed on the convex Kh (this is very simple because it is a projection on R -) [START_REF] Glowinski | Analyse numérique des inéquations variationnelles[END_REF]:

           u n+ 1 2 i = 1 a ii   f i - i-1 j=1 a ij u n+1 j - 2N h j=i+1 a ij u n j   u n+1 i = P K i (1 -ω)u n i + ωu n+ 1 2 i (68)
• if i ∈ I T (tangential component of a contact node), the absolute value is simply treated with two trials on the sign of u n+1 i [START_REF] Raous | Numerical methods for solving unilateral contact problems with friction[END_REF] :

                         u n+ 1 2 i = 1 a ii   f i - i-1 j=1 a ij u n+1 j - 2N h j=i+1 a ij u n j -sgn(u n+ 1 2 i ) g m i   if u n+ 1 2 i ∈ R * u n+ 1 2 i = 0 if not u n+1 i = (1 -ω)u n i + ωu n+ 1 2 i ( 69
)
As the residue (defect in the equilibrium) does not vanish in our case because of the constraints (it gives the contact forces), no convergence test on the residue is available. The convergence test is therefore written on the relative variation in norm L 2 of u i (with ε = .5 10 -5 ).

1 2N h 2N h j=1 (u n+1 j -u n j ) 2 2N h j=1 (u n+1 j ) 2 < ε (70)

b -Advantages and disadvantages

This method is very easy to implement. In a finite element code, the classical solver (Cholesky, Crout, ...) has to be replaced by the SSOR solver in which slight modifications [START_REF] Klarbring | Contact problems with friction by linear complementarity[END_REF] and [START_REF] Klarbring | A mathematical programming approach to contact problems with friction and varying surfaces[END_REF] have to be carried out.

As shown taking the examples at the end of the section, this method give an accurate determination of the contact conditions (as regards both displacements and forces). Moreover, the algoritm is robust : in [START_REF] Licht | Remarks on a numerical method for unilateral contact including friction[END_REF], examples with very large friction coefficients are presented.

However the number of iterates (fixed point iterates included) can be large : usually of the same order as the size of the finite element matrix, it can even be larger when the matrix is ill-conditioned. Nevertheless, the cost of one iterate is small, especiall when a sparse matrix storage is used (well adapted to evaluating sums). Computational times will be presented in the examples.

Another disadvantage is the need to estimate an optimum relaxation coefficient to ensure a rapid convergence, and no theoretical results are available for this estimate (mainly because of the projection). A dichotomic evaluation of this coefficient is worth while only in cases where several loading cases or an incremental problem are to be treated. To avoid this estimate, the Gauss-Seidel algorithm with an Aitken acceleration, which will be described in the next subsection, can be used.

d -Convergence results

The following theoretical convergence results have been proved.

a) -Frictionless contact

In this case, J is differentiable and the constraints are separate (K is a product set). Convergence has been proved by [START_REF] Glowinski | Analyse numérique des inéquations variationnelles[END_REF]].

b) -Bilateral contact with Tresca friction

Here, J is non-differentiable but is written i α i |u i |, α i ≥ 0 and there are no constraints. Because of the specific form of the nondifferentiable part, convergence is also established in [START_REF] Glowinski | Analyse numérique des inéquations variationnelles[END_REF]].

c) -Unilateral contact with Coulomb friction

The functional J is still non-differentiable and the minimization is conducted with separate constraints. Convergence is proved in Grégo's thesis [START_REF] Grégo | Méthodes multiniveaux pour des problèmes de contact unilatéral avec frottement[END_REF].

Gauss-Seidel with Aitken acceleration and projection a -The method

The Aitken acceleration procedure [START_REF] Aitken | On the iterative solution of a system of linear equations[END_REF] was developed to accelerate the Gauss-Seidel algorithm by Irons and Tuck [START_REF] Irons | A Version of the Aitken Accelerator for Computer Iterations[END_REF]. We have extended this method by dealing with the projection and the absolute value. A new implementation increases the convergence rate.

Given a convergent sequence of real numbers {x i }, the so-called Aitken-∆ 2 process transforms this sequence into another one, which in general converges faster than the original one. The {x ′ i } sequence is given by :

x ′ i = x i - (x i+1 -x i ) 2 x i+2 -2x i+1 + x i (71)
We generally do not want to compute the sequence {x i } to determine the sequence {x ′ i } ; we want to accelerate the convergence of the sequence {x i } by directly calculating a sequence {x * i } which differs from what sequence {x ′ i } would have been. Various alternatives can be used. In what follows, the sequence {x * i } will be denoted {x i } because we build only one sequence. In the sequel, T (x i ) will denote the term x i+1 of the Gauss-Seidel sequence.

Starting with an iterative process x i+1 = T (x i ), (i = 0, 1, 2 . . .), an accelerated one, proposed by Irons and Tuck [START_REF] Irons | A Version of the Aitken Accelerator for Computer Iterations[END_REF], can be written as a relaxation procedure as follows:

x i = T (x i-1 ) + R i (x i-1 -T (x i-1 )) ( 72
)
where R i is a relaxation factor that can be defined by the recurrence relation

R i = R i-1 + (R i-1 -1) (x i-1 -T (x i-1 )) (x i-2 -T (x i-2 )) -(x i-1 -T (x i-1 )) (73)
with the starting value R 0 = 0 (no acceleration). This version requires only one additional vector for storage. The vector case will be dealt with below.

After testing various forms of Aitken acceleration, the Irons and Tuck procedure was adopted as a basis, in the constrained minimization problem of the non-differentiable functional J(v). The Aitken procedure has turned out to be efficient provided a few precautions are taken. In this respect, let us now make a few comments about the implementation of the procedure for the non linear problem presented here.

First of all, when is the projection on the convex set K to be made? The two alternatives consist of either projecting after the Aitken acceleration or projecting twice : once during the Gauss-Seidel iteration and again after the Aitken procedure. After testing both possibilities, the second one was found to be preferable. The Gauss-Seidel procedure with projection is first applied to all degrees of freedom. The Aitken acceleration is then performed followed by an extra projection. In that way, the solution never leaves the convex set K.

When using this method, either chaotic behaviour was observed or the convergence tests were satisfied but the solution was incorrect. Irons and Tuck suggested a weighting procedure between R i and R i-1 in order to smooth the convergence process. We propose another technique which consists of applying Aitken's acceleration only once every k Gauss-Seidel iterations. The idea is to combine the smooth convergence of the Gauss-Seidel process (which is known to efficiently reduce the high error frequencies) and the acceleration properties of the Aitken procedure : the acceleration is applied only when the iteration counter i is a multiple of k.

The iteration is written in ( 14) and ( 15) for the N -dimensional case. The variable is the N dimensional component displacement vector V , operator T denotes the Gauss-Seidel procedure with Projection (GSP) and consequently T (V ) and ∆, the difference between V and T (V ), are also N dimensional vectors, and R is still a scalar :

∆ i-2 = V i-2 -T (V i-2 ) ( 74
) ∆ i-1 = V i-1 -T (V i-1 ) (75)
V i = T (V i-1 ) + R i ∆ i-1 (76)
R i = R i-1 + (R i-1 -1) (∆ i-2 -∆ i-1 ) T ∆ i-1 (∆ i-2 -∆ i-1 ) T (∆ i-2 -∆ i-1 ) (77)
As given in relations ( 76) and ( 77), the Aitken acceleration is applied using the previous values V i-2 and V i-1 . The number k must be larger than two, because for k = 1 or k = 2, the acceleration is computed using a combination of a previous Aitken accelerated iterate and a regular Gauss Seidel one. This is an important difference with respect to the Irons and Tuck algorithm : these authors calculated the acceleration using the previously accelerated iterates. Here we compute the acceleration every k th iteration of the regular Gauss Seidel iterates. In all the cases studied so far, the value k = 5 has proved to be a good choice. Because of the non-smooth convergence of the Aitken acceleration, the convergence test (70) can be satisfied although the solution is not correct. To avoid this, where the convergence tests are satisfied for the accelerated procedure, m additional Gauss-Seidel iterations are performed. The value of m is chosen so be equal to 2% or 3% of the total number of degrees of freedom.

b -Convergence

Despite the lack of convergence results for the Aitken procedure, we can that our iterates converge because the Aitken iterations can be viewed as just providing a very good starting vector for the final Gauss-Seidel iterations for which convergence has been proved. In practice, the projected Gauss-Seidel method with Aitken acceleration turns out to be as fast as the SSOR with projection previously presented, without having to estimate an optimum coefficient. This is a great advantage.

Projected Conjugate Gradient with Preconditioning [102]

Two difficulties arise when using conjugate gradient methods to solve Problem (P discret P discret P discret ) : how to treat the constraint (∈ K) and the non differentiable term (J(v). To implement an algoritm based on the conjugate gradient, a specic projection has been introduced to preserve the conjugate property between the successive descent directions, and the differentiability is obtained by a regularization as described in section 1.4.4. . Details of the method presented in this section can be found in Raous-Barbarin [START_REF] Raous | Conjugate Gradient for Frictional Contact[END_REF]. It is a variant of the methods proposed by May [START_REF]The conjugate gradient method for unilateral problems[END_REF], Dilintas et al [START_REF] Dilintas | A conjugate projected gradient method with preconditionning for unilateral problems[END_REF], and Marks-Salamon [START_REF] Marks | A projected conjugate gradient method for frictionless contact problems[END_REF] for frictionless problems, extended to the frictional case.

Preconditioning is essential to improve the efficiency of the conjugate gradient method. Four preconditioning procedures are tested here. Two of them are specific to friction problems and take the regularization into account ; this is an important point because the conditioning of the problem depends on the regularization parameters.

a -Treatment of the constraint v ∈ K. With a conjugate gradient method, it is not possible to use a simple projection as was done previously with the Gauss-Seidel method, because the conjugate property between two successive directions would be lost in that case. In order to preserve the fundamental property of the conjugate gradient, the solution will be sought for in the subspace R N -p , where p is the number of saturated constraints. This means that whenever a penetration occurs, the system will lose a degree of freedom so that the next directions and the conjugation will be defined relative to the new subspace. The convergence in a finite number of iterates is therefore preserved in each subspace, where the conjugation is now written. The condition of compressive normal tension is naturally tested, and a constraint can be relaxed if necessary. There exist various ways ot carrying out this process. May [START_REF]The conjugate gradient method for unilateral problems[END_REF] introduces evolutive boundary conditions into the matrix itself. The conjugation is still written with the current matrix, but these modifications are costly. Dilintas et al [START_REF] Dilintas | A conjugate projected gradient method with preconditionning for unilateral problems[END_REF] and Jeuzette-Sonzogni [START_REF] Jeuzette | A projected conjugate gradient method for structural stability analysis with linear constraints[END_REF] project the gradient direction. This means that, for a descent direction which would induce penetration, the solution is stopped on the boundary of the convex, and will keep going in the tangential direction, i.e. on the boundary of the convex. We use this process, and we check if the normal force is a tension, we relax the constraint.

The main steps in the algorithm are as follows :

-regular gradient conjugate iterate with evaluation of the optimal coefficientρ Opt for the descent, -for all the contact nodes, in the case of penetration, computation of the maximum value ρ i of the descent coefficient which ensures that the solution stays in the convex, -the final descent coefficient is the minimum between ρ opt and the ρ i , -when a new contact is obtained, thecorresponding normal component is eliminated. The next conjugate gradient iteration is then started with a regular gradient direction, -the dual condition F N ≤ 0 has to be checked over the whole the contact zone. If the condition is not fulfilled, the constraint u N = 0 is relaxed and the descent direction is computed again.

A diagonal process is developed : the condition on F N is checked after partial con-vergence of the conjugate gradient (in the convergence test, the final accuracy required is taken to be 10 -3 instead of .5 10 -5 ).

b -Treatment of the friction term

In order to compute the gradients, the non differentiable term i∈I T g m i |v i | of (59) has to be regularized. The four regularizations introduced in chapter 1 will be compared. The functional (59) becomes :

J ε (x) = 1 2 x T .A.x -f T .x + m i=1 g i ϕ ε (x num(i) ) ( 78
)
where the three regularizations, tested in the next section, are denoted by :

• Square root : ϕ 1 ε (x) = √ x 2 + ε 2 ,
• Hyperbolic :

ϕ 2 ε (x) = ε ln cosh x ε ,
• Piecewise polynomial :

ϕ 3 ε (x) =                -x - ε 2 , if x < -ε x 2 2ε , if -ε ≤ x ≤ ε x - ε 2 , if x > ε .

c -Conjugate Gradient for the non quadratic case

The classical conjugate gradient algorithm can be found in [START_REF] Lascaux | Analyse numérique matricielle appliquée à l'art de l'ingénieur[END_REF]. After regularization, the conjugation between the successive directions v k (k denotes here a conjugate gradient iterate) should be written relatively to the Hessian J ′′ ε of (78). To avoid the computation of the Hessian, J ′′ ε is approximated with a Taylor development of J ′ ε (x k ) = -r k . This gives the new conjugate direction (where P (.) denotes the projection previously defined):

V k+1 = P (r k+1 ) + γ k .V k with γ k = < r k+1 -r k , P (r k+1 ) > < r k , P (r k ) > (79)
The optimal descent coefficient is then computed by solving :

< J ′ ε (x k + ρ k .V k ), V k >= 0 ( 80
)
This non linear problem can be solved in two different ways :

-either by using a Taylor development of

J ′ ε (x k + ρ k .V k ) ρ k = < r k , V k > < J ′′ ε (x k ).V k , V k > (81)
-or by conducting a few iterations of the false position method to avoid having to computate J ′′ ε (with a prescribed acuracy test on the relative variation of ρ k equal to 10 -4 , the number of iterates is about twenty).

The two procedures have been tested and the first one is preferable because it requires no iterates and it yields nearly the same number of gradient iterates. The computational time is therefore shorter.

d -Preconditioning

It is well known that the convergence can be considerably improved by applying a preconditioning procedure (see [START_REF] Lascaux | Analyse numérique matricielle appliquée à l'art de l'ingénieur[END_REF]). In this case, a change of variable is made in the algorithm and an extra resolution is needed :

C(x k ).z k+1 = r k+1 (82)
where C(x k ) is the preconditioning matrix.

Let us decompose the finite element matrix A = D -(E + E t ) where D is the diagonal and E the opposite of the lower triangular part. Four preconditioning procedures have been tested. In the first two only the matrix A is used for the preconditioning, and the last two the Hessian J ′′ ε (x k ) is used, thus accounting for the terms resulting from the regularization (a few extra diagonal terms):

(J ′′ ε (x k )) ii = A ii + g i .φ" ε (x i k ) if i ∈ I (83) (J ′′ ε (x k )) ij = A ij if i = j or i = j ∈ I (84)
In this case, the preconditioning is related to the regularization parameter ε.

The four preconditioning procedures tested were :

• C = D = diag(A),
• Evans preconditioning (with ω = 1) :

C = [D + E] D -1 [D + E] t • C = D = diag(J ′′ ε (x k )),
• Evans preconditioning (with ω = 1) with

J ′′ ε (x k ) : C(x k ) = J ′′ ε (x k ) + E (J ′′ ε (x k )) -1 J ′′ ε (x k ) + E t
In conclusion, the various tests conducted in [Raous- [START_REF] Raous | Conjugate Gradient for Frictional Contact[END_REF]] have shown that :

-the piecewise polynomial regularization (or the hyperbolic one) is the most accurate using ε = 10 -7 (comparisons were made with the results obtained using methods with no regularization), -the Evans preconditioning methods applied to the Hessian J ′′ ε is the most efficient. In the numerical results presented, these choices will be used.

Mathematical programming -Lemke method

(see ( [START_REF] Lemke | A survey of complementarity theory[END_REF], [Klarbring, 1986[START_REF] Klarbring | Contact problems with friction by linear complementarity[END_REF][Klarbring, , 1988]], [START_REF] Raous | Numerical methods for solving unilateral contact problems with friction[END_REF]], [Chabrand-Dubois-Raous, 1995[START_REF] Martinet | Etude et mise en oeuvre de formulation ALE -Application au contact lubrifié en mise en forme[END_REF]])

In the present volume, A. Klarbring gives a wide presentation of this kind of method. We will only briefly recall these formulations and give the Lemke algorithm used in the next section, where we compare it with the other methods presented in this chapter. These mathematical programming methods are convenient for solving complementarity problems written as follows.

Problem (P 1 P 1 P 1 ) : Find w ∈ R p , z ∈ R p such that w -M z = q (85) w i ≥ 0, z i ≥ 0 i = 1 . . . p (86) w i z i = 0 i = 1 . . . p (87)

Complementarity formulation of the friction

Let us first write the frictional contact problem as a complementarity one. Unilateral contact conditions are already written as a complementarity problem. New variables have to be introduced in order to friction as a complementarity problem.

-unilateral contact

             u N ≤ 0 F N ≤ 0 u N F N = 0 (88)

-friction (introduction of an extra variable)

The tangential displacement u T is divided into a positive sliding component λ 1 and a negative sliding component λ 2 , and we set :

φ 1 = -F T + µF N (89) φ 2 = +F T + µF N ( 90
)
The friction problem can then be written :

                           u T = λ 2 -λ 1 F T = (φ 2 -φ 1 )/2 λ i ≤ 0 φ i ≤ 0 i = 1, 2
λ i φ i = 0 (91)

Important notes

-the matrix A * is full-sized, -the matrix A * is not symmetric and depends on the friction coefficient.

Condensation

In order to write the problem in the same form as problem (P 1 ), condensation is perfomed to reduce the problem to the contact variables alone. This consists of performing partial inversion on the linear equations of the system. This is possible without inversing any large matrix. In what follows, the subscript c denotes the degrees of freedom involved in the contact and l the other ones. The matrix A is written as follows.

  

A ll A lc

A cl A cc       u l u c    =    f l 0    +    0 F c    (92)
We can then write the problem as follows :

A * u c = F c + f * (93) u c ≤ 0 (94) F c ≤ 0 (95) u c F c = 0 (96) A ll u l = f l -A lc u c ( 97
)
where A * et f * are given by :

A * = A cc -A cl A -1 ll A lc (98) f * = -A -1 ll A lc f l ( 99
)
To calculate A -1 ll A lc , we solve:

A ll X i = Y i ( 100
)
where X i denotes the i th column of A -1 ll A lc and Y i the i th column of A lc .

The reduced discrete problem (the size of which is only twice the number of contact nodes) can then be written :

A * u c = f * + F c (101)

The Lemke method

Problem (P 1 P 1 P 1 ) can now be written as Problem (P 2 P 2 P 2 ) by taking :

z = -u c , w = -F c , q = f * ( 102
)
The method is a pivoting technique based on the following remark :

if q i ≥ 0 ∀ i = 1 . . . p q i ≥ 0 ∀ i = 1 . . . p q i ≥ 0 ∀ i = 1 . . . p, the solution is z i = 0 z i = 0 z i = 0 et w i = q i w i = q i w i = q i .
When a subscript "j" is such that q j < 0, a new real variable z 0 is introduced. Problem (P 1 P 1 P 1 ) becomes :

Problem (P 2 P 2 P 2 ) :                    Find w ∈ R p , z ∈ R p et z 0 ∈ R such that Iw -M z -1Iz 0 = q z 0 ≥ 0 w i ≥ 0, z i ≥ 0, w i z i = 0 i = 1 . . . p
where I is the identity matrix in R p , 1I T = (1,1,1, ...,1,1) .

Let A 0 = [ I, -M, 1I ] be the matrix with p lines and 2p+1 columns, X and q 0 the vectors of R 2p+1 and R p defined by :

X =        w z z 0        et q 0 = q
Introducing the variable z 0 enables us to get a positive second member, but z 0 is then introduced into the basis (see presentation of the algorithm (i)). The goal of the algorithm is to eliminate z 0 from the basis by pivoting, preserving at the same time the positivity of the second member.

Presentation of the algorithm (see [START_REF] Bazarra | Nonlinear programming-Theory and algorithms[END_REF] [54] for details) (i) Starting : Let s be the subscript such that :

q s = min i=1...p q i
The matrix A 1 and the second member q 1 are then defined by :

                   A 1 ij = A 0 ij -A 0 sj j = 1 . . . 2p = 1, i = s A 1 sj = -A 0 sj j = 1 . . . 2p + 1 q 1 i = q 0 i -q 0 s pour i = s q 0 s = -q 0 s
And the system to be solved becomes :

A 1 X = q 1
with the complementarity conditions of Problem (P 2 P 2 P 2 ).

(ii) Step k+1 : Let t be the column subscript of the variable which is complementary to the one eliminated from the basis during the previous step and r the line subscript such that :

q k r A k rt = min i=1...p q k i A k it , A k it > 0
The matrix A k+1 and the second member q k+1 are then :

                             A k+1 ij = A k ij - A k it A k rt A k rj j = 1 . . . 2p + 1, i = r A k+1 rj = A k rj A k rt j = 1 . . . p q k+1 i = q k i - A k it A k rt q k r i = r q k+1 r = q k r A k rt
X t therefore replaces the r th basis variable in the basis. If r = s, z 0 is eliminated from the basis and the solution is reached : the p variables of the basis will then be equal respectively to the p components of q k+1 and their complementary variables and z 0 will be zero.

(iii) Stop : Go to (ii) until z 0 is eliminated from the basis.

Convergence -theory

Matrix M has to be copositive (it depends on the friction coefficient) : it is not easy to give an explicit condition on µ.

-in practice :

. fails when the number of contact degrees of freedom is very large (full-sized matrix, accumulated round-off errors when the numbers of pivoting is very large, ...) ; but in this case, we can still use the Lemke method by introducing multigrid methods (see [START_REF] Raous | Convergence et efficacité de méthodes multigrilles pour un problème de contact avec frottement[END_REF]] [Grego, 1995]) because the solver works then only on the coarse grid.

. sometimes fails on medium-sized problems, if zeros occur on the diagonal (line permutations can eliminate this problem, but this increases the computational time).

Presentation of the examples used in the next section

In the last part of this chapter, we will describe some tests performed on the above methods, comparing of the computational times, the accuracy of the contact solutions and making some comments about the problem of how to identify the friction coefficients when dealing with industrial cases. Let us first present the various examples.

Dovetail assembly

Details of this dovetail assembly can be found in [START_REF] Raous | Numerical methods in mechanics of contact involving friction[END_REF]. The geometry is given in Fig. 15 and the mesh in Fig. 16 (4345 nodes, 49 contact nodes). The topology of the mesh is related to multigrid tests which have been conducted on this example.

A block sliding on a plane

This example (see Fig. 17) is a simulation of an experiment conducted by [START_REF] Zeghloul | Stress waves in a sliding contact -Part 1 : Experimental study[END_REF]] : a polyurethane parallelepipedic block is pressed on an araldite plate on which a tangential displacement is prescribed. Studies on stability analysis have been applied to this example : details can be found in [START_REF] Raous | Friction induced instabilities and sound generation[END_REF] in [START_REF] Raous | Numerical methods for solving unilateral contact problems with friction[END_REF]] and [START_REF] Lebon | Multibody contact problem including friction in structure assembly[END_REF]]. The geometry of the closure of the pressure vessel is given in Fig. 19 where there are three contact zones with three different friction coefficients. The mesh is given in Fig. 19 (674 nodes, 54 contact nodes).

Central plate

High pressure screw-press

Details about this example can be found in [START_REF] Grégo | Méthodes multiniveaux pour des problèmes de contact unilatéral avec frottement[END_REF]. The geometry and the mesh are given in Fig. 20 (11933 nodes, 250 contact nodes). Details about the following two examples can be found in [START_REF] Chabrand | Various numerical methods for solving unilateral contact problems with friction[END_REF]]. The mesh (using Q4/P0 elements) is given in Fig. 21 (441 nodes, 21 contact nodes).

Figure 21: Indentation test

Extrusion of an aluminium cylinder

The mesh (with Q4/P0 elements) is given in Fig. 22 (105 nodes, 21 contact nodes).

Figure 22: Extrusion of an aluminium cylinder

Comparison between the numerical methods

In this section, the numerical methods previously introduced in this chapter are tested and compared on the examples presented in section 3.4 in terms of the number of iterates and the computational times required. These computational times are essentially given to be compared to one another. They are quiet long because these tests were conducted by using a CISM processor (VAX 6310). The CPU times are actually about one hundred times smaller with a CISC processor (DEC 4100). Table 4: Comparison between the CPU times for the screw-press

The dovetail assembly

The plastic indentation test

The following two problems are dealt with in the context of finite plastic deformations (see [START_REF] Chabrand | Various numerical methods for solving unilateral contact problems with friction[END_REF]] [START_REF] Dubois | Contact, frottement, grandes déformations élastoplastiques[END_REF]. The SSOR with projection and the Lemke methods have been extended to this non linear behaviour (they operate on the tangent matrix), and they are compared here with augmented Lagrangian and penalty methods. The penalty method corresponds to the augmented Lagrangian with only one augmentation, but the solution is less accurate.

Accuracy of the contact solution

Solving a Signorini problem (no regularization of the unilateral conditions) with a non regularized Coulomb friction law (except for the conjugate gradient) lead to an accurate determination of the contact condition (in terms of both the displacements and the stresses). In order to illustrate this, let us take the evolution of the displacements and the forces on the contact zone in the case of the block sliding on a plate when the prescribed tangential displacement of the plate increases.

Table 7 gives the contact forces and the contact displacements at a given step during the loading. Three zones can be clearly identified :

-a separate part where u N < 0 and F N = F T = 0 with good accuracy (10 -6 ), -a sliding part where u N = 0, u T = 0, F T = µF N (with an accuracy of less than 0.1 per cent, µ = 0.9 here), -and an adhesive part where u N = 0, u T = 0, and F T = µF N .

In Fig. 23, the changes in the contact forces during the loading are presented. The curves

Some comments about the validation of the model and the identification of the parameters

A few comments need to be made about the validation of the computational code, the identification of the parameters and the validation of the model.

It is clear that the validation of a numerical method is not conducted by making comparisons with experimental results. Although it is not very easy to validate a computational code for non linear problems (and especially for unilateral problems with friction), there exist various ways, which consist of comparing the numerical results with :

-exact analytical solutions (Hertz solution of the indentation of an infinite plane by a cylinder, for example), -other numerical results obtained on classical Benchmarks (a few of which can be found in [START_REF] Raous | Numerical methods in mechanics of contact involving friction[END_REF]), -results obtained with a linear procedure for limit cases (frictionless problem with persistent contact for example).

On the contrary, a model has to be validated by comparing experimental results with the numerical simulation, once the computational code has been validated, and once the parameters of the model have been identified. It turns out that in the case of friction problems, the identification of the friction coefficients is somewhat difficult. On the one hand, tables giving the friction coefficients of various materials can be found in technical handbooks, and on the other hand, specific measurements can be carried out on samples of the constitutive material of the structure. Several difficulties arise, however :

-a contact involves two bodies and the friction coefficient is not a characteristic of one of the materials but of both, -the measurement of a friction coefficient (tribometer, Inland test, etc.) is very delicate : parallelism of the tools, influence of the experimental conditions (velocity, normal force, temperature, hygrometry, presence of a third body such as oil residue or microscopic fragments, etc.), -the surface condition of the sample may be different from that of the same material in the manufactured structure.

Moreover, depending on the geometry, the solution can be very sensitive to a slight variation in the friction coefficients, which means that an accurate determination is necessary.

It is recommended, when identifying the friction coefficients, to be as close as possible to the real situation. One procedure consists of conducting the identification on a significant step of the loading by comparing the results of the experiment with those of the simulation. Then, keeping these parameters constant, the comparison between theory and experiment on the other steps of the loading will either validate the model or not, depending on the agreement between these results.

We have illustrated this in [START_REF] Lebon | Multibody contact problem including friction in structure assembly[END_REF], taking the example of a scale structure of the closure of a pressure vessel (which was presented in section 3.4.3). In this structure, the three contact parts have been treated differently (different coatings, different polishings). The three friction coefficients were identified by comparing the experimental results with those of the simulation during the first step of the loading (bolt tightening). Solutions obtained with various friction coefficients can be found in Tab. 8. The comparison involves four quantities, denoted : deflection (at the center of the plate), rotation (at the level of the closure), upper part/plate (sliding length) and thrust (sliding length). It can be noted that the value of the various friction coefficients greatly affects the solution (experimental values are given at the bottom of the table). The selected friction coefficients where µ 1 = 0.2, µ 2 = 0.1, µ 3 = 0.1. For these values, the differences between experimental and theoritical results are given in Tab. 9. Then, with these coefficients, the simulations of the next two loading steps were computed (lower part pressure and upper part pressure). Conclusions as to the validity of the model can then be drawn, based on comparisons between theory and experiment as regards the new loading (see Tab. 11). The large difference observed between the "deflection" values was due to another strong approximation concerning the geometry of the plate, which also has some influence on the "rotation". Taking this fact into account and in view of the sensitivity of the system to the value of the parameters and the experimental accuracy, the Coulomb model has been considered to be an appropriate model here. Note that in the case of some other friction coefficients, differences of 300% were observed in Tab. 8.

µ
4 Model coupling adhesion to unilateral contact and friction

Introduction

We now consider a problem where the contact is initially adhesive, which means that it can admit normal traction (F N > 0) and shear force with neither loss of contact nor sliding. When the forces are large enough, the adhesion can be progressively broken and classical unilateral contact with friction will have to be considered.

A model ensuring a continuous transition from a total adhesive condition to a unilateral frictional problem is proposed. It is based on the formulation of adhesion, introduced by M. Frémond (see [START_REF] Frémond | Adhérence des solides[END_REF], [START_REF] Truong Dinh | Contact avec adhérence[END_REF], [START_REF] Frémond | Contact with adhesion[END_REF]), using an intensity of adhesion β on the interface (which is similar to a damage function). This approach is described in details by M. Frémond in the present volume. Here we will focus on the coupling to the friction. An extensive presentation of the model can be found in Raous-Cangémi-Cocu [START_REF] Raous | to appear, A consistent model coupling adhesion, friction and unilateral contact[END_REF] [28] and in other works [START_REF] Raous | Un modèle couplant adhérence et frotement pour le contact entre deux solides déformables[END_REF] [15] [16] [START_REF] Cangémi | Frottement et adhérence : modèle, traitement numérique et application à l'interface fibre/matrice[END_REF].

In the framework of continuum thermodynamics, the contact zone is taken to be a material boundary and the local constitutive laws are obtained by choosing two specific surface potentials : the free energy and the dissipation potential. Since these potentials are not regular, convex analysis is used to obtain the local behavior laws from the state and the complementary laws. The continuous transition from a total adhesive condition to a possible purely frictional one is enforced by using elasticity coupled with damage on the interface. The variational formulation for quasistatic problems is written as the coupling between an implicit variational inequality, a variational inequality, and a differential equation. The model will be used to simulate a micro-indentation experiment conducted in order to characterize the behavior of a fiber/matrix interface in a ceramic composite.

The model

We will first present the model and then specify the thermodynamic basis in the next section. Contact between two deformable solids is considered and sligtly different notations are introduced. Let Ω 1 and Ω 2 be two domains of R d (d = 2, 3) occupied by two continuous bodies, the boundary of each domain being separated into three separate parts :

∂Ω α = Γ α U ∪ Γ α φ ∪ Γ α C , α = 1, 2 . Let Γ 1 φ , Γ 2
φ denote the parts of the boundary where external forces φ 1 , φ 2 are respectively applied (see Fig. 24). Likewise, Γ 1 U and Γ 2

U C U Γ 1 U Γ φ 2 C Γ 2 φ Γ 1 n Γ 2 φ 1 φ 1 2 n f C Γ Ω 2 1 Γ 2 f Ω 1 1 2
Figure 24: Contact between two elastic bodies are the parts of ∂Ω 1 , ∂Ω 2 where the displacements are prescribed and f 1 and f 2 denote the imposed volume force densities. Let Γ 1 C , Γ 2 C denote the parts of the boundary where the two solids are initially in contact. Assuming small displacements hypothesis, we have n 1 = -n 2 , where n α denotes the outward normal unit vector to ∂Ω α , α = 1, 2.

Each particle of Γ 1

C is therefore coupled to a particle of Γ 2 C in a single valued correspondence. Consequently, the two material boundaries Γ α C ⊂ ∂Ω α (α = 1, 2) define at the beginning a common contact zone denoted Γ

C : Γ 1 C ≃ Γ 2 C ≃ Γ C .
The relative displacement between the two bodies is defined on Γ C by u = u 2 -u 1 where u 1 and u 2 are the displacements of two corresponding points. The vector n 1 is chosen (n = n 1 = -n 2 ) for the normal/tangential decomposition : u = u N n + u T with u N = u.n. According to this choice, the non penetration condition will be written u N ≥ 0 (which is convenient for a gap condition). When the second solid is rigid (u 2 = 0), the usual Signorini condition u N = u 1 .n = u 1 .n 1 ≤ 0 is obtained. The internal force on Γ C is denoted R. Using the principle of virtual power, we obtain : R = σ 1 n 1 = -σ 2 n 2 , according to previous choice of u. The contact force is separated into reversible and irreversible parts : R = R r + R ir .

The kinematic variables on the contact boundary are : • u, jump of displacement on an interface, • β, intensity of adhesion (surface damage variable) where β = 1 means total adhesion, 0 < β < 1 partial adhesion, and β = 0 breakdown. The interface law coupling adhesion, friction and unilateral contact is then written as follows.

If the loading remains constant (when the adhesion limit is overcome), the adhesion keeps decreasing (by relaxation). If the tangential loading is now backward, an opposite tangential displacement occurs only when the other side of the Coulomb cone is reached.

Thermodynamic framework

The general framework has been presented by M. Frémond. It is based on considering the contact area as a material boundary. Let us focus here on the specific aspects, namely the choice of the free energy and the dissipation potential.

Free energies and state laws

Onto Γ C , the surface density of free energy Ψ is chosen as follows :

Ψ(u N , u T , β) = C N 2 u 2 N β 2 + C T 2 u T 2 β 2 -w h(β) + I K (u N ) + I P (β) (106)
where K = {v / v ≥ 0} and P = {γ / 0 ≤ γ ≤ 1}.

Introducting of the indicator functions I K and I P imposes the unilateral condition u N ≥ 0 and the condition β ∈ [0, 1]. The interfacial forces induced by the adhesion are introduced under the form of a compliance law depending on the current state of adhesion β and characterized by the initial stiffness C N and C T . The term w h(β) is a general form of the energy of decohesion, which is assumed to be differentiable. With h(β) = β, w corresponds to the Dupré's energy. Other forms of h(β) could be adopted.

The pseudo-potential Ψ has a part (the first three terms) that is differentiable but not convex and a part (the last two terms) that is convex but not differentiable. To write the state laws, the two difficulties (lack of convexity and lack of differentiability) are overcome by using local or partial subdifferentiation (see [START_REF] Raous | to appear, A consistent model coupling adhesion, friction and unilateral contact[END_REF] [107] [17] [49]). The state laws can then be written as follows :

R r N ∈ ∂ u N Ψ(u N , u T , β) (107) R r T ∈ ∂ u T Ψ(u N , u T , β) (108) -G β ∈ ∂ β Ψ(u N , u T , β) (109)
where ∂ u and ∂ β denote the subdifferential with respect to the variables u and β respectively, and where R (contact force) and G β are the thermodynamic forces associated with the jump of contact displacement u and with the adhesion intensity β, respectively.

The following relation [START_REF] Simo | An augmented Lagrangian treatment of contact problems involving friction[END_REF] can be easily deduced from [START_REF] Raous | to appear, A consistent model coupling adhesion, friction and unilateral contact[END_REF]. By making the subdifferentials ∂ u N and ∂ β explicit in ( 107) and (109), we obtain the following relations on the normal components [START_REF] Strömberg | A generalized standard model for contact, friction and wear[END_REF] and on the thermodynamic force G β [START_REF] Tabor | Friction-The present state of our understanding[END_REF].

R r T = C T u T β 2 (110) u N ≥ 0 -R r N + C N u N β 2 ≥ 0 -R r N + C N u N β 2 u N = 0 (111)              G β ≥ w h ′ (β) if β = 0 G β = w h ′ (β) -(C N u 2 N + C T u T 2 ) β if β ∈ ]0, 1[ G β ≤ w h ′ (β) -(C N u 2 N + C T u T 2 ) if β = 1. (112)
The state laws show that :

-the reversible (elastic) part of the tangential force depends on the square of the adhesion β,

-the reversible parts of the normal components of R and u satisfy a generalized Signorini condition (unilateral contact),

-the thermodynamic force G β , if β ∈]0, 1[, is composed of the adhesive energy minus the elastic energy of the interface.

Dissipation and complementary laws

We choose the following form for the pseudo-potential of dissipation which can be shown to agree with the Clausius Duhem inequality :

Φ uT , β; χ N = µ R N -C N u N β 2 uT + b p + 1 β p+1 + I C -( β) (113)
with C -= {γ ∈ W/γ ≤ 0} and p ≤ 1. A power law is used to account for the evolution of the adhesion. The indicator function imposes that β ≤ 0 : the adhesion is allowed only to decrease and cannot be regenerated (it is not reversible). This is in agreement with our application, but other choices could be made for other situations.

The complementary laws are then written :

R ir N = 0 (114) (R ir T , G β ) ∈ ∂Φ( uT , β; χ N ) (115)
and ( 115) can be written as :

R ir T ∈ ∂ uT Φ( uT , β; χ N ) (116)
G β ∈ ∂ β Φ( uT , β; χ N ) (117)
The normal behavior has been assumed to be elastic. Making explicit the subdifferentials in [START_REF] Vola | to appear, Bifurcations and instabilities in some finite dimensional frictional contact problem[END_REF] and [START_REF] Vola | Frottement et instabilités en dynamique : bruits de crissement[END_REF], we obtain on Γ C :

R T -C T u T β 2 ≤ µ R N -C N u N β 2 (118)
with :

R T -C T u T β 2 < µ R N -C N u N β 2 ⇒ uT = 0 R T -C T u T β 2 = µ R N -C N u N β 2 ⇒ ∃ λ ≥ 0 , uT = λ(R T -C T u T β 2 ) β = -G - β /b 1/p , (119)
where G - β denotes the negative part of G β .

These relations give rise to the model given in section 4.2.

The variational formulation

The complete problem can be then written as follows :

Problem P Adh . Find the displacements u α , the stresses σ α (α = 1, 2), the strains ǫ, and the contact force R such that :

ǫ = grad s u α in Ω α (120) div σ α + f α = 0 in Ω α (121) σ α n α = φ α on Γ α φ (122) u α = 0 on Γ α U (123) u = u 2 -u 1 = u N n + u T on Γ C (124) σ 1 n 1 = -σ 2 n 2 = R N n + R T on Γ C (125) σ α = K α : ǫ in Ω α (126)
and on Γ C :

-R N + C N u N β 2 ≥ 0 , u N ≥ 0 , -R N + C N u N β 2 u N = 0 (127) R r T = C T u T β 2 R T -R r T ≤ µ R N -C N u N β 2 (128) R T -R r T < µ R N -C N u N β 2 ⇒ uT = 0 R T -R r T = µ R N -C N u N β 2 ⇒ ∃λ ≥ 0 , uT = λ(R T -R r T ) , β = -(1/b) w h ′ (β) -(C N u 2 N + C T u T 2 )β -1/p if β ∈ [0, 1[ , (129) β ≤ -(1/b) w h ′ (β) -(C N u 2 N + C T u T 2 )β -1/p if β = 1 .
In what follows, we make the simplest hypothesis of a linear dissipation for the adhesion evolution (case p = 1).

The variational formulation is similar to that presented in chapter 3 for the quasistatic problem. Nevertheless, the functional j associated with friction depends also on β, and the problem is now coupled with a differential equation. We obtain the coupling of an implicit variational inequality, a "variational inequality" and a differential equation. Details can be found in [START_REF] Raous | to appear, A consistent model coupling adhesion, friction and unilateral contact[END_REF].

We denote

V α = v α ∈ [H 1 (Ω α )] d ; v α = 0 a.e. on Γ α U , (α = 1, 2), V = V 1 × V 2 , and K = {v = (v 1 , v 2 ) ∈ V ; v N ≥ 0 a.e. on Γ C }.
We shall adopt the following notations :

• a : V × V -→ IR, a(u, v) = a 1 (u 1 , v 1 ) + a 2 (u 2 , v 2 ) ∀u = (u 1 , u 2 ), v = (v 1 , v 2 ) ∈ V, where a α (u α , v α ) = Ω α A α ijkl ǫ ij (u α )ǫ kl (v α )dx, α = 1, 2 . • j : H × V × V -→ IR , j(β, u, v) = Γ C µ|R N (P u 1 ) + C N β 2 u N | v T ds, • c N , c T : H × V × V -→ IR , c N (β, u, v) = Γ C C N β 2 u N v N ds and c T (β, u, v) = Γ C C T β 2 u T .v T ds. • y(β, u) = - 1 b w -(C N u 2 N + C T u T 2 )β -. • (F, v) = α=1,2 Ω α f α .v α dx + Γ α φ φ α .v α ds ∀v ∈ V.
The local problem P Adh then admits the following variational formulation (see [START_REF] Raous | to appear, A consistent model coupling adhesion, friction and unilateral contact[END_REF] [28] [107] [START_REF] Cangémi | Frottement et adhérence : modèle, traitement numérique et application à l'interface fibre/matrice[END_REF]).

Problem P V ar . Find (u, β) ∈ W 1,2 (0, T ; V ) × W 1,2 (0, T ; H) such that u(0) = u 0 , β(0) = β 0 and for almost all t ∈ [0, T ], u(t) ∈ K and

a(u, v -u) + j(β, u, v) -j(β, u, u) + c T (β, u, v -u) ≥ (F, v -u) + R N (u 1 ), v N -uN ∀ v ∈ V (130) R N (u 1 ), z N -u N + c N (β, u, z -u) ≥ 0 ∀ z ∈ K, (131) β
= y(β, u) a.e. on Γ C , (132)
where the initial conditions u 0 ∈ K, β 0 ∈ H, β 0 ∈ [0, 1[ a.e. on Γ C and satisfy the following compatibility condition :

a(u 0 , w -u 0 ) + j(β 0 , u 0 , w -u 0 ) + c T (β 0 , u 0 , w -u 0 ) ≥ (F (0), w -u 0 ) ∀w ∈ K.
(133)

Incremental formulation and mathematical results

As done in Chapter 3 for the quasistatic problem, an incremental formulation is obtained by operating a time discretization of problem P V ar , taking n ∈ IN * and setting ∆t = T /n, t i = i ∆t and F i = F (t i ) for i = 0, ..., n. For the differential equation, we use an implicit scheme. We obtain the following sequence of problems (P n i ), i = 0, ..., n-1, defined for a given (u 0 , β 0 ) ∈ K × H . Problem P n i . Find (u i+1 , β i+1 ) ∈ K × H such that : a(u i+1 , v -u i+1 ) + j(β i+1 , u i+1 , v -u i ) -j(β i+1 , u i+1 , u i+1 -u i ) +c(β i+1 , u i+1 , v -u i+1 ) ≥ (F i+1 , v -u i+1 ) ∀v ∈ K (134)

β i+1 -β i = ∆t y(β i+1 , u i+1 ) a.e. on Γ C , (135)
where c(•) = c N (•) + c T (•).

The equation( 135) is solved using a fixed point method on β i+1 where the problem (134) has to be solved at each step. The existence and uniqueness are established for the incremental solution if friction is small in [START_REF] Cocu | Approximation results for a class of quasistatic contact problem including adhesion and friction[END_REF]. The proof is based on the work on the quasistatic problem presented in Chapter 2 [START_REF] Cocu | Formulation and approximation of quasistatic frictional contact[END_REF]. Because of the use of an implicit integration of the differential equation, a fixed point has been introduced to treat equation (135) on β k+1 . The convergence of this fixed point method has been proved in [START_REF] Cocu | Approximation results for a class of quasistatic contact problem including adhesion and friction[END_REF].

As in the case of the quasistatic problem, it is possible to construct a piecewise polynomial function with the incremental solution which converges weakly towards a continuous solution. This gives the existence for the solution of the quasistatic problem on condition µ is small.

Numerical methods and algorithm

As seen before, the problem is very similar to the quasistatic problem treated in Chapter 2. The main difference concerns the implicit integration of the differential equation. An implicit Euler method was introduced in the previous section. The implicit or semiimplicit character of the integration is important because the intensity of adhesion β may decrease very fast during the loading. In the computations, a θ-method will be used, which improves the accuracy of the solution. For the sake of simplicity, only the Euler method is presented here.

In comparison with Problem (P inc1 P inc1 P inc1 ) in Section 2.3.2, the main differences are : -at each time step, successive approximation iterates have to be conducted on β k+1 , -extra terms have to be included in the stiffness matrix (because of the contact elasticity), -the friction term j(β, u, v) depends on β.

After introducing the fixed point on the friction threshold g, we obtain the following discrete minimization problem (instead of Problem (P discret P discret P discret ) of Section 3.2.3.

Problem P h i . Find u ∈ k such that :

J (u) ≤ J (v) ∀v ∈ k (136) with J (v) = 1 2 v T Av + G T v -u i h + 1 2 v T C(β)v -F i+1 h T v (137)
where :

-k = {ΠK i with Any of the methods presented in Chapter 3 can be used to solve this problem. Schematically, the algorithm involving the various procedures can be presented as follows.

K i = R + if i ∈ I N and K i = R if not } -I N is
                                                                
Loop : on time steps t i+1 : f ind (u i+1 , β i+1 u i+1 , β i+1 u i+1 , β i+1 )

                                                
Loop : on the fixed point : (β i+1 ) ν (β i+1 ) ν (β i+1 ) ν (implicit integration) :

(initial condition (β i+1 ) 0 = β i )

• Solve the unilateral contact problem with friction (depending on (β i+1 ) ν (β i+1 ) ν (β i+1 ) ν )using either the Lemke algorithm or the following :

             
Loop : on the friction threshold ((g i+1 ) ν ) α+1

• ((g i+1 ) ν ) α+1 = µR N (((g i+1 ) ν ) α )

• Resolution of ProblemP h i by using SSORP or Gauss-Seidel+Aitken or Conjugate Gradient. This exemple is part of an undergoing research on the interaction of a progressive crack in the matrix with a fiber/matrix interface for a SiC/SiC composite. Here, the model is used to simulate a microindentation experiment in order to identifiy the parameters of the model. These experiments are conducted in the ONERA (Office National des Etudes et de la Recherche Aérospatiale). The results given by this model and those given by some other recent interface models are compared in the case of an uniaxial modelisation of the microindentation experiment in [START_REF] Monerie | Comparaison de lois d'interface fibre/matrice sur la base d'un modèle simplifié d'expérience de microindentation[END_REF]].

The model of the experiment is given on Fig. 27 and the mesh on Fig. 28 (1419 nodes, 300 contact nodes, triangle P 1 ). The radius of the fiber is R f = 8.6µm. By analyzing the volume density of fibers (V f = 40%), the radius R m = 12.6µm is chosen to define a zone of pure matrix. The extra zone is defined as an equivalent homogeneous material equivalent to the composite one by using a mixture rule. The domain is defined by L = 200µm and h = 1600µm. The boundary conditions are given on Fig. 27. The elasticity coefficients are : E=200GPA, ν=0.25 for the fiber, E=350GPa, ν=0.2 for the matrix and E=290GPa; ν=0.22 for the equivalent composite. The characteristic length l is taken equal to the fiber radius : l = 8.6µm. A vertical prescribed displacement is applied on the indentor. Results of the identification of the model parameters are given on Table 12 and on Fig. 30. The same contact stiffness was chosen for the normal and the tangential interface behavior : C N = C T = C. Figure 30 shows the ability of the model to fit the force/displacement evolution of the push-in experiment. Mechanical considerations give a specified range of variations for each parameter. For the contact stiffness, the range is evaluated by considering the composition and the thickness of the layer of carbon and oxyde in the interface. On Fig. 30, we have also plotted the results of the simulation by making µ = 0 (no friction) and w = 0 (no adhesion) in the values of the parameters given on Table 12. This is helpful to better understand the experimental results : the change of the slope is strongly related to the lost of adhesion and the final slope to the residual friction between the fiber and the matrix when adhesion is broken.

On Fig. 31, the contact condition is presented for a given step of the loading. The val- Force on the indentor : F (Newton)

. experiment numerical : adhesion alone numerical : friction alone numerical : adhesion and friction ues of the ratio R T /R N (scale on the left), the adhesion intensity β, and the tangential sliding u T (scales on the right) are plotted along the interface. Three zones are clearly characterized :

• zone 1 : close to the indentor, adhesion is totally broken (β ≃ 0) and only friction remains active (R T /R N = µ),

• zone 2 : in this transition zone, adhesion is partial (0 < β < 1) and the elasticity acts (R T /R N > µ),

• zone 3 : on this part, the interface is still weakly affected, the adhesion is total (β = 1). The smooth evolution of the solution underlines the good property of the model, which gives a continuous transition from total adhesion to pure frictional contact.

X

Figure 2 :

 2 Figure 2: Asperity behaviour during compression

Figure 3 :

 3 Figure 3: Complexity of the composition of the material in the vicinity of the surface

2 Figure 4 :

 24 Figure 4: Contact with a rigid obstacle

Figure 5 :

 5 Figure 5: Graph of the Signorini law

Figure 6 :

 6 Figure 6: Graph of the compliance law

Figure 7 :

 7 Figure 7: Graph of the Coulomb law

-Figure 8 :

 8 Figure 8: Three regularizations of the Coulomb law

Figure 9 :

 9 Figure 9: Graphs of various friction laws for 2D problems

Figure 11 :

 11 Figure 11: Analytical model of Archard ([Archard, 1974])

Figure 12 :

 12 Figure 12: Computed normal force versus the the decrease in height of the asperity

Figure 13 :

 13 Figure 13: Shearing of two asperities

2. 3 . 1

 31 Variational formulation : coupling a variational inequality to an implicit variational inequality[25] [26] [27][START_REF] Telega | Quasi-static Signorini's contact problem with friction and duality[END_REF]

2. 4 . 4

 44 The static unilateral problem with Coulomb friction a -The static problem with local friction Since the contact forces are defined only as a linear form, no results have been obtained in the usual functional framework. Rather than the difficulty of defining | F N |, it is the lack of compactness which constitutes the main obstacle.By introducing non classical functional spaces, Necas-Jarusek-Haslinger have established some results on simple geometries in[Necas et al, 1980], and in a more general case in[Jarusek, 1983].

  Figure 15: Geometry

Figure 16 :

 16 Figure 16: Mesh of the dovetail (4345 nodes, 49 contact nodes)

  [START_REF] Raous | Stress waves in a sliding contact -Part 2 : Modelling[END_REF] [[START_REF] Barbarin | Instabilité et frottement en élasticité -Application à un problème d'ondes de contrainte[END_REF] [Martins et al, in preparation] [Vola et al, in press]. The mesh is given in Fig.18 (4193 nodes, 65 contact nodes). The problem is quasistatic and 200 time steps have been used to closely follow the evolution in the contact conditions.

Figure 17 :Figure 18 :

 1718 Figure 17: Experimental assembly scheme

Figure 19 :

 19 Figure 19: Geometry and mesh of the closure of the pressure vessel (674 nodes, 54 contact nodes)

F

Figure 20 :

 20 Figure 20: Geometry and mesh of the screw-press

Figure 23 :

 23 Figure 23: Evolution of the contact stresses and of the contact conditions along the contact zone.

  the set of the number of degrees of freedom concerning the normal components of the contact nodes, -A is the matrix of dimension N = dim(V ) : A ij = a(w i , w j ), -C is the diagonal matrix of dimension M (M is the number of contact nodes) :C kl = c(β h , w k , w l ),-G is the vector of dimension M : G j = Γ C g w j ds .

End• 4 . 7

 47 Implicit integration of differential equation on β β β :(β i+1 ) ν+1 = β i + ∆tf (β i+1 ) ν+1 , (u i+1 ) ν (u i+1 ) ν (u i+1 ) ν )EndEnd Application to the behaviour of the fiber/matrix interface for a composite

FibreFigure 27 :

 27 Figure 27: Model of the experiment.

Figure 30 :

 30 Figure 30: Simulation and experimental results.

Figure 31 :

 31 Figure 31: Contact forces, adhesion intensity β and tangential displacement u T along the interface for a given loading step.

Table 1 :

 1 Comparison of the CPU time for the dovetail 3.5.2 The block sliding on a plane

	Method	Nber of iterates CPU time
	Lemke	-	32s
	SSORP ω = 1.95	2543	45s
	AITKEN	1726	35s
	Conjugate Gradient	422	25s

Table 2 :

 2 Comparison of the CPU times for the block (200 time steps) 3.5.3 The pressure vessel assembly

	Method	Nber of iterates CPU time
	Lemke	-	3s
	SSORP ω = 1.85	342	8s
	GS+Aitken	373	8s
	Conjugate Gradient	425	9s

Table 3 :

 3 Comparison between the CPU times for the pressure vessel 3.5.4 The high pressure screw-press

	Method	Nber of iterates CPU time
	Lemke	-	-
	SSORP ω = 1.9	1273	86s
	GS+Aitken	998	76s
	Conjugate Gradient	541	109s

Table 5 :

 5 CPU times for the indentation test3.5.6 Extrusion of an aluminium cylinder

	Method	GS+Aitken Lemke Augm. Lag. Penalization
	CPU time	8mn30s	2mn	6mn30s	1mn45s

Table 6 :

 6 CPU times for the extrusion of the cylinder

Table 8 :

 8 Parameter identification

		Difference
	Deflection	1%
	Upper rotation	5%
	Sliding upper part/plate	5%
	Thrust displacement	14%

Table 9 :

 9 Theory/experiment difference for the 1 st loading step (identification)

F N and F T /µ are plotted. Three zones can be observed : no contact (F N = F T = 0), sliding (the zone where the two curves are superposed), and adhesion (the zone where the curves are separated). It is possible to see how the sliding begins : a separated part occurs on the right, and a central sliding region extends progressively until the steady sliding solution is reached. Unilateral conditions with adhesion

Coulomb friction with adhesion

Evolution of adhesion intensity

The contact variables are u N , u T , β, R N , R T . The six parameters involved in the model are :

-C N , C T , the initial normal and tangential stiffness of the interface if the adhesion is complete, -µ is the friction coefficient, -b the viscosity of the adhesion evolution, -w the limit of decohesion energy (Dupré energy), -p a power coefficient for β evolution, -h(β) is given (in the fiber/matrix example, h(β) is taken to be equal to β.

Let us now analyze the interface behavior for a 2D case. The initial conditions are assumed to be complete adhesion (β = 1) and zero displacement (u T = u N = 0). Dealing first with the normal behavior (see Fig. 25), under compressive action, the non penetration condition is strictly fulfilled (u N = 0). Under traction (u N ≥ 0), an adhesive resistance (R N = C N u N β 2 ) is active (elasticity with damage). The intensity of adhesion starts to decrease when the displacement is sufficiently large for the elastic energy to become larger than the limit of the adhesion energy w. The evolution of the adhesion is then governed by equation [START_REF] Raous | Friction induced instabilities and sound generation[END_REF]. When the adhesion is completely broken, the classical Signorini problem is obtained.

Looking now at the shear behavior (see Fig. 26), note first that friction acts only if a normal compression is applied ; if a normal traction is applied (u N > 0), the sliding limit (µ the sliding limit, no sliding occurs (u T = 0 as initial condition and uT = 0 in relation ( 104)). When the sliding limit is reached, an elastic tangential displacement occurs. The adhesion begins to decrease when the adhesive limit is reached and the evolution of β is then governed by [START_REF] Raous | Friction induced instabilities and sound generation[END_REF]. When adhesion is lost (β tends toward zero), the usual Coulomb friction conditions are obtained.

Zoom