Serge Dumont
email: dumont@lma.cnrs-mrs.fr

Frédéric Lebon
email: lebon@lma.cnrs-mrs.fr

Raffaella Rizzoni
email: rizzoni.raffaella@unife.it

An asymptotic approach

Introduction

Adhesive bonding is an assembly technique often used in structural mechanics. In bonded composite structures, the thickness of the glue is much smaller than the other dimensions. For example, in [START_REF] Goglio | Design of adhesive joints based on peak elastic stresses[END_REF], the glue thickness is 0.1 mm, whereas the dimension of the structure is about 150 mm, giving a typical dimensional ratio close to 1/1500. Thus, the glue thickness can be considered as a small parameter in the modeling process. Usually, the glue stiffness is taken as another small parameter when compared with the adherents stiffness (soft interface theory), as shown in [START_REF] Lebon | Analysis of non-linear soft thin interfaces[END_REF] and [START_REF] Lebon | Asymptotic study of a soft thin layer: the non convex case[END_REF]. For example, in [START_REF] Goglio | Design of adhesive joints based on peak elastic stresses[END_REF], two steel structures are bonded by a Loctite 300 glue and the ratio between the Young moduli of the materials is close to 1/230. Nevertheless, in the case of an epoxy based adhesive bonding of two aluminium structures, the ratio between the Young moduli is typically about 1/20 (see for example [START_REF] Cognard | Comparative study of the results of various experimental tests used for the analysis of the mechanical behaviour of an adhesive in a bonded joint[END_REF]). Thus, the glue stiffness cannot be considered as the smallest parameter (stiff interface theory). The aim of this paper is to analyze mathematically and numerically the asymptotic behavior of bonded structures in the case of only one small parameter: the thickness. In the following, the stiffness is not a small parameter, and the Young moduli of the glue and of the adherents are of the same order of magnitude.

The mechanical behavior of thin films between elastic adherents was studied by several authors: [START_REF] Abdelmoula | Comportement asymptotique d'une interface mince[END_REF], [START_REF] Benveniste | An O(h n ) interface model of a three-dimensional curved interphase in conduction phenomena[END_REF], Bertoldi et al. (2007a), Bertoldi et al. (2007b), [START_REF] Bigoni | Statics and dynamics of structural interfaces in elasticity[END_REF], [START_REF] Cognard | Comparative study of the results of various experimental tests used for the analysis of the mechanical behaviour of an adhesive in a bonded joint[END_REF], [START_REF] Duong | A layerwise finite element for multilayers with imperfect interfaces[END_REF], [START_REF] Goglio | Design of adhesive joints based on peak elastic stresses[END_REF], [START_REF] Hirschberge | Computational multiscale modelling of heterogeneous material layers[END_REF], [START_REF] Krasuki | Yield design of bonded joints[END_REF], [START_REF] Kumar | Methods for solving singular perturbation problems arising in science and engineering[END_REF], [START_REF] Lebon | Analysis of non-linear soft thin interfaces[END_REF], [START_REF] Lebon | Asymptotic study of a soft thin layer: the non convex case[END_REF], [START_REF] Lebon | Asymptotic analysis of a thin interface: the case involving similar rigidity[END_REF], Lebon and Rizzoni (2011a), Lebon and Rizzoni (2011), [START_REF] Lebon | First-order numerical analysis of linear thin layer[END_REF], [START_REF] Nguyen | Benchmark of a layerwise stress model for laminated and sandwich plates [Benchmark d'un modèle pour sandwiches et multicouches, de type layerwise en contrainte[END_REF], [START_REF] Rizzoni | Asymptotic analysis of an elastic thin interphase with mismatch strain[END_REF], and [START_REF] Sacco | A damage-friction interface model derived from micromechanical approach[END_REF]. The analysis was based on the classic idea that a thin adhesive film can be replaced by a contact law, like in [START_REF] Abdelmoula | Comportement asymptotique d'une interface mince[END_REF]. The contact law describes the asymptotic behavior of the film in the limit as its thickness goes to zero and it prescribes the jumps in the displacement (or in the displacement rate) and in the traction vector fields at the limit interface. The limit problem formulation involves the mechanical and the geometrical properties of the adhesive and the adherents, and in [START_REF] Lebon | Analysis of non-linear soft thin interfaces[END_REF], [START_REF] Lebon | Asymptotic study of a soft thin layer: the non convex case[END_REF], [START_REF] Lebon | Asymptotic analysis of a thin interface: the case involving similar rigidity[END_REF], Lebon and Rizzoni (2011a), Lebon and Rizzoni (2011), [START_REF] Rizzoni | Asymptotic analysis of an elastic thin interphase with mismatch strain[END_REF], and [START_REF] Lebon | First-order numerical analysis of linear thin layer[END_REF] several cases were considered: soft films [START_REF] Klarbring | Derivation of the adhesively bonded joints by the asymptotic expansion method[END_REF] and [START_REF] Lebon | Analysis of non-linear soft thin interfaces[END_REF]); adhesive films governed by a non-convex energy [START_REF] Lebon | Numerical study of soft adhesively bonded joints in finite elasticity[END_REF], Lebon and[START_REF] Lebon | Asymptotic study of a soft thin layer: the non convex case[END_REF]Licht and[START_REF] Licht | A modeling of elastic adhesive bonded joints[END_REF]; imperfect gluing 2002)); flat linear elastic films having stiffness comparable with that of the adherents and giving rise to imperfect adhesion between the films and the adherents [START_REF] Lebon | Asymptotic analysis of a thin interface: the case involving similar rigidity[END_REF] and Lebon and Rizzoni (2011a)); joints with mismatch strain between the adhesive and the adherents, see for example [START_REF] Rizzoni | Asymptotic analysis of an elastic thin interphase with mismatch strain[END_REF]. Several mathematical techniques can be used to perform the asymptotic analysis: -convergence, variational analysis, matched asymptotic expansions and numerical studies (see Lebon andRizzoni (2011) and[START_REF] Sánchez-Palencia | Non-homogeneous Media and Vibration Theory[END_REF] and references therein).

The first part of the paper extends the imperfect interface law given in Lebon and Rizzoni (2011a) to the case of a very thin interphase whose stiffness is of the same order of magnitude as that of the adherents, firstly when an elastic body is glued to a rigid base, and secondly in the plane strain case. In the second part of the paper, numerical methods adapted to solve the limit problems obtained in the first part are developed. In the case of the gluing of a deformable body with a rigid solid, the numerical scheme is very classical. On the contrary, the gluing of two deformable bodies leads to more complicated numerical strategies. The proposed method is based on an original method presented in [START_REF] Nitsche | Convergence of nonconforming methods (finite element solution for Dirichlet problem)[END_REF]. This kind of method is well known in the domain decomposition context. This method is implemented in a finite element software. In the third part, some numerical examples are presented and the numerical results are analyzed (in terms of mechanical interpretation, computed time, convergence, etc.) in order to quantify and justify the methodology.

Theoretical results for thin stiff films

Asymptotic analysis for an elastic body glued to a rigid base

Let us consider a linear elastic body ⊂ IR 3 of boundary ∂ . This structure is made of two parts (the adherents) perfectly bonded with a very thin third one (the glue or the interphase), see Fig. 1. Initially, one of the two adherents is considered as rigid. After introducing a small parameter ε > 0 denoting thickness of the glue, we define the following domains:

• B ε = {(x 1 , x 2 , x 3 ) ∈ : 0 < x 3 < ε} (the glue); • ε + = {(x 1 , x 2 , x 3 ) ∈ : x 3 > ε} (the deformable adherent); • S ε + = {(x 1 , x 2 , x 3 ) ∈ : x 3 = ε}; • = {(x 1 , x 2 , x 3 ) ∈ : x 3 = 0} (the interface); • B = {(x 1 , x 2 , x 3 ) ∈ : 0 < x 3 < 1}; • + = {(x 1 , x 2 , x 3 ) ∈ : x 3 > 1}; • S + = {(x 1 , x 2 , x 3 ) ∈ : x 3 = 1}; • 0 + = {(x 1 , x 2 , x 3 ) ∈ : x 3 > 0}.
On a part 1 of ∂ , an external load g is applied, and on 0 ⊂ ∂ , such that 0 ∩ 1 = ∅, a displacement u d is imposed. Moreover, we suppose that 0 ∩ B ε = ∅ and 1 ∩ B ε = ∅. A body force f is applied in ε + . We consider also that the interface is a plane normal to the third direction e 3 . We are interested in the equilibrium of such a structure.

