Ionel Nistor
email: ionel.nistor@edf.fr

Martin L E Guiton

Patrick Massin

Nicolas Moes

Samuel Géniaut

N Moës

An X-FEM approach for large sliding contact along discontinuities

Keywords: X-FEM, contact, large sliding

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Several recent advances in non-linear mechanics concern the investigation of kinematics with large displacement jumps along sliding discontinuities, in various application fields like biomechanics [START_REF] Heegard | The biomechanics of the human patella during passive knee flexion[END_REF], metal cutting process [START_REF] Liu | Finite element hydrodynamic friction model for metal forming[END_REF][START_REF] Marusich | Modelling and simulation of high-speed machining[END_REF] or structural geology [START_REF] Smart | Deformation behavior during blind thrust translation as a function of fault strength[END_REF][START_REF] Erickson | Mechanics of extension and inversion in the hanging walls of listric normal faults[END_REF][START_REF] Ellis | Comparisons between analogue and numerical models of thrust wedge development[END_REF][START_REF] Erickson | Sequential restoration and unstraining of structural cross sections: applications to extensional terraces[END_REF][START_REF] Maerten | Chronologic modeling of faulted and fractured reservoirs using geomechanically based restoration[END_REF][START_REF] Sanz | Mechanical aspects of thrust faulting driven by far-field compression and their implications for fold geometry[END_REF][START_REF] Zammali | 2D and 3D finite element restorations of geological structures with sliding contact along faults[END_REF], the latter motivating the present work. In these problems, the surfaces of discontinuity can result from either preexisting material discontinuities, or from a complex history of nucleation [START_REF] Bellec | A note on enrichment functions for modelling crack nucleation[END_REF] and coalescence of fractures [START_REF] Cardle | A method for multiple crack growth in brittle materials without remeshing[END_REF] or damage zones. Apart from the strong interest of the numerical engineering community to study the damage history in relation to the bulk constitutive law, it is commonly assumed that when the bandwidth of the damage zones is smaller than a characteristic scale of the representative elementary volume, the concept of displacement discontinuities (or 'strong discontinuities' in [START_REF] Oliver | Strong discontinuities and continuum plasticity models: the strong discontinuity approach[END_REF]) can be recovered. Doing so, the model topology is defined by a network of complex non smooth and sometimes intersecting surfaces at the limit of, or inside material bodies. Building a 3D finite element method (FEM) mesh that matches this network of discontinuities is still a time-consuming operation, which requires a strong user interaction (see for instance the work of Lepage [START_REF] Lepage | Génération de maillages tridimensionnels pour la simulation des phénomènes physiques en géosciences[END_REF] for geo-modeling). For similar reasons, updating such meshes when changing the discontinuity geometry is still not fully automatic to our knowledge. This limitation is of particular importance if one forecasts to constrain and analyze the uncertainty on the geometrical model, as it is the case for subsurface inter layers and faults of hydrocarbon reservoirs, which are really observed at given well locations only [START_REF] Thore | Structural uncertainties: determination, management and applications[END_REF][START_REF] Singh | Automatic seismic velocity inversion using multi-objective evolutionary algorithms[END_REF].

Mesh dependency on discontinuity is also one of the well-known reasons for FEM limitations when modeling the propagation of the discontinuities. Following the classical re-meshing techniques, a new mesh has then to be constructed on which all unknowns and internal variables of the model have to be projected, see for instance [START_REF] Yao | 3D shape optimal design and automatic finite element regridding[END_REF][START_REF] Jin | Two dimensional mesh generation; adaptive remeshing and refinement[END_REF]. The recent extended finite element method (X-FEM) developed in [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF] proposes an alternative to re-meshing, which relaxes the mesh dependency on discontinuity by enriching the displacement approximation, thanks to the partition of unity principle [START_REF] Babuška | The partition of unity method[END_REF]. It has been extensively applied to problems involving mode I opening or cohesive discontinuities in quasi static or dynamic frameworks (e.g. [START_REF] Moës | Extended finite element method for cohesive crack growth[END_REF][START_REF] Belytschko | Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment[END_REF][START_REF] Sukumar | Modelling quasi-static crack growth with the extended finite element method. Part I: computer implementation[END_REF][START_REF] Nistor | Numerical propagation of dynamic cracks using X-FEM[END_REF]). The use of X-FEM for modeling cracked structures involving large displacements remains also limited [START_REF] Areias | Analysis of fracture in thin shells by overlapping paired elements[END_REF]. Its application to a closed discontinuity requires the coupling of X-FEM with a contact algorithm, therefore adding a mechanical interaction between the two sides of the discontinuities. Several strategies have been developed to account for contact between deformable bodies (a state of the art is given in [START_REF] Wriggers | Computational Contact Mechanics[END_REF]). X-FEM has been coupled with various contact algorithms: the LATIN method in [START_REF] Dolbow | An extended finite element method for modeling crack growth with frictional contact[END_REF], a penalty method in [START_REF] Khoei | An enriched finite element algorithm for numerical computation of contact friction problems[END_REF], a mortar method in [START_REF] Kim | A mortared finite element method for frictional contact on arbitrary interfaces[END_REF] and the hybrid augmented lagrangian method of Ben Dhia and Zarrong [START_REF] Dhia | Hybrid frictional contact particle-in elements[END_REF] in [START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF].H o w e v e r ,a l l these works remained limited to the small perturbation hypothesis and to small sliding along discontinuities.