The equations of the problem are

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ div ε + f = 0 in ε + ∪ B ε ε n = g on 1 u ε = u d on 0 u ε = 0 on ε = A + e(u ε ) in ε + ε = Âe(u ε ) in B ε
(1) where e(u ε ) is the strain tensor (e ij (u ε ) = 1 2 (u i,j + u j,i ), i, j = 1, 2, 3) and A + , Â are the elasticity tensors of the deformable adherent and the adhesive, respectively. In the sequel, we consider that the glue is isotropic, with Lamé's coefficients equal to ˆ and ˆ in the interphase B ε . Let us emphasize that the Lamé's coefficients of the interphase do not depend on the thickness ε of the interphase (this will be referred as the case of a stiff interface hereinafter). Being thickness of the interphase very small, we seek the solution of problem (1) using asymptotic expansions with respect to the parameter ε:

u ε = u 0 + εu 1 + o(ε) ε = 0 + ε 1 + o(ε) (2)
In order to write the equations verified by u 0 , u 1 , 0 , 1 in + and on the interface , we consider the method developed by Lebon and Rizzoni (2011a) and based on the mechanical energy of the system:

J ε (u ε ) = 1 2 ε + A + e(u ε ) • e(u ε )dx - ε + f • u ε - 1 g • u ε ds + 1 2 B ε Âe(u ε ) • e(u ε )dx (3)
which is defined in the set of displacements

V ε = {u ∈ H( ; IR 3 ) : u = u d on 0 , u = 0 on }, (4)
where H( ; R 3 ) is the set of admissible displacements defined on (Fig. 2). At this level, the domain is rescaled using the following classical procedure:

• In the glue, we define the following change of variable:

(x 1 , x 2 , x 3 ) ∈ B ε → (z 1 , z 2 , z 3 ) ∈ B, with (z 1 , z 2 , z 3 ) = (x 1 , x 2 , x 3 ε )
and we denote ûε (z 1 , z 2 , z 3 ) = u ε (x 1 , x 2 , x 3 ). • In the adherent, we define the following change of variable:

(x 1 , x 2 , x 3 ) ∈ ε + → (z 1 , z 2 , z 3 ) ∈ + , with (z 1 , z 2 , z 3 ) = (x 1 , x 2 , x 3 + 1 -ε)
and we denote u ε (z 1 , z 2 , z 3 ) = u ε (x 1 , x 2 , x 3 ). We suppose that the external forces and the prescribed displacement u d are assumed to be independent of ε. As a consequence, we define

f (z 1 , z 2 , z 3 ) = f (x 1 , x 2 , x 3 ), g(z 1 , z 2 , z 3 ) = g(x 1 , x 2 , x 3 ) and u d (z 1 , z 2 , z 3 ) = u d (x 1 , x 2 , x 3 ).
Then, using these notations, the rescaled energy takes the form

J ε ( ûε , u ε ) = + 1 2 A + (e(u ε )) • e(u ε ) -f • u ε dz - 1 g • u ε dS + B 1 2 (ε -1 K33 ( ûε ,3 ) • ûε ,3 + 2 K˛3 ( ûε ,˛) • ûε ,3 + ε K˛ˇ( ûε ,˛) • ûε ,ˇ) dz (5)
where a comma is used to denote partial differentiation, ˛, ˇ∈ {1, 2} and K jl , j, l = 1, 2, 3, are the matrices whose components are defined by the relations

( Kjl ) ki := Âijkl . ( 6
)
In view of the symmetry properties of the elasticity tensor Â, the matrices Kjl have the property that Kjl = ( Klj ) T , j, l = 1, 2, 3. The rescaled equilibrium problem is formulated as follows: find the pair (u ε , ûε ) minimizing the energy (5) in the set of displacements

V = {(u ε , ûε ) ∈ H( + ; R 3 ) × H(B; R 3 ) : u ε = u d on 1 , u ε = ûε on S + , ûε = 0 on } (7)
where H( + ; R 3 ) and H(B ; R 3 ) are the sets of admissible displacements defined on + and B, respectively. We assume that the displacements minimizing J ε in V can be expressed as the sum of the series

ûε = û0 + ε û1 + ε 2 û2 + o(ε 2 ) , (8)
u ε = u 0 + εu 1 + ε 2 u 2 + o(ε 2 ) . (9)
Correspondingly, the rescaled energy ( 5) can be written as:

J ε ( ûε , u ε ) = 1 ε J -1 ( û0 ) + ε 0 J 0 ( û0 , u 0 , u 1 ) + εJ 1 ( û0 , u 0 , û1 , u 1 , û2 ) + ε 2 J 2 ( û0 , u 0 , û1 , u 1 , û2 , u 2 , û3 ) + o(ε 2 ), ( 10
)
where J -1 , J 0 , J 1 and J 2 will be written in the next paragraph. As proposed in [START_REF] Lebon | Asymptotic analysis of a thin interface: the case involving similar rigidity[END_REF], we now minimize successively the energies J -1 , J 0 , J 1 , and J 2 .

Minimization of J -1

The energy J -1 is minimized in the class of displacements

V -1 = {û 0 ∈ H(B; R 3 ) : û0 = 0 on }. ( 11
)
Because  is a positive definite tensor, the second order tensor K33 is also positive definite. Therefore, the energy J -1 is non-negative and the minimizers are such that û0

,3 = 0 a.e. in B . This result, together with the boundary condition in ( 11), yields û0 = 0, a.e. on S + .

(12)

Minimization of J 0

Based on ( 12), the energy J 0 turns out to depend only on u 0 and it takes the form:

J 0 = + 1 2 A + (e(u 0 )) • e(u 0 ) -f • u 0 dz - 1 g • u 0 dS. (13)
In view of (12) and of the continuity of the displacements at the surface Ŝ+ , we seek the energy minimizer in the class of displacements

V 0 = {u 0 ∈ H( + ; R 3 ) : u 0 = 0 on S + , u 0 = u d on 0 }. ( 14
)
Using standard arguments, we obtain the equilibrium equations div(A + (e(u 0 )))

+ f = 0in + , (15)
A + (e(u 0 ))n = g on 1 , (16)

A + (e(u 0 ))n = 0 on ∂ + \ ( 1 ∪ S + ).
(17)

Minimization of J 1

In view of ( 12), the energy J 1 simplifies as follows:

J 1 := + (A + (e(u 0 )) • e(u 1 ) -f • u 1 ) dz - 1 g • u 1 dS + B 1 2 K33 ( û1 ,3 ) • û1 ,3 dz. ( 18
)
We minimize this energy in the class of displacements

V 1 = {(u 1 , û1 ) ∈ H( + ; R 3 ) × H(B; R 3 ) : u 1 = û1 on S + , û1 = 0 on , u 1 = 0 on 0 }. (19)
Using (15-17), the Euler-Lagrange equations reduce to the following equation:

S+ (A + (e(u 0 ))n • Á 1 )dS + B ( K33 ( û1 ,3 ) • Á1 ,3 ) dz = 0, ( 20
)
where Á 1 , Á1 are perturbation of u 1 , û1 , respectively, and they are such that Á 1 = Á1 on S + . Integrating by parts the second integral, using the boundary conditions in ( 19) and the arbitrariness of Á 1 , Á1 , we obtain