In this contribution, we present a new approach generalizing the work of Géniaut et al. [START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF] by allowing large displacements on a closed frictionless discontinuity within the X-FEM framework. Basic principles of the contact mixed formulation developed in [START_REF] Dhia | Hybrid frictional contact particle-in elements[END_REF] are first introduced. Then we briefly recall the principles of X-FEM pointing particularities related to our implementation in which the discontinuity representation by the level set method [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF] allows now to account for large displacements along it. The weak formulation of the contact equations is then discretized in the X-FEM framework with details given on the numerical implementation for the construction of a new mixed contact element in case of large sliding and on the geometrical update algorithm. The methodology is finally tested on three simple 2D numerical applications with large sliding on a discontinuity associated with a small compression of elastic volumes. A comparison with a FEM solution is analyzed in the first problem. The second example illustrates the sensitivity on the choice of the pressure field interpolation in order to satisfy the LBB condition and reduce pressure oscillations for particular mesh intersection with the discontinuity. The third example finally accounts for large rotations and shows a convergence of the numerical solution toward an analytical one when the mesh is refined. All the corresponding numerical developments were conducted in the free finite element software Code Aster [http://www.codeaster.org].

CONTACT MIXED FORMULATION

A comprehensive description of computational approaches for contact problems can be found in [START_REF] Wriggers | Computational Contact Mechanics[END_REF]. We are interested in this paper in the quasi-static resolution of problems with large sliding contact along discontinuities. To provide more details on such boundary-value problem, let us consider a domain ⊂ R n dim where n dim is the number of spatial dimensions and loading conditions as represented in Figure 1(a) in the undeformed configuration. We consider in this section the case when contains one discontinuity without lack of generality if one wants to extend it to the case of several discontinuities. Note that in this paper tensorial quantities are represented by bold face characters.

The boundary of can be split into a set u with Dirichlet boundary condition, a set t with Neumann boundary condition and the discontinuity surface, which is decomposed into two opposite surfaces + c and - c as shown in Figure 1(a) on which contact traction r may be applied. The domain and the surfaces corresponding to , + c and - c in the deformed configuration are denoted by , + c and - c , respectively. We also use the upper indices + andon vectorial fields to denote their restriction on + c and - c , respectively. One should notice that in general it is not possible to predetermine exactly the domain on which contact efforts are applied as it evolves with the deformation. Instead, only potential contact surfaces can be defined initially, assuming at a given state that for the part of them which are not in contact r is set to the zero vector.

The corresponding strong form of the Lagrangian equilibrium problem completed by the boundary conditions is

Div P+f = 0 in P•N = t on t P•N + = r + on + c P•N -= r -on - c u = ū on u (1)
where P is the first Piola-Kirchhoff stress tensor, f is the body force and N is the outward unit normal vector to in the undeformed configuration. In a first step and for the sake of simplicity, we have restricted the contact in this paper to the frictionless case. Because of the evolving nature of contact, it is generally convenient to decompose the Lagrangian contact traction applied from one contact surface on the opposite in an Eulerian basis as follows:

r + = n - r -= n + (2)
where is the normal component (known also as the contact pressure) and n is the outward unit normal vector to in the deformed configuration. Considering V as the space of kinematical admissible field and V * the space of kinematical admissible field to zero, the principle of the Virtual Work can be written as:

find u ∈ V satisfying P: (∇ r u * ) T d - + c n -• u * + d - - c n + . u * -d = f•u * d + t t•u * d ∀u * ∈ V * (3)
where ∇ r u * is the gradient with respect to the reference coordinates of u * . At that point, several choices can be adopted for the surface on which the contact contribution will be integrated. In this paper, like in [START_REF] Dhia | Hybrid frictional contact particle-in elements[END_REF][START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF], we will adopt the biased master/slave approach for which we choose one of the contact surface labeled as 'slave' and check if it is penetrating or not the other which is labeled as 'master'. Consequently, the contact pressure multipliers are defined only on the slave surface. Note that this bias could significantly influence the resulting contact pressure in case of a marked contrast between the mesh refinement of the opposite contact surfaces [START_REF] Wriggers | Computational Contact Mechanics[END_REF]. Alternately, other approaches have been proposed to overcome this bias effect by considering a space of contact pressure multipliers on each of the two surfaces with a condition to enforce the continuity of the contact pressure field. This is achieved, for example, by means of a formulation in a sequence of two simultaneous sub-problems with an intermediate surface for frictionless contact in [START_REF] Papadopoulos | A novel finite element formulation for frictionless contact problems[END_REF], and by means of a mortar formulation with an intermediate surface for frictional contact in [START_REF] Mcdevitt | A mortar-finite element formulation for frictional contact problems[END_REF]. For the present paper, it is believed that the X-FEM described later on should help to minimize the constraints on the mesh design, and therefore to avoid pathological situations thanks to a rather uniform mesh refinement over the modeled domain.

In what follows, we consider + c to be the slave surface. The detection of contact generally involves the evaluation of the gap between + c and - c , see Figure 1(b), defined as:

d n = (x + -x+ )•n - (4)
The point x+ , called the master point, is the orthogonal projection on - c of a given slave point x + , understood in a generalized sense if - c is not regular at x+ . The gap and the contact pressure are then related through the following Hertz-Signorini-Moreau conditions:

0, d n 0, d n = 0 (5)
The first relation expresses that a closed gap leads to a compressive contact pressure, the second one expresses the non-penetration condition and the last one guarantees that the contact pressure vanishes when the gap is open.

We introduce now , the Heaviside function, as

(x) = 1i f x 0 0i f x<0 (6)
and g n , the augmented contact multiplier, as

g n = -n d n (7)
with n , the augmentation parameter, being strictly positive.

The presence of the augmented contact multiplier in the adopted mixed contact formulation offers some advantages compared with a penalty formulation: the solution is exact, not an approximation as in the penalization case and the augmentation parameter, n , allows to control the conditioning of the system [START_REF] Alart | A mixed formulation for frictional contact problems prone to Newton like solution methods[END_REF]. It also allows us to consider the Hertz-Signorini-Moreau conditions in a single equation more suitable for the weak formulation of the contact equation:

-(g n )g n = 0 (8)
In [START_REF] Maerten | Chronologic modeling of faulted and fractured reservoirs using geomechanically based restoration[END_REF] in case of contact g n 0s ot h a t = 1, and in case of absence of contact g n <0s ot h a t = 0. In the present work we use the contact algorithm introduced in [START_REF] Dhia | Hybrid frictional contact particle-in elements[END_REF], which was extended to the X-FEM framework for small perturbation in [START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF]. In this method, contact conditions are seen as an interface law and not as boundary conditions. It is related to the augmented lagrangian method of [START_REF] Laursen | A continuum element-based formulation for the implicit solution of multi-body, large deformation frictional contact problem[END_REF] but with the following particularities: contact non-differentiability is solved with an active set of constraint algorithm [START_REF] Dumont | Algorithme des contraintes actives et contact unilateral sans frottement[END_REF] and the non-differentiability related to the geometrical non-linearity is solved with a fixed point method.

From the principle of the virtual work (3), the equivalent form of the contact conditions (8) and the equilibrium of contact tractions, the weak frictionless contact formulation is now determined. Let H be the space of the contact pressure field. The quasi-static frictionless contact problem consists in finding (u, ) ∈ V × H , ∀(u * , * ) ∈ V * × H , satisfying the following equations:

P: (∇ r u * ) T d - + c (g n )g n n -•[[u * ]] d = f•u * d + t t•u * d (9) + c -1 n { -(g n )g n } * d = 0 ( 10
)
where

[[ u * ]] = u * + -u * -.
In the equilibrium equation ( 9), the contact pressure contribution is taken into account by the last term on the left side. Equation [START_REF] Zammali | 2D and 3D finite element restorations of geological structures with sliding contact along faults[END_REF] is the weak form of the Signorini contact law. In the following, this weak formulation is completed by considering the kinematics to allow large rotations but under a small deformation hypothesis. The material constitutive law is the Kirchhoff-St Venant linear elasticity, which relates the second Piola-Kirchhoff stress tensor to the Green-Lagrange strain tensor S = Tr(E)1+2 E [START_REF] Bellec | A note on enrichment functions for modelling crack nucleation[END_REF] where Tr(E) is the trace of the lagrangian strain tensor E and 1 is the second-order identity tensor. The relationship between the first and the second Piola-Kirchhoff stress tensors is given by

S = F -1 • P ( 12
)
where F is the deformation gradient tensor.

3. DISCRETIZATION WITHIN X-FEM FRAMEWORK

X-FEM and level sets in our approach

The main idea of the X-FEM, detailed in [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF], is to avoid the re-meshing technique by proposing instead the construction of an enriched approximation of the kinematic fields resulting from the intersection between the geometry of discontinuities and the mesh. A standard finite element approximation is enriched in the neighborhood of the discontinuity by local functions related to additional degrees of freedom dof. This enrichment is made possible through the partition of unity method [START_REF] Babuška | The partition of unity method[END_REF]. In the X-FEM framework, one introduces a discontinuous function to represent the jump in displacement across the discontinuity. One must mention that the enriched displacement approximation is global, whereas the support of these functions is local because they are multiplied by classical nodal shape functions. Note that we do not consider in this paper the case of cracks for which a special enrichment must be applied to account for the asymptotic displacement field near the crack tip [START_REF] Béchet | Improved implementation and robustness study of the X-FEM for stress analysis around cracks[END_REF][START_REF] Laborde | High order extended finite element method for cracked domains[END_REF].

For the sake of simplicity, we assume in what follows that the initial mesh for which we construct the X-FEM enrichment is constituted of 2D bilinear finite elements. We recall here, using a simple example shown in Figure 2, the X-FEM approach in order to build the enriched displacement approximation for a structure containing a discontinuity. For a given mesh, once the geometry of the discontinuity is known, numerical algorithms determine for which nodes the additional dof are necessary. Let us call as 'nodal support' the domain on which the interpolation functions associated with a given node are defined, considering all the elements to which it belongs (e.g. the gray domain for the bold circle node in Figure 2). For the nodal support entirely cut by the discontinuity, the classical shape functions are enriched with the generalized Heaviside function

H (d(X)) = -1i f d(X)<0 +1i f d(X)>0 ( 13
)
where d(X) is the distance to the discontinuity and its sign is related to which side of the discontinuity the node belongs to, depending on the location of X, see Figure 2. The corresponding nodes, circled in Figure 2, will have one additional vectorial dof accounting for the jump representation, which will be denoted here by a.

Finally, the X-FEM enriched displacement approximation is given at any location X by

u h (X) = I ∈N I (X)u I + I ∈N I (X)H (d(X))a I ( 14
)
where are the shape functions, u is the vector of the classical dof, N is the set of all nodes, N is the set of nodes enriched with the generalized Heaviside function.

For the numerical implementation of this method in our finite element code, level set method [START_REF] Sethian | Level Set Methods and Fast Marching Methods[END_REF] was used in order to localize the discontinuities. The main idea is to define a regular function (X) (Lipschitzian at least) such that (X) = 0 represents the interface describing the discontinuity. In the discrete representation, (X) is the normal level set function giving the distance from all nodes of the mesh to the discontinuity, d(X) in [START_REF] Oliver | Strong discontinuities and continuum plasticity models: the strong discontinuity approach[END_REF]. The sign of the normal level set is arbitrarily chosen, but is related to the sign of the generalized Heaviside function. Thus, for the nodes where (X)>0 and which are to be enriched, d(X)>0 in [START_REF] Oliver | Strong discontinuities and continuum plasticity models: the strong discontinuity approach[END_REF]. As we will see later, this gives us the possibility to change slave/master status of the contacting surface by only changing the sign of the level set function.

For the numerical implementation of the level set function, we have chosen to store its numerical values in a nodal field. To each node of the mesh is assigned the numerical value of (X), corresponding to the closest discontinuity in case of multiple ones. This field is approximated by the linear shape functions used for the displacement field approximation:

(X) = I ∈N I (X) I ( 15
)
In this way the contact segments (dotted lines on Figure 3) are determined by the intersection between the normal level set function and the edges of the elements having for its nodes opposite signs for the lsn values. These segments support the contact integration domains that we will introduce in the next subsection.

Discretization of the contact equations

The mixed contact formulation (Equations ( 9) and ( 10)) was discretized in the classical framework of FEM in [START_REF] Dhia | Hybrid frictional contact particle-in elements[END_REF], using mixed pressure-displacement contact elements. To point out the differences between the FEM and X-FEM in solving the contact problem using the above formulation, we present a simple FEM example of a hybrid contact element. In this example, the traces of the bilinear 2D finite elements on the contact boundaries are constituted of two-node segments. Let us consider two body candidates to be in contact as shown in Figure 4. The contact surfaces are declared as master and slave and the mixed contact elements are created by projecting given contact points belonging to the slave surface on the master surface. As shown in Figure 4, after the projection the master element is found (N1N2 in this case) and a hybrid contact element is created (a SEG2-SEG2 element, labeled N1N2N3N4 in our example). These elements have both geometrical dof (i.e. displacements) and contact ones (i.e. contact pressure and friction semimultipliers) stored at the node level. The geometrical dof allow for the gap evaluation, whereas the contact dof allow for contact pressure and frictional effort interpolation over the elementary domains using the suitable shape functions. The coupling of X-FEM with this contact formulation in case of small perturbations was achieved by Geniaut et al. [START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF]. The contact hybrid elements in this case are quite different from the classical ones presented in [START_REF] Dhia | Hybrid frictional contact particle-in elements[END_REF]. In fact, since in X-FEM a discontinuity surface is introduced in the mesh not by means of additional nodes but only with a set of intersection points, no surface elements can be created as it is the case of FEM. In this situation cut finite elements are used to create a new hybrid contact element just by adding the contact dof. This task was facilitated due to assumptions specifically related to the small perturbation framework:

• both sides of the interface are considered as a single surface, so that the master and slave distinction does not exist anymore, even if we need to choose a normal orientation to define the displacement gap; • there is no geometrical update;

• no contact search is necessary since the contact points and their projections are identical.

The generalization of the X-FEM contact formulation of [START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF] to large sliding contact is necessary for solving large kinematic problems, as presented in the introduction. A typical situation for large sliding along a discontinuity is shown in Figure 5. Obviously in this situation, the assumptions considered for small perturbations are not valid any more and treating this problem with X-FEM needs to find a way to communicate between a part of a cut element and parts belonging to other cut elements, in order to compute contact contributions. This limitation motivates the design of a new hybrid element adapted to large sliding. The creation of this X-FEM contact hybrid element and the changes in the discretized equations with respect to small perturbations is presented in the next section, which includes also the description of the numerical implementation.

The previous implementation allowing for the small perturbation analysis of structures under contact by X-FEM served as a starting point. At the present stage the approach has been applied to 2D models with a bilinear interpolation of the displacement (quadrilateral elements). In the method that we propose, the main stages for solving the contact problem involving large sliding are:

1. loop on the loading steps; 2. geometry update and contact search algorithm; 3. loop on active constraints (contact status checking); 4. Newton's iterations.

Using the simple example shown in Figure 5, we present now the geometrical update and contact search algorithm. As one can observe for the interface geometry, two sides are now considered denoted as master and slave (see Equation ( 4)). The slave or master status can be set by the user, based on the sign of the generalized Heaviside function. Two series of intersection points, p s i and p m i , are generated on the contact surfaces. Note that initially, a slave point and its associated master point are located at the same position, as for the small perturbation case. Updating the contact geometry means to compute before each loading step the coordinates of the intersection points as shown in Figure 5. To do so, different displacement fields are considered depending on which side of the discontinuity the displacement is interpolated. Using the enriched displacement approximation given in [START_REF] Lepage | Génération de maillages tridimensionnels pour la simulation des phénomènes physiques en géosciences[END_REF] in case of the considered example, we find for a slave intersection point k in the lower part of the discontinuity

u( p s k ) = 4 j=1 s j u j - 4 j=1 s j a j ( 16
)
and in a similar manner for a master intersection point l in the upper part of the discontinuity:

u( p m l ) = 4 j=1 m j u j + 4 j=1 m j a j ( 17
)
Once sliding begins, a contact search algorithm must provide the slave-master contact pairs in order to create the hybrid elements. Let us consider the situation shown in Figure 5: the contact point CP belonging to the slave contact segment of the element E7 is looking for a master contact segment. Looping on master contact segments of cut elements, one finds the nearest master intersection point as p m 2 belonging to E6a n dE7. The elements connected to this master intersection point are determined and projections are done on their contact segments. The shortest inside projection (PR in this example) provides the target element E6 to create the hybrid contact element. Note that to each contact point corresponds a hybrid contact element which may change from one time step to the next (i.e. the master element changes during the sliding). The flowchart of this algorithm is given in Algorithm 1.

With this information, hybrid contact elements are generated like the one shown in Figure 6. For its construction we adopted a similar 'phantom nodes' approach as the one introduced in [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF] and used in the X-FEM framework by [START_REF] Song | A method for dynamic crack and shear band propagation with phantom nodes[END_REF]. One uses the 'useful' parts of the original elements coming in contact: the lower part of the slave element (the left one, defined by H =-1) and the upper part of the master element (the right one, defined by H =+1) for the corresponding displacement discretizations. The contact unknowns are stored only at the slave element nodes plus additional mid-edge nodes in order to define a linear continuous interpolation (see Figure 6). A discussion about this choice for storing X-FEM contact unknowns in our code is given in [START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF]. The master element is a standard bilinear element that has only geometrical dof. The complete dof map for the hybrid contact element shown in Figure 6 are given in Table I. Note that the asymmetry of the contact dof distribution comes from the biased (master/slave) nature of our contact method.

The numerical integration of the contact contributions presented in Equations ( 9) and ( 10) needs the evaluation of the gap between the contact point and its projection (Equation ( 4)) and the approximation of the contact unknowns over the contact domain. Both evaluations can be performed on the hybrid contact element using three elementary interpolation domains: the slave one, the master one and the slave contact segment one. More clearly, the gap calculus formula for the hybrid contact element shown in Figure 6 is: where n CP are the reference coordinates for the contact point CP on the eight-node slave element, n PR are the reference coordinates on the useful upper part of the four-node master element for the projection of CP on the contact segment [P m 1 -P m 2 ], s i are the shape functions on the parent domain of the slave element, m j are the shape functions on the parent domain of the master element, n PR is the outward normal at the projection point belonging to the useful upper part of the master element.

d n = 4 i=1 s i (n CP )u i - 4 i=1 s i (n CP )a i - 12 j=9 m j (n PR )u j + 12 j=9 m j (n PR )a j •n PR (18)
If the master and slave elements are identical (in case of small perturbations), Equation ( 18) is simplified and we find the well-known relation of the jump inside a cut X-FEM element:

d n = -2 4 i=1 s i (n CP )a i •n PR ( 19
)
The contact unknowns approximation is achieved over the parent domain related to the contact segment [P s 1 -P s 2 ] by

= 2 k=1 k (f CP ) k ( 20
)
where f CP are the reference coordinates for the contact point on the two-node contact element, k are the shape functions on the parent domain of the 1D contact element.

In fact, we consider a fictitious 1D element P s 1 P s 2 to perform this approximation, but there are no dof associated to this element. Only the shape function values are evaluated according to the contact point position on this segment, but the contact dof are stored on the slave nodes from N 1 to N 8, as mentioned before. This fictitious element gives also the numerical integration domain for the contact contributions since the contact points inside this segment are the integration points.

The integration rule on the contact segments is chosen to be trapezoidal with integration points located at the position of nodes for the following reason. Consider a contact segment to be free (i.e. not in contact) adjacent to a segment in contact. If choosing a standard Gauss integration, the integration points inside the free segment are both with a free status. Because of the continuous linear interpolation of the contact pressure, the contact pressure unknown at the node common to the two segments is forced to be zero. To compensate, the contact pressure at the other node of the segment in contact must be higher than expected. This results in an oscillation artefact of the contact pressure. Such oscillation is avoided by a kind of smoothing when considering the trapezoidal rule, because the integration point is at the node location which is common to both contacting and non-contacting segments. In this way, the integration point is the point where we look at the results and avoids interpolation.

Therefore, in Equations ( 9) and ( 10) three sets of shape functions (slave and master for displacements and the third for the pressure) appear when computing the contact terms of the elementary matrix and external force vector. Generalizing this approach to the 3D case is easy since only a new hybrid contact element needs to be built, made of 3D slave and master elements with a contact surface partitioned into several triangles.

Before we proceed with the linearization of ( 9) and ( 10), let us identify the non-linearities that are involved. As already mentioned, frictionless contact is considered and the geometrical non-linearity that can modify the contact normal vector n and all related quantities are treated by a fixed-point algorithm. Note also that we do not consider here non-linear bulk constitutive laws. However, the kinematics can present finite rotations and justify to employ the Kirchhoff-St Venant elasticity. Thus, there is still a geometrical non-linearity that appears in the Lagrangian strain tensor requiring to linearize Equations ( 9) and ( 10) in order to produce a suitable form for a Newton-Raphson strategy. Note that n remains fixed during the Newton-Raphson iterations. Let us define the directional Gateaux derivatives of an integral term G(u, u * , , * ) of Equations ( 9) and ( 10) as follows:

D u [G(u, u * , , * )]= d d ⏐ ⏐ ⏐ ⏐ =0 G(u+ u, u * , , * ) D [G(u, u * , , * )]= d d ⏐ ⏐ ⏐ ⏐ =0 G(u, u * , + , * ) (21)
Considering a k Newton iteration, the displacement and contact unknowns at the end of the iteration are:

u k+1 = u k + u k+1 k k+1 = k k + k k+1 (22)
For the sake of brevity, in what follows n will denote n -. Using definition [START_REF] Moës | Extended finite element method for cohesive crack growth[END_REF] and relations [START_REF] Belytschko | Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment[END_REF], the linearized problem resulting from a Newton-Raphson strategy applied to Equations ( 9) and ( 10), at a given iteration k,i s :

D u P: (∇ r u * ) T d -D u + c (g n )g n n•[[u * ]] d • u k+1 -D k + c (g n )g n n•[[u * ]] d • k k+1 =- P: (∇ r u * ) T d + + c (g n )g n n•[[u * ]] d + f•u * d + t t•u * d (23) D u + c -1 n { k -(g n )g n } * d • u k+1 + D k + c -1 n { k -(g n )g n } * d • k k+1 =- + c -1 n { k -(g n )g n } * d (24)
We replace the augmented lagrangian contact multiplier by its Equation ( 7) and after some analytic calculations one obtains the following linear system:

*P *u • u k+1 : (∇ r u * ) T d - + c k+1 n•[[u * ]] d + + c n n•[[Du k+1 ]] n•[[u * ]] d = f•u * d + t t•u * d + + c ( k -n [[ u k ]] • n)[[ u * ]] n d + P(u k ): (∇ r u * ) T d (25) - + c 1- n k+1 * d - + c [[ u k+1 ]] • n * d = + c 1- n k * + [[ u k ]] • n * d (26)
The linearized formulation ( 25) and ( 26) can be expressed in a matrix form as follows, omitting the iteration indexes for the sake of clarity:

K meca +A u A T AC u = L meca +L 1 cont L 2 cont ( 27
)
where K meca is the mechanical stiffness matrix, A u is the augmented stiffness matrix due to the contact, A is the matrix linking the displacement terms to the contact ones, C is the matrix allowing to obtain the contact pressure in case of no-contact, L meca is the second member related to the internal forces, L 1 cont and L 2 cont are the second members due to the contact. In order to express the elementary terms of ( 27) related to the contact contributions, a unified notation was adopted for the shape functions of master and slave domains

I (x 1 ) = s I (x 1 ) 1<I <m s -m I -m s ( x1 ) m s <I <m s +m m ( 28
)
where m s and m m denote the number of slave and master dof, respectively. If we consider also n c the number of contact dof, and introduce in ( 25) and ( 26) the general approximations for the gap calculus [START_REF] Jin | Two dimensional mesh generation; adaptive remeshing and refinement[END_REF] and contact pressure [START_REF] Babuška | The partition of unity method[END_REF], we obtain

A u I,J = + c n I J (n⊗n) d (29)
with I, J = 1, m s +m m A I,J =-

+ c I J n d (30)
with

I = 1, n c and J = 1, m s +m m C I,J =- + c 1- n I J d (31)
with

I, J = 1, n c L 1 cont I = + c ( -n [[ u ]] • n) I n d ( 32
)
with I = 1, m s +m m and

L 2 cont I = + c 1- n + [[ u ]] • n I d (33)
with I = 1, n c .

LBB stabilization algorithm

It is well known that the existence and uniqueness of the solution to a discrete mixed formulation of a contact problem require to cautiously define the interpolation space of the Lagrange multipliers, in accordance to that of the displacement, in order to satisfy the inf-sup LBB condition. In particular, the satisfaction to this condition is related to the ratio

r h = dim W h dim Q h ( 34
)
where W h is the discrete space of the displacement dof on the contact surface and Q h is the discrete space of Lagrange multipliers associated with the contact constraints. For an FEM approach, several configurations for these discrete spaces are tested in [START_REF] Solberg | An analysis of dual formulations for the finite element solution of two-body contact problems[END_REF] and for the ones which satisfy the LBB condition for a two-body contact problem, r h >1(r h = 1 for the Signorini contact problem). For the X-FEM approach presented in this paper, an additional difficulty comes from the difference between the geometric support of the displacements (nodes) and the contact constraints (intersection points between the mesh and the level set). In particular, the ratio r h cannot be established on a unit element collection, but requires to consider the whole boundary-value problem. To overcome such limitations, Géniaut et al. [START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF] proposes an algorithm in the case of small perturbation hypothesis which increases r by imposing global linear relationships between the Lagrange multipliers. Details about the implementation of this solution and time succeeding references of its development are given below.

In [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF] an algorithm is proposed to reduce the Lagrange multiplier space by imposing Dirichlet conditions on an X-FEM interface. The superfluous multipliers are linked to the other multipliers, imposing equality or linear relations between the Lagrange multipliers. This algorithm tends to impose more equality relations than linear relations, so that the linear approximation of the multipliers [START_REF] Babuška | The partition of unity method[END_REF] is degraded. This algorithm was improved in [START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF] in order to make it more efficient and reverse this trend by imposing more linear relations in order to satisfy most of the time equation [START_REF] Babuška | The partition of unity method[END_REF]. First, one defines S e the set of edges that are strictly cut by the interface. The interface being represented by a level set function , an edge is strictly cut by the interface if the product of the values of at both vertices of the edge is strictly lower than 0. Note that if the interface coincides with one of the vertices of an edge, this edge is not in S e . One defines also two sets of nodes: N + e and N - e which include the vertices of the edges in S e , 'below' and 'above' the interface. Therefore, N + e (resp. N - e ) contains the nodes n connected by the elements of S e with (n)>0 (resp. (n)<0). Then the algorithm searches for S ve , a minimal sub-set of S e , which permits to connect nodes in N + e to nodes in N - e . This set is not necessarily unique. Each node is affected by a weight corresponding to the number of edges connected to it. For each edge we compute the minimum score of the two summit nodes connected to it. If this score is equal to 1, the edge is vital. If a node is connected to several edges, it should have at least one vital edge. We subtract from the set of nodes all nodes related together to the vital edges identified before. With this construction node 10 (see Figure 7) is left for which we choose the vital edge as the shortest edge. This choice reduces the area of P0 (i.e. constant) multiplier approximation with respect to previous algorithm [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF]. Edges in S ve are called 'vital edges', because if one of these edges is missing, at least one node in N + e ∪ N - e will be isolated on one side. It is important to notice that S e is composed by some disconnected edges and some connected edges. Those groups of connected vital edges are extracted from S e . Note that in such a group, all the edges are connected by a unique node (See Figure 7). Finally, one imposes the relations between the multipliers. All the multipliers linked to edges within a group are imposed to be equal. No relation is imposed on the Lagrange multipliers linked to single vital edges. These multipliers are free. The other multipliers are on non-vital edges. They are not essential for the contact pressure approximation. Therefore, one imposes them to be a linear combination of multipliers on vital edges. These linear combinations are determined by the following procedure. Let e , e ∈ S e \ S e be the Lagrange multiplier, which lies on a non-vital edge e, at the intersection between e and the level set function . For each vertex of the edge e, one searches the closest k lying on an edge connected to e.I f k1 and k2 are these two values, then a linear relation is imposed between e , k1 and k2 taking the following form: Figure 7 shows a 2D example for this algorithm.

The reason why the algorithm used in small displacements is still valid for large displacements comes from the fact that the displacement approximation space used for the gap expression is at least as rich as the displacement approximation space used for the displacement on the slave surface. As a consequence, the LBB satisfaction that has been observed (numerical inf-sup test) with the algorithm of [START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF] can be extended here to large sliding by considering the Lagrange multipliers on the slave contact surface only. We will show in particular in the second numerical example of the next section that this algorithm is still reducing pressure oscillation.

EXAMPLES

In this section we present three 2D numerical applications of our new approach with frictionless large sliding contact under the assumption of plane-strain elasticity. For all three problems now detailed, the material behavior is assumed to be Kirchhoff-St Venant linear elasticity with Young's modulus set to 100 MPa and Poisson's ratio set to 0.3. The 2D elements are bilinear with four nodes. The corresponding contact elements are two-node segments. The integration rule for the contact contributions is trapezoidal with integration points at the location of contact segment nodes.

Frictionless sliding along horizontal interfaces under small compression

This problem model is represented in Figure 8 crossing the whole plate (Figure 8(a)), symmetrically placed with respect to the middle horizontal and distanced by b = 4 m. We define the contact surfaces of the middle block as masters and the opposite contact surfaces belonging to the lateral blocks as slaves. Displacement boundary conditions are applied on the upper and lower faces of the plate, Depl Y = 10 -3 m, to ensure the closure of the contact interfaces and the horizontal displacements on the right edges of the top and bottom blocks are restricted for the well-posedness of the problem. The middle block is loaded by imposing a large displacement boundary condition on the left side, Depl X = 2m, so that it slides along the contact surfaces. The numerical resolution of this problem is realized with a coarse discretization of nine by four elements along the length and width, respectively. . Evolution of the relative error for the pressure magnitude of the X-FEM solution compared with that of the FEM solution; X is the distance from the left extremity of the surface; P xfem and P fem refer to the contact pressure for X-FEM and FEM model, respectively.

The final displacements and yy stress fields are shown in Figure 8(b) where the deformed mesh is obtained from the computed kinematics after a triangulation on each side of the discontinuity. Note that the vertical displacement magnitude is multiplied 500 times in order to point out the deformed shapes of the contact surfaces. As expected, compressive stress is found in contact zones. In the free zones (at the rear of the moving block) the numerical values of yy are small, but not exactly equal to zero because the finite element framework only approximately satisfies tractionfree boundary conditions. Numerical values for the contact pressure on the interface surfaces are compared with the ones obtained by modeling the same test with standard FEM (i.e. with contact interfaces matching the mesh) and using the same mesh density in Figure 9. Each curve represents the pressure distribution along one of the slave contact surfaces for a given linear increment of Depl X . Note that all the Dirichlet conditions are following the same linear loading relationship from 0 to their final values. The resulting pressure distribution is, as expected, vanishing to zero behind the middle block as this one advances in the right direction. The relative error between the FEM and X-FEM approaches is less than 7% (Figure 10) being related to the different design of the new hybrid contact element compared with the FEM one, especially because of a lack in the application of the Dirichlet conditions on the boundaries of X-FEM enriched elements.

Frictionless sliding oblique to mesh under small compression

In the second example we analyze the algorithm efficiency to ensure the stability of the contact formulation with respect to the LBB condition (also called inf-sup condition), see [START_REF] Bathe | The inf-sup condition and its evaluation for mixed finite element methods[END_REF][START_REF] El-Abbasi | Stability and patch test performance of contact discretization and a new solution algorithm[END_REF].T h e geometry of this problem is unchanged with respect to the first one, except for the discontinuities that are now oblique contact interfaces in order to cross the mesh over several rows. This configuration is favorable to triggering contact pressure oscillations if a naive algorithm keeping all the contact multipliers is employed [START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF]. Such pathology is well known to occur when contact imposed on too many points is over-constraining the solution and is related to the non-respect of the LBB condition. The LBB fulfilling algorithm was already implemented to ensure the LBB condition respect for the small perturbation X-FEM contact formulation [START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF], derived and improved from the algorithm introduced by Moës et al. [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF]. The boundary conditions (Figure 11(a)) are similar to that of the first problem. The numerical solution is computed with a finer mesh. The final displacement and yy stress fields are shown in Figure 11(b).

The problem is first solved with a naive algorithm and then with the LBB fulfilling algorithm. The contact pressure distributions are represented for the first case in Figure 12(a) and for the second case in Figure 12(b). One can observe large oscillations for the contact pressure corresponding to the contact dof placed in the positions where the discontinuity goes from one row of elements to another one. In this case the contact multiplier space is too rich and the LBB fulfilling algorithm, when activated, imposes some linear relations in between these contact dof leading to a more regular contact pressure distribution (Figure 12(b)).

Fric tionless large rotation under small compression

The last numerical example presented here tests the efficiency of our new approach in case of non-planar contact interfaces, when the outward normal vector defining the contact is varying with the sliding. For this example, the gap d n between the bodies in contact (see Equations ( 4) and [START_REF] Jin | Two dimensional mesh generation; adaptive remeshing and refinement[END_REF]) is forced to be zero. This condition enables to deactivate the loop on active constraints, the third stage for solving the contact problem as mentioned in Section 3.2.

The boundary-value problem is shown in Figure 13. An elastic plate is delimited by an inner and outer radii R 1 and R 2 , respectively. This plate is surrounding a concentric inner rigid disk of radius R 1 . We apply a large rotation on the mesh of the inner disk, while the outer boundary of the plate is submitted to a homogeneous pressure of low magnitude. Considering the problem in small strains and stresses, we can compare the solution to one of a classical elastic boundaryvalue problem in small deformation hypothesis as follows (see for example reference solutions in [START_REF] Salenc ¸on | Handbook of Continuum Mechanics[END_REF]). Let us denote (r, , z) the cylindrical coordinates with r the radius from the origin, the polar angle and z the coordinate along the out of plane axis. From the equilibrium and because of the problem symmetry, the radial displacement in the deformable zone around the rigid inclusion is

u r = ar + b r ( 36
)
Figure 13. Boundary conditions for the third example; denotes the imposed rotation on the inner disk and is an angular coordinate of a point on the outer slave surface with respect to point 1 measured anticlockwise.

where a and b are unknowns. The elastic non-zero Cauchy stress components are

rr = A - B r 2 (37)
= A + B r 2 (38) zz = 2 a ( 39
)
with

A = 2( + )a (40)
B = 2 b (41)
where ( , ) are the Lamé coefficients. From the above equations, the boundary condition can be set as follows:

aR 1 + b R 1 = 0 (42) A - B R 2 2 =-p (43)
In the sequel, p is set to 1 MPa (i.e. 1% of Young's modulus), R 1 is set to 3 m and R 2 is set to 20 m. The numerical model consists of a 2D squared plate with an edge dimension set to 40 m and centered as the circle and the plate. This choice allows us to control the refinement of the finite element model since a regular mesh can be used. These meshes are composed by regular grids of n by n bilinear elements, where n is set to either 20, 40 or 80 to analyze the influence of the mesh refinement. The radial displacement of Equation ( 36) is applied on the outer boundary of the square. We set a concentric discontinuity defined initially by a circle of radius R set to 8 m. A rigid rotation of 90 • is progressively imposed on the inner disk. With a convenient choice of the level set function defining the discontinuity, the slave surface is the outer side and because of frictionless contact conditions it should not be influenced by the rotation of the inner disk. Solving Equations ( 40)- [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF] for the numerical values of the problem, the radial Cauchy stress is found to be -1.047 MPa at a radius corresponding to that of the discontinuity and will be denoted as * rr in what follows. A quite similar model, called 'journal bearing' was considered in [START_REF] Belytschko | Arbitrary discontinuities in finite elements[END_REF] for capturing the lowest eigenvalue and the eigenvector corresponding to rigid-body rotation by X-FEM.

The main advantage of this simple model is to enable the comparison of the numerical and analytical solutions, in particular, because it avoids boundary heterogeneous conditions at the edge of contact surfaces. The final displacement is shown in Figure 14(b) and is dominated as expected by the rigid rotation of the inner disk. From Equations ( 39) to [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF], one can find that the trace of the Cauchy stress tensor is homogeneous and equal to p. It is represented for the numerical solution in Figure 14(a). The differences between the numerical and the analytical solutions are mainly located in the vicinity of the discontinuity and can be explained as follows. Because the interpolation of the displacement field is bilinear, a finite rotation increment results in mesh incompatibility. The 'sliding contact' conditions then introduce an artifact in the deformation by pushing (resp. pulling) the slave body when the master body is penetrating it (resp. when it is not). This effect is clearly responsible for the oscillations in traction and compression of the contact pressure. To alleviate these artifacts, it would be preferable to use a higher-order spatial interpolation in order to minimize the mesh incompatibility. The relative error of the contact pressure compared with * rr is given at the final state in Figure 15, at the contact points distributed along the slave surface. On this figure, one can see that when the mesh is refined, the mean estimate is improved and the deviation is reduced. Thus, the X-FEM solution converges toward the analytical one. Note the 90 • period profile on Figure 15 which results from the symmetries of the squared shape of the original mesh. For a given configuration obtained with a value of the imposed rotation , the error in contact pressure is measured along n p sample points of the slave surface as:

err = 1 n p i=1,n p i * rr -1 2 (44)
It is represented on the log-log diagram of Figure 16 with respect to the three mesh refinements: 20×20, 40×40 and 80×80 elements. The differences of the error magnitude with the different configurations given by are related to the already mentioned artefacts due to mesh incompatibility. However, a linear convergence of the X-FEM solution toward the analytical one is obtained with the mesh refinement.

In Figure 17 an evaluation of the Newton iteration convergence is given for this last test presented here. Note that we use a generalized Newton tangent operator. Using logarithmic scales and for the final loading step, as an example, the convergence slope for the absolute residual error is about 1.3 instead of a quadratical convergence as one would expect. The explanation for this convergence, more linear than quadratic, comes from the lack of the normal vector variation update after each Newton iteration. In fact, as mentioned before in Section 3.2, at this state of the numerical implementation in our code, the normal update is computed only during the geometrical update leading to a slower Newton iteration convergence compared with the one obtained for the first two numerical examples presented in this paper where the normal direction does not change and the Newton iteration procedure converges quadratically.

CONCLUSION AND PROSPECTS

A new X-FEM approach for frictionless large sliding contact along a discontinuity has been presented which the authors would like to consider as a preliminary step to many new applications taking advantage of X-FEM in the context of large deformations. The central development is a dedicated hybrid contact element that gathers the mechanical interaction between two contacting enriched elements. This contact pairing from which the kinematic gap between the contacting bodies is updated with the relative displacement of the bodies. Three 2D plane-strain elasticity examples with enriched four-node bilinear elements, using this methodology, have been satisfactorily confronted to an analytical solution or to a standard FEM solution. Of particular importance is the result that the LBB fulfilling algorithm developed in [START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF], in the case of small perturbations, is still efficient to avoid oscillations in the contact pressure field with large sliding. For the last example, in spite of obvious limitations due to the use of enriched bilinear elements with a constant normal per element to represent a curved discontinuity, the convergence of the solution with mesh refinement has also been verified.

Several perspectives and improvements are considered for the future. First, the extension to 3D should be realized without changing the principles presented in this paper. In addition, as illustrated by the third example, the application of this method to complex kinematics with geometrical nonlinearities should require the extension to quadratic or other higher-order displacement interpolation of elements in order to reduce drastic perturbations of the solution because of mesh incompatibilities during sliding. One must also account for the variation of the outward normal in the tangent operator during the Newton iterations, so that its update must be achieved before each Newton iteration to ensure quadratic convergence. Such difficulty is expected in particular for sliding along non-smooth surfaces like geological faults in sedimentary basin. Accounting for a rheology on the interface is also of interest; more work is needed on the use of a frictional Coulomb model with a focused attention on the stability of the model, namely the satisfaction of the LBB condition (see details in [START_REF] Géniaut | A stable 3D contact formulation for cracks using X-FEM[END_REF]), in the transition between stick and slip zones. To complete the interface rheology, non-linear bulk constitutive law should be implemented to represent the behavior of sedimentary rocks at geological time scale, like for example the three-invariant elasto-plasticity model of [START_REF] Matsuoka | Stress deformation and strength characteristics of soil under three different principal stresses[END_REF] used in [START_REF] Sanz | Mechanical aspects of thrust faulting driven by far-field compression and their implications for fold geometry[END_REF], or the elasto-plasticity model in [START_REF] Guiton | Activation of diffuse discontinuities and folding of sedimentary layers[END_REF] with yield surfaces corresponding to the activation in sliding and opening mode of pervasive fracture sets, or the poroelastoplastic model of [START_REF] Bernaud | Evolution of elastic properties in finite poroplasticity and finite element analysis[END_REF] which derives from a micromechanical approach.

Last but not the least, we would like to recall the interest of this X-FEM approach if one wants to change the model topology, either to update the data (e.g. the geometry of discontinuity surfaces) to better fit observations or if the problem of interest involves the propagation of the discontinuity. For the latter case, it is necessary to extend this work by considering appropriate enrichment functions for the asymptotic displacement field near the crack tip. To apply this method on a 3D model of geological faults, it is also necessary to design a new X-FEM enrichment based on an appropriate implicit representation of the discontinuity geometry and of the material distribution, in particular to handle arbitrarily full or partial intersections between the discontinuities.

Figure 1 .

 1 Figure 1. (a) Domain and boundary notations for a contact problem and (b) evaluation of the contact gap between the slave point x + and its projection x+ .

Figure 2 .

 2 Figure 2. X-FEM enrichment of nodes near a discontinuity; if a nodal support, e.g. the gray domain around the bold circle node, is fully cut, the enrichment function is of Heaviside type with opposite signs on each side of the discontinuity.

Figure 3 .Figure 4 .

 34 Figure 3. Level sets description for X-FEM.

Figure 5 .

 5 Figure 5. Illustration of the X-FEM contact search algorithm with large sliding along a discontinuity; circles are nodes of the enriched elements, squares are intersection points between the mesh and the discontinuity and crosses are integration points on the slave surface (case of Gauss integration rule).

Algorithm 1 Figure 6 .

 16 Figure 6. X-FEM hybrid contact element with reference mappings for the kinematic on the left and for the contact unknowns on the right.

35 )Figure 7 .

 357 Figure 7. Example of edges cut by an interface. Dotted edges are non-vital. Groups of connected vital edges are circled. The final multiplier approximation (P0o rP1) is plotted on the interface.

  Figure 8. (a) Boundary conditions and (b) deformed mesh with isocontours of the yy stress field; vertical displacement magnitude is multiplied by 500 times.

Figure 9 .

 9 Figure 9. Contact pressure distribution along one of the slave contact surfaces for FEM: (a) X-FEM and (b) analysis; X is the distance from the left extremity of the surface.

Figure 11 .

 11 Figure 11. (a) Boundary conditions and (b) deformed mesh with isocontours of the yy stress field.

Figure 12 .

 12 Figure 12. Contact pressure distribution along one of the slave contact surfaces with a naive algorithm (a) with the LBB fulfilling algorithm and (b) X is the distance from the left extremity of the surface.

Figure 14 Figure 15 .

 1415 Figure 14. (a) Isocontours of the trace of the Cauchy stress tensor at the final state for 80×80 elements and (b) deformed mesh for = 45 • with isolines of the displacement norm for 80×80 elements.

igure 16 .Figure 17 .

 1617 Figure 17. Convergence of the absolute residual error for the Newton-Raphson iterations; e i denotes the absolute error for the dof i, n being the total number of dof.

Table I .

 I dof table for a X-FEM hybrid contact element.

	Node	Standard dof	S tandard dof	Enriched dof	Enriched dof	Contact dof
	N1	DX	DY	HX	HY
	N2	DX	DY	HX	HY
	N3	DX	DY	HX	HY
	N4	DX	DY	HX	HY
	N5
	N6
	N7
	N8
	N9	DX	DY	HX	HY
	N10	DX	DY	HX	HY
	N11	DX	DY	HX	HY
	N12	DX	DY	HX	HY

ACKNOWLEDGEMENTS

This work has been realized during a one-year post-doctorate of the first author at the Geology-Geophysics-Geochemistry Division of IFP and was financed and supported by the Scientific Direction of IFP. Developments have been realized with the software Code Aster (version 9.0) which is a general code for the study of mechanical behavior of structures distributed by EDF R&D under a GNU GPL license on the web site http://www.code-aster.org. The authors are grateful to Jacques Pellet, Damijan Markovic and Mickael Abbas of EDF R&D for their support in the numerical implementation of the method and would like to acknowledge Jean-Franc ¸ois Barthélémy, Jean-Marc Daniel and Gilles Perrin of IFP for their constructive remarks and helpful discussions.

The final deformations of the numerical examples were illustrated with GNU GPL software Gmsh available at http://www.geuz.org.

Contract/grant sponsor: Scientific Direction of IFP