û1 ,33 = 0 in B, (21)
û1 ,3 = -( K33 ) -1 (A + (e(u 0 ))n on S + , (22)
which, together with the boundary condition on , give

û1 = -(( K33 ) -1 (A + (e(u 0 ))n)z 3 in B. (23) 2.1.4. Minimization of J 2
Using the results obtained so far, the energy J 2 simplifies as follows:

J 2 = + 1 2 A + (e(u 1 )) • e(u 1 dz + B K˛3 ( û1 ,˛) • û1 ,3 dz (24)
In view of ( 23) the second integral is a constant term and it can be dropped in the minimization procedure. Thus, we minimize the remaining term in the energy in the class of displacements

V 2 = {u 1 ∈ H( + ; R 3 ) : u 1 = -(( K33 ) -1 (A + (e(u 0 ))n) on S + }. ( 25
)
and we obtain the equilibrium equations:

div(A + (e(u 1 ))) = 0 in + , (26)
A + (e(u 1 ))n = 0 on ∂ + \ S + .

(27)

Limit equilibrium problems

Due to the continuity of the displacement across the surface S + , we have

ûε (z 1 , z 2 , 1 -) = u ε (z 1 , z 2 , 1 + ) = u ε (x 1 , x 2 , ε). ( 28
)
Note that the same condition is obtained for the stress field along S + . Using an asymptotic expansion, we have

u ε (x 1 , x 2 , ε) = u ε (x 1 , x 2 , 0 + ) + εu ε ,3 (x 1 , x 2 , 0 + ) + o(ε).
Using relations (2) and the last two equations, we obtain

u 0 (x 1 , 0 + ) = u 0 (x 1 , 1 + ), ( 29
)
u 0 ,3 (x 1 , 0 + ) + u 1 (x 1 , 0 + ) = u 1 (x 1 , 1 + ). ( 30
)
With the help of the above relations, we can rewrite the interface conditions obtained in the asymptotic analysis in terms of the displacement in the deformable adherent. In summary, we have obtained the following "problem at the order zero":

(P 0 ) ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ div(A + (e(u 0 ))) + f = 0 in 0 + , A + (e(u 0 ))n = g on 1 , A + (e(u 0 ))n = 0 on ∂ 0 + \ ( 1 ∪ ) u 0 = 0 on , (31)
and the following "problem at the first order":

(P 1 ) ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ div(A + (e(u 1 ))) = 0 in 0 + , A + (e(u 1 ))n = 0 on ∂ 0 + \ u 1 = -(( K33 ) -1 (A + (e(u 0 ))n) -u 0 ,3
on .

(32)

Plane strain problem

We consider the case of plane strain in the plane (x 1 , x 2 ), where the interface between the glue and the adhesive is a line orthogonal to the direction e 2 . After adopting natural notations and implementing the same kind of technique used in the previous section, we reobtain problem P 0 as in (31) and problem P 1 modified as follows:

(P 1 ) ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ div(A + (e(u 1 ))) = 0 in 0 + , A + (e(u 1 ))n = 0 on ∂ 0 + \ u 1 = -(( K22 ) -1 (A + (e(u 0 ))n) -u 0 ,2
on .

(33)

Case of the gluing of two elastic bodies

Let us now consider the adhesive bonding of two linear elastic bodies satisfying the plane strain hypothesis. We extend the notation used before, and we define the following domains:

• B ε = {(x 1 , x 2 ) ∈ : |x 2 | < ε/2} (the glue); • ε ± = {(x 1 , x 2 ) ∈ : ±x 2 > ε/2}; • S ε ± = {(x 1 , x 2 ) ∈ : x 2 = ±ε/2}; • = {(x 1 , x 2 ) ∈ : x 2 = 0} (the interface); • B = {(x 1 , x 2 ) ∈ : |x 2 | < 1/2}; • ± = {(x 1 , x 2 ) ∈ : ± x 2 > 1/2}; • S ± = {(x 1 , x 2 ) ∈ : x ± = ±1/2}; • 0 ± = {(x 1 , x 2 ) ∈ : ±x 2 > 0}.
The methodology and the notations used here are similar to the ones used in the previous section. The main differences are:

• The introduction of a jump of the stress vector at order 0 and 1.

• The displacement along the interface is replaced by a jump of the displacement across the interface between the two bodies.

• The minimization of J 1 leads to concentrated forces at the edges of the interface.

More precisely, the problem at order 0 becomes

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ div 0 + f = 0 in 0 ± 0 n = g on 1 u 0 = u d on 0 0 = A ± e(u 0 ) in 0 ± [u 0 ] = 0 on [ 0 n] = 0 on (34)
where [f](x 1 ) = f(x 1 , 0 + )f(x 1 , 0 -). The problem at order 1 becomes (see [START_REF] Lebon | Asymptotic analysis of a thin interface: the case involving similar rigidity[END_REF] p. 479, here simplified for the 2D case with the plane strain hypothesis).

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ div 1 = 0 in 0 ± 1 n = 0 on 1 u 1 = 0 on 0 1 = A ± e(u 1 ) in 0 ± [u 1 ] = C 1 ( 0 n) + C 2 (u 0 ) ,1 - 1 2 (u 0 (x 1 , 0 + ) + u 0 (x 1 , 0 -)) on [ 1 n] = C 3 ( 0 n) ,1 + C 4 (u 0 ) ,11 on 1 e 1 = F on ∂ (35)
where

C 1 = ⎛ ⎜ ⎝ 1 ˆ 0 0 1 ˆ + 2 ˆ ⎞ ⎟ ⎠ , C 2 = ⎛ ⎝ 0 -1 - ˆ ˆ + 2 ˆ 0 ⎞ ⎠ , C 3 = ⎛ ⎝ 0 - ˆ ˆ + 2 ˆ -1 0 ⎞ ⎠ , C 4 = ⎛ ⎝ -4 ˆ ˆ + ˆ ˆ + 2 ˆ 0 0 0 ⎞ ⎠ ,
and the localized forces which appear on the lateral boundary of the thin layer are given by

F = C 3 ( 0 e 1 ) + C 4 (u 0 ) ,1 (36)
Let us notice that this term appears naturally in this method, and has been first observed by one another technique in [START_REF] Zaittouni | Etude théorique et numérique du comportement d'un assemblage de plaques[END_REF] (see also [START_REF] Lebon | First-order numerical analysis of linear thin layer[END_REF] and [START_REF] Lebon | Asymptotic modelling of interface taking into account contact conditions: asymptotic expansions and numerical implementation[END_REF]).

Numerical method

In this paragraph, we present the numerical method developed to solved problem (35). The generic problem associated to this problem can be written

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ div (u) = 0 in 0 ± (u)n = 0 on 1 u = 0 on 0 = A ± e(u) in 0 ± [u] = D on [ (u)n] = G on (37)
where D and G are given functions, provided by the solutions u 0 and 0 of problem (34) at order 0. Note the solution of problem at order 0 (34) is very classic and can be solved using a classical finite element method. In the following, we will denote the restriction of u on 0 + (resp. 0 -) by u + (resp. u -). In order to write a discretized approximation of this equation, we use a Discontinuous Galerkin approach. For that purpose, we first write the variational formulation of the four first equations of (37) both in 0 -and 0 + , that leads, after an integration by parts, to

0 ± A ± e(u ± ) • e(v ± )dx - ∂ 0 ± (A ± e(u ± ))n ± • v ± dS = 0 (38) for v ± ∈ {v ∈ H 1 ( ) : v = 0 on ∂ 0 }. Then, introducing the boundary conditions, we obtain 0 ± A ± e(u ± ) • e(v ± )dx - (A + e(u + ))n + • v + dS -(A -e(u -))n -• v -dS = 0. ( 39
)
We now choose the normal n equal to the outward normal of 0 -(n = n -=n + ), and we denote

I = (A + e(u + ))n • v + dS -(A -e(u -))n • v -dS.
(40)

Then, using the jump condition (A + e(u + ))n = (A -e(u -))n + G, we have

I = (A -e(u -))n • (v + -v -)dS + Gv + dS. ( 41
)
Similarly, writing the jump condition as (A -e(u -))n = (A + e(u + ))n -G, we also have

I = (A + e(u + ))n • (v + -v -)dS + Gv -dS. ( 42
)
In order to have a symmetric variational formulation, numerically more efficient, we consider the half sum of ( 41) and ( 42):

I = 1 2 (A + e(u + ))n + (A -e(u -))n • (v + -v -)dS + G v + + v - 2 dS. (43)
Again, in order to have a fully symmetric formulation, we need to add in the left-hand side of Eq. ( 38) the term

(u + -u -) • 1 2 (A + e(v + ))n + (A -e(v -))n dS (44)
and we use the fact that u +u -= D on . Finally, we have the weak formulation 0

+ ∪ 0 - A ± e(u ± ) • e(v ± )dx + ( Ae(u)n • [v] + [u] • Ae(v)n ) dS = - G • v dS + D • Ae(v)n dS, (45)
for all v ∈ {H 1 ( ) : (v) = 0 on ∂ \ }, where • denotes the average of the value of the function on the both sides of the interface :

f = 1 2 (f + + f -)
. This formulation, known as the Nitsche's method [START_REF] Nitsche | Convergence of nonconforming methods (finite element solution for Dirichlet problem)[END_REF] is not stable and the discrete operator can be non-invertible after a discretization. It is then necessary to add a stabilization term such as

ȟ [u] • [v]dS = ȟ D • [v]dS,
where h is the size of the smallest element of the finite element discretization of 0 ± considered, and ˇ> 0 is a fixed number that must be sufficiently large to ensure the stability of the method. It can be shown that this formulation is equivalent to (37). In particular, solutions of ( 45) are weak solutions of (37) (see [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF], [START_REF] Dumont | Meshfree methods and boundary conditions[END_REF][START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF] for the complete study of this method and for a priori and a posteriori error estimates in the case D = 0).

Let us notice that this method is formally equivalent to the use of Lagrange multipliers to enforce the jump conditions (see [START_REF] Barbosa | Boundary Lagrange multipliers in finite element methods: error analysis in natural norms[END_REF] and [START_REF] Becker | A finite element method for domain decomposition with non-matching grids[END_REF]), but it takes its advantage on the fact that the Nitsche's method does not increase the number of unknowns. Unfortunately, this method does not work properly to solve the problem (37) as soon as D / = 0. To overcome this difficulty, we split the problem (37) into two parts. More precisely, we write u ± = w ± + z ± where the unknowns z ± and w ± solve the problems

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ div (z ± ) = 0 in 0 ± (z ± )n = 0 on 1 z ± = 0 on 0 (z ± ) = A ± e(z ± ) in 0 ± z ± = ± 1 2 D on ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ div (w ± ) = 0 in 0 ± (w ± )n = 0 on 1 w ± = 0 on 0 (w ± ) = A ± e(w ± ) in 0 ± [w] = 0 on [ (w)n] = G -[ (z)n] on (46) since [w] = w + -w -= [u] -z + + z -= (1 -(1/2) -(1/2))D = 0.
The two first problems defined in the left of both in 0 + and 0 -are standard and can be solved simultaneously using a standard finite element method. The problem on the right of ( 46) is solved using the Nitsche's method developed above.

Numerical results for an elastic body glued to a rigid base

In this paragraph, we consider a 2D solid composed of an aluminum adherent and an epoxy resin interphase, that glues the structure to a rigid base. The mechanical coefficients of the materials are the following:

• In the adhesive (epoxy resin): Ê = 4 GPa, ˆ = 0.33. • In the adherent (aluminium): E = 70 GPa, = 0.33. This application was intentionally selected in order to introduce a significative difference between the elastic moduli of the adherent and those of the adhesive materials.

The geometry of the problem is provided in Fig. 3. The meshes are realized using the GMSH software developed by [START_REF] Geuzaine | Gmsh: a three-dimensional finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. The finite element computations are made with the MATLAB ® software.

The computations are first realized for the interphase problem with various values of the thickness of the interphase. The values of the jumps of the displacement and of the stress components across the interphase are then computed. Independently, the interface problem at order 0 (Eqs. ( 31)) is solved numerically. Then, the jumps across the interface at order 1 are computed using the corresponding equations in (32), and they are compared with the jumps across the interphase. We can notice that, to compute the jump in the using (32), one needs to numerically compute the derivative of the displacement on the boundary. In order to make the computation in a suitable way, it is necessary to use at least quadratic finite elements. In the numerical experiments below, we use the quadratic T6 finite elements.

Jump [u] across the interface

In Fig. 4, we present comparative plots of displacement amplitudes u 1 and u 2 across the interface for various values of the thickness ε. More precisely, since the displacement of the rigid base is vanishing, we compare the displacement u ε (x 1 = ε, x 2 ), denoted [u ε i ], i = 1, 2 in Fig. 4 and computed using the real geometry of the adhesive, with the displacement u(x 1 = 0 + , x 2 ) = u 0 (x 1 = 0 + , x 2 ) + εu 1 (x 1 = 0 + , x 2 ), denoted [u i ], i = 1, 2 and computed using Eqs. ( 31) and (32).

In Table 1, the relative error is computed using the L 2 norm: e = u εu L 2 ( ) / u L 2 ( ) , written as a percentage. In Fig. 4, we can observe that, as expected, when ε = 0.1 m the difference between the jumps across the interphase and the jumps across the interface is relatively large (with a relative error of about 12%). For smaller values of the interphase thickness ε, the difference between the results obtained using the interphase problem and those obtained using the interface approximation is negligible.

In Table 1, one can observe that the relative error depends also on the ratios of the rigidity of the glue and the adherent. When the ratio becomes larger, the approximated model is efficient only for smaller values of the thickness.

Jump [ 12 ] across the interface

In this paragraph, we present a comparison between the traction amplitude 12 at the top of the interphase/interface and at the bottom of the interphase/interface computed for the original interphase problem and for the approximated (at order 1) problem with the interface. They are respectively referred as up and bottom in Figs. ( 5) and ( 6). In the comparison, the thickness ε ranges from 0.01 m to 0.001 m.

Table 1

Percentage of the relative error in the jump in the displacement with respect to the relative rigidity Ê/E of the glue and the adherent and the thickness of the glue ε. We also present the traction amplitude computed at order 0 for the approximated problem. Since the traction is continuous across the interface, i.e.

[ 12 ] = 0, the traction amplitude takes the same value at the top and at the bottom of the interface.

The case ε = 0.1 m is not presented here because the difference between the case with the interphase and the case with the interface is large.

We can observe that the stress 12 computed using the interphase problem numerically converges to the stress 12 computed using the interface problem when ε tends to 0. One can also observe that the traction amplitude calculated at order 0 converges much slower than the traction amplitude calculated at order 1.

Time computing

In this paragraph, we present the time computing necessary to obtain the solutions of the problems considered in the previous section (Tables 2 and3).

Even if only linear finite elements are necessary for the computations for the interphase problem, we can notice that the CPU times necessary to obtain the solution quickly increases as the thickness of the interphase tends to 0. The reason is that the mesh has to be sufficiently fine inside the interphase, at least four nodes along the thickness. Therefore, in order to keep a reasonable condition number for the rigidity matrix, the mesh has to be fine also in a large zone around the interphase. This necessity significantly increases the number of degrees of freedom as the thickness of the interphase tends to zero. For the interface problem, the computation is relatively independent of the parameter ε. As a consequence, the meshes and the CPU times of the computations increase very slowly as the thickness tends to zero. So, for small values of the thickness, the use of the imperfect interface model is very convenient.

Numerical results for two elastic bodies glued

In this section, we present an example of two elastic structures glued together. This example which is academic is however inspired by [START_REF] Goglio | Design of adhesive joints based on peak elastic stresses[END_REF] and is composed of two T-form elastic bodies (aluminium) glued with an epoxy resin (see Fig. 7 for the geometry). More precisely, the mechanical coefficients of the materials are the same as before.

We present results for ε = 0.01 m. The results show that the interface law at order 1 is able to reproduce the phenomena that occur in the interphase much more accurately that the interface law at order 0. For example, we observe that 12 is in good agreement with ε 12 (see Fig. 8). The model at order 1 is able to reproduce the displacement in the jump and the displacement with a very small error (see Fig. 9 for u 1 and Fig. 10 for u 2 ).

Conclusion

In this paper, we have presented an interface law at order 1 when the Lamé's coefficients of the adhesive do not rescale with the thickness of the interphase. Based on the proposed interface law, some numerical experiments were also presented that show the accuracy of the method when the interphase thickness becomes smaller and smaller. For this purpose, we have developed an original numerical scheme based on the Nitsche's method to simulate the adhesion between two elastic materials.

This model is very efficient (relative error less than 2% for a ratio of thickness smaller than 1%) and the interface law is able to reproduce the mechanical behavior of the real interface. On the other hand, the numerical model developed in this paper is less expensive than the solution of the real problem. More precisely, the method is independent of the thickness of the interphase and it becomes more and more efficient as the thickness decreases. For example, in the first numerical test proposed above, the CPU times of the interphase problem and the asymptotic interface problem are equivalent when ε = 0.1, but their ratio is lower than 1/10 when ε = 0.001.

Fig. 1 .

 1 Fig. 1. (a) Initial, (b) rescaled, and (c) limit configuration of a solid glued to a rigid base.

Fig. 2 .

 2 Fig. 2. Geometry of the interphase/interface problem. Left: the initial problem with an interphase of thickness ε. Middle: the rescaled problem with interphase height equal to 1. Right: the limit interface problem.

Fig. 3 .

 3 Fig. 3. Geometry of the problem (ε = 0 for the interface problem).

Fig. 4 .

 4 Fig. 4. Example 1 -Jump in the displacement [u](x1, x2 = 0) (m) along the interface: ε = 0.1 m (left top), ε = 0.01 m (right top), ε = 0.005 m (left bottom) and ε = 0.001 m (right bottom).

Fig. 5 .

 5 Fig. 5. Example 1 (ε = 0.01 m) -Stress 12 (x1, x2 = 0) (GPa) at the bottom of the elastic adherent (up) and on the rigid base (bottom) computed with the various approximations (zoom on the right).

Fig. 6 .

 6 Fig. 6. Example 1 (ε = 0.001 m) -Stress 12 (x1, x2 = 0) (GPa) at the bottom of the elastic adherent (up) and on the rigid base (bottom) computed with the various approximations (zoom on the right).

Fig. 7 .

 7 Fig. 7. Example 2: geometry of the problem (ε = 0 for the interface problem).

Fig. 8 .

 8 Fig. 8. Example 2 (ε = 0.01 m) -Stress 12 (x1 = 2.5, x2) (GPa) on a vertical cut.

Fig. 9 .

 9 Fig. 9. Example 2 (ε = 0.01 m) -Displacement u1(x1 = 2.5, x2) (m) on a vertical cut (zoom on the right).

Fig. 10 .

 10 Fig. 10. Example 2 (ε = 0.01 m) -Displacement u2(x1 = 2.5, x2) (m) on a vertical cut (zoom on the right).

Table 2

 2 Mesh properties (T3 finite elements) and times computing for the interphase problem and various values of ε.

	Thickness ε (m)	Number of nodes	Number of elements	Number of degrees of freedom	CPU time (s)
	0.1	2582	4966	4962	6.4
	0.01	28,172	55,166	54,342	48
	0.005	39,486	77,432	76,304	66
	0.001	150,054	294,831	290,106	274

Table 3

 3 Mesh properties (T6 finite elements) and times computing for the interface problem and various values of ε.

	Thickness ε (m)	Number of nodes	Number of elements	Number of degrees of freedom	CPU time (s)
	0.1	731	1363	1410	5.2
	0.01	10,887	5338	21,372	22
	0.005	12,511	6138	24,620	24
	0.001	13,417	6586	26,332	25

Acknowledgement

RR thanks the financial support of the Italian Ministry of Education, University and Research (MIUR) through the PRIN project "Multi-scale modeling of materials and structures" (code 2009XWLFKW).


Lori Badea
email: lori.badea@imar.ro

Frédéric Lebon
email: lebon@lma.cnrs-mrs.fr

Hector Ramirez

Schwarz method for dual contact problems

Keywords: Contact problems, Dual formulation, Domain decomposition methods, Schwarz method, Subspace correction methods, Variational inequalities Mathematics Subject Classification 65N55, 65K15, 74M10, 74M15

est destinée au dépôt

Schwarz method for dual contact problems

Introduction

Traditionally, frictional contact problems are formulated in terms of the displacements (primal variables) [START_REF] Cocu | Existence of solutions of Signorini problems with friction[END_REF][START_REF] Raous | Numerical methods for frictional contact problems and applications[END_REF][START_REF] Cocou | Formulation and approximation of quasistatic frictional contact[END_REF][START_REF] Lebon | Contact problems with friction: models and simulations[END_REF] or of forces-displacements (mixed variables) [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF][START_REF] Alart | Solution of frictional contact problems using ILU and coarse/fine preconditioners[END_REF]Wriggers 2002). There also exists a third kind of formulation in terms of stresses (dual variables) [START_REF] Telega | Quasi-static Signorini's contact problem with friction and duality[END_REF]. This formulation has given numerical developments in the last few years [START_REF] Bisegna | D-PANA : a convergent block-relaxation solution method for the discretized dual formulation of the Signorini-Coulomb contact problem[END_REF]). This approach is little used but it has a lot of advantages in point of view of numerical computation. The first advantage is that the stress tensor is the variable of interest in the mechanical design and is directly obtained by this formulation. Another one is the capability of deriving an a posteriori error estimator [START_REF] Kuss | Error estimation and mesh adaptation for SignoriniCoulomb problems using E-FEM[END_REF]. Also, in the dual case, the problem is formulated as a variational inequality (of the first kind) which is more suitable for the numerical point of view than the variational inequality of the second kind in the primal case which contains a non-differentiable term.In addition, the numerical experiments have shown that the linear systems in dual variable are better conditioned than those in primal variables [START_REF] Bisegna | Relaxation procedures for solving Signorini-Coulomb contact problems[END_REF]). On the other hand, this approach is less used than the primal or mixed ones because of numerical difficulties related to the finite elements which have to be used [START_REF] Kuss | Stress based finite element methods for solving contact problems: comparisons between various solution methods[END_REF].

In the case of the Schwarz method, one of the difficulties is even the formulation of the method. In the primal case, the displacements are thought in H 1 and they have traces on the boundary of the domain, and the Schwarz method transfers, during the iteration, the solution values on a subdomain toward the solution on the neighboring subdomains. But, the dual approach utilizes the space H div where only the normal trace of the stresses exists. This difficulty has been avoided by interpreting the Schwarz method as a subspace correction method. In this case, we have to find corrections in appropriate subspaces which may not have any connection with the notion of trace.

We prove in this paper that, for the weak formulation of the dual contact problem with Tresca friction, the Schwarz method converges and also gives an error estimate depending on the domain decomposition parameter. It is well known that the solution of the contact problem with Coulomb friction can be obtained by a fixed point procedure in which we solve a contact problem with Tresca friction at each iteration step [START_REF] Licht | Remarks on a numerical method for unilateral contact including friction[END_REF]. The study of the method for the discretized problem by finite elements, the one-level or multi-level methods, will be treated in a subsequent paper.

The paper is organized as follows. In Sect. 2, the mechanical and the weak formulations of the problem in dual variables are introduced. Section 3 is dedicated to the study of the Schwarz algorithm. As we said above, using the primal variables, the problem is formulated as a variational inequality of the second kind which contains a non-differentiable term. Oneand two-level Schwarz methods for such formulations of the frictional contact problems have been given in [START_REF] Badea | One-and two-level Schwarz methods for inequalities of the second kind and their application to frictional contact[END_REF]. When the problem is formulated in dual variables, we get a simply variational inequality (of the first kind) and we shall use the results in [START_REF] Badea | Convergence rate of a multiplicative Schwarz method for strongly nonlinear inequalities[END_REF] to prove the convergence of the method. For an abstract theory of Schwarz methods we recommend [START_REF] Toselli | Domain decomposition methods-algorithms and theory[END_REF]. The algorithm is introduced as a subspace correction method and we rewrite a general convergence theorem from [START_REF] Badea | Convergence rate of a multiplicative Schwarz method for strongly nonlinear inequalities[END_REF] in terms of the considered framework. We prove that the assumption made in the convergence theorem is verified and derive the constants in the error estimation. Finally, summing up the previous results we show that the method is geometrically convergent with a convergence rate which depends on the domain decomposition parameter.

Statement of the problem

We consider an elastic body occupying the open bounded set ⊂ R d , (d = 2, 3) , with sufficiently smooth boundary. The boundary is divided in three open parts D , F and C such that ∂ = ¯ D ∪ ¯ F ∪ ¯ C and we assume that the measure of D is positive, meas ( D ) > 0. We denote u = (u i ) 1≤i≤d the displacement field, ε = (ε i j (u)) 1≤i, j≤d , ε i j (u) = 1 2 (u i, j + u j,i ), the linearized strain tensor and σ = (σ i j (u)) 1≤i, j≤d , σ i j (u) = d k,l=1 A i jkl ε kl (u) the Cauchy stress tensor. Operator A is bounded and verifies usual symmetry and positivity conditions. We denote S = A -1 the compliance tensor. This tensor has the following properties of symmetry, S i jkl = S kli j = S jilk , and of positivity, i.e. there exists α > 0 such that

d i, j=1 d k,l=1 S i jkl τ i j τ kl ≥ α d i, j=1 τ i j τ i j .
The body is in receding contact with a rigid obstacle on the part of the boundary C . The contact is with friction and is represented by the unilateral contact law of Tresca.

On the contact boundary, displacement and stress vector are decomposed introducing the external normal unit vector n to :

u N = u • n, u T = u -u N • n σ N = (σ n) • n, σ T = σ n -σ N n (2.1)
The body is submitted to a body force density f in and surface force density g on

F . We suppose f ∈ (L 2 ( )) d and g ∈ (L 2 ( F )) d .
On D the displacement is given.

Classically, the mechanical problem of an elastic body in Coulomb frictional contact with a rigid obstacle is given by: Problem (P) Find the displacement field u and the stress field σ such that

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ σ = Aε(u) in div σ + f = 0 i n u = u 0 on D σ • n = g on F u N ≤ 0, σ N ≤ 0, u N σ N = 0 on C |σ T | ≤ -μσ * N on C and |σ T | < -μσ * N ⇒ u T = 0 |σ T |=-μσ * N ⇒ ∃λ ≥ 0, u T = -λσ T (2.2)
where μ is the friction coefficient (bounded and positive) and * is the normal smoothness operator. This operator * is linear and continuous from

H -1 2 ( C ) to L 2 ( C ). Usually, it is taken as a convolution, τ * = ω * τ , τ ∈ H -1/2 ( C ), where ω ∈ D(-δ, δ), δ -δ ω = 1, with δ ∈ R, δ > 0,
but other choices can be made (see Badea and Krause 2012, for instance). We suppose that u 0 ∈ (H 1 2 ( D )) d , and define the sets of statically admissible fields:

H = H div ( ) = {τ = (τ i j ) 1≤i, j≤d : τ i j = τ ji , 1 ≤ i, j ≤ d, τ ∈ (L 2 ( )) d×d , divτ ∈ (L 2 ( )) d } H f,g = {σ ∈ H ; div σ + f = 0 in , σ n = g on F } τ = {σ ∈ H f,g ; σ N ≤ 0 , |σ T | ≤ -μτ N on C } (2.3) where τ ∈ (L 2 ( C )) d .
The case of the Tresca friction is a particular case of the last equation in (2.2), where the normal stress σ * N is replaced by a given η ∈ (L 2 ( C )) d . Note that the solution of the contact problem with Coulomb friction can be obtained by a fixed point procedure in which we solve a contact problem with Tresca friction at each iteration [START_REF] Licht | Remarks on a numerical method for unilateral contact including friction[END_REF]. The stress field variational formulation of (P) with Tresca friction law is (see Telega 1991, for instance):

Problem (P s ). Find the stress field σ : → H f,g such that σ ∈ η s(σ, τ -σ ) ≥ l(τ -σ ) for any τ ∈ η (2.4) with s(σ, τ ) = Sσ : τ dx and l(τ ) = D u 0 • τ nds.
From the properties of the operator S, we can conclude that the bilinear form s is symmetric and there exist two constants α, β > 0 such that, for any σ, τ ∈ H f,g , we have

α τ -σ 2 H div ( ) ≤ s(τ -σ, τ -σ ) and s(σ, τ ) ≤ β||σ || H div ( ) ||τ || H div ( ) (2.5)
where

τ 2 H div ( ) = (τ : τ + divτ • divτ ).
Also, there exists M > 0 such that

l(τ ) ≤ M τ H div ( ) (2.6)
for any τ ∈ H . In view of the above properties (2.5) and (2.6) of the bilinear form s and of the functional l, we notice that problem (2.4) has a unique solution σ ∈ η which also minimizes the functional I (τ ) = 1 2 s(τ, τ )l(τ ) over all τ ∈ η . These properties of s and l are also needed in the convergence proof of the Schwarz method.

Schwarz algorithm

Let us consider a domain decomposition of ,

= ∪ m i=1 i (3.7)
with an overlapping parameter δ, and we assume that it has the property

meas(∂( i ∩ (∪ m j=i+1 j )) ∩ D ) > 0 for any i = 1, . . . , m -1 (3.8)
This is an existence condition for the solutions of some auxiliary problems used in the proof of Lemma 3.1 to prove that Assumption 3.1 below is satisfied. We associate to this decomposition the spaces

H i 0,0 = {σ ∈ H div ( ) : div σ = 0 in i , σ = 0 on \ i σ | i n = 0 on ( ∩ ∂ i ) ∪ ( F ∩ ∂ i )}, i = 1, . . . , m
We point out that the elements of H i 0,0 are elements of H 0,0 ( ) which vanish on \ i .We introduce the Schwarz algorithm, written as subspace correction method, for the solution of problem (2.4) as Algorithm 3.1 We start the algorithm with an arbitrary σ 0 ∈ η . At iteration n + 1, having σ n ∈ η , n ≥ 0, we compute sequentially for i = 1, . . . , m, the local corrections σ n+1 i as the solution of the variational inequalities

σ n+1 i ∈ H i 0,0 , σ n+ i-1 m + σ n+1 i ∈ η : s(σ n+ i-1 m + σ n+1 i , σ i -σ n+1 i ) ≥ l(σ i -σ n+1 i ) for any σ i ∈ H i 0,0 , σ n+ i-1 m + σ i ∈ η (3.9)
and then we update

σ n+ i m = σ n+ i-1 m + σ n+1 i .
4

In this algorithm, we can take l(σ i -

σ n+1 i ) = D u 0 • (σ i -σ n+1 i
)nds in the place of

D ∩∂ i u 0 • (σ i -σ n+1 i )nds because σ i -σ n+1 i
vanishes outside i . The convergence of the above algorithm is obtained from a more general convergence result in [START_REF] Badea | Convergence rate of a multiplicative Schwarz method for strongly nonlinear inequalities[END_REF] given for algorithms in a reflexive Banach space. In this general framework, the convergence is proved by making an assumption and we have to show that it holds true in order to prove the convergence of a certain particular algorithm. In the case of our algorithm, this assumption is written below and the main result of the paper is to show that it is satisfied for the spaces introduced above. We point out that this assumption has only a theoretical value, the fact that it is satisfied is a sufficient condition for the convergence of the algorithm, but it is not used in the practical application of the algorithm. Assumption 3.1 There exists a constant C 0 such that for any σ, τ ∈ η and

σ i ∈ H i 0,0 , σ + i j=1 σ j ∈ η , i = 1, . . . , m, there exist τ i ∈ H i 0,0 , i = 1, . . . , m, satisfying σ + i-1 j=1 σ j + τ i ∈ η for i = 1, . . . , m, (3.10) τ -σ = m i=1 τ i , (3.11) m i=1 τ i 2 H div ( ) ≤ C 0 τ -σ 2 H div ( ) + m i=1 σ i 2 H div ( ) .
(3.12)

We have the following result concerning the convergence of the algorithm.

Theorem 3.1 We suppose that Assumption 3.1 holds. Then, if σ is the solution of problem (2.4) and σ n , n ≥ 0, are its approximations obtained from Algorithm 3.1, we have the following error estimation

||σ n -σ || 2 H div ( ) ≤ C C C + 1 n
where C and C are two positive constants, C depends on σ , σ 0 and α, and C depends on the constants α, β, the number m of subdomains, and is an increasing function of the constant C 0 introduced in Assumption 3.1.

Proof Since the bilinear form s and the linear functional l satisfy (2.5) and (2.6), we can apply Theorem 2 in [START_REF] Badea | Convergence rate of a multiplicative Schwarz method for strongly nonlinear inequalities[END_REF] to get the convergence of Algorithm 3.1 and the error estimation. Constant C can be taken as

C = 2βm α(1 -κ) 1 + 2C 0 + 2βm ακ C 2 0 , (3.13)
where 0 < κ < 1 is chosen such that C has a minimum value. Also, we can take

C = 2 α 1 2 s(σ 0 , σ 0 ) + l(σ 0 ) - 1 2 s(σ, σ ) -l(σ ) . (3.14)
where σ 0 is the starting approximation in algorithms and σ is the solution of problem (2.4).

Constants α and β come from problem (2.4) and are given in (2.5). Also, the number m of the subspaces can be associated with the number of colors needed to mark the subdomains such that the subdomains with the same color do not intersect with each other. Since, in general, this number of colors depends on the dimension of the Euclidean space where the domain lies (see [START_REF] Toselli | Domain decomposition methods-algorithms and theory[END_REF], for instance), we can conclude that the convergence rate of the algorithms essentially depends on the constant C 0 . Moreover, C and the convergence rate, too, are increasing functions of C 0 . Indeed, by a simple calculus, we get from (3.13) that

C = 2βm α 2βm α C 2 0 + 2C 0 + 1 + 2βm α C 2 0 2 .
In the following, we prove that Assumption 3.1 holds and estimate this constant C 0 .

To prove that this assumption holds, we associate to the decomposition (3.7) of , some functions

θ i ∈ H 1 ( ), i = 1, . . . , m, such that 0 ≤ θ i ≤ 1 on , θ i = 0 on ∪ m j=i+1 j \ i , and θ i = 1 on i \ ∪ m j=i+1 j . (3.15)
Since the overlapping size of the domain decomposition is δ, the above functions θ i can chosen to satisfy

|∂ x k θ i | ≤ C/δ, a.e. in , for any k = 1, . . . , d, (3.16)
where C is a constant independent of the domain decomposition. Now, let us consider σ ∈ η , σ i ∈ H i 0,0 such that σ + i j=1 σ j ∈ η , i = 1, . . . , m, and let τ be another element in η . Using the functions θ i , i = 1, . . . , m defined in (3.15) and (3.16), we first define τ i , i = 1, . . . , m.

For i = 1, we consider

u 1 ∈ H 1 ( 1 ∩ (∪ m j=2 j )), u 1 = 0 on ∂( 1 ∩ (∪ m j=2 j )) ∩ D , the solution of the problem 1 ∩(∪ m j=2 j ) Aε(u 1 ) : ε(v 1 ) + ∂( 1 ∩(∪ m j=2 j ))∩ C [θ 1 (τ -σ ) + (1 -θ 1 )σ 1 ]n • v 1 + ∂(∪ m j=2 j )∩ 1 (τ -σ )n • v 1 = 0 (3.17) for any v 1 ∈ H 1 ( 1 ∩ (∪ m j=2 j )), v 1 = 0 on ∂( 1 ∩ (∪ m j=2 j )) ∩ D ,

and define

τ 1 = ⎧ ⎨ ⎩ τ -σ on 1 \(∪ m j=2 j ) Aε(u 1 ) on 1 ∩ (∪ m j=2 j ) 0 o n(∪ m j=2 j )\ 1 (3.18)
For i = 2, we make the same construction as above by taking τ -σ -τ 1 in the place of τ -σ , 2 in the place of 1 and ∪ m j=2 j in the place of = ∪ m j=1 j . In this way, we look for the solution

u 2 ∈ H 1 ( 2 ∩ (∪ m j=3 j )), u 2 = 0 on ∂( 2 ∩ (∪ m j=3 j )) ∩ D , the solution of the problem 2 ∩(∪ m j=3 j ) Aε(u 2 ) : ε(v 2 ) + ∂( 2 ∩(∪ m j=3 j ))∩ C [θ 2 (τ -σ -τ 1 ) + (1 -θ 2 )σ 1 ]n • v 2 + ∂(∪ m j=3 j )∩ 2 (τ -σ -τ 1 )n • v 2 = 0 (3.19) for any v 2 ∈ H 1 ( 2 ∩ (∪ m j=3 j )), v 2 = 0 on ∂( 2 ∩ (∪ m j=3 j )) ∩ D ,

and define

τ 2 = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 0 o n 1 \(∪ m j=2 j ) τ -σ -τ 1 on 2 \(∪ m j=3 j ) Aε(u 2 ) on 2 ∩ (∪ m j=3 j ) 0 o n (∪ m j=3 j )\ 2 (3.20) For a certain 3 ≤ i ≤ m -1, we consider τ -σ -τ 1 -• • • -τ i-1 , the domains i and ∪ m j=i+1 j and u i ∈ H 1 ( i ∩ (∪ m j=i+1 j )), u i = 0 on ∂( i ∩ (∪ m j=i+1 j )) ∩ D , the solution of the problem i ∩(∪ m j=i+1 j ) Aε(u i ) : ε(v i ) + ∂( i ∩(∪ m j=i+1 j ))∩ C [θ i (τ -σ -τ 1 -• • • -τ i-1 ) +(1 -θ i )σ i ]n • v i + ∂(∪ m j=i+1 j )∩ i (τ -σ -τ 1 -• • • -τ i-1 )n • v i = 0 (3.21) for any v i ∈ H 1 ( i ∩ (∪ m j=i+1 j )), v i = 0 on ∂( i ∩ (∪ m j=i+1 j )) ∩ D ,

and define

τ i = ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ 0 o n (∪ i-1 j=1 j )\(∪ m j=i j ) τ -σ -τ 1 -• • • -τ i-1 on i \(∪ m j=i+1 j ) Aε(u i ) on i ∩ (∪ m j=i+1 j ) 0 o n (∪ m j=i+1 j )\ i (3.22)
Finally, we define

τ m = τ -σ - m-1 j=1 τ j (3.23)
The construction of these τ i , i = 1, . . . , m, concerns only the verification of the assumptions, but not the application of the algorithm. We have the following result concerning the above defined τ 1 , . . . , τ m . Lemma 3.1 Let τ 1 , . . . , τ m be defined in (3.17)-(3.23). Then,

τ i ∈H i 0,0 for i = 1, . . . , m (3.24) τ -σ - i j=1 τ j = 0 on \(∪ m j=i+1 j ) for i = 1, . . . , m -1 (3.25) σ + τ 1 ∈ η and σ + i-1 j=1 σ j + τ i ∈ η for i = 2, . . . , m (3.26) τ - i j=1 τ j + i j=1 σ j ∈ η for i = 1, . . . , m -1 (3.27)
Proof We notice that, in view of the condition (3.8) on the domain decomposition, Eqs.

(3.17), (3.19) and (3.21) have unique solutions. The proof of lemma is made by induction.

(a) For i = 1, from (3.17), we get divAε(u 1 ) = 0 on 1 ∩ (∪ m j=2 j ) and Aε(u 1 )n = (τ -σ )n on ∂(∪ m j=2 j ) ∩ 1 . Since div(τ -σ ) = 0 on 1 \(∪ m j=2 j ), we get from (3.18) that τ 1| 1 ∈ H div ( 1 ) and div(τ 1| 1 ) = 0 on 1 Since (τ -σ )n = 0 on F and, in view of (3.17), we have Aε(u

1 )n = 0 on ∂( 1 ∩ (∪ m j=2 j )) ∩ F , we get τ 1| 1 n = 0 on ∂ 1 ∩ F Also, we get from (3.17) that Aε(u 1 )n = 0 on ∂ 1 ∩ (∪ m j=2 j ), or τ 1| 1 n = 0 on ∂ 1 ∩
and therefore, the extension of τ 1| 1 with zero in \ 1 is an element of H 0,0 ( ). Consequently, (3.24) holds for τ 1 defined in (3.18). From the definition of τ 1 we get that τ -σ -τ 1 = 0 on 1 \(∪ m j=2 j ), and consequently, (3.25) holds for i = 1. Since τ 1 ∈ H 1 0,0 and τ, σ ∈ H f,g , we get σ + τ 1 ∈ H f,g .

On ∂( 1 ∩ (∪ m j=2 j )) ∩ C , using again (3.17), we have

(σ + τ 1 )n = (σ + Aε(u 1 ))n = [θ 1 τ + (1 -θ 1 )(σ + σ 1 )]n
Since τ, σ, σ + σ 1 ∈ η , in view of the above equation, we have

(σ + τ 1 ) N ≤ 0 and |(σ + τ 1 ) T | ≤ -μη N on C
In this way, we get that (3.26) holds for i = 1. In a similar way, we prove that (3.27) holds for i = 1. First, we see that

τ -τ 1 + σ 1 ∈ H f,g Also, on ∂( 1 ∩ (∪ m j=2 j )) ∩ C , we have (τ -τ 1 + σ 1 )n = (τ -Aε(u 1 ) + σ 1 )n = [θ 1 (σ + σ 1 ) + (1 -θ 1 )τ ]n
and, from here and the fact that τ, σ, σ + σ 1 ∈ η , we get

(τ -τ 1 + σ 1 ) N ≤ 0 and |(τ -τ 1 + σ 1 ) T | ≤ -μη N on C
Consequently, we get that (3.27) holds for i = 1. (b) For i = 2, let us write

σ = σ + σ 1 and τ = τ -τ 1 + σ 1 (3.28)
First, from (3.25) for i = 1, we get that

τ -σ = τ -σ -τ 1 = 0 on \(∪ m i=2 i ) (3.29)
Then, we reduce the case of i = 2 to that of i = 1 by considering the decomposition 2 , . . . , m of ∪ m j=2 j instead of the decomposition 1 , . . . , m of ∪ m j=1 j , and σ and τ in the place of σ and τ , respectively. From (3.29) and the definition of τ 2 in (3.20), we get that if we prove that τ 2|∪ m j=2 j belongs to H 2 0,0 corresponding to 2 in the decomposition 2 , . . . , m of ∪ m j=2 j , then τ 2 belongs to H 2 0,0 corresponding to 2 in the decomposition 1 , . . . , m of ∪ m j=1 j , i.e. (3.24) is satisfied for i = 2. We see that by replacing σ and τ in the definition of τ 2 given in (3.19) and (3.20), we get that τ 2|∪ m j=2 j is identically defined as τ 1 in (3.17) and (3.18) where σ and τ are replaced by σ and τ , respectively. Besides this, in view of the conditions in Assumption 3.1, we have σ , σ + σ 2 ∈ η . Also, from (3.27) with i = 1, we have τ ∈ η . Consequently, repeating the reasoning from i = 1 in the new context, we get that, for the decomposition 2 , . . . , m of ∪ m j=2 j , τ 2|∪ m j=2 j satisfies (3.24)-(3.27) from i = 1 with σ and τ replaced by σ and τ , respectively. Now, replacing σ and τ in these equations with their values given in (3.28) and taking into account (3.29), we get that (3.24)-(3.27)

hold for i = 2. (c) For a given 3 ≤ i ≤ m -1 we define σ = σ + i-1 j=1 σ j and τ = τ - i-1 j=1 τ j + i-1 j=1 σ j
and, with the above reasoning, we get that (3.24)-(3.27) for that i. Proof Let us consider τ, σ ∈ η , σ i ∈ H i 0,0 such that σ + i j=1 σ j ∈ η , and let τ 1 , . . . , τ m be defined in (3.17)-(3.23). We shall use Lemma 3.1 to prove that (3.10) and (3.11) hold. First, (3.24) shows that τ i ∈H i 0,0 for i = 1, . . . , m. Then, (3.26) proves (3.10), and finally, (3.23) is (3.11).

(d) For i = m, we consider τ m defined in (3.23). Since τ , σ , σ 1 , . . . , σ m-1 ∈ H div ( ), it follows that τ m ∈ H div ( ). Moreover, from (3.25) for i = m -1, we get that τ m = 0 on \ m . Besides that, since τ -σ ∈ H 0,0 and, for i = 1, . . . , m -1, τ i ∈ H i 0,0 , we get τ m ∈ H m 0,
We now prove that (3.12) holds with C 0 given in (3.30), and we first evaluate the norms of τ i , i = 1, . . . , m.

(a) For i = 1, by replacing v 1 with u 1 in (3.17), and using the trace theorems in H 1 and H div (see Theorem 2.5, page 27, in [START_REF] Girault | Finite element methods for Navier-Stokes equations[END_REF] for the vectorial case, for instance), we get

1 ∩(∪ m j=2 j ) Aε(u 1 ) : ε(u 1 ) = - ∂( 1 ∩(∪ m j=2 j ))∩ C [θ 1 (τ -σ ) + (1 -θ 1 )σ 1 ]n • u 1 - ∂(∪ m j=2 j )∩ 1 (τ -σ )n • u 1 = - ∂( 1 ∩(∪ m j=2 j ))∩( C ∪ 1 ) [θ 1 (τ -σ )+(1-θ 1 )σ 1 ]n • u 1 ≤ ||θ 1 (τ -σ )n + (1 -θ 1 )σ 1 n|| H -1/2 (∂( 1 ∩(∪ m j=2 j ))) • ||u 1 || H 1/2 (∂( 1 ∩(∪ m j=2 j ))) ≤ C||θ 1 (τ -σ ) + (1 -θ 1 )σ 1 || H div ( 1 ∩(∪ m j=2 j )) ||u 1 || H 1 ( 1 ∩(∪ m j=2 j ))
Above, we have used the fact that

θ 1 = 1 on ∂(∪ m j=2 j )∩ 1 = ∂( 1 ∩(∪ m j=2 j ))∩ 1 .
Constant C can be taken depending only of , independently of the domain decomposition. From now on, we denote by C a generic constant with this property. Finally, we have written, for instance,

|| • || H 1 ( 1 ∩(∪ m j=2 j )) in the place of || • || H 1 ( 1 ∩(∪ m j=2 j )) d .
The superscripts d or d × d in the notation of the norms will also be omitted in the following. In view of (3.8), we have (see Théorème 3.3, p. 115, in Duvaut and Lions 1972, for instance)

||u 1 || 2 H 1 ( 1 ∩(∪ m j=2 j )) ≤ C 1 ∩(∪ m j=2 j )
Aε(u 1 ) : ε(u 1 ) Also, using the properties of S, we have

γ ||Aε(u 1 )|| 2 L 2 ( 1 ∩(∪ m j=2 j )) ≤ (S(Aε(u 1 )), Aε(u 1 )) L 2 ( 1 ∩(∪ m j=2 j )) × 1 ∩(∪ m j=2 j ) Aε(u 1 ) : ε(u 1 )
From the above three equations, we get

||Aε(u 1 )|| L 2 ( 1 ∩(∪ m j=2 j )) ≤ C||θ 1 (τ -σ ) + (1 -θ 1 )σ 1 || H div ( 1 ∩(∪ m j=2 j ))
and since div(Aε(u 1 )) = 0, we have

||Aε(u 1 )|| H div ( 1 ∩(∪ m j=2 j )) ≤ C||θ 1 (τ -σ ) + (1 -θ 1 )σ 1 || H div ( 1 ∩(∪ m j=2 j ))
Moreover, in view of (3.18) and since θ 1 = 1 on 1 \(∪ m j=2 j )), we get

||τ 1 || H div ( ) ≤ C||θ 1 (τ -σ ) + (1 -θ 1 )σ 1 || H div ( 1 )
or, in view of (3.16) and taking into account that div(τ -σ ) = divσ 1 = 0 on , we get

||τ 1 || H div ( ) ≤ C 1 + d δ 2 1 2 [||τ -σ || L 2 ( 1 ) + ||σ 1 || L 2 ( 1 ) ] = C 1 + d δ 2 1 2 [||τ -σ || H div ( 1 ) + ||σ 1 || H div ( 1 ) ] ≤ C 1 + d δ 2 1 2 [||τ -σ || H div ( ) + ||σ 1 || H div ( ) ]
Consequently, we can write

||τ 1 || 2 H div ( ) ≤ C mδ [||τ -σ || 2 H div ( ) + ||σ 1 || 2 H div ( ) ] (3.32)
where C mδ is given in (3.31). (b) For 2 ≤ i ≤ m -1, using the same arguments as above, we get

||τ i || H div ( ) ≤ C||θ i τ -σ - i-1 j=1 τ j + (1 -θ i )σ i || H div ( i ) and ||τ i || H div ( ) ≤ C 1 + d δ 2 1 2 ⎡ ⎣ ||τ -σ || H div ( i ) + i-1 j=1 ||τ j || H div ( i ) + ||σ i || H div ( i ) ⎤ ⎦ ≤ C 1 + d δ 2 1 2 [||τ -σ || H div ( ) + i-1 j=1 ||τ j || H div ( ) + ||σ i || H div ( ) ]
Also, we can write

||τ i || 2 H div ( ) ≤ C mδ ⎡ ⎣ ||τ -σ || 2 H div ( ) + i-1 j=1 ||τ j || 2 H div ( ) + ||σ i || 2 H div ( ) ⎤ ⎦ (3.33) (c)
In view of (3.32) and (3.33), we can prove by induction that i j=1

||τ j || 2 H div ( ) ≤ i j=1 C j mδ [||τ -σ || 2 H div ( ) + ||σ i+1-j || 2 H div ( ) ] for i = 1, . . . , m -1. Now, from (3.23), we have ||τ m || 2 H div ( ) ≤ m ⎡ ⎣ ||τ -σ || 2 H div ( ) + m-1 j=1 ||τ j || 2 H div ( ) ⎤ ⎦
From these two last equations we get

m i=1 ||τ i || 2 H div ( ) ≤ m m-1 j=0 C j mδ ||τ -σ || 2 H div ( ) + m m-1 j=1 C j mδ ||σ m-j || 2 H div ( )
and therefore, we can take constant C 0 as in (3.30).

Summing up the results in Theorem 3.1 and Proposition 3.

τ i = θ i (τ -σ - i-1 j=1 τ j ) + (1 -θ i )σ .
In this case, using property (3.16), we obtain, without using the recursivity in the proof of Proposition 3.1, that

||τ i || 2 H div ( ) ≤ C mδ ⎡ ⎣ ||τ -σ || 2 H div ( ) + i j=1 ||σ j || 2 H div ( ) ⎤ ⎦
in the place of (3.33). The introduction of Aε(u i ) in the definition of τ i imposes the use of the recursivity, and consequently, C 0 in (3.30) depends on the powers of C mδ . However, as we already said, the number m of the subdomains can be considered as being the minimal number of colors needed to mark the subdomains, and therefore, C 0 does not depend on the actual number of subdomains in the domain decomposition.

Conclusions

We have studied the convergence of the Schwarz method for contact problems with Tresca friction formulated in stress variables. In this dual variable, the problem is written as a variational inequality in the space H div ( ), being the domain of the problem. The method is introduced as a subspace correction algorithm and we have used a general result in [START_REF] Badea | Convergence rate of a multiplicative Schwarz method for strongly nonlinear inequalities[END_REF] which proves the convergence provided that an assumption is satisfied. The main result of this paper was to prove that the assumption made in the general convergence theorem is verified for our particular variational inequality. As in the case of the classical problems formulated in primal variables, we show that the method is geometrically convergent with a convergence rate which depends on the overlapping parameter of the domain decomposition.

  0 , i.e. (3.24) holds for i = m. Finally, using (3.27) for i = m -1, we get independent of the domain decomposition.

  1, we get Corollary 3.1 We assume that the domain decomposition (3.7) satisfies (3.8). Then, if σ is the solution of problem (2.4) and σ n , n ≥ 0, are its approximations obtained from Algorithm 3.1, we have the following error estimation ||σ n -σ || 2 H div ( ) ≤ C and C are given in (3.13) and (3.14), respectively, and constant C 0 in C is given in (3.30). For problems in primal variables (see Badea 2006, for instance), the construction of the functions similar to τ i in (3.17)-(3.23) is directly done by writing

	C
	C + 1
	Remark 3.1

n where C

Acknowledgments

The authors acknowledge the support of the Laboratoire Européen Associé CNRS Franco-Roumain de Mathématiques et Modélisation, LEA Math-Mode, for this paper.